2004-Fez conference on Differential Equations and Mechanics Electronic Journal of Differential Equations, Conference 11, 2004, pp. 33–40. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

MIXED FINITE ELEMENT DISCRETIZATION OF SOME VARIATIONAL INEQUALITIES ARISING IN ELASTICITY PROBLEMS IN DOMAINS WITH CRACKS

ZAKARIA BELHACHMI, SOUAD TAHIR

Abstract. We consider some mixed variational formulations of elasticity system in domains with cracks. Inequality type conditions are prescribed at the crack faces which results in a model of unilateral contact. Relying in a new variational formulation of these problems in the smooth domain, we study and implement various mixed finite elements methods. We derive convergence rates and optimal error estimates.

1. Introduction

The numerical approximation of solutions of partial differential equations in non smooth domains is often a very difficult and challenging task. Among them, unilateral contact crack problems in elasticity, arises some specific difficulties, both in theoretical and approximation grounds. Such problems are characterized by inequality type and nonlinear boundary conditions prescribed on non smooth part of the boundary [15], describing the mutual non penetration between the crack faces.

The variational formulations associated to these problems lead to solving variational inequalities arising in contact mechanics. We refer the reader to [10, 18, 17] for mathematical foundations and [15] for crack problems. A new approach to crack theory for linear elastic bodies proposed in [16] allows us to solve the problem in the entire domain (including the cracks), reducing the difficulties which occur when dealing with the numerical simulation of such problems. In the framework of this formulation, the discretization with finite elements method leads to similar results, from the accuracy point of view and for developing efficient algorithms, as for classical variational inequalities in unilateral contact problems [3]. For the discretization of variational inequalities, we refer the reader to [7, 12, 13, 14, 17, 23], and for more recent developments to [6, 9, 20, 4].

In [2, 3], the approach called smooth domain method is considered in the case of an elastic membrane. Using a regularization technique, the discretization by various mixed finite element methods is considered, the numerical analysis of the method is carried out and numerical simulations are performed.

2000 Mathematics Subject Classification. 65N30, 74M15, 35J85.

Key words and phrases. Crack problems; variational inequalities; smooth domain method;

augmented Lagrangian formulation; mixed finite elements; a priori estimates.

\begin{matrix}c \\ ◯ \end{matrix}2004 Texas State University - San Marcos.

Published October 15, 2004.

We extend in this work the formulation in the smooth domain to the case of the elasticity system. However, our discretization by finite elements method is based on an augmented Lagrangian formulation [11] since this approach yields many efficient numerical algorithms.

2. The continuous problem

Problem formulation. Let Ω be a bounded domain in ℝ^{2} with smooth boundary Γ, and Γ_{c} ⊂ Ω be a smooth curve without selfintersections. We assume that Γ_{c} can be extended to a closed smooth curve Σ ⊂ Ω, with Σ of class C^{1,1}, and Ω = Ω^{1} ∪ Σ ∪ Ω^{2} divided into two sub-domains Ω^{1}, Ω^{2}. In this case, Σ = ∂Ω^{1} is the boundary of Ω^{1} and Σ ∪ Γ = ∂Ω^{2} is the boundary of Ω^{2}. Let Ω_{c} be the domain Ω \ \overline{Γ_{c}}, then Γ_{c} is called a crack in the elastic body of the reference configuration Ω_{c} (see Figure 1).

The equilibrium problem for a linear elastic body occupying the domain Ω_{c} with the interior crack Γ_{c} can be formulated as follows: find u = (u_{1}, u_{2}), and σ = (σ_{ij}), i, j = 1, 2, such that

−div σ = f in Ω_{c},  (2.1)
Cσ − ϵ(u) = 0 in Ω_{c},  (2.2)
u = 0 on Γ,  (2.3)
[u] ν ≥ 0, [σ ν] = 0, σ_{ν}[u] = 0 on Γ_{C},  (2.4)
σ_{ν} ≤ 0, σ_{t} = 0 on Γ_{C}^{±}, .  (2.5)

