
Citation: Daniele Cerroni, Florin Adrian Radu, Paolo Zunino. Numerical solvers for a poromechanic problem with a moving boundary[J]. Mathematics in Engineering, 2019, 1(4): 824-848. doi: 10.3934/mine.2019.4.824
[1] | Marguerite Champion, Miguel A. Fernández, Céline Grandmont, Fabien Vergnet, Marina Vidrascu . On the analysis of a mechanically consistent model of fluid-structure-contact interaction. Mathematics in Engineering, 2024, 6(3): 425-467. doi: 10.3934/mine.2024018 |
[2] | Andrea Manzoni, Alfio Quarteroni, Sandro Salsa . A saddle point approach to an optimal boundary control problem for steady Navier-Stokes equations. Mathematics in Engineering, 2019, 1(2): 252-280. doi: 10.3934/mine.2019.2.252 |
[3] | Sudarshan Tiwari, Axel Klar, Giovanni Russo . Interaction of rigid body motion and rarefied gas dynamics based on the BGK model. Mathematics in Engineering, 2020, 2(2): 203-229. doi: 10.3934/mine.2020010 |
[4] | Stefano Berrone, Stefano Scialò, Gioana Teora . The mixed virtual element discretization for highly-anisotropic problems: the role of the boundary degrees of freedom. Mathematics in Engineering, 2023, 5(6): 1-32. doi: 10.3934/mine.2023099 |
[5] | Francesca Bonizzoni, Davide Pradovera, Michele Ruggeri . Rational-approximation-based model order reduction of Helmholtz frequency response problems with adaptive finite element snapshots. Mathematics in Engineering, 2023, 5(4): 1-38. doi: 10.3934/mine.2023074 |
[6] | Yves Achdou, Ziad Kobeissi . Mean field games of controls: Finite difference approximations. Mathematics in Engineering, 2021, 3(3): 1-35. doi: 10.3934/mine.2021024 |
[7] | Juan Pablo Borthagaray, Wenbo Li, Ricardo H. Nochetto . Finite element algorithms for nonlocal minimal graphs. Mathematics in Engineering, 2022, 4(2): 1-29. doi: 10.3934/mine.2022016 |
[8] | Antonello Gerbi, Luca Dedè, Alfio Quarteroni . A monolithic algorithm for the simulation of cardiac electromechanics in the human left ventricle. Mathematics in Engineering, 2019, 1(1): 1-37. doi: 10.3934/Mine.2018.1.1 |
[9] | Tommaso Tassi, Alberto Zingaro, Luca Dede' . A machine learning approach to enhance the SUPG stabilization method for advection-dominated differential problems. Mathematics in Engineering, 2023, 5(2): 1-26. doi: 10.3934/mine.2023032 |
[10] | Daniela De Silva, Giorgio Tortone . Improvement of flatness for vector valued free boundary problems. Mathematics in Engineering, 2020, 2(4): 598-613. doi: 10.3934/mine.2020027 |
Reconstructing the stress and deformation history of a sedimentary basin is a challenging and important problem in the geosciences [32]. The mechanical response of a sedimentary basin is the consequence of complex multi-physics processes involving mechanical, geochemical, thermal and hydrological aspects. Among these, the inclusion of erosion or deposition of material plays a key role in the reconstruction of a sedimentary basin history.
A sedimentary basin is usually represented by a porous medium filled with water. To model surface erosion or deposition, we let the top surface of the basin evolve according to a given profile. The dynamics of erosion/deposition is affected by many complex phenomena, which are only partially accounted for in this study. Since the scope of this work is more focused on the methodology, precisely the development of a numerical solver, we address a simplified problem configuration which features a slab of porous, deformable material that is homogeneous, completely saturated with water and affected by a prescribed mixed erosion deposition profile on the top surface. In this way, we combine in a sigle model the mechanical, hydrological and erosion/deposition effects.
In this work we tackle the evolution of the top surface using the cut finite element method (CutFEM). This method is combined with a fixed-stress approach for solving the coupled flow and mechanics models. To deal with a moving boundary, the physical domain, is embedded in a larger computational grid. The evolution of top surface of the slab is described by the zero value of a time dependent level set function. The CutFEM approximation allows to handle the evolution of the top layer of the slab avoiding re-meshing techniques. However, depending on the intersection between the computational grid and the level set function, it can worsen the condition number of the problem. To overcome the ill-conditioning due to the CutFEM approximation, we propose the use of a stabilization terms and of a preconditioner applied to the fixed-stress splitting scheme.
On the basis of seminal works using the finite element method for the approximation of elliptic problems featuring unfitted boundaries and interfaces [9,15,20,21], the CutFEM method was developed and formalized as a general framework for simultaneously discretizing the geometry of the domain and the related partial differential equations [8]. Despite its many advantages, it is well known that this approach significantly worsens the conditioning of the discretized operators, because the finite element discretization is set on cut elements that may become arbitrarily small and skewed, according to the configuration of the unfitted boundary or interface with respect to the computational mesh. To override this issue, several techniques have been proposed. Stabilization techniques acting on the problem formulation have been successfully applied to this purpose [7,9,10,25,29]. It has been also shown that the ill-conditioning of the CutFEM discretization can be addressed by means of suitable preconditioners, [24,28,34]. Such techniques make CutFEM become a stable, efficient and robust method for a wide variety of problems, including second order elliptic equations [11,15,20], incompressible flow problems such as Stokes [19,22,25,28] and Navier-Stokes [29] and linear elasticity models [21]. However, the application of CutFEM to multi-physics problems involving the coupled CutFEM approximation of multiple variables is still rather unexplored. The scope of the present study is to shed light on some of the difficulties that appear in this context and to show how the available theory of CutFEM can be used to overcome these issues. In particular, we address a poromechanics problem modeled by the Biot equations. We believe that this is a particularly interesting case because it couples the equations for incompressible flow in porous media with a linear elasticity model for the deformations of the material. Indeed, it puts together the problems for which CutFEM has been thoroughly studied.
The work is organized as follows. In Section 2 we introduce the mathematical formulation and the discretization of the problem. In Section 3 we discuss the relation between the condition number of the discrete problem and the size of cut elements. Then, we develop the new numerical solver for Biot equations, based on the fixed stress approach combined with stabilization and preconditioners for CutFEM. Illustrative numerical examples are performed to test the proposed algorithm in Section 4, followed by concluding remarks in Section 5.
In this section we present the mathematical setting of the poromechanics problem. Then, we introduce the numerical discretization, obtained by casting the problem into the CutFEM framework.
We deal with the mechanical deformation of a porous medium. We consider for simplicity an isotropic material (named the skeleton) filled with an isothermal single-phase fluid. We also assume that the Hypothesis of Small Perturbations (HSP) holds true, so that the linearity of the problem is ensured. In particular HSP implies several consequences that are described in the following. The material particles of the skeleton are subject to small displacements. This hypothesis allows us to identify the initial spatial configuration and the current configuration of the system. Small variations of the porosity and the fluid density are also considered, namely
|ϕ−ϕ0ϕ0|≃0,|ρf−ρf0ρf0|≃0, | (2.1) |
where the index (0) refers to the reference/initial configuration. This assumption allows us to replace ϕ and ρf with ϕ0 and ρf0 whenever required. We remark that the poromechanics problem is solved only in the physical part of the domain, namely Ω(t). A sketch of a domain example is show in Figure 1.
In particular, it is shown that the upper boundary is denoted by Γ(t) and it is free to move. The lower boudary is denoted with Σ and it is fixed. The lateral sides of the domain are ∂Ω(t)∖(Γ(t)∪Σ). Under the small perturbations hypothesis, the poromechanics problem is modeled through the following system of equations:
−∇⋅σ(u)+α∇p=fin Ω(t), | (2.2) |
∂t(pM+α∇⋅u)−∇⋅K∇p=0in Ω(t), | (2.3) |
where u is the solid matrix displacement vector with respect to the unloaded reference configuration, represented in Figure 1. The function p is the pore pressure and the forcing term f is the volumetric force due to the mass of the porous medium, namely
f(t,x)=−((1−ϕ)ρs+ϕρf)R(t)gez | (2.4) |
being ρs,ρf the densities of the solid and fluid phases, respectively, g the standard gravity, ez the unit vector in the vertical direction and R(t) a ramp function R(t)=η2max[(t+δtgs+|t+δtgs|),1] where η is a suitable slope coefficient to be defined later on and δtgs is the temporal interval needed to reach the geostatic conditions.
We remind that ∂t denotes the standard partial derivative with respect to time in the Eulerian framework. The parameters α, M and K are the Biot number, the Biot modulus and the the hydraulic conductivity, respectively. We recall that the hydraulic conductivity is defined as the ratio between the permeability (ks) and the dynamic viscosity of the fluid (μf), namely K=ks/μf. To complete the definition of the problem we also assume the linear elasticity behavior for the skeleton. This implies that the stress tensor σ, appearing in (2.2), is defined by
σ(u):=2με(u)+λ∇⋅u, | (2.5) |
where μ and λ are the Lamé coefficients and ε(u) is the symmetric gradient of the skeleton displacement, defined as
ε(u):=12(∇u+∇ut). | (2.6) |
For further details on poromechanics, the interested reader is referred to e.g. [13,14].
For a well-posed problem we must complement the previous governing equations with appropriate boundary and initial conditions. We consider the following,
p=0,σ(u)⋅n=0,on Γ(t), | (2.7) |
K∇p⋅n=Kρfg,u=0,on Σ, | (2.8) |
K∇p⋅n=0,σ(u)⋅n=0,on∂Ω(t)∖(Γ(t)∪Σ), | (2.9) |
where n is the unit outward normal to the boundary. Concerning the initial condition we consider the basin to be in the reference unloaded state so that u=0, p=0 enforced at the beginning of simulations t=−δtgs. We point out that, as illustrated in Figure 2, by the time t=0 the basin has reached the geostatic state characterized by u=ugs, p=pgs.
By taking into account the conditions (2.7), (2.8), (2.9), the weak formulation of the problem described in (2.2)–(2.3) reads: for each t∈(0,T], find (u,p)∈V(t)×Q(t) such that
{a(u,v)−α(p,∇⋅v)=(f,v)∀v∈V(t),1M(∂tp,q)+α(∇⋅∂tu,q)+c(p,q)−∫Γ(t)K∇p⋅nq=∫ΣKρfgq∀q∈Q(t), | (2.10) |
where V(t):={u∈(H1(Ω(t)))n∣u∣Σ=0} and Q(t):={p∈H1(Ω(t))∣p∣Γ(t)=0}. Here, H1(Ω) denotes the Hilbert subspace of L2(Ω(t)) of functions with first weak derivatives in L2(Ω(t)), (⋅,⋅) is the standard inner product in the space L2(Ω(t)). Equation (2.10) is characterized by the following bilinear forms
a(u,v):=2∫Ω(t)με(u):ε(v)dΩ+∫Ωλ(∇⋅u)(∇⋅v)dΩ,c(p,q):=∫Ω(t)K∇p⋅∇qdΩ. |
The main difficulty of problem (2.10) consists of the time-dependence of the domain and of the functional spaces. We aim at removing this dependence at the level of numerical discretization by casting the discrete problem into the CutFEM framework.
In the cut finite element method, briefly CutFEM, the boundary of the physical domain is represented on a background grid using a level set function. The background grid is also used to approximate the solution of the governing problem. In this framework the boundary conditions are built into the discrete variational formulation, leading to a method that can handle a time varying domain avoiding mesh moving algorithms or re-meshing techniques. A comprehensive review which also cover the implementation aspects of the CutFEM method can be found in [8,19].
Let ΩT be a bounded computational domain, constant in time, with an interface Γ(t) dividing ΩT into two non-overlapping subdomains Ω(t) and Ωout(t), so that ¯ΩT=¯Ω(t)∪¯Ωout(t), as show in Figure 1. The active zone, Ω(t), represents the physical domain, while Ωout is a fictitious region.
Let Th:={T} denote a triangulation of ΩT, not necessarily conformal to the interface Γ(t). The collection of elements that are intersected by Γ is denoted with ΩΓ:=T∈Th:T∩Γ≠∅.
Let us introduce the discrete spaces
Vrh:={vh∈H1(ΩT):vh|T∈Pr(T),∀T∈Th}d,Qsh:={qh∈H1(ΩT):qh|T∈Ps(T),∀T∈Th}, |
where Ps denotes the space of scalar piecewise polynomials of order s.
We consider a decomposition of the finite element spaces Vrh,Qsh as proposed in [28] and also used in [34]. For simplicity, we address the decomposition for Vrh solely, the one for Qsh follows similarly. We separate the nodes related to each finite element space into the ones neighboring the interface IΓ(Vrh) from the remaining ones in the active part of the domain IΩ(Vrh) and those in the fictitious domain IΩout(Vrh). Precisely, let I(Vrh) be the set of all degrees of freedom of Vrh and let ϕk be the corresponding nodal basis functions. The degrees of freedom related to basis functions entirely supported in Ω and Ωout are IΩ(Vrh):={k: supp(ϕk)⊂Ω} and IΩout(Vrh):={k: supp(ϕk)⊂Ωout} respectively. Their complementary set is IΓ(Vrh)=I(Vrh)∖(IΩ(Vrh)∪IΩout(Vrh). Then, we use those sets to construct the following subspaces,
VΓh:=span{ϕk:k∈IΓ(Vrh)},VΩh:=span{ϕk:k∈IΩ(Vrh)},VΩouth:=span{ϕk:k∈IΩout(Vrh)}, |
and we decompose the finite element space as Vrh=VΓh⊕VΩh⊕VΩouth, so that any test function vh∈Vrh can be uniquely represented as as vh=vΓh+vΩh+vΩouth.
The semi-discrete scheme can be written as follows: find (uh,ph)∈Vrh×Qsh such that:
{a(uh,vh)−α(ph,∇⋅vh)=(f,vh)∀vh∈Vrh,1M(∂tph,qh)+α(∇⋅∂tuh,qh)+c(ph,qh)−∫Γ(t)K∇ph⋅nqh=∫ΣKρfgqh∀qh∈Qsh. | (2.11) |
Since all the integrals are defined either on Ω or on Γ, it would be equivalent to use test functions defined on VΓh⊕VΩh solely. In practice, we find more convenient to use formulation (2.11) and to set the degrees of freedom of IΩout(Vrh) and IΩout(Qsh) to zero.
For the imposition of the pressure constraint in the internal unfitted interface Γ(t) in the second equation of (2.11), we use Nitsche's method following the approaches proposed in the seminal work [20], see also [6,22] for applications of the method to Stokes problem and poromechanics. This technique allows to weakly enforce interface conditions at the discrete level by adding appropriate penalization terms to the variational formulation of the problem. This results in the modified problem formulation,
1M(∂tph,qh)+(α∇⋅∂tuh,qh)+c(ph,qh)−∫Γ(t)K∇ph⋅nqh−∫Γ(t)K∇qh⋅nph+γKh∫Γ(t)phqh=∫ΣKρfgqh∀qh∈Qsh, | (2.12) |
where h is the characteristic size of the quasi-uniform computational mesh and γ>0 denotes a penalization parameter independent on h and on the parameters of the problem. Finally, we consider a backward-Euler time discretization scheme. Denoting with τ the computational time step, the fully discretized problem at time tn, n=1,2,...,N reads as follows: find (unh,pnh)∈Vrh×Qsh such that
{a(unh,vh)−α(pnh,∇⋅vh)=(f,vh),1M(pnh,qnh)+α(∇⋅unh,qh)+τc(pnh,qh)−τ∫Γ(t)K∇pnh⋅nqh−τ∫Γ(t)K∇qh⋅npnh+γKhτ∫Γ(t)pnhqh=∫ΣKρfgqh+1M(pn−1h,qh)+α(∇⋅un−1h,qh) | (2.13) |
for all vh∈Vrh and qh∈Qsh. We notice that the finite element space of pn−1h may not coincide with the one of qh, relative to the current time tn. To compute the integrals on the right hand side, we project pn−1h onto the space QΓh⊕QΩh, such that the degrees of freedom of the projected pn−1h and qh are conformal. If IΩ,n−1(Qsh)∪IΓ,n−1(Qsh)⊂IΩ,n(Qsh)∪IΓ,n(Qsh) the projection is just a simple extension to zero of the new active degrees of freedom. Otherwise, we just discard at time n the degrees of freedom that were active at time n−1. We finally remark that bilinear forms of problem (2.13) are defined either on Ω or Γ. As a result the degrees of freedom IΩout(Vrh),IΩout(Qsh) are not properly constrained. In practice, all the matrix blocks related to these degrees of freedom are set to the identity, with a corresponding null right hand side, such that all the degrees of freedom on Ωout are constrained to vanish.
We focus the attention on the solution of the system (2.13) through the fixed stress iterative scheme as presented e.g. in [4,23,26,30] and used in several studies, see for example [1,3,17,33] as a solver and preconditioner for poromechanics problems.
In this algorithm the flow is solved first, followed by the solution of the mechanical problem. At each time step, the procedure is iterated until the solution converges within an acceptable tolerance.
The first step reads as follows: given (un,kh,pn,kh)∈Vrh×Qsh find pn,k+1h∈Qsh such that
(1M+β)(pn,k+1h,qh)+τc(pn,k+1h,qh)−τ∫ΓK∇pn,k+1h⋅nqhdΓ−τ∫ΓK∇qh⋅npn,k+1hdΓ+γKhτ∫Γpn,k+1hqhdΓ=∫ΣKρfgqh+1M(pn−1h,qh)+(α∇⋅un−1h,qh)+β(pn,kh,qh)−(α∇⋅un,kh,qh) ∀qh∈Qsh, | (2.14) |
where (⋅k) denotes the iteration index. The second step of the algorithm reads: given pn,k+1h∈Qsh find un,k+1h∈Vrh such that
a(un,k+1h,vh)=(f,vh)+α(pn,k+1h,∇⋅vh)∀vh∈Vrh. | (2.15) |
The iterative steps are performed until the following convergence criterion is fulfilled
pn,k+1h−pn,khpn,0h<ηp,un,k+1h−un,khun,0h<ηu, | (2.16) |
where ηp and ηu are the desired tolerances (here chosen equal to 10−7).
For notational convenience, we reformulate the fixed-stress steps in a more compact form as
step 1Bppn,k+1h=Lp(un−1h,un,kh,pn,kh,pn−1h),step 2Buun,k+1h=Lu(pn,k+1h). | (2.17) |
We will refer at Bp and Bu as the system matrices of the first (pressure) and the second (displacement) step of the fixed-stress algorithm.
The real number β≥0 is a stabilization parameter and has to be chosen such that the convergence of the scheme is ensured. We refer to e.g. [23,30] for the choice of β and to [1,4,5,26,31] for a mathematical interpretation and analysis of it. In particular, the optimization of the parameter β is performed in [31]. We remark that all the works mentioned above are referring to classical discretizations of the Biot equations, for the CutFEM method presented in this paper the convergence analysis of the iterative method is thoroughly addressed in Section 3.6.
In this section we analyze the performance of the original fixed stress algorithm, namely (2.17), in order to investigate the effects of the CutFEM discretization in the solution of the problem. To this purpose we use two computational test cases arising from different choices of the finite element space for the discretization of the displacement. In the first test case we take (uh,ph)∈V1h×Q1h while in the second case we choose (uh,ph)∈V2h×Q1h. We consider a cubic computational domain characterized by an edge of 4 Km. The physical domain is described by the negative part of the following level set function
f(x,y,z)=(z−(0.7+14((y/l)2−1)((x/l)2−1))l)(1−t+1.26×1011)+(z−0.844l)t+1.26×1011, | (2.18) |
where t+ is the physical time in the temporal window [0,4Ky] (Ky stands for kilo-years) and l is the characteristic length of the domain. A sketch of the initial geometry of the physical domain is shown in the top left part of Figure 2, while the material properties considered are reported in Table 1. According to (2.18), Γ(t) is such that in some points [x,y] the elevation of ground surface decreases as time evolves, while it increases in some other points. With this choice, we aim at demonstrating the capability of the approach to model erosion as well as deposition of material. On a sufficiently large spatial scale, the two phenomena must necessarily take place together, because of mass conservation.
Young Modulus | E | 1010 | Pa |
Rock density | ρs | 2.2×103 | Kg/m3 |
Rock permeability | ks | 10−13 | m2 |
Water density | ρl | 103 | Kg/m3 |
Water viscosity | μl | 10−3 | Pas |
Dimension | l | 4×103 | m |
To reach the geostatic condition we consider a temporal interval of 3.99Ky in wich the ramp function R(t), introduced in Eq. (2.4), reads R(t)=η2max[(t+δtgs+|t+δtgs|),1] with δtgs=3.99Ky and η=(tr+δtgs+|tr+δtgs|)−1 with tr≈−1Ky. A schematic of the simulation set up is reported in the timeline in the bottom part of Figure 2. We notice that the simulation is split in three parts: before t1=0Ky, between t1=0Ky and t3=3.99Ky and after t3=3.99Ky. The first part from t0=−δtgs to t1 is used to reach the geostatic condition that are reached at t1=0Ky. The second part of the simulation, namely from t1=0Ky to t3=3.99Ky is when the erosion/deposition take place, while between t3 and tend≈8Ky there isn't any change in the basin configuration. We use a time step size of 3.15⋅109s≈0.1Ky so that the simulation requires 40 time steps for the first, the second and the third part. The level set function is updated just for 40 time steps after t1=0 so that remains steady for t>t3=3.99Ky.
The performance of the numerical solver is measured by three indices:
a=npnmax,p,b=nunmax,u,c=nfs, | (2.19) |
where np and nu are the average number of GMRES iterations needed for the solution of the pressure and the displacement subproblem, namely step 1 and step 2 of (2.17), nmax,p and nmax,u are the maximum number of iterations for the pressure and the displacement problem respectively, while the number of fixed stress iterations per time step is nfs.
The behavior of a and b over time for the case (uh,ph)∈V1h×Q1h is reported in the left part of Figure 3. We notice that the average number of iterations needed for the solution of the pressure and the displacement sub-problem (a and b indices) is highly variable and a few peaks can be observed. We remark that the indices a and b are related to the convergence rate of GMRES and, as a consequence, to the conditioning of the two sub-problems. Depending on how the level set intersects the computational domain, the two sub-problems may become ill-conditioned and a larger number of iterations is needed for the fixed stress algorithm to converge. In the first test case, we notice that the indices a and b do not reach their maximum (unit) value, meaning that the sub-problems are properly solved at every time step. In such condition, as we see from the right panel of Figure 3 (c index), the number of fixed stress iterations at each time step remains almost constant.
The behavior of the three indices for the test case (uh,ph)∈V2h×Q1h is reported in Figure 4. We notice that the variability of the indices a and b is greater that in the previous case and furthermore, at certain time steps, the b index reaches its maximum. More precisely, the second step of the fixed stress algorithm does not reach convergence. The index c of Figure 4, shows that the ill-conditioning of the sub-problems affects also the global convergence of the fixed stress algorithm. In particular, we notice that the c index follows the same trend of b and severe peaks occur when the solver of the subproblem b does not converge.
From these preliminary test cases we infer that as long as the sub-problems are accurately solved, the global convergence of the fixed stress algorithm is not significantly affected, while is highly deteriorated when one of the two steps is not properly resolved. Then, the aim of this work is exploring different numerical techniques to overcome the dependence of the sub-problems form the size of cut-elements and, as a consequence, of the fixed stress algorithm from the position of the unfitted interface.
In this section we investigate the spectral properties of the systems associated to the two steps of the fixed stress algorithm. We show that the CutFEM approximation may cause the ill-conditioning of the sub-problems, depending on how the level set function intersects the computational grid. In the second part of the section, we introduce a stabilization term and a preconditioner that cure the ill-conditioning of the CutFEM approximation method.
In the CutFEM framework interface conditions are enforced weakly and numerical integration over sub-elements is required. The integrals over very small sub-elements lead to very small element matrix entries and therefore to ill-conditioned systems. Since the introduction of CutFEM various authors have proposed different techniques to avoid this type of ill-conditioning [7,9,10,24,25,28,29,34]. Among this non exhaustive list, we refer to [7], where the ghost penalty stabilization is introduced.
The idea of [7] is to add a penalty term in the interface zone that extends the coercivity of the (physical) operator also in the region of the computational domain where the solution has no physical significance. Let us introduce the following norms for any qh∈Qsh (for simplicity we use the scalar notation, but these definitions can be straightforwardly extended to vector valued functions)
||qh||0,Ω:=(qh,qh)1/2Ω, ||qh||0,Γ:=(qh,qh)1/2Γ, ||qh||1,h,Ω:=||∇qh||0,Ω+h−1/2||qh||0,Γ. |
Following [7] we introduce the following discrete inequalities that will be used in the analysis of the problem:
Poincaré-Friedrichs inequality||qh||0,Ω≤CP||qh||1,h,Ω, | (3.1) |
Inverse inequality||qh||1,h,Ω≤CIh−1||qh||0,ΩT, | (3.2) |
where the constants CP,CI are independent of the size of cut elements. Furthermore, we need estimates of the discrete functions in terms of their degrees of freedom. More precisely, we denote the degrees of freedom of qh∈Qsh with qh:=∑Ni=1qiϕh,i as q:={qi}Ni=1∈RN. Then, there exist constants cmin,0 and cmax,0 such that
c2min,0hd|q|2≤||qh||20,ΩT≤c2max,0hd|q|2. | (3.3) |
When properties of discrete functions are considered on meshes that are cut by the boundary, the smallest size of cut elements may affect the validity such properties. To define the smallest cut element, we assume there exist a node xn such that the basis function related to such node have small support in the physical domain, i.e. measd(∪T:xn∈TT∩Ω)<<hd. Let xh be the nodal basis function related to that node. Let zn be the node and zh be the basis function such that measd(∪T:zn∈TT∩Ω)/measd(T) is minimal among the whole mesh. Then we define the smallest (relative) support of a basis function defined on cut element as
ν:=measd(∪T:zn∈TT∩Ω)measd(∪T:zn∈TT). | (3.4) |
In the case of a mesh cut by the boundary, following [34] and in particular Lemma 3.3 and Lemma 3.4 respectively, for piecewise linear finite elements, namely r=1,s=1, and for the special function zh related to definition (3.4), we have that there exist constants cmin,1,cmin,2,cmax,1,cmax,2 such that
c2min,1hdν2/d+1|z|2≤||zh||20,Ω≤c2max,1hdν2/d+1|z|2, | (3.5) |
c2min,2hd−2ν|z|2≤||zh||21,h,Ω≤c2max,2hd−2ν|z|2. | (3.6) |
Finally, in Lemma 3.7 of [34] it is shown that there exists a constant CT, independent of the mesh cuts, such that for any qh∈Q1h it holds,
hν1/d||qh||20,Γ≤CT||qh||20,Ω, ∀qh∈Q1h. | (3.7) |
On the basis of these results, in the following we present the analysis of the two sub-problems of the fixed-stress approach, showing how the interface position affects their conditioning in both the standard and the stabilized formulation.
Equation (2.14), with a generic right hand side, can be written in a more compact form as
mh(ph,qh)+τch(ph,qh)=(g,qh), | (3.8) |
where
mh(ph,qh):=(1M+β)(ph,qh)Ω,ch(ph,qh):=(∇ph,∇qh)−(∇ph⋅n,qh)Γ−(∇qh⋅n,ph)Γ+γKh(ph,qh)Γ. |
Before proceeding, we observe that using standard arguments based on the norms and the inequalities defined in the previous Section, it is straightforward to prove the following properties of the pressure problem, with constants mc and Mc uniformly independent to the mesh characteristic size and the cut elements,
Coercivitymc||ph||21,h,Ω≤ch(ph,ph), | (3.9) |
Continuity|ch(ph,qh)|≤Mc||ph||1,h,Ω||qh||1,h,Ω. | (3.10) |
To study the conditioning of the system matrix associated with (3.8), we address the spectral properties of the matrix Bp of (2.17). It is the sum of a mass ad a stiffness matrix, namely Bp=M+C. The eigenvalues/eigenfunctions of (3.8), namely ϑ,yh are the solution of the following problem
mh(yh,qh)+τch(yh,qh)=ϑ(yh,qh)Ω. | (3.11) |
Introducing a suitable parameter CB, equation (3.11) can be written as
τch(yh,qh)=(ϑ−CB)(yh,qh)Ω=ζ(yh,qh)Ω, | (3.12) |
and, as a consequence, the eigenvalues ϑ are related to the eigenvalues ζ associated to the operator ch. For this reason, we study the spectral properties of the matrix associated to (3.12).
Let k(C):=||C||||C−1|| be the condition number of a generic matrix C, where ||C||:=supx∈RN|Cx||x| and |⋅| denotes the euclidean norm on RN. Let us consider x∈RN, associated to xh∈Qsh such that supp(xh)⊂Ω, and therefore ||xh||0,Ω=||xh||0,ΩT. It can be shown, using the coercivity relation (3.9), the Poincaré inequality (3.1) and (3.3) that
||C||=supw∈RN|Cw||w|≥(Cx,x)|x|2≥ch(xh,xh)|x|2≥mc||xh||21,h,Ω|x|2≥mcC−2P||xh||20,Ω|x|2≥mcC−2Pc2min,0hd. |
We focus now on ||C−1||. By definition we have
||C−1||=supw∈RN|C−1w||w|=supw∈RN|w||Cw|≥|z||Cz|. |
Let us now take the special function zh of (3.4),
|Cz|=supw∈RN(Cz,w)|w|=supw∈RNch(zh,wh)|w|. |
Using the continuity of the bilinear form (3.10) combined with (3.6) applied to zh, (3.3) applied to wh, and reminding that |z|=1 we obtain that
ch(zh,wh)≤Mc||zh||1,h,Ω||wh||1,h,Ω≤Mccmax,2h(d−2)/2ν12|z|||wh||1,h,Ω≤Mccmax,2h(d−2)/2ν12CIh−1|z|||wh||0,ΩT≤Mccmax,2h(d−2)/2ν12CIcmax,0h−1hd/2|z||w|. | (3.13) |
In conclusion we have that
k(C)=||C||||C−1||≥mcC−2Pc2min,0McCIcmax,0cmax,2h2ν−1/2 |
meaning that the conditioning of the pressure sub-problem tends to infinity when ν→0.
Let EB be the set of element edges in the boundary zone, more precisely,
EB:={K∩K′, where either K∈ΩΓ or K′∈ΩΓ}. | (3.14) |
Following [11], for any ph and qh∈Q1h the ghost penalty stabilization operator j(⋅,⋅) is given by
jp(ph,qh)=∑E∈EB(γgh∇ph⋅nE,∇qh⋅nE)E, | (3.15) |
and Equation (3.8) is modified as follows,
mh(ph,qh)+τch(ph,qh)+τjp(ph,qh)=(g,qh). | (3.16) |
Then, the bilinear form of the stabilized pressure problem is
Ch(ph,qh):=ch(ph,qh)+jp(ph,qh), | (3.17) |
and we denote by G the related matrix.
In [7] it is shown that the stabilization operator does not affect the continuity and the coercivity of the problem, such that there exist constants mC,MC
mC|||ph|||21,h,ΩT≤Ch(ph,ph),|Ch(ph,qh)|≤MC|||ph|||1,h,ΩT|||qh|||1,h,ΩT. |
Using (3.3) and the standard inverse inequality we obtain,
||G||=supp∈RN|Gp||p|=supp∈RNsupw∈RNCh(ph,wh)|w||p|≤MCC2Ic2max,0hd−2. |
For ||G−1|| using the coercivity, the inverse inequality and the bound (3.3) we have,
mC|||ph|||21,h,ΩT≤Ch(ph,ph)≤|Gp||p|≤|Gp|c−1min,0h−d/2||ph||0,ΩT. |
Using the Poincaré type inequality (3.1) and the bound (3.3)
|||ph|||21,h,ΩT≥C−2P||ph||20,ΩT≥C−2Pcmin,0hd/2||ph||0,ΩT|p|. |
As a result we have
|Gp||p|≥mCC−2Pc2min,0hd,∀p∈RN. |
Using the previous inequality we estimate ||G−1|| as
||G−1||=supw∈RN|G−1w||w|=supp∈RN|p||Gp|≤1mCC−2Pc2min,0hd. |
It follows that k(G)=||G||||G−1||≃h−2 is not dependent by the smallest cut element. The condition number of G indeed scales with h as expected for the standard FEM method.
We focus on the operator a(u,v) defined as
a(uh,vh):=2∫Ωinμε(uh):ε(vh)dΩ+∫Ωinλ(∇⋅uh)(∇⋅vh)dΩ. | (3.18) |
For the upper bound of the bilinear form it can be shown that the following inequalities hold true,
|a(uh,vh)|≤2μ||ε(uh)||0,Ω||ε(vh)||0,Ω+λ||∇⋅uh||||∇⋅vh||≤(2μ+λd)||ε(uh)||0,Ω||ε(vh)||0,Ω≤Ma||uh||1,Ω||vh||1,Ω. | (3.19) |
For the lower bound of the bilinear form we have straightforwardly that a(uh,uh)≥ma||∇uh||20,Ω. Using Korn's inequality, in particular ||uh||0,Ω≤CK||∇uh||0,Ω,∀uh∈[H10(Ω)]d, (see for example [16]) we obtain that
a(uh,uh)≥maC−2K||uh||20,Ω,∀uh∈[H10(Ω)]d. |
We now estimate the condition number of the system matrix A as k(A)=||A||||A−1||. Let us consider xh∈V1h such that supp(xh)⊂Ω. Therefore we have ||xh||0,Ω=||xh||0,ΩT. Using this function in the coercivity of a(⋅,⋅) and exploiting the bound (3.3) we obtain
||A||=supw∈RdN|Aw||w|≥(Ax,x)|x|2≥a(xh,xh)|x|2≥mac2min,0hd. | (3.20) |
As in the previous case we evaluate ||A−1|| as
||A−1||=supw∈RdN|A−1w||w|=supw∈RdN|w||Aw|≥|z||Az|, | (3.21) |
where z is the vector of degrees of freedom of the function zh related to (3.5), (3.6). Then, there holds
|Az|=supw∈RdNa(zh,wh)|w|≤Ma||zh||1,Ω||wh||1,Ω|w|, | (3.22) |
and then following (3.13) we obtain
|Az|≤MaCIcmax,0cmax,2hd−2ν12|z|. |
The condition number is therefore bounded from below by the constant
k(A)=||A||||A−1||≥maC−2Kc2min,0MaCIcmax,0cmax,2h2ν−1/2. |
We notice that, also in the displacement subproblem, the conditioning of the system matrix depends on the smallest cut element.
We address the second step of the fixed-stress approach (2.17) associated to the linear elasticity problem that reads,
Buun,k+1h=Lu(pn,k+1h). | (3.23) |
On the basis of the decomposition Vrh=VΓh⊕VΩh⊕VΩouth, after filtering out the redundant degrees of freedom of VΩouth, the matrix Bu admits the following block representation,
Bu=[BΩuBΩΓu(BΩΓu)tBΓu], |
where
(uΩh)tBΩuvΩh:=a(uΩh,vΩh),(uΓh)tBΓuvΓh:=a(uΓh,vΓh),(uΩh)tBΩΓuvΓh:=a(uΩh,vΓh). |
The system (3.23) is naturally ill-conditioned because it consists of an almost incompressible elasticity problem. In addition, the corresponding matrix is assembled on possibly small cut-elements. We simultaneously reduce the stiffness of the system and cure the instability due to the small cut-elements by means of a preconditioner and a stabilization term. Following [34] we assume that the only genuinely stiff block of Bu is BΩu. As a result we propose the following preconditioner
P:=[BΩu00diag(BΓu)]. | (3.24) |
We solve the problem of Bu by means of preconditioned Krylov iteration that requires the inversion of P. The solution of the diagonal lower block sub-system diag(BΓu) is trivial, while the inversion of the upper block diag(BΩu) is computationally demanding. We observe that the latter operation is the solution of a standard linear elasticity problem. We achieve this task by means of an Algebraic Multigrid (AMG) algorithm, which has been successfully used for this class of operators, as discussed in [18].
We also adopt a stabilization operator for the elasticity problem, to ease the ill-conditioning related to the presence of small cut-elements. Following the approach described in section 3.3 we add to equation (2.15) the following ghost penalty operator
ju(uh,vh)=∑E∈EBr∑i=1(γgh2i−1∂iuh⋅nE,∂ivh⋅nE)E∀uh,vh∈Vrh, | (3.25) |
where ∂i denotes the i−th order derivative and EB is the set of element edges in the boundary zone defined in Eq. (3.14). This is a more general form of the ghost-penalty stabilization operator introduced for the pressure in section 3.3, which can be applied either to linear or quadratic finite elements, namely the case uh∈Vrh with r=1 and 2. For a comprehensive review of the ghost penalty stabilization approach we remand the interested reader to [12]. In conclusion, the stabilized form of the displacement sub-problem reads as,
a(uh,vh)+ju(uh,vh)=(f,vh)+α(ph,∇⋅vh)∀vh∈Vrh. | (3.26) |
In this section we investigate the convergence of the fixed stress algorithm in relation with ill-conditioning of the two sub-problems, due to the presence of small cut-elements. For this reason, we address the pressure and elasticity problems of (2.17) in absence of the ghost penalty stabilizations. Anyway, it is easy to see that the same conclusions of this section apply to the stabilized formulation including jp(⋅,⋅) and ju(⋅,⋅).
Let us introduce the splitting error of the fixed stress approach as
ekp=ph−pkhandeku=uh−ukh. | (3.27) |
We subtract the fixed-stress equations form the monolithic problem formulation. The pressure equation becomes
(1M+β)(ekp,qh)+τc(∇ekp,∇qh)+τγKh(ekp,qh)Γ+τ(K∇ekp⋅n,qh)Γ−τ(K∇qh⋅n,ekp)Γ=β(ek−1p,qh)−α(∇⋅ek−1u,qh),∀qh∈Qsh. | (3.28) |
The displacement error equation reads
2μ(ε(eku),ε(vh))+λ(∇⋅eku,∇⋅vh)=α(ekp,∇vh),∀vh∈Vrh. | (3.29) |
First, we choose the test function as qh=ekp and vh=ek−1u. Equation (3.28) becomes now
(1M+β)||ekp||20,Ω+τK||∇ekp||20,Ω+τγKh||ekp||20,Γ−2Kτ(∇ekp⋅n,ekp)Γ=β(ek−1p,ekp)−α(∇⋅ek−1u,ekp), | (3.30) |
while (3.29) gives
2μ(ε(eku),ε(ek−1u))+λ(∇⋅eku,∇⋅ek−1u)=α(ekp,∇ek−1u). | (3.31) |
Using the polarization identity a⋅b=1/4(a+b)2−1/4(a−b)2 we get
(ε(eku),ε(ek−1u))=14||ε(eku)+ε(ek−1u)||20,Ω−14||ε(eku)−ε(ek−1u)||20,Ω, | (3.32) |
(∇⋅eku,∇⋅ek−1u)=14‖∇⋅eku‖20,Ω−14‖∇⋅ek−1u‖20,Ω, | (3.33) |
and by the the binomial identity (a−b,a)=12a2+12(a−b)2−12b2 we obtain
β||ekp||0,Ω−β(ek−1p,ekp)=β(ekp−ek−1p,ekp)=12||ekp||20,Ω+12||ekp−ek−1p||20,Ω−12||ek−1p||20,Ω. | (3.34) |
Regarding the terms arising from the CutFEM approximation, using (3.7) we obtain,
τK(∇ekp⋅n,ekp)Γ≤τK2δh−1||ekp||2Γ+δ2τKh||∇ekp⋅n||2Γ≤τK2δh−1||ekp||20,Γ+12δτKCTν−1/d||∇ekp||20,Ω. |
By choosing δ=ν1/dCT, 1−δ2CTν1/d=12 and γ=CTν1/d we obtain
τK(1−δ2CTν1/d)||∇ekp||20,Ω+τKh(γ−12δ)||ekp||20,Γ=τK2(||∇ekp||20,Ω+CThν1/d||ekp||20,Γ). | (3.35) |
Using (3.30), (3.31) (3.32), (3.33), (3.34) and (3.35) we get
2M||ekp||20,Ω+τK||∇ekp||20,Ω+τKCThν1/d||ekp||2Γ+β||ekp||20,Ω+β||ekp−ek−1p||20,Ω+μ||ε(eku)+ε(ek−1u)||20,Ω+λ2||∇⋅eku+∇⋅ek−1u||20,Ω−μ||ε(eku)−ε(ek−1u)||20,Ω−λ2||∇⋅eku−∇⋅ek−1u||20,Ω≤β||ek−1p||20,Ω. | (3.36) |
Second, we address the terms arising from the splitting error of the displacement problem. From equation (3.31) taken at iteration k and k−1 and using the test function vh=eku−ek−1u, we have
2μ‖ε(eku−ek−1u)‖20,Ω+λ‖∇⋅(eku−ek−1u)‖20,Ω=α(ekp−ek−1p,∇⋅(eku−ek−1u)). | (3.37) |
We recall that
||∇⋅(eku−ek−1u)||20,Ω≤d||ε(eku−ek−1u)||20,Ω. |
Using Cauchy-Schwarz and Young's inequality we obtain,
α(ekp−ek−1p,∇⋅(eku−ek−1u))≤α2212μd+λ‖ekp−ek−1p‖20,Ω+12(2μd+λ)‖∇⋅(eku−ek−1u)‖20,Ω≤α2212μd+λ‖ekp−ek−1p‖20,Ω+μ||ε(eku−ek−1u)||20,Ω+λ2‖∇⋅(eku−ek−1u)‖20,Ω. |
As a result we have,
μ||ε(eku−ek−1u)||20,Ω+λ2||∇⋅(eku−ek−1u)||20,Ω≤α22(2μd+λ)||ekp−ek−1p||20,Ω. | (3.38) |
Equation (3.38) shows the convergence of the displacement is dominated by the convergence of the pressure. Then, combining inequalities (3.36) and (3.38) we obtain,
2M||ekp||20,Ω+τK||∇ekp||20,Ω+τKCThν1/d||ekp||2Γ+β||ekp||20,Ω+(β−α22(2μd+λ))||ekp−ek−1p||20,Ω+μ||ε(eku)+ε(ek−1u)||20,Ω+λ2||∇⋅eku+∇⋅ek−1u||20,Ω≤β||ek−1p||20,Ω, | (3.39) |
which, on the basis of the Banach contraction theorem, proves that ||ek−1p||20,Ω→0, provided that
β≥α22(2μd+λ). |
Finally, owning to the Poincaré inequality, for the splitting error of the pressure we obtain,
||∇ekp||20,Ω+CThν1/d||ekp||20,Γ≥C−1Pmin(CTν1/d,1)||ekp||20,Ω, | (3.40) |
showing that the global convergence of the fixed stress algorithm is not affected by the position of the unfitted interface Γ, because when ν→0 the constant of (3.40) is uniformly bounded. We conclude that as long as the two sub-problems are well posed, in the case of piecewise linear finite elements for the pressure approximation, the global convergence of the fixed stress algorithm depends only on the material properties of the problem and we expect a constant number of fixed stress iteration. This theoretical result is in accordance on what is obtained in the preliminary test case where the c index remains almost constant as shown in the right part of Figure 3.
In this section we investigate the performance of the stabilization term introduced in Section 3.2 and of the preconditioner described in Section 3.5, by means of several numerical experiments. In Test 1 we consider a simple two dimensional problem, because in such simplified configuration the smallest (relative) size of the cut element (ν) is easily controlled. In the second test, named Test 2, we consider a more realistic three dimensional case.
We consider a squared domain as shown in the left part of Figure 5 with a characteristics length l=4 Km divided into the physical (Ω) and fictitious region (Ωout) by the curve Γ. The material properties are the ones reported in the right panel of Figure 1. Regarding the boundary conditions, we consider vanishing displacement and null fluid flow on the bottom surface y=0, while free stress condition and vanishing pressure are imposed on the remaining faces. The surface Γ is defined as the zero value of the level set function Γ(y)=y−H. Using this type of level set function and a uniform triangulation, we control the minimum dimension of the cut element as shown in the right panel of Figure 5.
Both the first and the second step of the fixed stress algorithm (2.17) (pressure and displacement step, respectively) are solved by means of GMRES. The number of iterations performed by each step is used to investigate the performance of the pressure stabilization and the preconditioner for the displacement. In Figure 6 (panel A) we show the conditioning of the matrices associated to the fixed-stress sub-problems with respect to the minimum cut-element size. From that figure we notice that the numerical experiments are consistent with the lower bound (dashed line) derived in Section 3. However, we notice that the theoretical bound is not sharp, because when ν tends to zero the conditioning of both matrices explode with a rate higher than ν−12. When the pressure stabilization and the preconditioning of the displacement are applied, the condition number of the pressure and the displacement systems become insensitive to the parameter ν, as shown in Figure 6 (panel B). The effect of the pressure stabilization term can be seen in Figure 6 (panel C), where we show the number of GMRES iterations needed by the pressure step of the fixed stress algorithm. We notice that the introduction of the stabilization term makes the number of iterations become insensitive to the parameter ν. Concerning the solver for the displacement, the application of the preconditioner (3.24) requires the inversion of the matrix
P:=[BΩu00diag(BΓu)]. |
In alternative to P, we also consider the diagonal matrix P1 defined as
P1:=[diag(BΩu)00diag(BΓu)]. |
The number of GMRES iterations for the displacement system using these preconditioners, is shown in Figure 6 (Panel D). We notice that both preconditioners optimally cure the stiffness of the matrix due to presence of small cut-elements, because in both cases the convergence rate of GMRES seems to be independent of ν. However, the full preconditioner P is more effective, because it highly reduces the number of GMRES iterations.
In this section we show and comment the results for a more realistic test case. We consider a tridimensional slab subject to mechanical compaction. The physical domain is described by the intersection of a cube, of edge l=4 Km, and the negative part of the level set function defined in (2.18).
A sketch of the physical domain is shown in the top part of Figure 2 at different time steps. We model the slab as a water saturated porous material, characterized by the properties reported in the Table 1. Such parameters are in the range used for modeling sedimentary basins. We remand the interested reader to [2,13,14,27] for further details.
Concerning the boundary conditions we assume that the base of the slab is impermeable and fixed. The lateral surfaces are impermeable, while a free stress condition is imposed for the displacement problem. On the unfitted boundary at the top, we prescribe an homogeneous Dirichlet and an homogeneous Neuman conditions for the pressure and the displacement, respectively. We study the evolution of the slab in a temporal interval as described in Section 2.4. In particular for the second part of the simulation, after the geostatic condition is reached, we investigate the evolution of the slab for 8Ky using a time step of approximately 0.1Ky. In the top part of Figure 7 we show the displacement field at t1=0, t2=1.99 and t3=3.99Ky. The deformed domains (in blue) are obtained warping the current configuration (in gray) by the displacement field magnified by a factor of 60 so that they can be appreciated. We notice that the displacement field is highly driven by the shape of the physical domain. In particular when the top surface is not planar, the largest deformation occurs in correspondence of the thickest zones of the slab while, in the final part of the simulation when the level set is planar, the displacement field becomes almost constant in the cross section of the slab. Similar considerations also apply to the pressure field, as confirmed by the lower part of Figure 7 where we show the geostatic pressure on the left and the difference between the current and the geostatic pressure at t2 and t3.
We now investigate the performance of the numerical solver based on the fixed stress algorithm applied to Test 2. In particular we show that the pressure stabilization together with the preconditioning of the elasticity system are effective to reduce the computational cost of the problem. We run numerical experiments for two discretization choices, (uh,ph)∈V1h×Q1h and (uh,ph)∈V2h×Q1h.
In Figure 8 we show the three indices introduced in Section 2.4, obtained for the case (r,s)=(1,1). We compare the performance of the fixed stress algorithm in the original set up (black curves) and using the stabilization term for the pressure problem, namely jp(⋅,⋅) and the preconditioner P for the elasticity problem, as described in (3.24) (red curves).
From the index a we notice that, using the stabilization term the number of iterations for the pressure problem is reduced. Precisely, the peak that occurs near the time step 30 for the black curve, is no longer present for the red one. Focusing on the index b, we notice that the preconditioner P effectively reduces the number of iterations needed for the solution of the displacement problem. Finally, from the index c we see that the number of fixed stress iterations is almost constant both in the original and the stabilized formulation of the problem. This observation supports the conclusions of Section 3.6, confirming that the global convergence of the fixed stress does not depend on the position of the level set, resulting in an almost constant number of iterations.
The results obtained in the case of ((r,s)=(1,2)) are shown in Figure 9. In particular we compare the following cases:
(nsnp) non stabilized, non preconditioned (black line); it refers to the original formulation of the problem.
(sp) pressure-stabilized and preconditioned (red line); it refers to the stabilized formulation of the pressure problem and the preconditioner P for the elasticity problem.
(ssp) pressure and displacement stabilized and preconditioned (blue line); these data are are obtained using the ghost penalty stabilization ju(⋅,⋅) in the elasticity problem.
Again, we see that the number of iterations for each subproblem of the fixed stress algorithm decreases for both the sp and ssp cases. Concerning the index b, the overall number of iterations for the displacement significantly lowers in the sp and ssp cases, despite some peaks that are still present. However, the major advantage of using the preconditioner and the stabilization term of the displacement is appreciated from the index c. More precisely, we notice that the peaks of the fixed stress iterations drastically decrease in the sp case and almost vanish in the ssp cases, resulting in a great reduction in the computational cost of the whole coupled system.
In this work we have considered a free-boundary poromechanic problem, relevant for studying the erosion/deposition of a sedimentary basin. The mathematical model for this problem consists of Biot's equations with a moving boundary. We have proposed a discretization method based on CutFEM and the fixed-stress splitting scheme. The use of CutFEM avoids re-meshing the computational domain along with the evolution of the boundary. However, it worsens the conditioning of the matrices in the presence of cut-elements of small size. To overcome this issue, we have introduced a stabilisation term into the problem formulation and used a preconditioner for the numerical solver. Illustrative numerical examples have shown that these numerical techniques are effective and enable the application of the combined CutFEM/fixed stress approach to realistic applications in geomechanics.
The authors D.C. and P.Z. acknowledge financial support from ENI S.p.A. The authors F.A.R. and P.Z. thank Equinor for the travel grant through the Akademia program.
The authors declare no conflict of interest in this paper.
[1] |
Bause M, Radu FA and Koecher U (2017) Space-time finite element approximation of the Biot poroelasticity system with iterative coupling. Comput Method Appl M 320: 745–768. doi: 10.1016/j.cma.2017.03.017
![]() |
[2] | Bense VF and Person MA (2008) Transient hydrodynamics within intercratonic sedimentary basins during glacial cycles. during glacial cycles. J Geophys Res-Earth 113. |
[3] |
Borregales M, Kumar K, Radu FA, et al. (2019) A partially parallel-in-time fixed-stress splitting method for Biot's consolidation model. Comput Math Appl 77: 1466–1478. doi: 10.1016/j.camwa.2018.09.005
![]() |
[4] |
Both JW, Borregales M, Nordbotten JM, et al. (2017) Robust fixed stress splitting for Biot's equations in heterogeneous media. Appl Math Lett 68: 101–108. doi: 10.1016/j.aml.2016.12.019
![]() |
[5] | Both JW, Kumar K, Nordbotten JM, et al. (2018) Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media. Comput Math Appl 77: 1479–1502. |
[6] |
Bukač M, Yotov I, Zakerzadeh R, et al. (2015) Partitioning strategies for the interaction of a fluid with a poroelastic material based on a Nitsche's coupling approach. Comput Method Appl M 292: 138–170. doi: 10.1016/j.cma.2014.10.047
![]() |
[7] | Burman E (2010) Ghost penalty. CR Math 348: 1217–1220. |
[8] |
Burman E, Claus S, Hansbo P, et al. (2015) Cutfem: discretizing geometry and partial differential equations. Int J Numer Meth Eng 104: 472–501. doi: 10.1002/nme.4823
![]() |
[9] | Burman E and Hansbo P (2012) Fictitious domain finite element methods using cut elements: II. a stabilized Nitsche method. Appl Numer Math 62: 328–341. |
[10] | Burman E and Hansbo P (2014) Fictitious domain methods using cut elements: III. a stabilized Nitsche method for Stokes' problem. ESAIM-Math Model Num 48: 859–874. |
[11] | Burman E, Hansbo P and Larson M (2018) A cut finite element method with boundary value correction. Math Comput 87: 633–657. |
[12] | Burman E and Zunino P (2012) Numerical Approximation of Large Contrast Problems with the Unfitted Nitsche Method. In: Blowey J and Jensen M (Eds.) Frontiers in Numerical Analysis - Durham 2010. (pp. 227-282). Springer Berlin Heidelberg: Berlin, Germany. |
[13] | Cheng AHD (2016) Poroelasticity. volume 27, Springer. |
[14] | Coussy O (2004) Poromechanics. John Wiley & Sons. |
[15] |
Dolbow J and Harari I (2009) An efficient finite element method for embedded interface problems. Int J Numer Meth Eng 78: 229–252. doi: 10.1002/nme.2486
![]() |
[16] | Ern A and Guermond JL (2013) Theory and practice of finite elements. Springer. |
[17] |
Gaspar FJ and Rodrigo C (2017) On the fixed-stress split scheme as smoother in multigrid methods for coupling flow and geomechanics. Comput Method Appl M 326: 526–540. doi: 10.1016/j.cma.2017.08.025
![]() |
[18] |
Griebel M, Oeltz D and Schweitzer MA (2003) An algebraic multigrid method for linear elasticity. SIAM J Sci Comput 25: 385–407. doi: 10.1137/S1064827502407810
![]() |
[19] | Gross S and Reusken A (2011) Numerical methods for two-phase incompressible flows. volume 40 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin. |
[20] |
Hansbo A and Hansbo P (2002) An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. Comput Method Appl M 191: 5537–5552. doi: 10.1016/S0045-7825(02)00524-8
![]() |
[21] |
Hansbo A and Hansbo P (2004) A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput Method Appl M 193: 3523–3540. doi: 10.1016/j.cma.2003.12.041
![]() |
[22] |
Hansbo P, Larson MG and Zahedi S (2014) A cut finite element method for a Stokes interface problem. Appl Numer Math 85: 90–114. doi: 10.1016/j.apnum.2014.06.009
![]() |
[23] |
Kim J, Tchelepi HA and Juanes R (2011) Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits. Comput Method Appl M 200: 1591–1606. doi: 10.1016/j.cma.2010.12.022
![]() |
[24] |
Lehrenfeld C and Reusken A (2017) Optimal preconditioners for Nitsche-xfem discretizations of interface problems. Numer Math 135: 313–332. doi: 10.1007/s00211-016-0801-6
![]() |
[25] |
Massing A, Larson MG, Logg A, et al. (2014) A stabilized Nitsche fictitious domain method for the Stokes problem. J Sci Comput 61: 604–628. doi: 10.1007/s10915-014-9838-9
![]() |
[26] |
Mikeli´ c A and Wheeler MF (2013) Convergence of iterative coupling for coupled flow and geomechanics. Computat Geosci 17: 455–461. doi: 10.1007/s10596-012-9318-y
![]() |
[27] |
Nasir O, Fall M, Nguyen ST, et al. (2013) Modeling of the thermo-hydro-mechanical–chemical response of sedimentary rocks to past glaciations. Int J Rock Mech Min 64: 160–174. doi: 10.1016/j.ijrmms.2013.08.002
![]() |
[28] |
Reusken A (2008) Analysis of an extended pressure finite element space for two-phase incompressible flows. Computing and Visualization in Science 11: 293–305. doi: 10.1007/s00791-008-0099-8
![]() |
[29] |
Schott B and Wall WA (2014) A new face-oriented stabilized xfem approach for 2d and 3d incompressible Navier-Stokes equations. Comput Methods Appl M 276: 233–265. doi: 10.1016/j.cma.2014.02.014
![]() |
[30] |
Settari A and Mourits FM (1998) A coupled reservoir and geomechanical simulation system. SPE J 3: 219–226. doi: 10.2118/50939-PA
![]() |
[31] |
Storvik E, Both JW, Kumar K, et al. (2019) On the optimization of the fixed-stress splitting for Biot's equations. Int J Numer Meth Eng 120: 179–194.. doi: 10.1002/nme.6130
![]() |
[32] |
Tuncay K, Park A and Ortoleva P (2000) Sedimentary basin deformation: an incremental stress approach. Tectonophysics 323: 77–104. doi: 10.1016/S0040-1951(00)00095-0
![]() |
[33] |
White JA, Castelletto N and Tchelepi HA (2016) Block-partitioned solvers for coupled poromechanics: A unified framework. Comput Method Appl M 303: 55–74. doi: 10.1016/j.cma.2016.01.008
![]() |
[34] |
Zunino P, Cattaneo L and Colciago CM (2011) An unfitted interface penalty method for the numerical approximation of contrast problems. Appl Numer Math 61: 1059–1076. doi: 10.1016/j.apnum.2011.06.005
![]() |
1. | Daniele Cerroni, Mattia Penati, Giovanni Porta, Edie Miglio, Paolo Zunino, Paolo Ruffo, Multiscale Modeling of Glacial Loading by a 3D Thermo-Hydro-Mechanical Approach Including Erosion and Isostasy, 2019, 9, 2076-3263, 465, 10.3390/geosciences9110465 | |
2. | Mathias Anselmann, Markus Bause, Higher order Galerkin–collocation time discretization with Nitsche’s method for the Navier–Stokes equations, 2020, 03784754, 10.1016/j.matcom.2020.10.027 | |
3. | Zhijun Liu, Yimin Zhang, Yao Jiang, Han Yang, Yongtao Yang, Unfitted finite element method for fully coupled poroelasticity with stabilization, 2022, 397, 00457825, 115132, 10.1016/j.cma.2022.115132 | |
4. | Zhijun Liu, Yuxin Tong, Yimin Zhang, Hong Zheng, Fanyu Zhang, Stabilized Unfitted Finite Element Method for Poroelasticity With Weak Discontinuity, 2025, 130, 2169-9313, 10.1029/2024JB030155 |
Young Modulus | E | 1010 | Pa |
Rock density | ρs | 2.2×103 | Kg/m3 |
Rock permeability | ks | 10−13 | m2 |
Water density | ρl | 103 | Kg/m3 |
Water viscosity | μl | 10−3 | Pas |
Dimension | l | 4×103 | m |