Numer. Math. (2002) 91: 649–673 Digital Object Identifier (DOI) 10.1007/s002110100347
Numerische Mathematik
From compressible to incompressible materials via an asymptotic expansion Philippe Guillaume, Mohamed Masmoudi, Anis Zeglaoui UMR MIP 5640, D´epartement de Math´ematiques, INSA, Complexe Scientifique de Rangueil, 31077 Toulouse Cedex, France; e-mail: {guillaum,zeglaoui}@gmm.insa-tlse.fr,
[email protected] Received April 21, 2000 / Revised version received February 28, 2001 / c Springer-Verlag 2001 Published online October 17, 2001 –
Summary. In linear elasticity problems, the pressure is usually introduced for computing the incompressible state. In this paper is presented a technique which is based on a power series expansion of the displacement with respect to the inverse of Lam´e’s coefficient λ. It does not require to introduce the pressure as an auxiliary unknown. Moreover, low degree finite elements can be used. The same technique can be applied to Stokes or Navier-Stokes equations, and can be extended to more general parameterized partial differential equations. Discretization error and convergence are analyzed and illustrated by some numerical results. Mathematics Subject Classification (1991): 74S05, 74B05, 35Q72, 35Q30, 30E15
1. Introduction In solid and structural Mechanics, the displacement field u is the solution to the equilibrium equation − div σ(u) = f in Ω, u = 0 on ΓD , (1) σ(u)n = g on ΓN , where f and g are respectively volume and boundary forces and n is the unit outward normal to the boundary Γ = ΓD ∪ ΓN of a domain Ω. In the linear Correspondence to: Ph. Guillaume
650
Ph. Guillaume et al.
homogeneous isotropic elasticity model, the stress tensor σ(u) is related to the strain tensor ε(u) = (Du + DuT )/2 by σ(u) = 2µε(u) + λdiv uI, where λ and µ > 0 are the Lam´e’s coefficients and I denotes the identity matrix. Then, Eq. (1) can also be written −2µ div ε(u) − λ∇div u = f in Ω, u = 0 on ΓD , (2) (2µ ε(u) + λdiv u I)n = g on ΓN . Coefficient λ is related to Poisson’s ratio ν by 2µν 1 − 2ν and it is not defined for ν = 0.5, which corresponds to the case of an incompressible material. For ν close to 0.5, it is well-known that locking may appear when solving Eq. (2) with a displacement-based finite element method, particularly when low degree finite elements are used [2, 3]. One way for overcoming this difficulty is to use a mixed finite element method involving an additional unknown, the pressure [7, 18, 3]. The drawback is that the matrix resulting from the discretization of the problem becomes larger and is not anymore positive definite, making the linear system more difficult to solve numerically. Moreover, the finite element spaces for the displacement and the pressure must be carefully chosen in order to satisfy the Babuˇska-Brezzi inf-sup condition, otherwise locking or instability may still occur [7]. Another possibility, which does not necessarily need introducing pressure, is to use p or hp finite element methods [19, 20], which consists in increasing the degree of approximation polynomials. They are not easy to implement, and the error analysis requires the solution to be very regular. Methods which do not require such a regularity are the non conforming finite elements methods, see [9, 12, 5] and [8] where the reader can find an extended review on locking phenomena. Interesting and practical also are the popular reduced integration techniques [21], although their error analysis seems to be incomplete. Another method which was recently introduced by engineers [1] using parameterization techniques [15], consists in the following. With the change of variable α = 1/λ, problem (2) reads −2αµ div ε(u) − ∇div u = αf in Ω, u = 0 on ΓD , P (α) (3) (2αµ ε(u) + div u I)n = α g on ΓN . λ =
Naturally, the difficulty for solving P (α) with α close to 0 is the same than for solving Eq. (2) with λ close to infinity, there is no uniform coercivity of
From compressible to incompressible materials via an asymptotic expansion
651
the operator with respect to α. However, the solution u(α) to Eq. (3) is not singular at α = 0 : it will be shown in Sect. 2 that u(α) is holomorphic with respect to α on an open subset of the complex plane, containing the half space {z ∈ C; real(z) > −1/L}, where L > 0 is a stability constant which appears in the mixed formulation associated to the incompressible case ν = 0.5 (see Lemma 2.1). Besides, it means that Eq. (3) can also be defined for a negative bulk modulus κ = E/(3 − 6ν) corresponding to a Poisson’s ratio ν = 1/(2µα + 2) with α > −1/L. It follows that for a fixed α0 > 0 (see Sect. 4 for the choice of α0 ), the displacement u(α) can be expanded into a power series around α0 (see also [10] for an expansion around α0 = 0), (4)
u(α) =
∞
un (α − α0 )n ,
n=0
whose radius of convergence ρu satisfies ρu ≥ α0 + 1/L. The coefficient u0 is the solution to P0 (α0 ) := P (α0 ) (3), and for n ≥ 1, each coefficient un is the solution to −2α0 µ div ε(un ) − ∇div un = 2µ div ε(un−1 ) + δn f in Ω, on ΓD , Pn (α0 ) un = 0 (2α0 µ ε(un ) + div un I)n = − 2µ ε(un−1 )n + δn g on ΓN , where δ1 = 1 and δn = 0 for n ≥ 2. These systems are obtained by differentiating Eq. (3) with respect to α up to the order n (see also Eq. (17)). The key point is that, for any given α0 > 0 such that Pn (α0 ) can be numerically solved without difficulties, series (4) converges at the point α = 0 (i.e., for ν = 0.5), and u(0) solves the equations associated to an incompressible material, p being the pressure: −2µ div ε(u) + ∇p = f in Ω, div u = 0 in Ω, Q(0) u = 0 on ΓD , (2µ ε(u) − pI)n = g on ΓN . Hence, instead of computing directly the solution u(0) to problem Q(0), one can start by evaluating the first coefficients un in Eq. (4) by solving the associated systems Pn (α0 ), and next obtain an approximation of u(0) by using a Taylor expansion. An important property shared by systems Pn (α0 ), n ≥ 0, is that they all involve the same operator, only the right-hand sides are modified. Thus,
652
Ph. Guillaume et al.
when solving the corresponding discretized problem, a single factorization of the associated matrix is needed for computing the finite element approximation uh,n of un . This technique allows to solve the incompressible (or nearly incompressible) problem without involving the pressure, and for a computational cost comparable to the cost of the solution of the compressible problem. The pressure can then be recovered from the displacement field by using also a power series expansion. Another advantage of this technique is that the resulting matrix remains symmetric positive definite, and smaller than in the case of a mixed finite element method. Moreover, low degree finite elements can be used as for example linear elements, and do not produce any locking phenomenon when α0 is chosen reasonably far from 0, that is, when ν is not too close from 0.5. The method is very easy to implement in a numerical program which is designed only for compressible material. It can also be applied to the solution of Stokes equations, and could be extended to Navier-Stokes equations. In Sect. 3 is studied the discretization error on each coefficient of power series (4), together with the convergence of a particular sequence of approximate solutions: it is shown that for a given α0 > 0, there exist some integers N (h) such that the error on the finite element approximation Taylor expansion defined by N (h)
(5)
Uh,N (h) :=
uh,n (−α0 )n ,
n=0
satisfies (6)
u(0) − Uh,N (h) ≤ chξk , 1
where h is the usual notation for the mesh size, k is the degree of the finite element, c > 0 and 0 < ξ < 1 both depend on α0 (but are independent on h). It is noticeable that N (h) increases slowly: we will see that N (h) behaves like log(1/h) when h → 0. The usual error bound in finite element problems is in O(hk ). Error bound (6) is weaker because ξ < 1. The reason of this loss is that the discrete problem does not behave exactly like the continuous problem: although the finite element approximation uh (α) of u(α) may be holomorphic on some half space {z ∈ C; real(z) > −1/Lh }, the bounds 1/Lh can have an accumulation point at 0 when h → 0 (if not, then one could obtain a uniform discrete inf-sup condition with respect to ν close to 0.5, which is usually not the case). However, we will see that ξ → 1 when α0 → 0, but unfortunately with c = O(1/α0 ) in Eq. 6. Hence the best error bound is obtained for some intermediate value of α0 (see Remark 3.2). This situation reminds the behavior of penalty methods, where the divergence
From compressible to incompressible materials via an asymptotic expansion
653
is forced to be small by adding a term ∇div u/$ (see e.g., [3, 17]) with $ being small enough, but not too small or else locking appears. Finally, some numerical results are reported in Sect. 4, which confirm the error analysis. 2. Holomorphy of the displacement Let Ω be an open and bounded subset of R3 , with boundary Γ = ΓD ∪ ΓN , ΓD ∩ ΓN = ∅. It is supposed throughout this paper that both ΓD and ΓN have a nonnegative Lebesgue measure. The case of a pure homogeneous Dirichlet boundary condition (ΓD = Γ ) can be treated similarly, with the difference that a condition Ω p dx = 0 must be prescribed in Eq. (10). Functional spaces V and W are defined by V = {u ∈ H 1 (Ω)3 ; u = 0 on ΓD }, W = V × L2 (Ω), and their dual spaces are denoted respectively V and W . For notational simplicity we do not mention C in the involved spaces, for example L2 (Ω, C) is simply written L2 (Ω). We need to consider also the fractional Sobolev space 1/2
H00 (ΓN )3 = {g ∈ H 1/2 (Γ )3 ;
g = 0 in Γ \ ΓN }.
1/2
The space H00 (ΓN )3 (the restrictions to ΓN of its elements) is a strict subspace of H 1/2 (ΓN )3 , and consists in functions who vanish at the boundary 1/2 of ΓN in some special sense (cf. [14, 16]). The dual space of H00 (ΓN )3 is −1/2 denoted H00 (ΓN )3 . −1/2 For a complex number α = / 0, f ∈ L2 (Ω)3 and g ∈ H00 (ΓN )3 , the variational formulation of problem (3) is the following: find u(α) ∈ V such that (7)
a(α, u(α), v) = αl(v),
∀v ∈ V ,
with (8)
a(α, u, v) = αa1 (u, v) + a2 (u, v) , a1 (u, v) = 2µ ε (u) : ε (v) dx , Ω a2 (u, v) = div u div v dx , Ω f.v dx + g.v dγ(x). l(v) = Ω
ΓN
654
Ph. Guillaume et al.
Here x.y denotes the usual dot product xi yi , ε (u) : ε (v) = εij (u) εij (v) i
i,j
and dγ(x) is the Lebesgue measure on the boundary. It follows from meas(ΓD ) > 0, µ > 0 and Korn’s inequality that a1 (u, u)1/2 is a norm on the space V which is equivalent to the usual norm induced by H 1 (Ω)3 . This allows to misuse the notation ·1 by writing u1 := a1 (u, u)1/2 . The space W is equipped with the inner product (and the associated norm) defined by (u, p), (v, q)W := a1 (u, v) + pq dx. Ω
For real s = / 1, the usual norm on H s (Ω)3 is denoted ·s , and similarly the −1/2 1/2 norms on H00 (ΓN )3 , H00 (ΓN )3 or H0s (ΓN )3 for s > 1/2 are denoted · s,ΓN with s = −1/2, 1/2 or s > 1/2. When solving Eq. (3) for ν close to 0.5, it is usual to introduce the unknown pressure given by p = −κ div u. It is here more convenient to consider a “pseudo” pressure p(α), homogeneous to the pressure and defined by 3 1 p(α) = − div u(α) = (9) (−κdiv u(α)). α 3 + 4µα Although there is here a slight misuse of notation, observe that for α = 0 the usual pressure is recovered. With this change of variable, problem ( 3) can be written −2µ div ε(u) + ∇p = f in Ω, div u + α p = 0 in Ω, (10) Q(α) u = 0 on ΓD , σ(u, p)n = g on ΓN , with σ(u, p) = 2µε(u) − pI. Clearly, problems (2), P (α) and Q(α) are equivalent for ν = / 0.5 (i.e., α = / 0). The next lemma is given without proof, and follows from the standard inf-sup condition [13] if α = 0 and from the coercivity of a(α, ·, ·) if α 0. The set of negative real numbers is denoted R∗− . Lemma 2.1. For α ∈ C\R∗− and b ∈ W , the problem: find (u, p) ∈ W such that 2µε (u) : ε (v) − p div v + div u q + αpq dx = b(v, q) ∀(v, q) ∈ W, Ω
(11)
From compressible to incompressible materials via an asymptotic expansion
655
has a unique solution (u(α), p(α)) ∈ W, and there exists L(α) > 0 such that (u(α), p(α))W ≤ L(α)bW . The main result of this section is the following. Theorem 2.2. Let L := L(0). The functions u(α), p(α) defined by Eq. (11) are holomorphic on C\R∗− , and can be extended into holomorphic functions on the half space {z ∈ C; real(z) > −1/L}. Proof. Let L(W) denote the Banach algebra whose elements are the endomorphisms on W. It follows from Riesz’s theorem that for all complex number α, there exists A(α) ∈ L(W) and F ∈ W such that A(α)(u, p), (v, q)W = 2µε (u) : ε (v) − p div v + div u q + αpq dx, Ω
b(v, q) = F, (v, q)W
∀(v, q) ∈ W, ∀(v, q) ∈ W.
One can observe that function A, A : C → L(W) α → A(α) is affine and continuous, thus holomorphic. Let Is(W) denote the set of automorphisms of L(W). The set Is(W) is open in L(W), thus O := A−1 (Is(W)) is open in C. With this notation, Eq. (11) reads A(α)(u, p) = F. Hence we have (u(α), p(α)) = (A(α))−1 F, and, by composition of holomorphic functions, α → u(α) and α → p(α) are holomorphic on O. We know from Lemma 2.1 that O contains C\R∗− . In particular, the open set O contains the origin, and u(α) and p(α) can be expanded into a power series at the origin: (12) (13)
u(α) = p(α) =
∞ n=0 ∞ n=0
un α n , pn α n ,
656
Ph. Guillaume et al.
where (u0 , p0 ) = (u(0), p(0)) is the solution to Eq. (11) with α = 0, and for n ≥ 1, (un , pn ) is the solution to the following problem: 2µε (un ) : ε (v) − pn div v + div un q dx Ω pn−1 q dx ∀(v, q) ∈ W. =− Ω
This equation is obtained by substituting (12)-(13) for u, p, in Eq. (11) and by identifying the coefficients of αn in each term. It follows from Lemma 2.1 that (u0 , p0 )W ≤ LbW , (un , pn )W ≤ Lpn−1 0 ≤ L (un−1 , pn−1 )W ,
∀n ≥ 1,
which by induction yields (un , pn )W ≤ Ln+1 bW ,
∀n ≥ 0.
This implies that the convergence radius of series (12) and (13) are bounded below by 1/L. Thus O contains C\R∗− ∪ B (0, 1/L), and in particular O contains the half space {z = x + iy ∈ C; real(z) > −1/L}.
A straightforward consequence of Theorem 2.2 follows. Corollary 2.3. For α > −1/L and b(v, q) = l(v) in Eq. (11), u(α) and p(α) defined by Eq. (11) are solution to problem Q(α). For a given α0 > −1/L, u(α) and p(α) can be expanded into a power series around α0 : (14) (15)
u(α) = p(α) =
∞ n=0 ∞
un (α0 )(α − α0 )n , pn (α0 )(α − α0 )n ,
n=0
with un (α0 ) ∈ V and pn (α0 ) ∈ L2 (Ω) for all n ≥ 0. Let ρu and ρp denote respectively the radius of convergence (which depend on α0 ) of these two series. Then 1 1 ρu ≥ α0 + , ρp ≥ α0 + . (16) L L Particularly, these two series converge at α = 0 to the solution to problem Q(0).
From compressible to incompressible materials via an asymptotic expansion
657
ρ (α0) u
1/L
0
α0
Fig. 1. The convergence radius ρu
For α0 > 0, the pressure p(α) and its coefficients pn (α0 ) are not necessary to recover the coefficients un (α0 ). The latter can be obtained by substituting (14) into (7), which yields ∞
n (α − α0 + α0 )a1 un (α0 )(α − α0 ) , v + a2
n=0 ∞
n
un (α0 )(α − α0 ) , v
n=0
= (α − α0 + α0 )l(v), ∀v ∈ V, and by identification of the coefficient of (α − α0 )n , it leads to the following problems Pn (α0 ) : a(α0 , u0 (α0 ), v) = α0 l(v), ∀v ∈ V , a(α0 , u1 (α0 ), v) = −a1 (u0 (α0 ), v) + l(v), ∀v ∈ V , (17) Pn (α0 ) a(α0 , un (α0 ), v) = −a1 (un−1 (α0 ), v), ∀v ∈ V, ∀n ≥ 2. Observe that all these systems involve the same operator, only the right-hand side are modified. 3. Finite element approximation Let Th be a mesh of Ω in the general meaning of usual approximations by a finite element method [10], and let Vh ⊂ V be a finite-element space of degree k ≥ 1. As usual for simplicity, we do not take into account the consistency error coming from the approximation of the geometry. However, this aspect could be treated by usual techniques and Bernardi’s general
658
Ph. Guillaume et al.
results [4]. It is assumed that the family (Vh )h>0 is regular [10], that is, there exists C > 0 such that (18)
v − Πh v1 ≤ Chk vk+1
∀v ∈ V ,
where Πh denotes the usual Lagrange interpolation operator. The finite element approximation of u(α) in Eq. (7) is: find uh (α) ∈ Vh such that (19)
a(α, uh (α), vh ) = αl(vh ),
∀vh ∈ V .
Similarly, the approximation of un (α0 ) in Eq. (17) is: find uh,n (α0 ) ∈ Vh such that a(α0 , uh,0 (α0 ), vh ) = α0 l(vh ), ∀vh ∈ Vh , a(α0 , uh,1 (α0 ), vh ) = −a1 (uh,0 (α0 ), vh ) + l(vh ), ∀vh ∈ Vh , Ph,n (α0 ) a(α0 , uh,n (α0 ), vh ) = −a1 (uh,n−1 (α0 ), vh ), ∀vh ∈ Vh , ∀n ≥ 2. (20) The functions uh,n (α0 ) are nothing else than the coefficients of the power series expansion of uh (α) with respect to α at α0 . Although the functions uh (α) are clearly holomorphic on the right half plane {z ∈ C; real(z) > 0}, nothing can be said on a possible (uniform in h) extension to a part of the left half plane, because the arguments used in the continuous case, based on the well-posedness of the incompressible problem, cannot be applied to the discrete case. Henceforth, the functions un (α0 ) and uh,n (α0 ) in Eq. (17) and (20) are respectively denoted un and uh,n , and the solution u(0) to the incompressible problem Q(0) is simply denoted u. Unlike u(α) which is holomorphic around α = 0, uh (α) may be singular at (or near) α = 0. Hence, the series ∞ n n=0 uh,n (α − α0 ) may not (or hardly) converge at the point α = 0, whereas series (4) was convergent. However, we will see that a good approximation can still be recovered from this series. For a given N > 0, we aim to compare the truncated series Uh,N :=
N
uh,n (−α0 )n ,
n=0
with the exact solution u. First we estimate the error un − uh,n 1 . Next, for a particular choice N = N (h), we show that the error between the approximate solution Uh,N (h) and the exact solution u behaves like O(hξk ) with 0 < ξ < 1.
From compressible to incompressible materials via an asymptotic expansion
659
3.1. Error estimate on the coefficient of the series Let αmax > 0 be a fixed number, and w ∈ V. For α ∈]0, αmax ] consider the solution y(α) to (21)
a(α, y(α), v) = αa1 (w, v) + αl(v),
∀v ∈ V .
The function v → a1 (w, v) can be considered an element of V , and a direct consequence of Theorem 2.2 is that y(α) can be extended into an holomorphic function on {z ∈ C; real(z) > −1/L}, still denoted y(α). Particularly, for (22)
0 < η < 1/L,
the function (α, f, g, w) → y(α)1 is continuous on −1/2
{z ∈ C; real(z) ≥ η − 1/L} × V × H00
(ΓN ) × V.
Thus there exists R0 > 0 such that the following uniform a priori estimate holds: (23)
y(α)1 ≤ R0 (f V + g−1/2,ΓN + w1 ), 1 ∀α ∈ B(αmax , αmax + − η), L
where as usual B(z, r) denotes the closed ball centered at z with radius r. For the error analysis, we need to consider the following stronger regularity assumption which depends on Ω, ΓN and ΓD : there exists a number R > 0 such that (24)
y(α)k+1 ≤ R(f k−1,Ω + gk−1/2,ΓN + wk+1 ), 1 ∀α ∈ B(αmax , αmax + − η). L
It is well-known that in the case of mixed Dirichlet-Neumann boundary conditions, assumption (24) is not fulfilled for arbitrary Ω, ΓN and ΓD . An example for which this assumption is satisfied is when ΓN and ΓD are two regular and separate (i.e., dist(ΓN , ΓD ) > 0) components of the boundary Γ of Ω. In the general case, this assumption could be weakened by replacing the integer k ≥ 1 by a real number s > 0, involving fractional Sobolev spaces for which the analysis becomes more complicate. Throughout this section, the number α0 is chosen in the interval ]0, αmax ] and η defined by (22) is fixed once for all.
660
Ph. Guillaume et al.
Lemma 3.1. If assumption (24) is satisfied, then n R un k+1 ≤ (1 + R) (f k−1,Ω + gk−1/2,ΓN ), α0
∀n ≥ 0,
where un solves Eq. (17). Moreover, one has necessarily R ≥ α0 /ρu , where ρu was defined in Eq. (16). Proof. We proceed by induction on n. The term in (17) corresponding to w in (21) is 0 for n = 0 and −u0 /α0 for n = 1, and l is substituted for α0 l for n = 1. It follows from (24) that u0 k+1 ≤ R(f k−1,Ω + gk−1/2,ΓN ), R (f k−1,Ω + gk−1/2,ΓN + u0 k+1 ) u1 k+1 ≤ α0 R (1 + R)(f k−1,Ω + gk−1/2,ΓN ). ≤ α0 Hence the result is true for n = 0 and n = 1. Let the result be true for n ≥ 1. The term in (17) corresponding to w in (21) is −un−1 /α0 . It follows from (24) that R un k+1 α0 n+1 R ≤ (1 + R) (f k−1,Ω + gk−1/2,ΓN ). α0 This proves that the radius of convergence of series un z n is bounded k+1 k+1 1 below in H (Ω) by α0 /R. As H (Ω) ⊂ H (Ω), it is bounded above by the radius of convergence ρu of the same series in H 1 (Ω). Thus un+1 k+1 ≤
α0 /R ≤ ρu .
The bilinear form a(α, ·, ·) is continuous (see Eq. (8) and recall u1 = a1 (u, u)1/2 ), with
1 u1 v1 , ∀u, v ∈ V, ∀α > 0, (25) |a(α, u, v)| ≤ α + 2µ and α-coercive, (26)
a(α, u, u) ≥ αu21 ,
∀u ∈ V,
∀α > 0.
Combining C´ea’s lemma, Strangs’s lemma and these two bounds yields the following error estimate.
From compressible to incompressible materials via an asymptotic expansion
661
Lemma 3.2. If assumption (24) is satisfied, then the discretization error on the coefficients is bounded by
n 1 n−m+1 1 un − uh,n 1 ≤ (α0 + um k+1 , ∀n ≥ 0, )C hk 2µ α0 m=0
where C was defined in Eq. (18). Proof. We proceed by induction on n. The equations to compare are (17 ) and (20), which define respectively un and uh,n . Let M = α0 + 1/(2µ). For n = 0 in Eq. (20), it follows from C´ea’s lemma [10] with (25)-(26) and (18) that M inf u0 − vh 1 α0 vh ∈Vh MC k ≤ h u0 k+1 . α0
u0 − uh,0 1 ≤
Let the result be true for a given n − 1 ≥ 0. For n ≥ 1 in Eq. (20), it follows from Strang’s lemma [10] with (25)-(26) and (18) that M 1 inf un − vh 1 + a1 (un−1 , ·) α0 vh ∈Vh α0 −a1 (uh,n−1 , ·)V MC k 1 ≤ h un k+1 + un−1 − uh,n−1 1 . α0 α0
un − uh,n 1 ≤
Hence, using the induction hypothesis we obtain n−1 1 n−m MC k 1 k un − uh,n 1 ≤ h un k+1 + M Ch um k+1 α0 α0 α0 m=0
n−m+1 n 1 = M Chk um k+1 , ∀n ≥ 0. α0 m=0
Lemmas 3.1 and 3.2 are summarized in the next proposition. Proposition 3.3. If assumption (24) is satisfied, then there exist constants c > 0 and β ≥ 1 (independent on α0 and h) such that
n+1 β 1 (27) un − uh,n 1 ≤ c α0 + 2µ α0 ×(f k−1,Ω + gk−1/2,ΓN ) hk
∀n ≥ 0.
Moreover, (27) holds with β = max(1, R) if R = / 1 and β = (1 + τ ) if R = 1, with τ > 0 arbitrarily close to 0.
662
Ph. Guillaume et al.
Proof. Let M = α0 + 1/(2µ). It follows from Lemmas 3.2 and 3.1 that
n 1 n−m+1 un − uh,n 1 ≤ M Ch um k+1 α0 m=0 n+1 1 k ≤ M (1 + R)Ch α0 n × Rm (f k−1,Ω + gk−1/2,ΓN ) k
m=0 n+1 1 k M (1 + R)Ch α0 n+1 × 1−R (f k−1,Ω + gk−1/2,ΓN ) 1−R = n+1 1 k M (1 + R)Ch α0 ×(n + 1)(f k−1,Ω + gk−1/2,ΓN )
if
R= / 1,
if
R = 1.
Remark 3.1. The “best” situation occurs when β = 1, that is, when R < 1. However, we know from Lemma 3.1 that necessarily R ≥ α0 /ρu , and α0 /ρu may stand quite close to 1, especially when α0 = (1 − 2ν)/(2µν) is large, that is, when the corresponding Poisson’s ratio ν is chosen close to 0. Hence it is likely to happen that R ≥ 1, in which case β > 1. 3.2. Error estimate on the approximate solution The main result of this section is the following. Theorem 3.4. Let u be the solution to the incompressible problem Q(0) (10). Let α0 be a given positive number, and for n ≥ 0, let uh,n be the solution to problem Ph,n (α0 ) (20). If assumption (24) is satisfied, then for all $ > 0, there exist N ≥ 0 and h0 > 0 such that for all h ≤ h0 , one has N (28) uh,n (−α0 )n ≤ $. u − n=0
1
Moreover, for 0 < η < 1/L, define (29)
ρ = α0 + 1/L − η, 0 0, α0 log(ρ/α0 )
From compressible to incompressible materials via an asymptotic expansion
663
If β > 1, there exists a constant c (independent of α0 , N and h) such that for − log(γhk log β) N (h) = −1 , log(ρβ/α0 ) where [x] = round(x), one has N (h) n (30) u − uh,n (−α0 ) ≤ c log(ρβ/α0 ) n=0 1
×γ ξ (f V + g−1/2,ΓN ) hξk .
If β = 1, there exists a constant c (independent of α0 , N and h) such that for − log(γhk ) −1 , N (h) = log(ρ/α0 ) one has N (h) n u − uh,n (−α0 ) ≤ c (f V + g−1/2,ΓN ) γhk log(γhk ). n=0 1
Proof. Let N ≥ 0. We consider only the case β > 1, the case β = 1 can be treated in the same way by substituting N + 1 for β N +1 in Eq. ( 33). We have N N (31) u − uh,n (−α0 )n ≤ u − un (−α0 )n n=0 n=0 1 1 N n + (un − uh,n )(−α0 ) . n=0
1
We bound separately each term of the right-hand side. Let Γ+ denote the standard oriented circle, centered at α0 with radius ρ. Due to ρ < α0 + 1/L, we have Γ+ ⊂ {z ∈ C; real(z) > −1/L} where u(α) is well-defined, and Cauchy’s formula yields u(α) 1 dα. un = 2iπ Γ+ (α − α0 )n+1 It follows from Γ+ ⊂ B(αmax , αmax + un 1 ≤
1 L
− η) and (23) with w = 0 that
R0 (f V + g−1/2,ΓN ). ρn
664
Ph. Guillaume et al.
Hence, using 0 < α0 /ρ < 1 we obtain N un (−α0 )n u − n=0 1 ∞ n = un (−α0 ) n=N +1
(32)
1
n ∞ α0 ≤ R0 (f V + g−1/2,ΓN ) ρ n=N +1 N +1 ρ α0 = R0 (f V + g−1/2,ΓN ) ρ − α0 ρ
N +1 Lα0 α0 (f V + g−1/2,ΓN ) = R0 1 + 1 − Lη ρ N +1 α0 ≤ c1 (f V + g−1/2,ΓN ) , ρ
with c1 = R0 (1 + Lαmax /(1 − Lη)). For the second term in right hand side of (31), it follows from Proposition 3.3 with M = α0 + 1/(2µ) and assumption β > 1 that N n (u − u )(−α ) n 0 h,n n=0 1
N β n+1 n k ≤ M ch (f k−1,Ω + gk−1/2,ΓN ) α0 α0 n=0
= M chk (f k−1,Ω + gk−1/2,ΓN )
β β N +1 − 1 . α0 β − 1
c2 k h (f k−1,Ω + gk−1/2,ΓN )β N +1 , α0 with c2 = cβ/(β − 1). Hence, gathering (31), (32) and (33) yields with c3 = max(c1 , c2 ): N n uh,n (−α0 ) u − n=0 1 α0 N +1 hk N +1 . ≤ c3 (f V + g−1/2,ΓN ) +M β ρ α0 (33)
≤M
Since 0 < α0 /ρ < 1, the first part of the theorem is obtained by choosing first N, and then h. The best error bound is obtained by considering the minimum of the function
From compressible to incompressible materials via an asymptotic expansion
j(x) =
α0 ρ
x
+M
665
hk x β , α0
which is given by 1 log x = log(ρβ/α0 ) ∗
j(x∗ ) =
1 k γh log β
log(ρβ/α0 ) (γ log β)ξ hkξ , log β
,
γ=
0 0, α0 log(ρ/α0 ) log(ρ/α0 ) < 1. log(ρβ/α0 )
It follows that the optimal N, for which the error bound is minimum, is given by : − log(hk γ log β) −1 . N= log(ρβ/α0 )
The term (log β)ξ / log β which remains bounded with respect to α0 (see Remark below) has been put into constant c in Eq. (30). Remark 3.2. One can observe that ξ given by (29) satisfies (recall ρ = α0 + 1/L − η) lim ξ = lim
α0 →0
α0 →0
log(ρ/α0 ) = 1, log(ρβ/α0 )
hence the smaller α0 , the better the behavior of hξk with respect to h. However, the constant which appears in front of h in (30) is (using α0−ξ (α0 β)−1 for small α0 ): c log(ρβ/α0 )
(α0 + 1/(2µ)) α0 log(ρ/α0 )
ξ
c . 2µβα0
Thus what is gained on one side is loosed on the other side, and in fact, the best bound is obtained for some intermediate values of α0 corresponding to ν0 somewhere around 0.3 as illustrated by Fig. 2 (recall α = (1 − 2ν)/(2µν)). This figure plots some curves ν0 −→ x∗ and ν0 −→ j(x∗ ) for several values of h. The constants which are involved in Theorem 3.4 were chosen as follows: k = 1, µ = 0.5, L = 1, η = 0.1, β = 1.2, c = 1, h = 2−1 /10, 2−2 /10, 2−3 /10, 2−4 /10, 2−5 /10. One can observe that the optimal ν0 is not very sensitive to the mesh size h.
666
Ph. Guillaume et al.
Fig. 2. Optimal degree of Taylor expansion and optimal error bounds for different h’s
3.3. Numerical implementation In this section is presented the algorithm for computing the displacement and the stress tensor. Let (ϕj )1≤j≤d be a basis of space Vh . The functions ϕj vanish on ΓD , that is, we are considering what is usually called the reduced problem. In this basis, the displacement uh (α) has the form uh (α) =
d
xi (α)ϕi .
i=1
The power series expansion of xi (α) reads i
x (α) =
∞ n=0
xin (α − α0 )n .
With this notation, we have uh,n =
d
xin ϕi ,
∀n ≥ 0.
i=1
The vector (xin )1≤i≤d is denoted xn , and (xi (α))1≤i≤d is denoted x(α). The d × d matrices A and B and vector b are defined by Aij = a1 (ϕj , ϕi ), Bij = a2 (ϕj , ϕi ), bi = l(ϕi ). Thus, Eq. (19) reads (αA + B)x(α) = αb,
From compressible to incompressible materials via an asymptotic expansion
667
and Eq. (20) are equivalent to the following systems: (α0 A + B)x0 = α0 b, (α0 A + B)x1 = −Ax0 + b, (34) Sn (α0 ) (α0 A + B)xn = −Axn−1 , ∀n ≥ 2. Then, the displacement u = u(0) at the incompressible state is approximated by the following truncated series: Uh,N :=
N
uh,n (−α0 )n .
n=0
In the analysis of solids, it is also important to determine the stress tensor. Clearly it cannot be evaluated from σ(u) = 2µε(u) + div u I/α by putting α = 0. It must be evaluated from σ(u) = 2µε(u) − pI, that is, although not necessary for computing u, the pressure p must be computed for retrieving the stress. It is obtained from u in the following way at the post-processing stage. Likewise the displacement, the “pressure” p(α) = −
div u(α) α
can be expanded into a power series around α0 (recall that p(0) coincides with the exact pressure). The power series expansion of div u(α) around α0 is ∞ (35) div un (α − α0 )n . div u(α) = n=0
The power series expansion of −1/α around α0 is ∞
1 − = α
n=0
−1 α0
n+1
(α − α0 )n .
Thus, the power series expansion of p(α) around α0 is p(α) =
∞
pn (α − α0 )n ,
n=0
pn =
−1 p+1 div uq . α 0 p+q=n
The resulting algorithm is the following.
668
Ph. Guillaume et al.
Algorithm 1. Choose the order of approximation N. 2. Solve (34) for n = 1, 2, . . . , N. 3. Compute the approximate displacement Uh,N =
N d n=0 i=1
(−α0 )n xin ϕi .
4. Compute the approximate strain tensor 1 T ). ε(Uh,N ) = (DUh,N + DUh,N 2 5. Compute the coefficients of series (35) for n = 1, 2, . . . , N, div uh,n =
d
xin div ϕi .
i=1
6. Compute the approximate pressure ph =
N n=0
n
ph,n (−α0 ) , ph,n
−1 p+1 = div uh,q . α0 p+q=n
7. Compute the approximate stress tensor σ(Uh,N ) = 2µε(Uh,N ) − ph I. The choice of N is not straightforward. In complement to the informations given by Theorem 3.4 and Remark 3.2, a guideline for chosing N is the following. Denote respectively bynρ and ρh the radius ∞of convergence of n the series u(α) = ∞ u (α − α ) and u (α) = 0 h n=0 n n=0 uh,n (α − α0 ) . Hadamard’s formula yields (36)
lim sup un 1/n = 1/ρ, lim sup uh,n 1/n = 1/ρh . n→∞
n→∞
We consider here the most simple situation where un = cρ−n , which is sufficient for illustrating the guideline. Due to limh→0 uh,n − un = 0, we have for small n uh,n 1/n un 1/n = c1/n /ρ, instead of (36) for large n. The place where this change of rate (ρ to ρh ) occurs in the decay of the coefficients uh,n seems to be a reasonable choice for N.
From compressible to incompressible materials via an asymptotic expansion
0
oh=2 +h=1 * h = 0.5
Poisson coefficient = 0.5 ν0 = 0.1
0
−1
−1
−1.5
−1.5
−2
−2 5
0
oh=2 +h=1 * h = 0.5
10
15
number of terms
20
25
5
Poisson coefficient = 0.5 ν0 = 0.3
0
oh=2 +h=1 * h = 0.5
10
15
number of terms
20
25
Poisson coefficient = 0.5 ν0 = 0.35
−0.5
log10(error)
−0.5
log10(error)
Poisson coefficient = 0.5 ν0 = 0.2
−0.5
log10(error)
log10(error)
−0.5
oh=2 +h=1 * h = 0.5
669
−1
−1.5
−1
−1.5
−2
−2 5
10
15
number of terms
20
25
5
10
15
number of terms
20
25
Fig. 3. The relative error log10 (uh − uref 0 / uref 0 ) with respect to N
4. Numerical results We consider here the example of a 2D cantilever beam with a Dirichlet boundary condition along the left side and loaded along the right side. The domain is a rectangle Ω = [0, L] × [−c, c] and the Neumann condition is g(x, c) = g(x, −c) = (0, 0)T , x ∈ [0, L], y ∈ [−c, c]. g(L, y) = (0, P )T , The discrete problem has been solved with MATLAB, using a triangular mesh and the linear finite element P1. The exact solution of the incompressible state is not known, however a reference displacement uref and pressure pref have been obtained from ABAQUS with a fine mesh (h = 0.1) and a quadrilateral element Q1 . The data used in the computation are P = −0.001, L = 16, c = 2 and E = 1, where E is Young’s modulus. In Figs. 3, 4, 5 and 6, the displacement and the pressure at the incompressible state (Poisson coefficient ν = 0.5) are approximated respectively
670 0
Ph. Guillaume et al. ν0 ν0 ν0 * ν0 − o +
Poisson coefficient = 0.5
= 0.35 = 0.3 = 0.3 = 0.1
0
ν0 ν0 ν0 * ν0
= 0.35 = 0.3 = 0.3 = 0.1
2
4
Poisson coefficient = 0.5
−0.5 log10(error)
log10(error)
−0.5
− o +
−1
−1
−1.5
−1.5
−2
−2 2
4
6
8
10
12
14
16
18
20
h=1
6
8
10
12
14
16
18
20
h = 0.5
Fig. 4. The relative error log10 (uh − uref 0 / uref 0 ) with respect to N
by (37)
N n=0
uh,n (α − α0 )n ,
N
ph,n (α − α0 )n
n=0
with α = 0. The computation are performed for different values of N, h and α0 , and are compared to the displacement uref and the pressure pref at the incompressible state computed by ABAQUS. The value ν0 = α0 /(2µ+2α) corresponding to the chosen α0 where the Taylor expansion is computed is indicated on each figure. Figure 3 shows the relative error on the displacement with respect to the order of approximation N . Each plot concerns a particular value of ν0 and 3 values of h. In Fig. 4, ν0 and h are exchanged. Figures 5 and 6 show the relative errors on the pressure. One can observe that a good precision is obtained for a small number N of terms in the Taylor expansion, and that the optimal N increases approximately linearly with respect to log 1/h, which agrees with Theorem 3.4. Also, the smaller is ν0 , the larger is this optimal N and the less sensitive to the choice of N is the error. The numerically observed values of ξ are approximately 2/3 for ν0 = 0.35, which agrees quite well with the theoritical analysis, in spite of the loss of regularity of the solution at x = 0, y = ±c. All these tests were performed for computing the incompressible state corresponding to the Poisson coefficient ν = 0.5. Of course, the method works also for other values of ν < 0.5, in which case the problem is easier to solve. In Fig. 7 are shown the results for ν = 0.1 and ν = 0.4 when the Taylor expansion is computed at the point corresponding to ν0 = 0.3, that is by putting α = (1 − 2ν)/(2µν) and α0 = (1 − 2ν0 )/(2µν0 ) in (37).
From compressible to incompressible materials via an asymptotic expansion Poisson coefficient = 0.5 ν0 = 0.1
0
−0.2
−0.1 −0.2
−0.4
log10(error)
−0.3
−0.4
log10(error)
−0.3
−0.5
−0.5
−0.6
−0.6
−0.7
−0.8
−0.9
−0.9 2
4
6
8
10
12
14
number of terms
16
18
−1
20
Poisson coefficient = 0.5 ν0 = 0.3
0
−0.2
−0.3
−0.3
−0.4
−0.4
−0.5
4
6
8
10
12
14
number of terms
16
18
20
Poisson coefficient = 0.5 ν0 = 0.35
log10(error)
−0.1
−0.2
log10(error)
2
0
−0.1
−0.5
−0.6
−0.6
oh=1 + h = 0.5 * h = 0.35
−0.7 −0.8 −0.9 −1
oh=1 + h = 0.5 * h = 0.35
−0.7
−0.8
−1
Poisson coefficient = 0.5 ν0 = 0.2
0
oh=1 + h = 0.5 * h = 0.35
−0.1
671
oh=1 + h = 0.5 * h = 0.35
−0.7 −0.8 −0.9
2
4
6
8
10
12
14
number of terms
16
18
−1
20
2
4
6
8
10
12
14
number of terms
16
18
20
Fig. 5. The relative error log10 (ph − pref 0 / pref 0 ) with respect to N
Poisson coefficient = 0.5
= 0.35 = 0.3 = 0.3 = 0.1
0.2
0
log10(error)
0
ν0 ν0 ν0 * ν0 − o +
− o +
ν0 ν0 ν0 * ν0
= 0.35 = 0.3 = 0.3 = 0.1
2
4
Poisson coefficient = 0.5
log10(error)
0.2
−0.2
−0.2
−0.6
−0.6
−0.8
−0.8 2
4
6
8
10
12
14
number of terms
h=1
16
18
20
6
8
10
12
14
number of terms
16
h = 0.5
Fig. 6. The relative error log10 (ph − pref 0 / pref 0 ) with respect to N
18
20
672
Ph. Guillaume et al.
Poisson coefficient = 0.1 ν0 = 0.3
−0.5 −1
oh=1 + h = 0.5
−2
−1
oh=1 + h = 0.5
−2 log10(error)
log10(error)
−1.5
Poisson coefficient = 0.4 ν0 = 0.3
0
−3
−2.5
−4
−3 −5
−3.5
−6
−4 −4.5 0
2
4
6
8 10 12 14 number of terms
16
18
20
−7
0
2
4
6
8 10 12 14 number of terms
16
18
20
Fig. 7. The relative error log10 (uh − uref 0 / uref 0 ) with respect to N
It can be observed that the Taylor series of the discrete problem converges, and it is sufficiant to take a large N to obtain a good precision. For these values of ν < 0.5, it was not necessary to use ABAQUS for computing a reference solution, which was simply computed by using MATLAB. Acknowledgement. The authors are grateful to the referees for there thorough reading and there suggestions concerning the presentation.
References 1. P. Aubert, J.D. Beley, P. Guillaume, A. Zeglaoui: An asymptotic expansion for incompressibility. Eccomas 2000, Barcelone, sept. 2000, Book of abstracts, 1012 2. I. Babuˇs ka, M. Suri: Locking effects in the finite element approximation of elasticity problems. Numerische Mathematik 62, 439–463 (1992) 3. K. J. Bathe: Finite Element Procedures, Prentice-Hall, 1996 4. C. Bernardi: Optimal finite-element interpolation on curved domains. SIAM J. on Num. Anal. 26, 1212–1240 (1989) 5. S. Brenner, L. Sung: Linear finite element methods for planar elasticity. Math. of Comp. 59, 321–338 (1992) 6. H. Brezis: Analyse fonctionnelle. Masson, 1993 7. F. Brezzi, M. Fortin: Mixed and Hybrid Finite Element Methods. Springer Verlag, New York, 1991 8. D. Capatina-Papaghiuc: Contribution a` la pr´e vention de ph´e nom`e nes de verrouillage num´e rique. Th`e se de l’Universit´e de Pau, 1998 9. M. Crouzeix, P.A. Raviart: Conforming and non-conforming finite element methods for solving the stationary Stokes equations. RAIRO 7 3, 33–76 (1973) 10. P.G. Ciarlet: The finite element method for elleptic problems, North Holland, 1978 11. G. Duvaut, J.L. Lions: Les in´e quations en M´e canique et en Physique, Dunod, Paris, 1972 12. M. Fortin, M. Soulie: A non conforming piecewise quadratic finite element on triangles. Internat. J. Numer. Methods Engrg. 19, 505–520 (1983)
From compressible to incompressible materials via an asymptotic expansion
673
13. V. Girault, P.A. Raviart: Finite Element Methods for Navier-stokes Equations. Theory and Algorithms, Springer, 1981 14. P. Grisvard: Caract´e risation de quelques espaces d’interpolation. Arch. Rational Mech. Anal. 25, 40–63 (1967) 15. Ph. Guillaume and M. Masmoudi: Solution to the time-harmonic Maxwell’s equations in a waveguide, use of higher order derivatives for solving the discrete problem, SIAM J. on Num. Anal. 34(4), 1306–1330 (1997) 16. J.L. Lions and E. Magenes, Non-Homogeneous Boundary Value Problems and Applications, Springer-Verlag, Berlin, Heidelberg and New-York, 1972 17. A. Quarteroni, A. Valli: Numerical Approximation of Partial Differential Equations. Springer Verlag, 1994 18. J.E. Roberts, J.M. Thomas: Mixed and Hybrid Methods. in Handbook of Numerical Analysis, vol. II: Finite Element Methods, P.G. Ciarlet and J.L. Lions editors, NorthHolland, Amsterdam, 523–639, 1991 19. R. Sternberg, M. Suri: Mixed hp finite element methods for problems in elasticity and Stokes flow. Numerische Mathematik 72, 367–389 (1996) 20. M. Vogelius: An analysis of the p-version of the finite element method for nearly incompressible materials. Numerische Mathematik 41, 39–53 (1983) 21. O.C. Zienkiewicz, R.L. Taylor, J.M. Too: Reduced integration techniques in general analysis of plates and shells. Internat. J. Numer. Methods Engrg. 5, 275–290 (1971)