Here [u] = u^{+} − u^{−} denotes the jump of the displacement field across Γ_{c}, and the signs ± indicate the positive and negative directions of the normal ν. f = (f_{1}, f_{2}) ∈ L^{2}(Ω)^{2} is a given external force acting on the body. We have used the following standard notation:

σ_{ν} = σ_{ij}ν_{j}ν_{i}, σ_{t} = σ ν − σ_{ν} ν = {σ_{t}^{i}}_{i=1}^{2} , σν = {σ_{ij}ν_{j}}_{i=1}^{2} ,

ϵ_{ij}(u) = \frac{1}{2}(u_{i,j} + u_{j,i}), i, j = 1, 2, ϵ(u) = (ϵ_{ij})_{i,j=1}^{2},

{Cσ}_{ij} = c_{ijkℓ}σ_{kℓ}, c_{ijkℓ} = c_{jikℓ} = c_{kℓij}, c_{ijkℓ} ∈ L^{∞}(Ω).

The tensor C satisfies the ellipticity condition

c_{ijkℓ}ξ_{ji}ξ_{kℓ} ≥ c_{0}|ξ|^{2}, ∀ξ_{ji} = ξ_{ij}, c_{0} > 0.  (2.6)

We use the summation convention over repeated indices.

The smooth domain method. The smooth domain method is based on a mixed variational formulation of problem (2.1)-(2.5). We consider the space

X = {σ = {σ_{ij}} ∈ (L^{2}(Ω_{c}))^{4} : div σ ∈ (L^{2}(Ω_{c}))^{2}} ,

equipped with the norm

∥σ∥_{X} = (∥σ∥_{(L^{2}(Ω_{c}))^{4}}^{2} + ∥div σ∥_{(L^{2}(Ω_{c}))^{2}}^{2})^{1/2} ,

and we define the closed convex set

K = {σ ∈ X : [σ.ν] = 0, on Γ_{c}, σ_{ν} ≤ 0, σ_{t} = 0 on Γ_{c}^{±}}

Figure 1. A crack in the reference configuration.

where [·] denotes the jump across Γ_{c}.

Remark 2.1. The contact condition on Γ_{c} is to be understood in the weak sense of traces by using fractional Sobolev spaces H_{00}^{1/2}(Γ_{c}) and its dual space [19].

The mixed formulation for problem (2.1)-(2.5) reads: Find u = (u_{1}, u_{2}) ∈ (L^{2}(Ω_{c}))^{2}, σ = {σ_{ij}} ∈ K, such that

(Cσ, τ − σ)_{Ω_{c}} + (u, div τ − divσ)_{Ω_{c}} ≥ 0, ∀\overline{σ} ∈ K,

div σ = −f, in Ω_{c}. (2.7)

The well-posedness of this problem follows from the Brezzi-Babuska theory [8]. A proof based on a regularization argument is given in [16] (see also [22] and [20]); we state only the result.

Proposition 2.2. There exists a unique solution to problem (2.7).

The formulation in the smooth domain consists of extending σ and u to the whole domain Ω = Ω_{c} ∪ \overline{Γ}_{c}, and it results in the (close) formulation obtained by replacing Ω_{c} with Ω in (2.7), with the obvious modification of the spaces that we denote also X(Ω) and K(Ω) (see. [16] for details). This formulation allows us to solve the problem in the smooth geometric domain while the constraint at the crack faces is included in the functional space.

The extended problem is also well-posed and it is readily checked that the restriction of the solution (u, σ) to Ω_{c} is the solution of problem (2.7). The converse is (trivially) true if additional smoothness holds for solutions of problem (2.7).

The Discretization of problem (2.7) by various finite elements of low order based on affine finite elements and Raviart-Thomas elements, is considered in [22]. In particular, the study of approximation properties as well as numerical simulations of the discrete problem is performed. In what follows, we consider a new formulation based on an augmented Lagrangian technique.

3. Augmented Langrangian formulation

