Modeling, with a unified level-set representation, of the expansion of a hollow in the ground under different physical phenomena

Modeling, with a unified level-set representation, of the expansion of a hollow in the ground under different physical phenomena Régis Cottereau ·Pedro Díez ·Antonio Huerta Abstract This paper builds on the flexibility of the level-set representation to model in a unified manner the expansion of a hollow in the ground under different physical phenomena. In particular, the dissolving action of a flow of water in a saturated soil, and that of a jet of particles of water in a nonsaturated one, are represented in a common framework. In that manner, the complex geometrical evolutions of the hollow can be followed without the need for remeshing and this approach allows for a smooth transition between saturated and non-saturated models of the soil. Implementation and numerical difficulties are discussed and two applications of industrial interest are considered. The first one describes the modeling of the piping phenomenon, and the second one the evolution of an excavation created by a leaking duct. Keywords Level-set method ·Porous media ·Piping · Internal erosion The support of Gas Natural Distribución SDG S.A. and Ministerio de Educación y Ciencia de España, through grants DPI2007-62395 and BIA2007-66965, is gratefully acknowledged. R. Cottereau International Center for Numerical Methods in Engineering (CIMNE), E.T.S. d’Enginyers de Camins, Canals i Ports de Barcelona, Universitat Politècnica de Catalunya, Barcelona, Spain R. Cottereau Laboratoire MSSMat, CNRS UMR 8579, École Centrale Paris, Châtenay-Malabry, France e-mail: [email protected] P. Díez ·A. Huerta (B ) Laboratori de Càlcul Numèric, Escola Tècnica Superior d’Enginyers de Camins, Canals i Ports de Barcelona, Universitat Politècnica de Catalunya, Barcelona, Spain e-mail: [email protected] 1 Introduction The requirement to capture the position and evolution of an interface is central to many areas of engineering and science, including metal forging, oceanography, imaging, flame modeling, melting of materials, and more generally, the modeling of heterogeneous or multi-phases materials. Along the years, several methods have been devised for such problems [1]. Among the most used ones are front tracking methods, in which the mesh is refined or deformed to follow the displacement of the interface [2], and the marker-in-cell method [3,4], in which a large set of markers follows the material in a lagrangian way, hence describing the position of each phase. However, both methods induce high computational costs. Indeed, to keep the numerical accuracy to an appropriate level, it is necessary, in the former class of methods, to re-mesh the domain when the deformations increase. Likewise, the number of markers necessary to follow appropriatelyaninterfaceisshowntobeprohibitiveforcomputational implementation [5]. Level-set functions [6,7] provide a very efficient and elegant alternative to these techniques. In the simplest setting, they allow to discriminate between two areas of a domain, with no explicit parameterization of the actual interface. Level-sets are functions defined on the entire domain, whose sign indicates the belonging to one or the other of the two areas. Usually, their absolute value represents the smallest distance to the interface, which is hence indicated by the cancellation of the level-set. Conceptually, they are constructed in a space of higher dimension than the interface they intend to represent, with an improved mathematical behavior that allows for an easier manipulation. Since the first use of level-set functions in the description of dynamical two-phase fluid systems [8,9], their power has been acknowledged for the parameterization of 1 complex evolving phases. In particular, their ability to deal with changes in topology without any remeshing has been recognized [10,11]. They have been used in several fields of geophysics and geomechanics, including modeling of two-phase flows and permeability estimation in reservoir simulations [11,12], tectonic plates subduction [13,14], seismic waves travel time computation [15], and, generally, for inverse problems and optimal design [16–18]. Apart from this potential in considering complex geometries, we want here to advocate the capabilities of the level set representation in integrating different physical models in a common framework. To our knowledge, this aspect has not been emphasized before. In particular, we consider a leaking duct, generating a hollow by dissolution and tearing of the soil. Depending on the height of the groundwater table, two very different physical phenomena occur. When the groundwater table is low, the hollow is full of air, and its expansion is mainly due to the impact of the jet of water coming out of the leak. When the groundwater table is high, the hollow is full of water, and the main phenomenon controlling the evolution of the excavation is the dissolution and transport of soil particles with the water flow. In this paper, both these phenomena will be considered in the same framework of the Hamilton–Jacobi equation. There, the choice of the model of erosion (jet or flow) is entirely contained in that of the velocity of advancement of the interface. Conceptually and practically, this change of focus with respect to the classical approach allows to shift very easily from one model of the advancement of the interface to the other. In real life (e.g., for the rise and fall of the water table in the soil), this may be a very common situation. Note that this change of focus means that the questions of whether the models of normal velocity are actually useful for real-life applications, or whether the numerical schemes used to compute them are accurate, are of secondary importance for this paper. We will concentrate on the integration of both models in a common framework, and on the numerical accuracy of the level set representation itself. The latter will be described on an example representing the piping phenomenon, with only one model involved. Apart from its interest in showing the numerical accuracy of the level set representation, this example is also, to the knowledge of the authors, the first application of level sets to this problem, for which it seems appealing. Indeed, piping is one the main causes of failure of dams and embankments, while its modeling is still based on semi-analytical [19] or experimental [20,21] approaches. In particular, the semi-analytical approaches are heavily dependentonsymmetryhypotheses that are probably not very relevant in practice. Using a level-set parameterization for the interface between the soil and the hollow allows to follow any geometry of the piping hole. The paper is constructed as follows. In a first part (Sect. 2), we recall the main definitions and results for level set descriptions, and introduce two examples of models for the normal velocity of the interface (Sects. 2.3,2.4). These two models are the ones that are used in the examples of the last section. We then describe some numerical issues related to the use of level-set descriptions (Sect. 3), in particular concerning the size of the mesh and re-distancing of the level-set function. Finally, we present the two applications of interest: a level-set description of the piping phenomenon (Sect. 4.2), and the expansion of a hollow in the ground under two different physical phenomena (Sect. 4.3). 2 The level-set description The goal of this section is to provide a general description of the level-set approach for the representation of a multi-phase medium. The original ideas were presented in [8]. Although we only discuss here the case where there are two phases, the ideas can be extended to multi-phases problems [22]. 2.1 Basic principles Let us consider a domain separated into exactly two subdomains 1and 2. The interface between these two subdomainsis denotedandevolvesintimewith anormal velocity vn(x,t). Note that the definition of this velocity field implies the choice of an “interior” and “exterior”. Two examples of such situations, with closed and open interfaces, are shown in the left column of Fig. 1. The main idea of the level-set approach is to consider a function, defined on the entire domain of interest, with values that depend on the position in one phase or the other. The most common approach consists in defining the levelset function as the distance function to the interface, with the sign indicating the subdomain. In mathematical terms, it is defined by ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ |(x)|=(x) (x)>0x∈1 (x)<0x∈2 (x)=0x∈ ,(1) where (x)indicates the distance from point xto the closest point of curve =0. Hence, the level-set function is defined in a space with a higher dimension than the interface it attempts to parameterize. However, the additional dimension allows for an increased smoothness of the function, which can hence be more easily manipulated. Further, the level-set description of the interface fits properly in a Finite Element(FE)contextbecause the functionmaybe interpolated using the same mesh. This allows describing arbitrary interfaces, not restricted to contain the mesh nodes, with an accuracy and regularity up to the resolution of the FE 2 2.3 Erosion due to a jet The first physical situation we consider is that of the erosion by a jet of water of an air-filled cavity within the earth. The jet of water is modeled as a set of ejected particles (small volumes of water, drops), launched at different times and from different points of a source. The volumes of the different particles are denoted by Vp,p=1,2,..., such that the particles move through space as point masses mp=ρwVp, where ρwis the density of water. The source corresponds geometrically to a small part of the boundary. Each particle undergoes a free flight, only submitted to gravity, starting at time t0 pfrom a point x0 pof the source and with an initial velocity v0 p. Thus, the position xp(t)and velocity vp(t)of a particle pat time tare given by xp(t)=x0 p+v0 p(t−t0 p)+g1 2(t−t0 p)2 vp(t)=v0 p+g(t−t0 p)(5) where g=−gizis the acceleration of gravity, with g=9.81 m/s2, and izpoints upwards. For a given time interval [t−t/2,t+t/2], the flow is modeled using a given number of particles N, to each of which is assigned the same mass and velocity (norm), so as to fit the macroscopic jet values, in particular the flow, Qin m3/s, and the velocity, in m/s. Both the location x0 pand the direction (angle) of the initial velocity v0 pare selected randomlyusingaproperassumptionintheirprobabilitydistribution (for instance, uniform density for the location, Gaussian distribution for the direction). Tuning the parameters of these probability distributions allows modeling wider or narrower jets, with different dispersions. The velocity of the particles is assumed to be large, in the sense that the flying time of the particles is much lower than the time interval t. Under this assumption, we concentrate the launching of all the particles at the mid time of the interval, namely t0 p=t,for p=1,2,...,N. Given the level-set parametrization (x)=0 of the interface between the earth and the air, we then calculate the point and velocity of impact of each particle with the interface as the first time ti p(after t0 p) when the particle verifies (xp(ti p)) =0, that is ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ti p=mint>t0 parg (xp(t)) =0 xi p=xp(ti p)=x0 p+v0 p(ti p−t0 p)+g(ti p−t0 p)2 2 vi p=vp(ti p)=v0 p+g(ti p−t0 p) .(6) When using a FE approach for the representation of the levelsetfunction,asdescribedinSect.3,findingxi pandvi preduces to finding the intersection of the parabolic trajectory of the particle pwith the polygonal interface described by the FE approximation to (x)(this polygonal interface corresponds to linear elements, and is replaced by a piecewise quadratic curve for elements of degree two). In the general case, this intersection can be located very easily by simply sampling the level set function along the 1D, curved, trajectories of the particles, and detecting the first change of sign. We then choose to model the erosion by taking the normal velocity of the interface as proportional to an averaged normal impacting momentum. Thus, in order to assess this normal velocity, at a given point xof the interface, and at time t,vn(x,t), we consider the particles impacting during a short time interval ta small surface (segment in 2D) dSx with normal n, centered around x. The average is taken as: vn(x,t)=V0 dSt P n·vi pmp,(7) where Pis the set of particles that hit the interface on surface dSxand in the time span [t−t/2,t+t/2], that is P={p|xi p∈dSx,ti p∈[t−t/2,t+t/2]},V0is a parameter that depends on the cohesion of the material, and dS =|dSx|is the measure of the surface dSx. If one additionally supposes that all the particles have the same mass m=ρwQ(t)t/N, where Q(t)is the flow of water coming out of the jet at time t, this equation simplifies to vn(x,t)=ρwV0Q(t) NdS  P n·vi p.(8) This formulation is well adapted to a Monte Carlo resolution because the number Nof particles of water in which the jet has been separated appears explicitely. 2.4 Erosion due to a flow For this second model problem, we consider water flowing through a water-filled cavity and a soil. The flow of water is modeled using Darcy’s law [23], which, together with a mass balance, specifies the evolution of the hydraulic head h(x,t)=z+P(x,t)/ρ wg, with zthe altitude and P(x,t)the hydraulic pressure, ∇·(K(x,t)∇h(x,t))=0,(9) The permeability K(x)describes how easily the water flows through the medium. Note that we consider here a quasistatic problem, so that we neglected a possible influence of ∂h/∂t. However, it would be possible to include dynamic effects through the inclusion of a time derivative term [24]. The two domains of our problem, namely the soil and the cavity, are modeled using two different permeability values K(x∈s,t)=K(x∈s)=Ks K(x∈c,t)=K(x∈c)=1000Ks ,(10) where sand care respectively the domain occupied by the soil and the water-filled cavity, such that s∪c=. Note that we suppose that the permeability coefficients do 4 not explicitly depend on the time variable. However, as the definitions of the domains sand cdo, then the value of the permeability fields in a given point in space will change with time. In practice, this equation is solved using a FE method. Other methods, in particular the eXtended FE (XFE) method combinedwiththelevel-setrepresentationofthe interfacefor the discretization in space, could have been used to increase the accuracy of the approximation (see Sect. 3.3 for more details). We now describe two models of flow erosion, designed respectively for flows parallel and normal to the surface that is being eroded. In the literature, the former type of model appears much more often, due to its importance in the creation of gullies and rills, which is an important issue for agriculture in particular (see [25] for a general review). Models dealing with the erosion of a surface under a normal flow are more scarce, and related to problems of the oil industry [26,27]. Starting with models for perpendicular flows, the most widely used describes the rate of erosion per surface area ˙ε(sometimes denoted Dc) as proportional to the hydraulic shear stress τ[25]: ˙ε=Kc(|τ|−τcr ), (11) where Kcis the coefficient of erosion (with values ranging from 10−1to 10−6sm −1[20]) and τcr is a critical shear stress (with values ranging from 6 to 160 N m−2[21]), both material-dependent. Other models replace the shear stress by thestreampoweras thecontrollingparameter butexperimental evidence seems to indicate that the shear force model is appropriate, both for laboratory and in situ cases [28]. More refined models try to include the influence of the sediments detached from the soil (see for example [19,29]), or consider probabilistic models for the erosion factors [30], however it remains unclear whether the alternatives they provide are truly required with respect to the simpler model of Eq. (11). Following [20,25], an approximation can be used to derive the hydraulic shear stress as a function of the gradient of the hydraulic head, τ=ρwg0∇h(x,t), (12) where 0=1 m in 2D problems. This finally leads to a normal velocity for the cavity-soil interface vn(x,t)=˙ε ρs=Kc ρs(ρwg0|∇h(x,t)|−τcr ),(13) where ρsis the density of the soil material. In the case of an erosion taking place along the axis of the flow, it is not reasonable to use the shear stress as a controlling parameter for erosion. Following [26], and simplifying the model described there because we are only interested in the erosion at the interface between the cavity and the soil, we get a linear relation between the rate of erosion per surface area and the norm of the Darcy flow qn=Ks∇h, which leads to vn(x,t)=˙ε ρs=λ|∇h|,(14) where λis a proportionality coefficient, that would have to be determined experimentally. If we neglect τcr in Eq. (13), then the two equations for the normal velocity have the same form (although derived for different directions), but with different proportionality coefficients a priori. As the real flow going through an interface is probably in between that of a normal incidence and grazing incidence, we will consider in all cases a model for the normal velocity as in Eq. (14), with the additional hypothesis that the proportionality coefficient is the same for all incidences. A more refined model would require the identification of the coefficient as a function of the angle of incidence. 3 Numerical issues In this section, we discuss several numerical points that arise when dealing with level-set methods. In particular, the Sects. 3.1 and 3.2 describe the important issues of, respectively, the influence of the size of the mesh on the representationof the interface,and the requirementofregular renormalization of the distance function. 3.1 Numerical description of the level-set Tothispoint,thelevel-setfunction(x)wasdefinedasageneral field. However, in practice, the level-set function is represented as a FE approximation over a mesh, and described by the nodal values of . This means that the position of the interface =0 is never explicitly stored. Rather, it is recovered, when necessary, by interpolation of the values of (xn)at the nodes xnof the mesh. When using linear finite elements, the interface is therefore approximated inside each element it crosses by a straight line. Hence, strong curvatures of the interface may be smoothed out if a coarse mesh is used. Note that this representation is totally independent of any FE scheme that might be used for the computation of the normal velocity that appears in the Hamilton–Jacobi equation. In particular, when using the model for the erosion by a flow of water (Sect. 2.4), it is possible to use two different meshes: one for the computation of the normal velocity, and one for representing the level-set function. In that case, it is usually interesting to use a mesh much more refined for the representation of the level-set than for the resolution of the flow problem required to compute the front velocity. Due to the explicit character of the method introduced in Sect. 3.2, the additional cost of the level-set transport associated with 5 the mesh refinement is not significant. Moreover, the refined mesh can easily be obtained from a coarse one because there is no requirement of shape for the elements. As they are not used for integration purposes, they can be as deformed as required by the needs of the representation of the curvature. However, in practice, when a FE scheme is used for the computation of the normal velocity, the two meshes are often taken as the same. This is mainly dictated by the simplicity of thenumericalimplementation,and this skips an interpolation step that would be otherwise necessary to pass the information about the position of the interface from one mesh to the other. Note that, when using only one FE mesh, one should pay attention to the fact that the size of the elements of the mesh should be small enough with respect to the smoothness of both the interface and the FE solution that is being computed(thehydraulic head inthecase of the erosion byaflow). Note, also, that this technique was not implemented in the examples of Sect. 4, because the mesh used for the resolution of Darcy’s equation was fine enough for the representation of the curvature of the interface. 3.2 Transporting the level set: solving the Hamilton–Jacobi equation Once vn(xj,t)is obtained from the models discussed in Sects. 2.3 and 2.4, the interface has to be updated solving theHamilton–Jacobiequation(4).Recall that this is a nonlinearhyperbolicadvectionproblem.Asummarizeddescription of the general numerical approaches to tackle this problem using FE may be found in [31], and ad-hoc strategies for the particular case of Hamilton–Jacobi are discussed in detail in [6]. The latter are based on upwind schemes for regular grids. Here, we present a new technique keeping the upwind approach in a general FE context. It is usable with unstructured meshes and allows the use of adaptivity techniques. First, the Hamilton–Jacobi equation (4) is discretized along the time dimension using a first-order forward Finite Difference method. We consider that the level-set function at one point tof the time grid, =(x,t), is known, and we look for the value of at the next point in the grid, t+1. Thus, the Hamilton–Jacobi equation is approximated at the time tin the span [t,t+1]as +1− t+vn(x,t)|∇(x,t)|=0,(15) where t=t+1−t. Hence, selecting the sampling time tequal to t(forward Euler type explicit scheme), the semidiscretized Hamilton–Jacobi equation reads +1=−vn(x,t)t|∇|,(16) and +1can be solved for, given vn(x,t)and . As previously mentioned, the space discretization strategy is adapted to exploit the information contained in the FE description of . We are not using however a standard FE approach, in the sense that the evolution of the nodal values of is computed independently, in a node-by-node basis. This results in an explicit low-cost algorithm. In this paper, we use linear FE interpolation functions, so that |∇|is not defined at the nodes of the mesh, but rather on each element of that mesh. Therefore, the values at the nodes |∇|j, = |∇|(xj,t)are approximated by least-square fitting over the elements surrounding the node xj. It is worth noting that the least squares fitting strategy is a very flexible tool. For instance, it allows introducing the upwind by simply increasing the weight of the elements located upstream from the point of interest. This is an attractive feature of the methodology because as pointed out by many authors [6], upwinding is required to stabilize the transport schemes in the presence of sharp fronts. In the application examples we include in the next section, however, the solutions are pretty stable and upwinding did not introduce any significant improvement. The discretized form of the Hamilton–Jacobi equation is then j,+1=j, −vn(xj,t)t|∇|j,,(17) Note that, following the remark of the previous section, it might be very interesting here to use a mesh for the representation of the level-set function independent from the one used to compute the normal velocity field. This new mesh could be structured and therefore allow for the use of FD schemes for the discretization along the space dimension. This may provide for substantial savings of computational power. The derivation of the Hamilton–Jacobi equation stems from two basic hypotheses: (1) the space points for which the equation is defined are along the interface =0, and (2) the level-set function is a distance function. This means that the Hamilton–Jacobi equation should not be solved on the entire domain but rather only close to the interface. Indeed, by solving it on the entire domain, there is no certainty that the level-set obtained would remain a distance function to the interface. This is particularly critical when the controlling parameters on the two sides of the interface are very different (in this paper, the ratio of permeability coefficients is 1000). After having moved the interface (by modifying the value at the nodes close to that interface), the value of the level-set function at all other nodes should therefore be obtainedthrough the actual computationofthe distance function. Note that this is a computationally involved step for which specific methods can be implemented [32,33]. Note also that some authors [9,34] have devised iterative algorithms to solve the Hamilton–Jacobi equation on the entire domain, with the added objective of obtaining directly a distance function for the solution. 6