c
TU Delft, The Netherlands, 2006
ANALYSIS OF FINITE ELEMENT METHODS FOR
COUPLED HYDRO-MECHANICAL PROBLEMS
A. Ern∗, S. Meunier∗,† and G. Nicolas† ∗CERMICS, ENPC,
6-8 avenue Blaise Pascal, Cit´e Descartes
77455 Champs-sur-Marne, Marne la Vall´ee Cedex 2, France e-mail: ern@cermics.enpc.fr
Web page: http://cermics.enpc.fr/∼ern/home.html e-mail: meunier@cermics.enpc.fr
Web page: http://cermics.enpc.fr/∼meunier
†EDF R&D, 1 avenue du G´en´eral de Gaulle,
92141 Clamart Cedex, France e-mail: gerald.nicolas@edf.fr
Key words: Finite element methods, Hydro-Mechanical coupling, a posteriori error analysis, Error indicator, Adaptive mesh, Code Aster
Abstract. This work is concerned with the analysis of a non-stationary Hydro-Mechanical problem. The well-posedness of the continuous, time semi-discretized and fully discretized problems is addressed. An a posteriori error analysis is performed, leading to a family of space-time error indicators. The reliability of this family, i.e. a global upper bound of the error using the error indicators, is proved.
1 INTRODUCTION
Darcy’s law. The mathematical formulation of the fully dynamic Biot system is ∂t(ρ∂tu) − ∇·σ0(u) + b∇p = f, ∂t 1 Mp + b∇·u − ∇· (κ∇p) = g, with σ0(u) = λM 1 (∇·u) 1 + λM2 ∇u+ ∇ut ,
where the unknowns u and p respectively denote the displacement of the solid matrix and the pressure distribution, σ0(u) is the so-called effective stress tensor, ρ denotes the total density of the medium, b is the Biot-Willis constant and M the Biot modulus. The quantity κ denotes the hydraulic permeability, f is the volume-distributed external forces and g denotes any volume-distributed source density.
In this work, we restrict our considerations to the linear quasi-static case (ρ = 0, κ constant). The proofs of the results presented below are omitted ; we refer to [2] for full detail.
2 ANALYSIS OF THE CONTINUOUS PROBLEM
In this section, a stability result for the continuous problem is stated. 2.1 Problem formulation
Let T > 0 be a fixed time, let Ω be a bounded open set whose boundary is denoted by Γ. Consider two partitions of the boundary in the form Γ = ΓM
D ∩ ΓMN = ΓHD∩ ΓHN. We focus on the following problem
−∇·σ0(u) + b∇p = f, in Ω×]0, T [, ∂t M1 p + b∇·u − κ∆p = g, in Ω×]0, T [, u = 0, in ΓM D×]0, T [, (σ0(u) − βbpχ S1) ·n = σnor, in ΓMN×]0, T [, p = 0, in ΓH D×]0, T [, −b∂t((1 − β) u·n) χS+ κ∇p·n = φnor, in ΓHN×]0, T [, u(x, 0) = u0, in Ω, p(x, 0) = p0, in Ω, (1)
The quantity β designates the fraction of the matrix pores which are sealed along the boundary [3]. The complement consists of those for which the diffusive flow paths are directly exposed by the boundary. We denote by χS the characteristic function of ΓS = ΓM
For the sake of simplicity, we restrict ourselves to the case β = 0 without loss of generality. Furthermore, solutions are supposed to be sufficiently regular, namely
u∈ C1[0, T ];H1(Ω)3, p ∈ C1 [0, T ]; H1(Ω).
Existence and uniqueness of solutions in a weaker sense are analyzed, e.g., in [3]. Fur-thermore, a possible incompatibility between initial conditions and boundary conditions, that would result in regularity loss at t = 0+, is not discussed here. The weak formulation of problem (1) is the following :
Seek (u, p) ∈ C1 [0, T ], UM,0× C1 [0, T ], UH,0 such that ∀t ∈ [0, T ], ∀v ∈ UM,0, Z Ω σ0(u):(v) + Z Ω b∇p·v = Z Ω f·v + Z ΓM N σnor·v, ∀q ∈ UH,0, Z Ω ∂t 1 Mp + b∇·u q + Z Ω κ∇p·∇q − Z ΓS b∂t(u·n) q = Z Ω gq + Z ΓH N φnorq, (2) with the initial conditions
u(·, 0) = u0, p(·, 0) = p0, where UM = [H1(Ω)]3, UH= H1(Ω) and UM,0 = {v ∈ UM/v|ΓM D = 0}, UH,0 = {q ∈ UH/q|ΓH D = 0}. 2.2 Stability
Define the mechanical energy and the hydraulic power of the system respectively by EM(u) = Z Ω σ0(u):(u) = Z Ω λM1 (∇·u)2+ 2λM2 (u)2, EH(p) = Z Ω κ (∇p)2,
and for all t ∈ [0, T ], for x = (u, p), introduce the norm |||x|||(t) = 1 2EM(u)(t) + 1 2 Z Ω 1 Mp 2(·, t) + Z t 0 EH(p) (s)ds 1 2 .
• u0 ∈ [H1(Ω)] 3 , p0 ∈ L2(Ω), • σnor ∈ C1 [0, T ]; L2(ΓMN) , f ∈ C1([0, T ]; L2(Ω)), • φnor ∈ L2(]0, T [; L2(ΓHN)), g ∈ L2(]0, T [; L2(Ω)), • b constant, M bounded, ∂tλM1 = ∂tλM2 = ∂tM = 0.
Under the above assumptions, we have the following stability property [2] Theorem 2.1. For all t ∈ [0, T ],
|||x|||2(t) .1 2EM(u0) + 1 2 Z Ω 1 Mp 2 0+ kfk20,Ω(0) + kσnork20,ΓM N(0) + T Z t 0 k∂tfk20,Ω + T Z t 0 k∂tσnork20,ΓM N + Z t 0 kgk20,Ω+ Z t 0 kφnork20,ΓH N.
3 ANALYSIS OF THE TIME SEMI-DISCRETIZED PROBLEM
In this section, we investigate the time semi-discretization of problem (1) by an implicit Euler scheme. We state the well-posedness of the problem. Then, we perform an a posteriori-type error analysis.
3.1 The semi-discretized problem
We consider a partitioning of the time interval [0, T ] in the form [0, T ] =SNn=1[tn−1, tn] with 0 = t0 < t1 < · · · < tN = T . For all n ∈ [[1, N]], we denote τn = tn− tn−1 by the length of the time interval [tn−1, tn]. In this section, the notation f . g means that there exists c > 0 independent of T and of the time partitioning such that f ≤ cg.
We consider a regular time partitioning, namely for all n ∈ [[1, N]], 1 . τn
τn−1 .1.
Without loss of generality, we assume that max1≤n≤N(τn) . 1. For all q ∈ N? and for any function g from Ω × [0, T ] to Rq, we set
∀n ∈ [[0, N]], gn= g(·, t n).
We use the lower suffix τ to denote a continuous, time-affine function defined over Ω × [0, T ]. Since gn
τ = gτ(·, tn) owing to the above definition, we clearly have ∀t ∈ [tn−1, tn], gτ(·, t) = t − tn−1 τn gn τ + tn− t τn gn−1 τ . We define the solution spaces
The time semi-discretization of problem (2) issued from Euler’s implicit scheme is the following :
Seek (uτ, pτ) ∈ UM,0τ × UτH,0 such as for all n ∈ [[1, N]] ∀v ∈ UM,0, Z Ω σ0(unτ):(v) + Z Ω b∇pnτ·v = Z Ω fn·v + Z ΓM N σnnor·v, ∀q ∈ UH,0, Z Ω 1 M pn τ − pn−1τ τn + b∇·u n τ − ∇·un−1τ τn q + Z Ω κ∇pn τ·∇q − Z ΓS bu n τ − un−1τ τn ·nq = Z Ω gnq + Z ΓH N φnnorq, (3) with the initial conditions
u0τ = u0 et p0τ = p0.
We introduce the norm |||·|||τ(t), defined for any function y = (v, q) ∈ C0 [0, T ], UM,0 × C0 [0, T ], UH,0 as |||y|||τ(t) = 1 2EM(v)(t) + 1 2 Z Ω 1 Mq 2(·, t) + X m≥1,tm≤t τmEH(qm) !1 2 . The following theorem holds [2]
Theorem 3.1. The problem (3) is well-posed : • There exists one and only one solution
• For all n ≥ 1, we have the following stability property : for x = (uτ, pτ),
|||xτ|||2τ(tn) . 1 2EM(u0) + 1 2 Z Ω 1 Mp 2 0+ n X m=1 τmkgmk20,Ω+ τmkφmnork20,ΓH N+ kf m k20,Ω+ kσmnork20,ΓM N . 3.2 A posteriori-type error analysis
In this paragraph, an a posteriori-type error estimation for the time semi-discretized problem is derived. Before, let us introduce some notation. For all q ∈ N?, for any function g defined over Ω × [0, T ] and taking values in Rq, we denote by Π0
τg the projection of g onto time-piecewise constant functions and we denote by Π1
τg the interpolant of g by time-piecewise affine and continuous functions. More specifically, we set for all t ∈ ]tn−1, tn],
Π0 τg(·, t) = g n, Π1τg(·, t) = t − tn−1 τn gn+tn− t τn gn−1.
Theorem 3.2. For all t ∈ [0, T ], for x = (u, p) and xτ = (uτ, pτ), we have |||x − xτ|||2(t) . Z t 0 g − Π0 τg 2 0,Ω+ φnor − Π0τφnor 2 0,ΓH N + Z t 0 κ12∇ p τ − Π0τpτ 2 0,Ω + T Z t 0 ∂t f − Π1τf 2 0,Ω+ ∂t σnor− Π1τσnor 2 0,ΓM N .
4 ANALYSIS OF THE FULLY DISCRETIZED PROBLEM
This section deals with the space discretization of the time semi-discretized problem (3). We investigate the well-posedness of the problem and we present an a posteriori error analysis.
4.1 The fully discretized problem
The open set Ω is assumed to be a polyhedron in R3, so that the latter can be exactly covered by a standard finite element mesh. Let Tnhbe the mesh of Ω at time tn. The mesh depends on the n suffix because it may change at each time step. Let k and l two positive integers, k being the interpolation degree for displacements and l that for pressures. We introduce the solution spaces
UM,0nh = {v ∈C0 Ω3∩ UM,0/∀K ∈ Tnh, v|K ∈ [Pk]3}, UnhH,0 = {q ∈ C0 Ω∩ UH,0/∀K ∈ Tnh, q|K ∈ Pl}.
In this section, the notation f . g means that there exists c > 0 such that f ≤ cg, where c is independent of T , of the time subdivision (tn)0≤n≤N and of the meshes (Tnh)0≤n≤N.
The meshes are supposed to be shape-regular in the usual sense of Ciarlet. Hence, the spaces UM,0nh and UnhH,0 are endowed with the following local interpolation property : for all (v, q) ∈ UM,0× UH,0, there exists (v
h, qh) ∈ UM,0nh × UnhH,0 such that for all K ∈ Tnh and for all F ⊂ ∂K, h−1K kv − vhk0,K+ h −12 F kv − vhk0,F .kvk1,DK, h−1K kq − qhk0,K + h −12 F kq − qhk0,F .kqk1,DK,
where hK is the diameter of K, hF is the diameter of F and DK is the set of simplices that share at least one point with K. The fully discretized problem is the following :
Seek (uhτ, phτ) ∈ UM,0hτ × U H,0
with the initial conditions
u0hτ = ΠMhu0 et p0hτ = ΠHhp0, where ΠM
h and ΠHh stand for appropriate projection operators taking their values in U M,0 0h and U0hH,0 respectively, and where
UM,0hτ = {vhτ continuous time-affine function/∀n ∈ [[0, N]], vhτn ∈ U M,0 nh }, UH,0hτ = {qhτ continuous time-affine function/∀n ∈ [[0, N]], qnhτ ∈ U
H,0 nh }. The following theorem holds [2]
Theorem 4.1. The problem (4) is well-posed. • There exists an unique solution.
• For all n ≥ 1, we have the following stability property |||xhτ|||2τ(tn) . 1 2EM(u0) + 1 2 M− 1 2p 0 2 0,Ω + n X m=1 τmkghτmk20,Ω+ τm φm nor,hτ 2 0,ΓH N + kf m hτk20,Ω+ σm nor,hτ 2 0,ΓM N . 4.2 A posteriori error analysis
In this paragraph, we are interested in deriving a posteriori space-time error indicators for the fully discretized problem (4). The literature dealing with a posteriori error analysis for parabolic problems is fairly extensive. Early work focuses on time adaptation. For example, Johnson, Nie and Thom´ee [4] dealt with an a posteriori error analysis in a linear case for the time scheme. Nochetto, Savar´e and Verdi [5] presented an a posteriori study for non linear parabolic problems. Picasso [6] derived an a posteriori analysis for a space-time discretization. Estimations are shown under the somewhat restrictive hypothesis that meshes are nested in time, which prevents mesh coarsening. Verf¨urth [7] performed an a posteriori analysis for a non-stationary advection-diffusion-reaction problem discretized by a space-time finite element method. The robustness of the error indicators when advection effects dominate is addressed. This leads to global space error estimations. Bergam, Bernardi and Mghazli [8] derived a posteriori space-time indicators with reliability and optimality properties without requiring the meshes to be nested in time. The goal of the present work is to generalize the techniques presented in [8] to mixed elliptic-parabolic type problems such as non-stationary HM problems.
Our family of error indicators is composed of one error indicator related to time dis-cretization and of one error indicator related to space disdis-cretization. For all n ∈ [[1, N]], we define the time error indicator as
For all n ∈ [[1, N]], for each element K ∈ Tnh, we define the local mechanical and hydrauli-cal spatial error indicators respectively by
ηM n,K =hKkfn+ ∇·σ0(unhτ) − b∇p n hτk0,K+ 1 2 X F∈Fi K h 1 2 F k[[σ0(u n hτ)·n]]k0,F + X F∈F∂ K∩Γ M N h 1 2 F σn nor,τ − n·σ0(u n hτ) 0,F, ηHn,K =hK g n− 1 M pn hτ − p n−1 hτ τn − b∇· u n hτ − u n−1 hτ τn + ∇· (κ∇pn hτ) 0,K +1 2 X F∈Fi K h 1 2 Fk[[κ∇p n hτ·n]]k0,F + X F∈F∂ K∩Γ H N h 1 2 F φnnor− κ∇pnhτ·n + b unhτ − un−1hτ τn ·nχS 0,F , where Fi
K and FK∂ respectively stand for faces of K located in Ω and for the faces of K located on the boundary of Ω. For all φ ∈ [L2(Ω)]3
, for all F ∈ Fi
K, we denote by [[φ·n]] the jump of the normal component of φ across F .
Our goal is to bound the error |||x − xhτ|||(tn) (where x = (u, p) solves (1) and xhτ = (uhτ, phτ) solves (4)) by a combination of the time error indicators ηn and the space error indicators ηM
n,K and ηn,KH . This property ensures the reliability of the exhibited error indicators. To this purpose, we proceed as follows : we first bound the quantity |||x − xτ|||(tn), and then the quantity |||xτ − xhτ|||(tn). To conclude, we use a triangle inequality.
Lemma 4.2. For all n ∈ [[1, N]], |||x − xτ|||2(tn) . Z tn 0 g − Π0 τg 2 0,Ω+ φnor− Π0τφnor 2 0,ΓH N + T Z tn 0 ∂t f − Π1τf 2 0,Ω+ ∂t σnor− Π1τσnor 2 0,ΓM N + |||xτ − xhτ|||2(tn) + n X m=1 ηm2.
Lemma 4.2 is a straightforward consequence of Theorem 3.2. Let us now bound |||xτ − xhτ|||2(tn).
If the data (f, σnor, g and φnor) are not polynomials, the space error indicators cannot be computed exactly. Hence, we consider modified error indicators in which forces and loadings are discretized in space by polynomials over each element of the mesh Tnh. We denote by ζlh the discretization of a function ζ by l-degree polynomials. This discretiza-tion can be constructed by Lagrange finite element interpoladiscretiza-tion or by projecdiscretiza-tion over discontinuous finite element spaces. Let
b ηMn,K =hK fn l1h+ ∇·σ 0(un hτ) − b∇p n hτ 0,K+ 1 2 X F∈Fi K h 1 2 Fk[[σ 0(un hτ)·n]]k0,F + X F∈F∂ K∩Γ M N h 1 2 F σn nor,l2h− n·σ 0(un hτ) 0,F , b ηHn,K =hK g n l3h− 1 M pn hτ − p n−1 hτ τn − b∇· u n hτ − u n−1 hτ τn + κ∆pn hτ 0,K +1 2 X F∈Fi K h 1 2 Fk[[κ∇p n hτ·n]]k0,F + X F∈F∂ K∩Γ H N h 1 2 F φnnor,l4h− κ∇p n hτ·n + b unhτ − un−1hτ τn ·nχS 0,F .
Lemma 4.3 can be rewritten as Lemma 4.4. For all n ∈ [[1, N]],
|||xτ − xhτ|||2(tn) . n X m=1 τm X K∈Tmh b ηHm,K2+ h2K gm− gm l3h 2 0,K + X F∈F∂ K∩Γ H N hF φm nor− φmnor,l4h 2 0,F + n X m=1 X K∈Tmh bηM m,K 2 + h2K fm− fm l1h 2 0,K + X F∈F∂ K∩Γ M N hF σm nor− σmnor,l2h 2 0,F + p0τ− ΠHhp0 2 0,Ω+ EM u 0 τ − ΠMhu0 .
Collecting the results of Lemmas 4.2 and 4.4 leads to the following reliability theo-rem [2].
Theorem 4.5. For all n ∈ [[1, N]],
with the data error indicators δ2M=T ∂t f − Π1τf 2 L2(0,tn,L2(Ω))+ ∂t σnor − Π1τσnor 2 L2(0,tn,L2(ΓM N)) + n X m=1 X K∈Tmh h2 K fm − flm1h 2 0,K+ X F∈F∂ K∩Γ M N hF σm nor− σ m nor,l2h 2 0,F + EM u0τ − ΠMh u0 , δH2 = g − Π0τg 2 L2(0,t n,L2(Ω))+ φnor− Π0τφnor 2 L2(0,t n,L2(ΓHN)) + n X m=1 τm X K∈Tmh h2 K gm − glm3h 2 0,K+ X F∈F∂ K∩Γ H N hF φm nor− φ m nor,l4h 2 0,F + p0 τ − ΠHhp0 2 0,Ω. 5 CONCLUSION
This work provides the first a posteriori error analysis for a non-stationary Hydro-Mechanical coupled problem. Space-time error indicators are analyzed. A reliability property is obtained by proving a global upper bound of the error. Future work includes the derivation of lower bounds for the error and the extension of the a posteriori error analysis to (non-linear) THM problems. The a posteriori error indicators are currently being implemented in EDF’s finite element program Code Aster, which offers a full range of multiphysical analysis and modelling methods : seismic analysis, porous environments, acoustics, fatigue, stochastic dynamics, etc. Its modelling, algorithms and solvers are constantly undergoing work to improve them.
REFERENCES
[1] Official website of Code Aster, www.code-aster.org
[2] S. Meunier, “Analyse d’erreur a posteriori pour une mod´elisation HM instationnaire”, Technical Report, EDF R&D, Clamart, France, (2006).
[3] R. E. Showalter, “Diffusion in poro-elastic media”, J. Math. Anal. Appl., 251, 310– 340, (2000).
[4] C. Johnson, Y.- Y. Nie and V. Thom´ee, “An a posteriori error estimate and adaptive timestep control for a backward Euler discretization of a parabolic problem”, SIAM J. Numer. Anal., 27, 277–291, (1990).
[5] R.H. Nochetto, G. Savar´e and C. Verdi, “A posteriori estimates for variable time-step discretizations of nonlinear evolution equations”, Comm. Pure Appl. Math., 53, 525–589, (2000).
[7] R. Verf¨urth, “Robust a posteriori error estimates for non-stationary convection-diffusion equations”, SIAM J. Numer. Anal., 43, 1783–1802, (2005).