Mixed variational formulation. The main advantage of the augmented Lagrangian method, compared to the method used in [3], is that it converges without imposing to the penalization parameter to be too small, so the numerical resolution is more stable and efficient. The new formulation is described by the following problem: For r > 0, find u = (u_{1}, u_{2}) ∈ (L^{2}(Ω))^{2}, σ = {σ_{ij}} ∈ K, such that for all τ ∈ K, (Cσ, τ − σ) + (u, div τ − divσ) + r(div σ, div τ − divσ) ≥ r (−f, div τ − divσ),

div σ = −f, in Ω.

(3.1)

Clearly, this problem is equivalent to the minimization over L^{2}(Ω)^{2} × K of the functional

L(v, τ) = \frac{1}{2}(Cτ, τ) + (v, div τ + f) + \frac{r}{2}∥div τ + f∥_{L^{2}(Ω)^{2}}^{2}.

The bilinear form

a_{r}(σ, τ) = (Cσ, τ) + r(div σ, div τ),

is elliptic on K and the bilinear form b(., .) satisfies the usual Brezzi-Babuska condition ([8]). The well-posedness of problem (3.1) and the convergence of the solution (σ(r), u(r)) to the solution (σ, u) of problem (2.7) (see [22]), follows from standard elliptic variational inequalities [18, 12] by applying Stamppacchia’s theorem on the cone K × L^{2}(Ω). Thus, we have the following result.

Proposition 3.1. There exists a unique solution (σ(r), u(r)) to problem (3.1). Moreover, (σ(r), u(r)) converges to (σ, u), the solution of (2.7) when r goes to zero.

Mixed hybrid formulation. The algorithms used to solve the discrete problem corresponding to (3.1) are the steepest descent methods. Another strategy for solving the problem consists in solving a hybrid formulation where the constraint in K is taken into account with a Lagrange multiplier. The resulting algorithm yields to solve a quadratic programming problem of small size to compute the Lagrange multiplier and then a Large linear system to compute the other unknowns. This approach is based on the following mixed hybrid formulation. For the sake of brevity, we will denote by M the space H^{1/2}(Γ_{c})^{2}. We introduce also the following closed convex cone

M_{+} = {μ ∈ H^{1/2}(Γ_{c}); μ ≥ 0, a.e.} .

and we define the new functional L̃ over V × X × M × M_{+}

L̃(v, τ, μ_{t}, μ_{n}) = L(v, τ) + ⟨⟨τ_{t}, μ_{t}⟩⟩_{1/2,Γ_{c}} + ⟨τ_{ν}, μ_{n}⟩_{1/2,Γ_{c}}.

where ⟨·, ·⟩_{1/2,Γ_{c}} denotes the duality product between (H^{1/2}(Γ_{c})) and its dual space. ⟨⟨·, ·⟩⟩_{1/2,Γ_{c}} is the duality product defined by

⟨⟨τ_{t}, ϕ⟩⟩_{1/2,Γ_{c}} = ⟨τ_{t1}, ϕ_{1}⟩_{1/2,Γ_{c}} + ⟨τ_{t2}, ϕ_{2}⟩_{1/2,Γ_{c}}, ∀ϕ = (ϕ_{1}, ϕ_{2}) ∈ M, ϕ_{i}ν_{i} = 0.

Therefore, problem (3.1) can be written as: Find (u, σ, λ_{t}, λ_{n}) ∈ V × X × M × M_{+}, λ_{ti}ν_{i} = 0, such that

a_{r}(σ, τ) + b(u, τ) + ⟨⟨τ_{t}, λ_{t}⟩⟩_{1/2,Γ_{c}} + ⟨τ_{ν}, λ_{n}⟩_{1/2,Γ_{c}} = −r(f, div τ), ∀τ ∈ X,

b(v, σ) = −(f, v), ∀v ∈ V,
⟨⟨σ_{t}, μ_{t}⟩⟩_{1/2,Γ_{c}} = 0, ∀μ_{t} = (μ_{1}, μ_{2}) ∈ M, μ_{i}ν_{i} = 0,
⟨σ_{ν}, μ_{n} − λ_{n}⟩_{1/2,Γ_{c}} ≤ 0, ∀μ ∈ M_{+}.

(3.2)

It is readily checked that solution (u, σ, λ_{t}, λ_{n}) of problem (3.2) is a saddle-point of the Lagrangian L̃.

