A Meshfree Generalized Finite Difference Method for Solution Mining Processes
Isabel Michel, Tobias Seifarth, Joerg Kuhnert, Pratik Suchde
TThis is a pre-print of an article published in Computational Particle MechanicsThe final version can be found with the doi: 10.1007/s40571-020-00353-2August 26, 2020
A Meshfree Generalized Finite Difference Method for SolutionMining Processes
Isabel Michel ∗ , Tobias Seifarth , J¨org Kuhnert , Pratik Suchde Fraunhofer Institute for Industrial Mathematics ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany
SUMMARYExperimental and field investigations for solution mining processes have improved intensely in recent years. Due totoday’s computing capacities, three-dimensional simulations of potential salt solution caverns can further enhancethe understanding of these processes. They serve as a “virtual prototype” of a projected site and support planningin reasonable time. In this contribution, we present a meshfree Generalized Finite Difference Method (GFDM)based on a cloud of numerical points that is able to simulate solution mining processes on microscopic as wellas macroscopic scales, which differ significantly in both the spatial and temporal scale. Focusing on anticipatedindustrial requirements, Lagrangian and Eulerian formulations including an Arbitrary Lagrangian-Eulerian (ALE)approach are considered.KEY WORDS: Meshfree Methods ; Generalized Finite Difference Method ; Lagrangian Formulation ;Arbitrary Lagrangian-Eulerian Formulation ; Solution Mining
1. INTRODUCTIONThe basic motivation of this research is to provide a method that is able to simulate the long-termdevelopment of a salt cavern during a double-well solution mining process. Solution mining is used toextract underground water-soluble minerals such as salt and potash. A double-well convection processhas been a preferred choice for solution mining due to it’s large recovery rate [2, 33]. As the namesuggests, this involves the use of two boreholes or wells for the extraction process: an injection well,and a recovery or extraction well. For the extraction of salt, fresh water is pumped into a salt depositthrough the first “injection” well. Salt present in the cavern dissolves in the water to produce a saturatedbrine solution. This is then extracted at the second “extraction” well. A schematic of this process isshown in Fig. 1. The main direction of dissolution is vertical, which is controlled by alternate lifting ofthe injection and the extraction well.In this work, we focus on modeling of the fluid flow involved in such a double-well solution miningprocedure, including the formation of the salt-water solution. An essential aspect of this is to accuratelymodel the long-term geometrical evolution of the salt cavern. This is needed to steer the actual processof solution mining, in terms of, for example, determining when and at what rate the injection andextraction wells are raised. However, numerically modeling this is very challenging, as it is a highlydynamic three-dimensional process involving different spatial and temporal scales. Over the time scaleof several years, as salt in the cavern dissolves in the water, the cavern starts to erode, causing significantdeformations in it’s overall shape. However, the dissolution process relies on a smaller time scale ofseveral minutes. On the spatial scale, the former involves the modeling of the entire salt cavern, whilethe latter is more localized and is relevant near the cavern walls. In the present work, we model boththese processes in separate simulation setups. A macroscopic simulation is carried out to model theevolution of the cavern over many years. This is done on actual salt cavern geometries. The computation ∗ Correspondence to: E-mail: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] A ug Figure 1. Schematic of double-well solution mining (adapted from [25]). of the diffusion rate of the salt (and related minerals) to be used in these macroscopic simulations aredone in a separate simulation, in the so-called microscopic setup. This involves simulations over thesmaller time scale of a few minutes, and over representative geometries several orders of magnitudesmaller than the size of the salt cavern in the macroscopic simulations.Over the past few decades, meshfree simulation methods have emerged as an alternative to theconventionally used mesh-based simulation procedures, especially in the context of modeling fluidflow. The advantages of meshfree methods are most notable for applications with complex domains,or those with moving geometry parts, free surfaces, phase boundaries, or large deformations. Whilemodeling each of the latter cases with mesh-based methods, the highly dynamic nature of thesimulations often requires an expensive global remeshing procedure. On the other hand, movingLagrangian and semi-Lagrangian procedures fit in naturally with meshfree methods, making thesimulation of dynamic geometries or phase boundaries a lot easier. In the application at hand, themodeling of the changing domain during the long-term evolution of the salt cavern falls into thiscategory. We thus choose a meshfree approach.In this paper, we use a meshfree Generalized Finite Difference Method (GFDM) [5, 7, 12, 20]based on a cloud of numerical points. This method has already been successfully applied in variousCFD and continuum mechanics applications. Prominent examples include water crossing of cars, waterturbines, hydraulic valves, soil mechanics, metal cutting, and mold filling [10, 14, 21, 23, 31, 32]. Themethods presented here are part of the in-house developed software MESHFREE † , which combinesthe advantages of GFDM and the fast linear solvers of SAMG [22].We start by using a Lagrangian formulation where the discretizing point cloud moves according tothe flow velocity [10, 13]. This results in an accurate and natural transport of physical information.The basic physical model consists of the conservation equations for mass, momentum, and energy.For solution mining processes, we extend it by the standard k - ε turbulence model and equations forthe concentration of the occurring species (see Sect. 2). The GFDM specific numerics are presentedin Sect. 3 with special focus on the Lagrangian and Eulerian formulations. Microscopic simulationsare presented in Sect. 4, and are used to determine the necessary effective model parameters of themacroscopic problem which follows. For macroscopic simulations, the Lagrangian formulation leadsto a significant restriction of the time step size due to the explicit movement of the point cloud. Toenable simulations in reasonable time, an Eulerian formulation should be preferred in this context.Here, the point cloud is fixed and convective terms represent the transport of physical information.The movement of the boundary of the salt cavern is based on the solution rate of the salt in the flowingwater. To accurately handle this moving boundary, interior points close to the boundary are subject to an † The basic underlying physical model is given by the conservation equations of mass, momentum, andenergy in Lagrangian formulation. d ρ dt + ρ · ∇ T v = , (1) ddt ( ρ · v ) + ( ρ · v ) · ∇ T v =( ∇ T S ) T − ∇ p + ρ · g , ddt ( ρ · E ) + ( ρ · E ) · ∇ T v = ∇ T ( S · v ) − ∇ T ( p · v )+ ρ · g T · v + ∇ T ( λ · ∇ T ) , for density ρ , velocity v ∈ R , stress tensor S ∈ R × (deviatoric part, i.e. tr ( S ) = p ,body forces g ∈ R , total energy E = c v · T + · ( v T · v ) , heat capacity c v , temperature T , and heatconductivity λ . Further, ddt = ∂∂ t + v T ∇ denotes the material derivative, and ∇ = ( ∂∂ x , ∂∂ y , ∂∂ z ) T denotesthe nabla operator.In general, the stress tensor is split into its viscous and solid parts by S = S visc + S solid [10, 13]. Forthe present application, the stress tensor is purely viscous, S solid = . The viscous part is defined by S visc = ( η + η turb ) · (cid:18) ∇ v T + ( ∇ v T ) T − · ( ∇ T v ) · I (cid:19) , (2)where I ∈ R × is the identity.To incorporate turbulent effects, the standard k - ε turbulence model [18] is considered for turbulentkinetic energy k and turbulent dissipation ε dkdt = ρ · ∇ T (cid:18)(cid:18) η + η turb σ k (cid:19) · ∇ k (cid:19) − ε + ρ · ( P pr + P b ) , (3) d ε dt = ρ · ∇ T (cid:18)(cid:18) η + η turb σ ε (cid:19) · ∇ ε (cid:19) − C ε · ε k + ρ · C ε · ε k · ( P pr + C ε · P b ) , where η is the laminar viscosity, and η turb = ρ · C η · k ε is the turbulent viscosity. Fluctuating dilatationand source terms are omitted [18]. The turbulent production rate is defined by P pr = η turb ·(cid:107) ∇ v T (cid:107) withvon Mises matrix norm (cid:107)·(cid:107) M . The turbulent buoyancy is given by P b = − ρ · η turb Pr turb · ∂ρ∂ T · ( g · ∇ T ) . For thismodel, well-established values for the constants are used σ k = . σ ε = . C ε = . C ε = . C ε = − . C η = .
09, and turbulent Prandtl number Pr turb = .
85. Furthermore, a logarithmic wallfunction is used in the vicinity of walls.In order to simulate solution mining processes, the basic model above is extended by convection-diffusion equations to represent the different minerals or species present in the salt mixture. For the3concentration c i of species i = , . . . , N with effective diffusion coefficient D i , eff , we have dc i dt + c i · ∇ T v = ∇ T ( D i , eff · ∇ c i ) . (4)In the Eulerian formulation, the material derivative is replaced by its definition, i.e. ddt = ∂∂ t + v T ∇ . Consider the general form of the equation of state ρ = ρ ( T , c , . . . , c N ) , (5)where density depends on the temperature and each of the concentrations. Based on the formulation in[16, 17], the density of a solution of N species in water is given by ρ sol = (cid:32) w H O ρ H O + N ∑ i = w i ρ apparent , i (cid:33) − , (6)where w H O and w i are the mass fraction of water and species i , respectively. Additionally, w H O + ∑ Ni = w i = ρ H O (7) = ((((( A · T + A ) · T + A ) · T + A ) · T + A ) · T + A ) + A · T , with A , . . . , A and C ∗ , . . . , C ∗ defined according to [17]. The apparent density of species i is given by ρ apparent , i (8) = ( C ∗ · ( − w H O ) + C ∗ ) · exp (cid:0) . · ( T + C ∗ ) (cid:1) ( − w H O ) + C ∗ + C ∗ · T . The mass fractions w i are defined by the concentrations c i as w i = c i ∑ Ni = c i + ρ H O . (9)The viscosity of the solution, η sol ( T , c , . . . , c N ) , and its heat capacity c v , sol ( T , c , . . . , c N ) aremodeled in a similar manner. For the viscosity of a solution of N species in water, we use a modifiedversion of the Arrhenius equation η sol = ( η H O ) w H2O N ∏ i = ( η i ) w i , (10)where the viscosity of water depends on the temperature as η H O = T + ( . · T + . ) T + . . (11)Furthermore, the viscosity of species i is given by η i = exp (cid:18) V ∗ ( − w H2O ) V ∗ + V ∗ V ∗ T + (cid:19) V ∗ ( − w H O ) V ∗ + V ∗ , . . . , V ∗ according to [15]. 4 Table I. Model constants for the sodium chloride solution according to [15, 16, 17]. i A i C ∗ i V ∗ i B ∗ i − . × − . × − . . − . . × − . . − . . × − . . . . × − . . − . . . . . . . . N species in water c v , sol = w H O c v , H O + N ∑ i = w i c v , i . (13)Furthermore, the heat capacity of species i is modeled by c v , i = B ∗ exp ( a ) + B ∗ ( − w H O ) B ∗ , (14)where a = B ∗ T + B ∗ exp ( . · T ) + B ∗ ( − w H O ) and constants B ∗ , . . . , B ∗ are according to [16].We use quadratic interpolation (extrapolation) for the definition of the heat capacity of water. Assumegiven temperatures T , T , T , with T = T + ∆ T and T = T + ∆ T ( ∆ T > c v , T , c v , T , and c v , T are also assumed given. Then, the heat capacity of water atarbitrary temperature T is determined by c v , H O = c v , T + ( c v , T − c v , T ) T − T T − T (15) + c v , T − c v , T + c v , T T − T T − T (cid:18) T − T T − T − (cid:19) . T , T , T are chosen adaptively with ∆ T = ◦ C, depending on the value of T . The range of c v , T k values, k = , ,
3, are taken from [16], which provides the values between 0 ◦ C and 95 ◦ C.We restrict the study in this paper to sodium chloride as the species of interest. All the modelconstants mentioned in this section for sodium chloride are summarized in Table I.3. NUMERICS BASED ON GFDM
In the GFDM approach, the computational domain is discretized by a cloud of numerical points. Thepoint cloud is composed of NP = NP ( t ) number of points, which includes points in the interior of thedomain, and those at the boundary. The initial seeding of these point clouds is done by a meshfreeadvancing front technique, details of which can be found in [19, 24]. The density of the point cloud isgiven by a sufficiently smooth function h = h ( x , t ) , the so-called interaction radius or smoothing length.Thus, h prescribes the resolution of the point cloud. It is also used to define the neighborhood of eachpoint. For a point x j in the point cloud, all approximations are performed using only nearby pointswithin a distance h from it. This set of nearby points is referred to as the neighborhood or support of x j , and is denoted by S j = { x l : (cid:107) x l − x j (cid:107) ≤ h ( x j ) } .To ensure a sufficient quality of the point cloud, it is ensured that no two points are closer than r min h apart, and that every sphere of radius r max h in the domain has at least one point. Thus, the5inter-point distance between each point and its nearest neighbor lies in the range ( r min h , r max h ) . Wefollow conventionally used values of these parameters in Lagrangian meshfree GFDM literature, andset r min = . r max = . −
50 points in each interior neighborhood,with lesser at and near the boundary.For the Lagrangian and ALE formulations, the movement of (parts of) the point cloud with thefluid velocity can result in the minimum and maximum inter-point distance criteria being violated.This happens in form of accumulation or scattering of points which would reduce the quality of thenumerical results. To prevent this, points are added in holes containing insufficient points, and aremerged in regions of accumulation. This method of fixing distortion is entirely local, and much cheaperthan the remeshing done in mesh-based methods. Details about these procedures of adding and deletingpoints follow from [4, 13, 14, 24, 28].
GFDMs generalize classical finite differences to arbitrarily spaced point clouds, using a specializedweighted moving least squares approach. Consider a function φ defined on each point of the pointcloud. At each point x j , numerical derivatives of φ are defined as a linear combination of functionvalues in it’s neighborhood ∂ ∗ φ ( x j ) ≈ ˜ ∂ ∗ j φ = ∑ l ∈ S j c ∗ jl φ l , (16)where ∗ = x , y , z , ∆ , . . . denotes the derivative of interest, ∂ ∗ is the continuous differential operator, ˜ ∂ ∗ j is the numerical differential operator at point x j , and φ l = φ ( x l ) . The numerical differential operatorsare thus given by the coefficients c ∗ jl , which are independent of the function being differentiated. Theyare computed by a norm minimization process that ensures that monomials up to a specified order aredifferentiated exactly. ∑ l ∈ S j c ∗ jl m l = ∂ ∗ j m , ∀ m ∈ M , (17)min ∑ l ∈ S j (cid:18) c ∗ jl W jl (cid:19) , where M is the set of monomials being differentiated exactly. To compute the Laplacian, themonomials are complemented by the delta function to control the central stencil value c ∆ j j , whichimproves stability in the pressure Poisson equations [26]. In the present work, we consider monomialsup to the order of 2. The weighting function W is defined such that neighboring points with the smallestdistance to the considered point obtain the highest weight. In the present work, we use a truncatedGaussian weighting function W jl = exp (cid:18) − c W (cid:107) x j − x l (cid:107) h j + h l (cid:19) , if x l ∈ S j , elsewhere (18)for a constant c W >
0. We note that the same differential operators as defined above can also beequivalently derived by minimizing errors in Taylor expansions [26].Using this procedure, we compute numerical gradient operators, and a numerical Laplacian. Formore details on the computation of the differential operators, we refer to [4, 14, 26].
A strong form discretization of the physical model (Sect. 2) is doneusing the numerical differential operators defined above, and a chosen time integration scheme. Forsimplicity, the following considerations are based on a first order time integration.6Starting with the Lagrangian formulation, equations (1) can be rewritten as d ρ dt = − ρ · ∇ T v , (19) d v dt = ρ · ( ∇ T S ) T − ρ · ∇ p + g , ( ρ · c v ) · dTdt = ∇ T ( S · v ) − ( ∇ T S ) · v − p · ∇ T v + ∇ T ( λ · ∇ T ) . To improve readability, we henceforth use the shorthand ρ = ρ sol and c v = c v , sol .Together with equations (2)–(4), this is the starting point of the numerical discretization. Thecontinuous spatial derivatives are replaced by their least squares approximated counterparts describedin Sect. 3.2. We consider the superscript n + n for the current one,giving the time step size ∆ t = t n + − t n . Below, we explain each of the steps of the discretization in theLagrangian formulation for the microscopic scale simulations. Most of the steps are the same also forthe macroscopic scale simulations, and the few differences are explained in Sect. 5. Step 1. Point cloud movement
The discretization procedure begins by moving the point cloud according to a second order method[27] by x n + = x n + ∆ t · v n + v n − v n − ∆ t · ( ∆ t ) , (20)with previous time step size ∆ t = t n − t n − . Step 2. Temperature
A semi-implicit time integration is then carried out to compute the new temperature T n + by ( I T + D T ) · T n + = ( ρ n · c n v ) · T n + f T , (21)with I T = ρ n · c n v · I , (22) D T = − ∆ t · ˜ ∇ T ( λ · ˜ ∇ ) , f T = ∆ t · ( ˜ ∇ T ( S n · v n ) − ( ˜ ∇ T S n ) · v n − p n · ˜ ∇ T v n ) , where the overhead ∼ indicates the discrete differential operators.To simplify notation, the index of the points has been omitted. Equation (21) forms a sparse linearsystem of equations with unknowns T n + at each point of the point cloud. All sparse implicit linearsystems arising in this and the coming steps are solved with a BiCGSTAB solver, without the use of apreconditioner. Step 3. Concentrations
A similar procedure as that done for the temperature is carried out for the concentrations. We use asemi-implicit time integration for the concentration of each species c n + i , i = , . . . , N , ( I c i + D c i ) · c n + i = c ni , (23)with I c i = ( I + ∆ t · ˜ ∇ T v n ) , (24) D c i = − ∆ t · ˜ ∇ T ( D i , eff · ˜ ∇ ) . Step 4. ρ , η , and c v The updated density ρ n + , viscosity η n + , as well as heat capacity c n + are then determinedaccording to the definitions in Sect. 2.2.Using the updated solution viscosity, a preliminary viscosity for the momentum equation iscomputed as ˆ η n + = η n + + η n turb . Step 5. Hydrostatic pressure
The pressure is split into its hydrostatic (body forces) and dynamic parts (movement of the fluid) as p = p hyd + p dyn . (25)First the updated hydrostatic pressure p n + is computed˜ ∇ T (cid:18) ρ n + · ˜ ∇ p n + (cid:19) = ˜ ∇ T g . (26)Using the updated hydrostatic pressure, a pressure guess ˆ p is computed which will be used whilecomputing the new velocity ˆ p = p n + + p n dyn . (27) Step 6. Coupled velocity-pressure
Time integration of the first equation in (19) provides the targeted divergence of velocity ˜ ∇ T v n + .To solve for v n + and p n + in an implicit time integration scheme, we use the penalty formulationintroduced in [10, 13]. Using the pressure guess defined above, we obtain the following coupledvelocity-pressure-system for preliminary velocity ˆ v n + and correction pressure p n + : (cid:18) I − ∆ t ρ n + · ˜ ψ n + η n + (cid:19) · ˆ v n + + ∆ t ρ n + · ˜ ∇ p n + (28) = v n − ∆ t ρ n + · ˜ ∇ ˆ p + ∆ t · g , ˜ ∇ T (cid:18) ∆ t virt ρ n + · ˜ ∇ p n + (cid:19) = ˜ ∇ T ˆ v n + − ˜ ∇ T v n + , with ( ˜ ψ n + η n + ) T = ˜ ∇ T ( ˆ η n + · ˜ ∇ )( ˆ v n + ) T (29) + ( ˜ ∇ ˆ η n + ) T · ( ˜ ∇ ( ˆ v n + ) T ) T + ˆ η n + · ( ˜ ∇ ( ˜ ∇ T ˆ v n + )) T − · ( ˜ ∇ T ˆ v n + ) · ( ˜ ∇ ˆ η n + ) T , and ∆ t virt = A virt · ∆ t , 0 ≤ A virt ≤
1. If A virt =
1, the scheme corresponds to an implicit Chorin projection,see [3]. Theoretically, choosing A virt = . ≤ A virt ≤ .
1, conditioning of the linearsystem is sufficiently good. Furthermore, the resulting preliminary velocity features a divergence whichis very close to the targeted one. We note that in equations (28) and (29), the stress tensor S n + wasdetermined according to equation (2). 8 Step 7. Update velocity and pressure
The updates of velocity and dynamic pressure are given by v n + = ˆ v n + − ∆ t virt ρ n + · ˜ ∇ p n + , (30) p n + = p n dyn + p n + . Step 8. Turbulence
For the k - ε turbulence model, we derive a singularity formulation from equation (3): ddt (cid:18) k ε (cid:19) =( C ε − ) + C η · ( − C ε ) · (cid:107) ˜ ∇ v T (cid:107) · (cid:18) k ε (cid:19) + C η · ( C ε · C ε − ) ρ · Pr turb · ∂ ρ∂ T · ( g · ˜ ∇ T ) · (cid:18) k ε (cid:19) + ρ · ˜ ∆ η ∗ (cid:18) k ε (cid:19) , (31) ddt (cid:16) ε k (cid:17) =( − C ε ) · (cid:16) ε k (cid:17) + C η · ( C ε − ) · (cid:107) ˜ ∇ v T (cid:107) + C η · ( − C ε · C ε ) ρ · Pr turb · ∂ ρ∂ T · ( g · ˜ ∇ T )+ ρ · ˜ ∆ η ∗ (cid:16) ε k (cid:17) , where ˜ ∆ η ∗ (cid:18) k ε (cid:19) = ε · ˜ ∆ η k k − k · ˜ ∆ η ε εε , (32)˜ ∆ η ∗ (cid:16) ε k (cid:17) = k · ˜ ∆ η ε ε − ε · ˜ ∆ η k kk with ˜ ∆ η k = ˜ ∇ T (cid:18)(cid:18) η + η turb σ k (cid:19) · ˜ ∇ (cid:19) , (33)˜ ∆ η ε = ˜ ∇ T (cid:18)(cid:18) η + η turb σ ε (cid:19) · ˜ ∇ (cid:19) . If k , ε > t n ≤ t ≤ t n + , numerical mean values can be determined from (31): k ε (cid:12)(cid:12)(cid:12)(cid:12) m = ∆ t (cid:90) t n + t n ddt (cid:18) k ε (cid:19) dt , (34) ε k (cid:12)(cid:12)(cid:12) m = ∆ t (cid:90) t n + t n ddt (cid:16) ε k (cid:17) dt . We use the mean values to avoid singularities in the discretized k - ε turbulence model. dkdt = ˜ ∆ η k k ρ − ε k (cid:12)(cid:12)(cid:12) m · k + C η · P prb , k · k ε (cid:12)(cid:12)(cid:12)(cid:12) m · k , (35) d ε dt = ˜ ∆ η ε ερ − C ε · ε k (cid:12)(cid:12)(cid:12) m · ε + C ε · C η · P prb , ε · k ε (cid:12)(cid:12)(cid:12)(cid:12) m · ε , P prb , k = (cid:107) ˜ ∇ v T (cid:107) − ρ · Pr turb · ∂ ρ∂ T · ( g · ˜ ∇ T ) , (36) P prb , ε = (cid:107) ˜ ∇ v T (cid:107) − C ε ρ · Pr turb · ∂ ρ∂ T · ( g · ˜ ∇ T ) . A fully implicit time integration scheme for the turbulent kinetic energy k n + can now be developedas k n + − ∆ t · ˜ ∆ η k k n + ρ + ∆ t · ε k (cid:12)(cid:12)(cid:12) m · k n + (37) − ∆ t · C η · P n + , k · k ε (cid:12)(cid:12)(cid:12)(cid:12) m · k n + = k n . A similar procedure is used to compute the updated turbulent dissipation ε n + − ∆ t · ˜ ∆ η ε ε n + ρ + ∆ t · C ε · ε k (cid:12)(cid:12)(cid:12) m · ε n + (38) − ∆ t · C ε · C η · P n + , ε · k ε (cid:12)(cid:12)(cid:12)(cid:12) m · ε n + = ε n . The mean values are determined analytically. This is illustrated in detail for k ε (cid:12)(cid:12) m . Assuming that thediffusion term ρ · ˜ ∆ η ∗ (cid:0) k ε (cid:1) is negligible as well as defining x = k ε , a = C ε − , b = C η · ( C ε − ) · (cid:107) ˜ ∇ v T (cid:107) + C η · ( − C ε · C ε ) ρ · Pr turb · ∂ ρ∂ T · ( g · ˜ ∇ T ) , we can rewrite equation (31) as dxdt = a − b · x . (39)For x = (cid:112) ab , we obtain x n + (40) = x · tanh (cid:16) ∆ t · √ a · b + arctanh (cid:16) x n x (cid:17)(cid:17) , x n < x x , x n = x x · coth (cid:16) ∆ t · √ a · b + arccoth (cid:16) x n x (cid:17)(cid:17) , x n > x . Finally, the updated turbulent viscosity is determined by η n + = ρ n + · C η · ( k n + ) ε n + . (41) In case of the Eulerian formulation, [25] shows that a second order timeintegration scheme should be applied to numerically solve transport terms of the form v T ∇ in theGFDM context. For this purpose, the SDIRK2 method is proposed (see [1]), which features the samestability properties as an implicit Euler time integration scheme. Furthermore, an upwind discretizationby means of a MUSCL reconstruction with a Superbee limiter is used.101The majority of the steps are the same as those carried out in the Lagrangian formulation. Themovement step of the Lagrangian formulation is skipped here. And the coupled velocity-pressuresystem is modified to the following two-step procedure: (cid:18) I ˆ v n + α − α · ∆ t ρ n + α ˜ ψ n + α ˆ η n + α (cid:19) · ˆ v n + α + α · ∆ t ρ n + α · ˜ ∇ p n + α corr (42) = v n − α · ∆ t ρ n + α · ˜ ∇ ˆ p + α · ∆ t · g , ˜ ∇ T (cid:18) ∆ t virt ρ n + α · ˜ ∇ p n + α corr (cid:19) = ˜ ∇ T ˆ v n + α − ˜ ∇ T v n + α , with I ˆ v n + α =( I + α · ∆ t · ( (cid:103) v T ∇ ) ˆ v n + α ) , (43) ( ˜ ψ n + α ˆ η n + α ) T = ˜ ∇ T ( ˆ η n + α · ˜ ∇ )( ˆ v n + α ) T + ( ˜ ∇ ˆ η n + α ) T · ( ˜ ∇ ( ˆ v n + α ) T ) T + ˆ η n + α · ( ˜ ∇ ( ˜ ∇ T ˆ v n + α )) T − · ( ˜ ∇ T ˆ v n + α ) · ( ˜ ∇ ˆ η n + α ) T , and α = − √ . Density and viscosity for the intermediate step can for instance be determined bylinear interpolation between time levels n and n + v n + − ∆ t · α · V ( ˆ v n + , p n + ) (44) = v n + ∆ t · ( − α ) · V ( ˆ v n + α , p n + α corr ) , ˜ ∇ T (cid:18) ∆ t virt ρ n + · ˜ ∇ p n + (cid:19) = ˜ ∇ T ˆ v n + − ˜ ∇ T v n + with V ( ˆ v n + , p n + ) = − ρ n + · ( (cid:103) v T ∇ ) ˆ v n + + ρ n + · ˜ ψ n + η n + (45) − ρ n + · ˜ ∇ ˆ p n + − ρ n + · ˜ ∇ p n + + g , V ( ˆ v n + α , p n + α corr ) = ˆ v n + α − v n α · ∆ t . For more details on the Eulerian procedure, we refer to [25], and for similarGFDM Eulerian formulations, we refer to [29].For numerical validations of the velocity-pressure scheme used here, their implementations within aGFDM framework, and a comparison of GFDM results with other numerical methods on benchmarkproblems, we refer to our earlier work [4, 10, 13, 21, 24, 25, 30].4. MICROSCOPIC SCALETo study the smaller scale (both spatially and temporally) dissolution of the salt species in the water, weconsider representative geometries of the salt cavern in a so-called microscopic setup. In this section,we identify effective parameters of the dissolution process. Specifically, we compute the effectivediffusion coefficient and the effective transition coefficient between water and surrounding species.These will be used later, in Sect. 5, in the macroscopic procedure to simulate the overall evolution112of the salt cavern. The Lagrangian formulation is used here. The time integration of the underlyingequations is done as presented in Sect. 3.3.1.For the sake of brevity, we restrict the following description to sodium chloride as the species ofinterest. The same procedure can directly be transferred to any other species.
We consider a cylinder with diameter of 5m and height of 10m which is initially filled with pure water,i.e. c NaCl ( t = ) =
0. During the simulation, the temperature is fixed to T = ◦ C.The roof of the cylinder acts as an inexhaustible supply of sodium chloride which is modeled byapplying a Dirichlet condition with saturation concentration c sNaCl = c sNaCl ( T ) =
357 kgm . (46)For the hull of the cylinder, a homogeneous Neumann condition is applied. Aiming at a quasi-steadystate, the bottom of the cylinder models an outflow boundary. In the interior, we solve dc NaCl dt + c NaCl · ∇ T v = ∇ T ( D micro · ∇ c NaCl ) , (47)where D micro = D laminar + D turb . The laminar diffusion coefficient for sodium chloride is given by D laminar = . · − s (see [6]). For the turbulent part, we have D turb = C η · k ε . (48)Standard boundary conditions (Dirichlet and Neumann) are prescribed for velocity, pressure, and theturbulent quantities. The simulation runs until a quasi-steady state is reached, which will be explainedbelow. In order to determine the effective quantities, the cylinder is split in the axial direction (z-direction) intoequal sub-cylinders SC j , j = , . . . , J . These are used to estimate the mass flow. The planes between thesub-cylinders are denoted by help-planes HP j , j = , . . . , J − v = − ) for 5 consecutive time steps. The mass flow of sodium chloride is given by dmdt = − D NaCl , eff · ∂ ¯ c NaCl ∂ n , (49)where ¯ c NaCl is the mean concentration. The mass flow and the mean concentration in sub-cylinder SC j are determined by dmdt ( SC j ) = (cid:82) SC j c NaCl · v dV SC j (cid:82) SC j dV SC j , (50)¯ c NaCl ( SC j ) = (cid:82) SC j c NaCl dV SC j (cid:82) SC j dV SC j . (a) t =
10s (b) t =
55s (c) t = Figure 2. Evolution of concentration in the microscopic simulation for interaction radius h = .
18m (Lagrangianformulation).
Based on the mean concentration in a sub-cylinder SC j , we can approximate its normal derivativewith respect to the help plane HP j . This yields the effective diffusion coefficients in each sub-cylinder D NaCl , eff ( SC j | HP j ) = − dmdt ( SC j ) ∂ ¯ c NaCl ∂ n (cid:12)(cid:12)(cid:12) HP j , j = , . . . , J − . (51)Once a quasi-steady state is reached, an overall effective diffusion coefficient can be determined. Toaccommodate the “quasi-steady” character of the simulation, we use a time-averaged effective diffusioncoefficient, over a small time interval, and over each of the sub-cylinders. This value will later be usedin the macroscopic setup. The effective transition coefficient γ NaCl , eff is derived in a mannersimilar to that done for the effective diffusion coefficient D NaCl , eff above. γ NaCl , eff ( SC j ) = − dmdt ( SC j ) c sNaCl − ¯ c NaCl ( SC j ) , j = , . . . , J − . (52)Once again, the time-averaged values of the effective transition coefficient in each of the sub-cylindersat the quasi-steady state gives the overall effective transition coefficient which will be used in themacroscopic simulations in Sect. 5. With the help of γ NaCl , eff , we can define the solution rate of sodiumchloride for given temperature T by R NaCl ( T ) = γ NaCl , eff ( c sNaCl − c NaCl ) . (53) In the simulations carried out, we choose J =
10 to divide the cylinder domain considered into 10sub-cylinders of height 1m each. We consider several levels of resolution to study the convergence ofthe effective parameters being determined to resolution-independent values. The coarsest resolutionused is h = .
8m corresponding to 70400 points in the domain. h is consecutively halved till h = . t = h = .
18m is illustrated in Fig. 2 in the time interval [ , ] .As expected, the flow is characterized by viscous fingering.134 Figure 3. Convergence of the effective diffusion coefficient in the microscopic simulations. NP denotes the numberof the points in the initial domain. -5 -4 -3 Figure 4. Convergence of the effective transition coefficient in the microscopic simulations. NP denotes the numberof the points in the initial domain. The convergence of effective diffusion as well as transition coefficient with decreasing h is shownin Fig. 3 and Fig. 4, respectively. The values plotted are also tabulated in Table II, along with therelation between the interaction radius h and the number of points in the domain NP . The time stepsize is governed by ∆ t = CFL
Lag · h | v | , with CFL
Lag set to 0 .
2. The diffusion coefficient converges to D NaCl , eff = . m s , while the transition coefficient converges to γ NaCl , eff = . ms . Using equation(53), we obtain a maximum solution rate of R NaCl , max ( ◦ C ) = . kgm · s . Compared to the solutionrate of 0 . kgm · s for T = ◦ C determined in [11] at a crystal level, the estimated solution rate is ofthe correct order of magnitude. 5. MACROSCOPIC SCALEWe now model the overall evolution of the salt cavern during the double-well solution mining process.Both the Lagrangian as well as the Eulerian formulation are evaluated for this.The model equations and time integration procedures for the macroscopic scale simulations arethe same as those described in Sect. 2 and Sect. 3.3, respectively, with a few variations. Firstly, the145
Table II. Estimated effective diffusion D [ m s ], transition coefficient γ [ ms ], and maximum solution rate R [ kgm · s ],with varying interaction radius h [m] and corresponding initial number of points NP . h NP D NaCl , eff γ NaCl , eff R NaCl , max ( ◦ C ) .
80 70 400 0 . . × − . .
40 443 318 0 . . × − . .
30 974 918 0 . . × − . .
20 3 076 350 0 . . × − . .
18 4 139 040 0 . . × − . .
15 7 014 494 0 . . × − . .
10 20 892 871 0 . . × − . D NaCl , eff · ∂ c NaCl ∂ n = γ NaCl , eff · ( c sNaCl − c NaCl ) . (54)Here, the effective diffusion coefficient D NaCl , eff , as well as the effective transition coefficient γ NaCl , eff are the values determined in the microscopic simulation in Sect. 4, D NaCl , eff = . m s and γ NaCl , eff = . ms .A further difference in the time integration procedure comes in Steps 2 and 4 described in Sect. 3.3.Here, we fix the temperature to T = ◦ C and, subsequently, obtain the corresponding saturationconcentration c sNaCl = kgm . For simplicity, the following linearized relations for density and viscosityof the solution are used (see [25]) ρ ( c NaCl ) ≈ ( . · c NaCl + ) kgm , (55) η ( c NaCl ) ≈ ( . · − · c NaCl + − ) Pas . We are interested in the geometrical evolution of the double-well salt cavern. The initial geometry isgiven by a small cavern filled with pure water that is surrounded by sodium chloride, see Fig. 5. Thedimensions of the initial cavern are approximately: width of 90m, height of 50m, and depth of 26m.The sodium chloride deposit is limited to impermeable surrounding rock. The pipe on the left side actsas an inlet of fresh water with inflow velocity | v in | = ms , whereas the pipe on the right side acts as theoutlet.In reality, the maximum diameter of the pipes is of the order of 1m. Hence, the resolution of thepoint cloud close to the inlet and the outlet has to be of the order of 0 .
1m to ensure accurate results incase of the Lagrangian formulation. This would lead to an extremely small time step size compared tothe desired simulation time of several years/decades due to the CFL-condition ∆ t Lag ≤ CFL
Lag · h min | v | . (56)Numerically, we observe that stable results are achieved for CFL
Lag = .
15, which results in ∆ t Lag = O ( . ) . Performing long-term simulations over months of simulation time is not feasiblewith such a small time step. This results in the need for using the Eulerian formulation for the problemat hand.In order to allow for a comparison of Lagrangian and Eulerian formulation, we consider pipes ofdiameter 12m. This also decreases the required actual time being simulated. Numerically, we observe156 Figure 5. Macroscopic simulation setup – initial geometry, see [25]. that with the significantly larger diameter used here, the evolution of the cavern only requires a fewhours of physical time to be simulated, compared to the few months or years with the actual diameter.We note that despite this time reduction, this still corresponds to a time scale two orders of magnitudegreater than that used in the microscopic simulations in Sect. 4.The simulations are performed with a constant interaction radius of h = The movement of the boundary of the cavern can be defined by the Stefan condition ρ v (cid:63) = γ NaCl , eff ( c sNaCl − c NaCl ) , (57)see [9]. This yields v (cid:63) = γ NaCl , eff ρ · ( c sNaCl − c NaCl ) (58)and, consequently, a movement of the boundary in normal direction n with velocity v boundary = v (cid:63) · n .To speed up computation, a time lapse procedure can be applied [25]. Due to small flow velocitiesinside the salt cavern, an additional speed-up factor A can be introduced in the definition of v (cid:63) by v (cid:63) A = A · γ NaCl , eff ρ · ( c sNaCl − c NaCl ) . (59)For stability reasons, we require ∆ t · v (cid:63) A ≤ . · h . (60)The maximum movement velocity of the boundary occurs in case of pure water, i.e. c NaCl =
0, andis given by v (cid:63) A , max = A · γ NaCl , eff · c sNaCl ρ . (61)Given a maximum time step size ∆ t max , equation (60) leads to the constraint A ≤ . · h · ρ ∆ t max · γ NaCl , eff · c sNaCl . (62)167 Figure 6. Translational velocity in the macroscopic simulation, see [25] (Eulerian formulation including ALE atthe moving boundary).
Due to the movement of the boundary, interior points close to this boundary have to move in theEulerian formulation also. For this purpose, the ALE-approach presented in [8] is used. Based oncurrent and future position of an affected interior point, the translational velocity v trans = x n + − x n ∆ t (63)is determined, see Fig. 6. Due to the explicit movement of these points, the convection terms in thenumerical model in Eulerian form in Sect. 3.3 must refer to the relative velocity v − v trans instead of v .Furthermore, this introduces a CFL-condition of the form ∆ t ALE ≤ CFL
ALE · h min v (cid:63) . (64)This depends on the boundary velocity v (cid:63) = O ( . ms ) which is considerably smaller than the flowvelocity (the inflow velocity is 1 ms , while the maximum velocity in the domain is even bigger). h min issubject to the desired resolution at the moving boundary. At the inlet and the outlet, a coarse resolutionis sufficient in this case. Starting from the initial domain as shown in Fig. 5, the salt cavern expands till the outer domain isfilled. Physically, this outer domain can represent either a rock formation where the salt cavern ends,or prescribed limits of the region where the solution mining is to be carried out. Fig. 7 illustrates theevolution of the salt cavern in the Eulerian formulation according to C = c NaCl . An animation of thisprocess can be found in the Online Resource. This expansion is quantified by plotting the volume as afunction of time in Fig. 8. We note that once the entire outer domain is filled, the volume of the domainabout 27 . h =
4m which corresponds to NP =
90 744 points at the initial state. Thesimulations are run until t = = . (a) t = t = t = Figure 7. Evolution of the macroscopic simulation for interaction radius h =
4m – concentration, see [25] (Eulerianformulation including ALE at the moving boundary). Figure 8. Expansion in the volume of the computational domain of the macroscopic simulation as the simulationprogresses. boundary Q c ( t ) = (cid:90) t (cid:18) (cid:90) ∂ Ω out c NaCl v · n d A (cid:19) d τ , (65)where ∂ Ω out is the outflow boundary located at the top of the extraction well. Physically, this representsa measure of the concentration of salt being extracted. The time evolution of Q c is shown in Fig. 9. Itillustrates that both formulations produce very similar results.To emphasize the need of the ALE formulation for such a simulation, we compare the time stepsrequired in both the ALE and Lagrangian formulations for stability. Considering the simulation time of7200s, we observe that the Lagrangian formulation required at least 22915 time steps to obtain stableresults, which corresponds to an average time step size of ∆ t ≈ . ∆ t ≈ . Figure 9. Comparison of the concentration flux across the outflow pipe between the Eulerian (with movingboundaries) and Lagrangian simulations. on a macroscopic, as well as a microscopic scale. Both Lagrangian and Eulerian approaches wereconsidered.On the macroscopic scale, we considered the expansion of the salt cavern as a result of erosionoccurring as the salt dissolves in the water. In reality, this procedure occurs over the time span of severalmonths or years. A simplified geometry was considered here, which enabled a comparison betweenthe Eulerian and Lagrangian formulations. In this simplified macroscopic set-up, the expansion ofthe salt cavern occurred over the time scale of several hours. Since the dissolution of salt in wateroccurs on a much smaller time level we also considered a microscopic set-up over a duration of a fewminutes. This was used to determine effective parameters governing the dissolution process. Using theexample of sodium chloride as the species of interest, effective diffusion and transition coefficientswere determined in the microscopic simulations. These values were then used in the macroscopicsimulations to determine the evolution of the concentration inside the salt cavern and to specify thesolution rate of the salt species at the boundary, i.e. to model the geometrical evolution of the saltcavern.A comparison of the numerical results of the Lagrangian and Eulerian formulations (extended by anALE-approach) in the macroscopic case illustrates the advantages of the latter one due the possibilityof using much larger time step sizes. Aiming at a simulation time of several years, the forecastcomputation time for a simulation of a double-well solution mining process based on the Lagrangianformulation would be of the order of years. In contrast to that, the flexibility of the Eulerian formulationregarding the resolution of the point cloud (local refinement only at the moving boundary) enablesmeshfree simulations in reasonable time – especially in terms of real applications.CONFLICT OF INTERESTOn behalf of all authors, the corresponding author states that there is no conflict of interest.
REFERENCES1. R. Alexander. Diagonally implicit runge–kutta methods for stiff odes.
SIAM Journal on Numerical Analysis , 14(6):1006–1021, 1977.2. J. Chen, D. Lu, W. Liu, J. Fan, D. Jiang, L. Yi, and Y. Kang. Stability study and optimization design of small-spacingtwo-well (sstw) salt caverns for natural gas storages.
Journal of Energy Storage , 27:101131, 2020.3. A. J. Chorin. Numerical solution of the navier-stokes equations.
Mathematics of computation , 22(104):745–762, 1968.4. C. Drumm, S. Tiwari, J. Kuhnert, and H.-J. Bart. Finite pointset method for simulation of the liquid - liquid flow field inan extractor.
Computers & Chemical Engineering , 32(12):2946 – 2957, 2008.
5. C.-M. Fan, C.-N. Chu, B. ˇSarler, and T.-H. Li. Numerical solutions of waves-current interactions by generalized finitedifference method.
Engineering Analysis with Boundary Elements , 2018.6. M. Flury and T. Gimmi. Solute diffusion. In J. Dane and G. Topp, editors,
Methods of soil analysis. Part. 4. Physicalmethods. , pages 1323–1351, Madison, WI., 2002. SSSA Book Ser. 5.7. L. Gavete, F. Ure˜na, J. Benito, A. Garc´ıa, M. Ure˜na, and E. Salete. Solving second order non-linear elliptic partialdifferential equations using generalized finite difference method.
Journal of Computational and Applied Mathematics ,318:378 – 387, 2017. Computational and Mathematical Methods in Science and Engineering CMMSE-2015.8. C. W. Hirt, A. A. Amsden, and J. Cook. An arbitrary lagrangian-eulerian computing method for all flow speeds.
Journalof computational physics , 14(3):227–253, 1974.9. E. Javierre-Perez.
Literature Study: Numerical methods for solving Stefan problems . Delft University of Technology,2003.10. A. Jefferies, J. Kuhnert, L. Aschenbrenner, and U. Giffhorn. Finite pointset method for the simulation of a vehicle travellingthrough a body of water. In M. Griebel and A. M. Schweitzer, editors,
Meshfree Methods for Partial Differential EquationsVII , pages 205–221, Cham, 2015. Springer International Publishing.11. O. Karsten. L¨osungsgeschwindigkeit von natriumchlorid, kaliumchlorid und kieserit in wasser und in w¨asserigen l¨osungen.
Zeitschrift f¨ur anorganische und allgemeine Chemie , 276(5-6):247–266, 1954.12. A. Katz and A. Jameson. Meshless scheme based on alignment constraints.
AIAA journal , 48(11):2501–2511, 2010.13. J. Kuhnert. Meshfree numerical scheme for time dependent problems in fluid and continuum mechanics. In S. Sundar,editor,
Advances in PDE Modeling and Computation , pages 119–136, New Delhi, 2014. Anne Books.14. J. Kuhnert, I. Michel, and R. Mack. Fluid structure interaction (fsi) in the meshfree finite pointset method (fpm): Theoryand applications. In M. Griebel and A. M. Schweitzer, editors,
Meshfree Methods for Partial Differential Equations IX,IWMMPDE2017 , pages 73–92. Springer, 2019.15. M. Lalibert´e. Model for calculating the viscosity of aqueous solutions.
Journal of Chemical & Engineering Data ,52(2):321–335, 2007.16. M. Lalibert´e. A model for calculating the heat capacity of aqueous solutions, with updated density and viscosity data.
Journal of Chemical & Engineering Data , 54(6):1725–1760, 2009.17. M. Lalibert´e and W. E. Cooper. Model for calculating the density of aqueous electrolyte solutions.
Journal of Chemical& Engineering Data , 49(5):1141–1151, 2004.18. B. Launder and D. Spalding. The numerical computation of turbulent flows.
Computer Methods in Applied Mechanicsand Engineering , 3(2):269 – 289, 1974.19. R. L¨ohner and E. O˜nate. An advancing front point generation technique.
Communications in Numerical Methods inEngineering , 14(12):1097–1108, 1998.20. M. Luo, C. G. Koh, W. Bai, and M. Gao. A particle method for two-phase flows with compressible air pocket.
InternationalJournal for Numerical Methods in Engineering , 108:695–721, Nov. 2016.21. I. Michel, S. M. I. Bathaeian, J. Kuhnert, D. Kolymbas, C.-H. Chen, I. Polymerou, C. Vrettos, and A. Becker.Meshfree generalized finite difference methods in soil mechanics—part ii: numerical results.
International Journal onGeomathematics , 8(2):191–217, Nov 2017.22. F. Nick, H.-J. Plum, and J. Kuhnert. Parallel detection of subsystems in linear systems arising in the meshfree finitepointset method. In M. Griebel and M. A. Schweitzer, editors,
Meshfree Methods for Partial Differential Equations IX ,pages 93–115, Cham, 2019. Springer International Publishing.23. F. R. Saucedo-Zendejo, E. O. Res´endiz-Flores, and J. Kuhnert. Three-dimensional flow prediction in mould fillingprocesses using a gfdm.
Computational Particle Mechanics , 6(3):411–425, 2019.24. B. Seibold.
M-Matrices in Meshless Finite Difference Methods . PhD thesis, Kaiserslautern University, 2006.25. T. Seifarth.
Numerische Algortihmen f¨ur gitterfreie Methoden zur L¨osung von Transportproblemen . PhD thesis, Universityof Kassel, Kassel, 2017.26. P. Suchde.
Conservation and Accuracy in Meshfree Generalized Finite Difference Methods . PhD thesis, University ofKaiserslautern, Kaiserslautern, Germany, 2018.27. P. Suchde and J. Kuhnert. Point cloud movement for fully lagrangian meshfree methods.
Journal of Computational andApplied Mathematics , 340:89 – 100, 2018.28. P. Suchde and J. Kuhnert. A fully lagrangian meshfree framework for PDEs on evolving surfaces.
Journal ofComputational Physics , 395:38 – 59, 2019.29. P. Suchde and J. Kuhnert. A meshfree generalized finite difference method for surface PDEs.
Computers & Mathematicswith Applications , 78(8):2789 – 2805, 2019.30. P. Suchde, J. Kuhnert, S. Schr¨oder, and A. Klar. A flux conserving meshfree method for conservation laws.
InternationalJournal for Numerical Methods in Engineering , 112(3):238–256, 2017.31. E. Uhlmann, E. Barth, T. Seifarth, M. H¨ochel, J. Kuhnert, and A. Eisentr¨ager. Simulation of metal cutting with cuttingfluid using the finite-pointset-method, 2020. 9th CIRP Conference on High Performance Cutting. Submitted to ProcediaCIRP.32. E. Uhlmann, R. Gerstenberger, and J. Kuhnert. Cutting simulation with the meshfree finite pointset method.
ProcediaCIRP , 8:391 – 396, 2013.33. G. Zhang, Z. Wang, K. Zhang, Y. Li, Y. Wu, Y. Chen, and H. Zhang. Collapse mechanism of the overlying strata above asalt cavern by solution mining with double-well convection.
Environmental Earth Sciences , 77(16):588, 2018., 77(16):588, 2018.