4. Discrete mixed hybrid problem

We denote by T_{h} a triangulation of Ω made of elements which are triangles (or quadrilateres) with a maximum size h satisfying the usual admissibility assumption, i.e. the intersection of two different elements is either empty, a vertex, or a whole edge. In addition, T_{h} is assumed regular, i.e. the ratio of the diameter of any element T ∈ T_{h} to the diameter of its largest inscribed ball is bounded by a constant σ independent of T and h. We will assume that the endpoints of Γ_{c} are vertices of the triangulation. The set of nodes on Γ_{c} are denoted by c_{1} = x_{0}, x_{1}, …, x_{I−1}, x_{I} = c_{2} and we set t_{i} =]x_{i−1}, x_{i}[.

In [22] we study various discretizations based on the finite elements method for both problem (3.1) and problem (3.2). In particular, we perform the complete analysis and the numerical simulations with the (lower order) PEERS element and the BDMS (Brezzi-Douglas-Marini and Steinberg) elements. The use of the above elements is based on a formulation obtained by modifying the Hellinger-Reissner principle. In this approach, the symmetry of the stress tensor σ is relaxed and only imposed by means of the Lagrange multiplier.

We present here another example of finite elements also studied in [22] due to Arnold and Winther [1], based on the unaltered Hellinger-Reissner principle.

For any T ∈ T_{h}, and for any Banach space E, we denote by P_{k}(T, E) the space of polynomial functions over T with values in E, of degree less than k. As usual we denote by P_{k}(T) the space P_{k}(T, ℝ). We define the associated local space Σ_{T} as

Σ_{T} = P_{2}(T, ℝ_{sym}^{2×2}) + {τ ∈ P_{3}(T, ℝ_{sym}^{2×2}); div τ = 0}
= {τ ∈ P_{3}(T, ℝ_{sym}^{2×2}); div τ ∈ P_{1}(T)^{2}} .

We introduce the following discrete spaces: For h > 0,

X_{h} = {σ_{h} ∈ X; σ_{h|T} ∈ Σ_{T} , ∀T ∈ T_{h}} ,
V_{h} = \{v_{h} ∈ (C(\overline{Ω}))^{2}; v_{h|T} ∈ (P_{1}(T))^{2}, ∀T ∈ T_{h}\}.

Concerning the approximation of Lagrange multipliers, we introduce the space

W_{h}^{1}(Γ_{c}) = \{μ_{h} ∈ C(\overline{Γ}_{c}), μ_{h|t_{i}} ∈ P_{1}(t_{i}), 0 ≤ i ≤ I − 1\} .

We also denote by W_{h}^{1}(Γ_{c}) the space (W_{h}^{1}(Γ_{c}))^{2}.

Therefore, we can define various version of the discrete space M_{h} and the discrete convex cone M_{h+}. An example of such choices consists to take [22]

M_{h} = W_{h}^{1}(Γ_{c})
M_{h+} = {μ_{h} ∈ W_{h}^{1}(Γ_{c}), ∫_{Γ_{c}} μ_{h}ψ_{h} dΓ ≥ 0, ∀ψ_{h} ∈ W_{h}^{1}, ψ_{h} ≥ 0}.

We will denote for brevity M_{h} the convex cone M_{h} × M_{h+} and λ_{h} = (λ_{ht} = (λ_{ht1}, λ_{ht2}), λ_{hn}), λ_{hti}ν_{i} = 0. The discrete problem is now the same as problem (3.2) when replacing the unknowns and spaces by the finite dimensional analogues.

The following result is proved in [22].

Proposition 4.1. Assume that the set M_{h} and M_{h+} are given as above. Then, the discrete problem corresponding to (3.2) admits a unique solution.

The error analysis and the study of convergence rates is based on uniform inf-sup condition (with respect to h) (see [3]) we only states the main error estimate when the multiplier spaces are chosen as in Proposition 4.1 (see [22] for details). Let us define

Z = {σ = {σ_{ij}} ∈ (H^{1}(Ω))^{4} : div σ ∈ (H^{1}(Ω))^{2}} ,

and denotes by Z(Ω^{ℓ}), ℓ = 1, 2 the space of the restrictions of functions of Z to Ω^{ℓ}. We have the following result.

Theorem 4.2. Let (u, σ, λ = (λ_{t1}, λ_{t2}, λ_{n})) be the solution of problem (3.2). Suppose that u_{|Ω^{1}} ∈ H^{2}(Ω^{1}), u_{|Ω^{2}} ∈ H^{2}(Ω^{2}) and also σ_{|Ω^{1}} ∈ Z(Ω^{1}), σ_{|Ω^{2}} ∈ Z(Ω^{2}). Let (U_{h}, λ_{h} = (λ_{ht1}, λ_{ht2}, λ_{hn})) be the solution of the discrete problem associated to (3.2) with the specified choice of M_{h} and M_{h+} given above. Then the following estimate holds

\begin{matrix}∥u − u_{h}∥_{V} + ∥σ − σ_{h}∥_{X} + ∥λ_{t} − λ_{ht}∥_{(H^{1/2}(Γ_{c}))^{2}} + ∥λ_{n} − λ_{hn}∥_{H^{1/2}(Γ_{c})} \\ ≤ C(r, u, σ, λ)h^{3/4} . \end{matrix} (4.1)

The constant C(r, u, σ, λ) depends linearly on ∥u_{|Ω^{ℓ}}∥_{H^{2}(Ω^{ℓ})}, ∥σ_{|Ω^{ℓ}}∥_{H^{1}(Ω^{ℓ})^{4}}, and ∥div σ_{|Ω^{ℓ}}∥_{H^{1}(Ω^{ℓ})^{2}}, ℓ = 1, 2.

Remark 4.3. More technical arguments allow us to avoid the regularity assumption div σ ∈ (H^{1}(Ω^{ℓ}))^{2}. In this case σ is estimated in the L^{2}-norm.

5. Implementation details

To perform the computations for the mixed hybrid problem, the matrix formulation of discrete problem (3.2) is derived. Let V, U denote the vectors with the entries given by the nodal values of the functions (v_{h}, τ_{h}) and (u_{h}, σ_{h}), respectively. Let M and Λ be the vectors with the entries given by the nodal values of μ_{h} and λ_{h}, respectively, for the various choices of the space M_{h} and the convex set M_{h+}. Therefore, the saddle-point problem for the Lagrangian can be rewritten in finite dimensional setting: Find U = (u_{h}, σ_{h}) and Λ, defined by the following max-min condition

max_{SM≥0}(min_{V} \frac{1}{2}^{t}VKV −^{t}VF + (^{t}V L)SM),  (5.1)

where K denotes the stiffness matrix, F is the vector corresponding to the external loading and the matrix S expresses the sign conditions for the multipliers.

The solution (U, Λ) of (5.1) satisfies the saddle-point conditions and we have

U = K^{−1}(F − L SΛ).  (5.2)

Therefore, for Φ = SΛ, the saddle-point problem (5.1) can be rewritten as a quadratic programming problem

min_{Φ≥0}(\frac{1}{2}^{t}Φ^{t}LK^{−1}LΦ −^{t}Φ^{t}LK^{−1}F + \frac{1}{2}^{t}FK^{−1}F).  (5.3)

If \overline{Φ} is the solution of (5.3) then Λ = S^{−1}\overline{Φ}. The solution U is obtained by solving (5.2).

Several numerical simulations and experiments are performed in [22] for this augmented Lagrangian method applied to the new formulation of crack problems in the smooth domain. The computations with various mixed finite element discretizations, for both unsymmetric stress tensor formulation (PEERS, BDMS) and the symmetric one considered in this article, give satisfactory results confirming the theoretical estimates. The use of the augmented Lagrangian formulation (i.e. r ≠ 0) is unnecessary for practical computations when compatible pairs of finite elements for the displacements and stresses are chosen, however it allows using efficient algorithms to solve the discrete problem.

References

[1] D. N. Arnold, R. Winther; Mixed finite elements for elasticity, Numer. Math., 92 (2002), 401-419.

[2] Z. Belhachmi, J. M. Sac-Epée, J. Sokolowski; Approximation par la méthode des éléments finis de la formulation en domaine régulier de problèmes de fissures, C. R. Acad. Sci. Paris, Série I 338 (2004), 499-504.

[3] Z. Belhachmi, J.M. Sac-Epée, J. Sokolowski; Mixed finite element methods for a smooth domain formulation of a crack problem. submitted.

[4] Z. Belhachmi, F. Ben Belgacem; Quadratic finite element for Signorini problem, Math. Comp. 72, (2003), 83-104.

[5] F. Ben Belgacem, P. Hild, P. Laborde; Extention of the mortar finite element method to a variational inequality modeling unilateral contact, Math. Models Methods. Appl. Sci., 9, (1999), 287-303.

[6] F. Ben Belgacem, Y. Renard; Hybrid finite element methods for the Signorini problem, Math. Comput. 72, (2003), 1117-1145.

[7] F. Brezzi, W. W. Hager and P. A. Raviart; Error estimates for the finite element solution of variational inequalities, Numer. Math., 28, (1977), 431-443.

[8] F. Brezzi, M. Fortin; Mixed and hybrid finite element methods, Springer Verlag, New York, Springer Series in Computational Mathematics, 15, (1991).

[9] P. Coorevits, P. Hild, K. Lhalouani, T. Sassi; Mixed finite element methods for unilateral problems: convergence analysis and numerical studies. Math. Comp., 71, 237, (2001), 1-25.

[10] G. Duvaut, J.-L. Lions; Les inéquations en mécanique et en physique, Dunod, (1972).

[11] R. Glowinski, P. Le Tallec; Augmented Lagrangian and operator-splitting methods in nonlinear mechanics. Siam, Philadelphia.

[12] R. Glowinski; Lectures on numerical methods for nonlinear variational problems, Springer, Berlin, (1980).

[13] J. Haslinger and I. Hlaváček; Contact between Elastic Bodies -2.Finite Element Analysis, Aplikace Matematiky, 26, (1981), 263–290.

[14] J. Haslinger, I. Hlaváček, J. Nečas; Numerical Methods for Unilateral Problems in Solid Mechanics, in the Handbook of Numerical Analysis, Vol. IV, Part 2, P.G. Ciarlet & J.-L. Lions eds, North-Holland, (1996).

[15] A. M. Khludnev, V. A. Kovtunenko; Analysis of cracks in solids, Southampton-Boston, WIT press, (2000).

[16] A. M. Khludnev, J. Sokolowski; Smooth domain method for crack problems. to appear in Quarterly of Applied Mathematics.

[17] N. Kikuchi, J. Oden; Contact problems in elasticity: A study of variational inequalities and finite element methods, SIAM, 1988.

[18] D. Kinderlehrer, G. Stamppachia; An introduction to variational inequalities and their applications, Academic Press, (1980).

[19] J.-L. Lions, E. Magenes; Problèmes aux limites non homogènes, Dunod, (1968).

[20] L. Slimane, A. Bendali, P. Laborde; Mixed formulations for a class of variational inequalities, M2AN, 38, 1, (2004), 177-201.

[21] L. Slimane; Méthodes mixtes et traitement du verouillage numérique pour la résolution des inéquations variationnelles, Ph.D. Thesis, INSA Toulouse, France (2001).

[22] S. Tahir; Méthodes d’approximation et analyse a posteriori des inéquations variationnelles et applications à des problèmes de fissures en élasticité, Ph.D. Thesis, université de Metz, to be published.

[23] Z.-H. Zhong; Finite Element Procedures for Contact-Impact Problems, Oxford University Press, 1993.

Zakaria Belhachmi Laboratoire de Mathématiques (LMAM), Université de Metz, ISGMP, Batiment A, Ile du Saulcy, 57045 Metz, France.

E-mail address: belhach@poncelet.sciences.univ-metz.fr

Souad Tahir Laboratoire de Mathématiques (LMAM), Université de Metz, ISGMP, Batiment A, Ile du Saulcy, 57045 Metz, France.

E-mail address: tahir@poncelet.sciences.univ-metz.fr