- Email: [email protected]
Numerical treatment of a class of thermal contact problems
Available online at www.sciencedirect.com
ScienceDirect Mathematics and Computers in Simulation (
)
– www.elsevier.com/locate/matcom
Original articles
Numerical treatment of a class of thermal contact problems Oanh Chau, Rachid Oujja ∗ Universit´e de La R´eunion, Laboratoire PIMENT EA4518, 97715 Saint-Denis Messag cedex 9 La R´eunion, France Received 13 March 2014; received in revised form 9 December 2014; accepted 22 December 2014
Abstract A mathematical problem modeling dynamical behavior of viscoelastic materials is studied. The problem is related to highly nonlinear and non-smooth phenomena like contact, friction and a class of dynamic contact problems with friction and thermal effects. The applied approach leads to a system involving variational inequalities and a fully discrete scheme for numerical approximations and the analysis of error estimate order is provided. Numerical computations are given so as to illustrate the mathematical model. c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Viscoelastic materials; Sub-differential contact condition; Heat transfer; Friction; Numerical simulations
1. Introduction Due to their growing complexity in the engineering sciences, contact problems for viscoelastic materials still remain a serious challenge. For example, contact phenomena are omnipresent in human skeletal systems, see for example [18]. Otherwise, many food materials used in process engineering are viscoelastic, see e.g. [11], and consequently, mathematical models can be very helpful in understanding various problems related to the product development, packing, transport, shelf life testing, thermal effects and heat transfer. It is thus important to study mathematical models that can be used to describe the dynamical behavior of a given viscoelastic material subjected to various highly nonlinear and even non-smooth phenomena like contact, friction and thermal effects. The earlier studies relating to contact problems in the framework of variational formulation and infinitesimal deformations could be found in the reference works [6,12,19]. Thermoviscoelastic contact problems by taking into account the evolution of the temperature parameter could be found in [2,7,10]. Since considerable progress and development have been made in contact mechanics, as widely attested the pioneering works in [9,14,22,21], including mathematical modeling, mathematical analysis, with specific studies on the error estimates analysis, numerical approximation and numerical simulations. In these works the authors show the cross-fertilization between numerous frictional new models arising in contact mechanics and various types of abstract variational inequalities. Further extensions to non convex contact conditions with non-monotone and possible multi-valued constitutive laws led to the active domain of non-smooth mechanics within the framework of the so-called hemivariational ∗ Corresponding author.
E-mail addresses: [email protected] (O. Chau), [email protected] (R. Oujja). http://dx.doi.org/10.1016/j.matcom.2014.12.007 c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.
2
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
inequalities. For a mathematical as well as mechanical treatment we refer to [8,15,16,20]. A large deformation formulation for frictional contact problem has been treated in [5]. This work constitutes in some sense a continuation paper of the results obtained in [1]. The work in [1] has been devoted to a mathematical study of a class of dynamical long memory viscoelastic thermal problems where the contact is governed by a general sub-differential condition. Qualitative results like existence and uniqueness result of weak solutions on displacement and temperature fields have been proved but no numerical approximations have been performed. Here we follow the latter work and propose a numerical scheme for the approximation of the solution fields so as to elaborate a general numerical analysis of error estimates. The theoretical results are then illustrated by different numerical simulations. The paper is organized as follows. In Section 2 we give a short description of the mathematical model and recall the main existence and uniqueness result. In Section 3, we introduce a fully discrete approximation scheme, derive and prove an optimal order error estimate under certain solution regularity assumptions. Finally, in Section 4, we present a set of numerical simulations showing the evolution of the displacement field, the temperature field as well as the Von Mises stress norm. 2. The contact problems
A viscoelastic body occupies the domain Ω with surface Γ that is partitioned into three disjoint measurable parts Γ1 , Γ2 and Γ3 , such that meas (Γ1 ) > 0. Let [0, T ] be the time interval of interest, where T > 0. The body is clamped on Γ1 × (0, T ) and therefore the displacement field vanishes there. We also assume that a volume force of density f0 acts in Ω × (0, T ) and that surface traction of density f2 acts on Γ2 × (0, T ). The body may come in contact with an obstacle, the foundation, over the potential contact surface Γ3 . The model of the contact is specified by a general sub-differential boundary condition, where thermal effects may occur in the frictional contact with the basis. We are interested in the dynamic evolution of the body. We denote by Sd the space of second order symmetric tensors on Rd (d = 1, 2, 3) in practical case, while “ · ” and |·| will represent the inner product and the Euclidean norm on Sd and Rd . Let Ω ⊂ Rd be a bounded domain with a Lipschitz boundary Γ and let ν denote the unit outer normal on Γ . Everywhere in the sequel the indexes i and j run from 1 to d, summation over repeated indices is implied and the index that follows a comma represents the partial derivative with respect to the corresponding component of the independent variable. We also use the following notation: d H = L 2 (Ω ) , H = {σ = (σi j )|σi j = σ ji ∈ L 2 (Ω ), 1 ≤ i, j ≤ d}, H1 = {u ∈ H |ε(u) ∈ H},
H1 = {σ ∈ H|Div σ ∈ H }.
Here ε : H1 −→ H and Div : H1 −→ H are the deformation and the divergence operators, respectively, defined by: 1 (u i, j + u j,i ), Div σ = (σi j, j ). 2 The spaces H , H, H1 and H1 are real Hilbert spaces endowed with the canonical inner products given by: (u, v) H = u i vi d x, (σ , τ )H = σi j τi j d x, ε(u) = (εi j (u)),
Ω
εi j (u) =
(u, v) H1 = (u, v) H + (ε(u), ε(v))H ,
Ω
(σ , τ )H1 = (σ , τ )H + (Div σ , Div τ ) H .
We recall that C denotes the class of continuous functions and C m , m ∈ N∗ is the set of m times differentiable functions. Finally D(Ω ) denotes the set of infinitely differentiable real functions with compact support in Ω ; and W m, p ,
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
3
m ∈ N, 1 ≤ p ≤ +∞ for the classical Sobolev spaces; and H0m (Ω ) := {w ∈ W m,2 (Ω ), w = 0 on Γ1 },
m ≥ 1.
The contact mechanical problem is now formulated. Problem. Q: Find a displacement field u : Ω × [0, T ] −→ Rd and a stress field σ : Ω × [0, T ] −→ Sd and a temperature field θ : Ω × [0, T ] −→ R+ such that for a.e. t ∈ (0, T ): t σ (t) = Aε(u(t)) ˙ + Gε(u(t)) + B(t − s) ε(u(s)) ds − θ (t) Ce in Ω (2.1) 0
u(t) ¨ = Div σ (t) + f0 (t) u(t) = 0
in Ω
(2.2)
on Γ1
σ (t)ν = f2 (t)
(2.3)
on Γ2
(2.4)
ϕ(w) − ϕ(u(t)) ˙ ≥ −σ (t)ν · (w − u(t)) ˙ ∂ u˙ i θ˙ (t) − div(K c ∇θ (t)) = −ci j (t) + q(t) in Ω ∂ xj
u(t) ∈ U,
−ki j
∂θ (t) n i = ke (θ (t) − θ R ) ∂ xj
θ (t) = 0 θ (0) = θ0 u(0) = u0 ,
∀w ∈ U on Γ3
on Γ3
(2.5) (2.6) (2.7)
on Γ1 ∪ Γ2
(2.8)
in Ω
(2.9)
u(0) ˙ = v0
in Ω .
(2.10)
Here, (2.1) is the long memory thermo-visco-elastic constitutive law of the body, σ the stress tensor, A is the viscosity operator, G for the elastic operator, Ce := (ci j ) represents the thermal expansion tensor, and B is the so called tensor of relaxation which defines the long memory of the material, as an important particular case, when B ≡ 0, we find again the usual visco-elasticity of short memory. Eq. (2.2) is the equation of motion where the mass density ϱ ≡ 1. The equation in (2.3) is the clamped condition and in (2.4) is the traction condition. On the contact surface, the general relation (2.5) is a sub-differential boundary condition, verifying D(Ω )d ⊂ U ⊂ H1 , where U represents the set of contact admissible test functions and D(Ω )d is the distribution space. σ ν denotes the Cauchy stress vector on the contact boundary and ϕ : Γ3 × Rd −→ R is a given function. Various situations may be modeled by such a condition, and a concrete example will be given below. The differential equation (2.6) describes the evolution of the temperature field, where K c := (ki j ) represents the thermal conductivity tensor, q(t) the density of volume heat sources. The associated temperature boundary condition is given by (2.7)–(2.8), where θ R is the temperature of the foundation, and ke is the heat exchange coefficient between the body and the obstacle. Finally in (2.9)–(2.10), u0 , v0 , θ0 represents the initial displacement, velocity and temperature, respectively. To derive the variational formulation of the mechanical problems (2.1)–(2.10) we need additional notations. Thus, let V denote the closed subspace of H1 defined by D(Ω )d ⊂ V = {v ∈ H1 | v = 0 on Γ1 } ∩ U ; E = {η ∈ H 1 (Ω ), η = 0 on Γ1 ∪ Γ2 },
F = L 2 (Ω ).
Since meas Γ1 > 0, Korn’s inequality holds: there exists C K > 0 which depends only on Ω and Γ1 such that ∥ε(v)∥H ≥ C K ∥v∥ H1
∀v ∈ V.
A proof of Korn’s inequality may be found in [17, p. 79]. On V we consider the inner product given by (u, v)V = (ε(u), ε(v))H
∀ u, v ∈ V,
4
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
and let ∥ · ∥V be the associated norm, i.e. ∥v∥V = ∥ε(v)∥H
∀ v ∈ V.
It follows that ∥ · ∥ H1 and ∥ · ∥V are equivalent norms on V and therefore (V, ∥ · ∥V ) is a real Hilbert space. Moreover, by Sobolev’s trace theorem, we have a constant C0 > 0 depending only on Ω , Γ1 , and Γ3 such that ∥v∥ L 2 (Γ3 ) ≤ C0 ∥v∥V
∀ v ∈ V.
In the study of the mechanical problem (2.1)–(2.10), we assume the following assumptions (see e.g. [1,13]). The viscosity operator A : Ω × Sd −→ Sd , (x, τ ) −→ (ai jkh (x) τkh ) is linear on the second variable and satisfies the usual properties of ellipticity and symmetry, i.e. (i) ai jkh ∈ W 1,∞ (Ω ); (ii) Aσ · τ = σ · Aτ ∀σ , τ ∈ Sd , a.e. in Ω ; (iii) there exists m A > 0 such that Aτ · τ ≥ m A |τ |2 ∀τ ∈ Sd , a.e. in Ω .
(2.11)
The elasticity operator G : Ω × Sd −→ Sd , (x, τ ) −→ (bi jkh (x) τkh ) is linear on the second variable, with bi jkh ∈ L ∞ (Ω ).
(2.12)
The relaxation tensor B : [0, T ] × Ω × Sd −→ Sd , (t, x, τ ) −→ (Bi jkh (t, x) τkh ) satisfies (i) Bi jkh ∈ W 1,∞ (0, T ; L ∞ (Ω )); (ii) B(t)σ · τ = σ · B(t)τ ∀σ , τ ∈ Sd , a.e. t ∈ (0, T ), a.e. in Ω .
(2.13)
We suppose the body forces and surface tractions satisfy f0 ∈ W 1,2 (0, T ; H ),
f2 ∈ W 1,2 (0, T ; L 2 (Γ2 )d ).
(2.14)
For the thermal tensors and the heat sources density, we suppose that Ce = (ci j ),
ci j = c ji ∈ L ∞ (Ω ),
q ∈ W 1,2 (0, T ; L 2 (Ω )).
(2.15)
The boundary thermic data satisfy ke ∈ L ∞ (Ω ; R+ ),
θ R ∈ W 1,2 (0, T ; L 2 (Γ3 )).
(2.16)
The thermal conductivity tensor verifies the usual symmetry and ellipticity. For some m K > 0 and for all (ξi ) ∈ Rd , K c = (ki j ),
ki j = k ji ∈ L ∞ (Ω ),
ki j ξi ξ j ≥ m K ξi ξi .
(2.17)
We assume that the initial data satisfy the conditions v0 ∈ V ∩ H02 (Ω )d ,
u0 ∈ V,
θ0 ∈ E ∩ H02 (Ω ).
On the contact surface, the following frictional contact function ψ(w) := ϕ(w) da Γ3
(2.18)
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
verifies (i) ψ : V −→ R is well defined, continuous and convex; (ii) there exists a sequence of differentiable convex functions (ψ ) : V −→ R such that ∀w ∈ L 2 (0, T ; V ), Tn T ψn (w(t)) dt −→ ψ(w(t)) dt, for n −→ +∞; 0
0
(iii) for all sequences (wn ) and w in W 1,2 (0, T ; V ) such that wn ⇀ w, w′n ⇀ w′ weakly in L 2 (0, T ; V ), T T then lim inf ψ (w (t)) dt ≥ ψ(w(t)) dt; n n n−→+∞ 0 0 (iv) ∀w ∈ V, (w = 0 on Γ3 =⇒ ∀n ∈ N, ψn′ (w) = 0V ′ ).
5
(2.19)
Here ψn′ (v) denotes the Fr´echet derivative of ψn at v. By using Green’s formula, we obtain the variational formulation of the mechanical problem Q in the abstract form as follows. Problem. QV : Find u : [0, T ] → V , θ : [0, T ] → E satisfying a.e. t ∈ (0, T ): ⟨u(t) ¨ + A u(t) ˙ + B u(t) + C θ (t), w − u(t)⟩ ˙ V ′ ×V t + ψ(w) − ψ(u(t)) ˙ + B(t − s) ε(u(s)) ds, ε(w) − ε(u(t)) ˙ 0
≥ ⟨f (t), w − u(t)⟩ ˙ ∀w ∈ V ; V ′ ×V ˙ + Q(t) in E ′ ; θ˙ (t) + K θ (t) = R u(t) u(0) = u0 , u(0) ˙ = v0 , θ (0) = θ0 .
H
Here the operators and functions A, B : V −→ V ′ , C : E −→ V ′ , K : E −→ E ′ , R : V −→ E ′ , f : [0, T ] −→ V ′ , and Q : [0, T ] −→ E ′ are defined by ∀v ∈ V , ∀w ∈ V , ∀τ ∈ E, ∀η ∈ E: ⟨A v, w⟩V ′ ×V = (A(εv), εw)H ; ⟨B v, w⟩V ′ ×V = (G(εv), εw)H ; ⟨Cτ, w⟩V ′ ×V = −(τ Ce , εw)H ; ⟨f (t), w⟩V ′ ×V = (f0 (t), w) H + (f2 (t), w)(L 2 (Γ2 ))d ; ′ ⟨Q(t), η⟩ E ×E = ke θ R (t) η da + q(t) η d x; Γ3
⟨K τ, η⟩
E ′ ×E
=
i, j=1
⟨R v, η⟩
Ω
d
∂τ ∂η ke τ · η da; dx + ki j ∂ x j ∂ xi Γ3 Ω
E ′ ×E
ci j
=− Ω
∂vi η d x. ∂x j
Let us recall now the following main mathematical result (see for details [1]): Theorem 2.1. Assume that (2.11)–(2.19) hold. Then there exists a unique solution {u, θ} to problem QV with the regularity: u ∈ W 2,2 (0, T ; V ) ∩ W 2,∞ (0, T ; H ) θ ∈ W 1,2 (0, T ; E) ∩ W 1,∞ (0, T ; F). 3. A numerical scheme Let {u, θ} be the unique solution of the problem QV , and introduce the velocity variable v(t) = u(t), ˙
∀t ∈ [0, T ].
6
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
Then t
u(t) = u0 +
v(s) ds,
∀t ∈ [0, T ].
0
From Theorem 2.1 we see that {v, θ} verify for all t ∈ [0, T ]: ⟨˙v(t) + A v(t) + B u(t) + C θ (t), w − v(t)⟩V ′ ×V t B(t − s) ε(u(s)) ds, ε(w) − ε(u(t)) ˙ + + ψ(w) − ψ(v(t)) 0 H ≥ ⟨f (t), w − v(t)⟩V ′ ×V , ∀ w ∈ V. ⟨θ˙ (t), η⟩ F + ⟨K θ (t), η⟩ E ′ ×E = ⟨Rv(t), η⟩ E ′ ×E + ⟨Q(t), η⟩ E ′ ×E , u(0) = u0 ,
v(0) = v0 ,
(3.1)
∀ η ∈ E.
(3.2)
θ (0) = θ0 ,
(3.3)
with the regularity:
v ∈ W 1,2 (0, T ; V ) ∩ W 1,∞ (0, T ; H ) θ ∈ W 1,2 (0, T ; E) ∩ W 1,∞ (0, T ; F).
(3.4)
We make the following additional assumptions on the solution and contact function: v ∈ W 2,2 (0, T ; H );
(3.5)
θ∈W
(3.6)
2,2
(0, T ; F);
ψ is Lipschitz continuous on V.
(3.7)
Before we continue, it is important to remark that proposition (3.7) is verified in the typical example of sub-differential condition we treated (see Section 4). Now let V h ⊂ V and E h ⊂ E be a family of finite dimensional subspaces, with h > 0 a discretization parameter. We divide the time interval [0, T ] into N equal parts: tn = n k, n = 0, 1, . . . , N , with the time step k = T /N . For a continuous function w ∈ C([0, T ]; X ) with values in a space X , we use the notation wn = w(tn ) ∈ X . Then from (3.1) and (3.2) we introduce the following fully discrete scheme. N N Problem. P hk . Find vhk = {vnhk }n=0 ⊂ V h , θ hk = {θnhk }n=0 ⊂ E h such that
v0hk = v0h ,
θ0hk = θ0h
(3.8)
and for n = 1, . . . , N , vhk − vhk n n−1
hk , wh − vnhk ⟩V ′ ×V + ⟨A vnhk , wh − vnhk ⟩V ′ ×V + ⟨B un−1 n−1 hk h hk hk h hk ′ h hk B(tn − tm ) ε(um ), ε(w ) − ε(vn ) + ⟨C θn−1 , w − vn ⟩V ×V + ψ(w ) − ψ(vn ) + k
k
, wh − vnhk
H
m=0
≥ ⟨fn , w
h
θ hk − θ hk n n−1 k
− vnhk ⟩V ′ ×V , , ηh
F
H
∀w ∈ V . h
h
+ ⟨K θnhk , ηh ⟩ E ′ ×E = ⟨R vnhk , ηh ⟩ E ′ ×E + ⟨Q n , ηh ⟩ E ′ ×E ,
(3.9) ∀ ηh ∈ E h
(3.10)
where hk unhk = un−1 + k vnhk ,
u0hk = u0h .
Here u0h ∈ V h , v0h ∈ V h , θ0h ∈ E h are suitable approximations of the initial values u0 , v0 , θ0 .
(3.11)
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
7
–
The main result of this section is the following. Theorem 3.2. We keep the assumptions of Theorem 2.1. Under the additional assumption (3.5)–(3.6), then for the unique solution vhk ⊂ V h , θ hk ⊂ E h of the discrete problem P hk , we have the following error estimate max
1≤n≤N
∥vn − vnhk ∥2H
1≤n≤N
n=1
c ∥u0 − u0h ∥2V
≤
N N hk 2 hk 2 hk 2 + k ∥vn − vn ∥V + max ∥θn − θn ∥ F + k ∥θn − θn ∥ E + c ∥v0 − v0hk ∥2H
n=1
+ c ∥θ0 − θ0h ∥2F
(3.12)
+ c max ∥vn − wnh ∥ H + c max ∥θn − ηnh ∥2F + c k 1≤n≤N
+c
N −1
1≤n≤N
∥(v j − whj ) − (v j+1 − whj+1 )∥ H
∥v j − whj ∥2V + c k
j=1
2 +c
j=1
+ c k2 + c k
N
N −1
N
∥θ j − ηhj ∥2E
j=1
∥θ j − ηhj − (θ j+1 − ηhj+1 )∥ F
2
j=1 N
∥v j − whj ∥V ,
(3.13)
j=1
where for j = 1, . . . , N , whj ∈ V h , ηhj ∈ E h are arbitrary. Proof. Denote by en = vn − vnhk ,
εn = θn − θnhk ,
0 ≤ n ≤ N,
the numerical solution errors. We begin by the estimate of (εn ). Let fix n = 1, . . . , N . From (3.2) and (3.10), we have hk θnhk − θn−1 θ˙n − , ηh + ⟨K θn − K θnhk , ηh ⟩ E ′ ×E = ⟨R vn − R vnhk , ηh ⟩ E ′ ×E , F k
Writing then θ˙n −
hk θnhk − θn−1
k
= θ˙n −
θn − θn−1 εn − εn−1 + , k k
and replacing ηh by ηnh − θn + εn we obtain ε − ε n n−1 , εn + ⟨K θn − K θnhk , εn ⟩ E ′ ×E F k hk θnhk − θn−1 = ⟨R vn − R vnhk , ηnh − θn + εn ⟩ E ′ ×E − θ˙n − , ηnh − θn F k θn − θn−1 hk h − ⟨K θn − K θn , ηn − θn ⟩ E ′ ×E − θ˙n − , εn . F k Consider the quantity for n = 1, . . . , N , ε − ε n n−1 Ξn := , εn + ⟨K θn − K θnhk , εn ⟩ E ′ ×E . F k We have Ξn ≥
1 ∥εn ∥2F − ∥εn−1 ∥2F + m K ∥εn ∥2E . 2k
∀ ηh ∈ E h .
8
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
Summing then Ξ j from j = 1 to j = n, and after some manipulation, we obtain n n n θ j − θ j−1 1 ˙ 2 ∥ε j ∥2E ≤ c ∥θ j − ηhj ∥2E ∥εn ∥2F − ∥ε0 ∥2F + θ j − +c F 2k k j=1 j=1 j=1 +c
n
∥e j ∥2V + c
j=1
n ε j − ε j−1 , θ j − ηhj . F k j=1
Denote now by Mε := max ∥εn ∥ F , 1≤n≤N
A0 := max ∥θn − ηnh ∥ F , 1≤n≤N
A1 :=
A2 :=
N θ j − θ j−1 ˙ 2 θ j − , F k j=1 N
∥θ j − ηhj ∥2E ,
j=1
A3 :=
N
∥θ j − ηhj − (θ j+1 − ηhj+1 )∥ F .
j=1
We deduce from the last inequality that for some constant c > 0 and for n = 1, . . . , N , ∥εn ∥2F + k
n
∥ε j ∥2E ≤ c ∥ε0 ∥2F + c A20 + c k (A1 + A2 ) + c A3 Mε + c k
j=1
n
∥e j ∥2V .
(3.14)
j=1
We now turn to the estimate of (en ). Let fix n = 1, . . . , N . Using (3.1) where we put t = tn and w = vnhk , and adding to (3.9) where wh = wnh , we have hk hk vn − vn−1 , wnh − vnhk + ⟨A vn , vnhk − vn ⟩V ′ ×V + ⟨A vnhk , wnh − vnhk ⟩V ′ ×V (˙vn , vnhk − vn ) H + k H hk hk + ⟨B un , vnhk − vn ⟩V ′ ×V + ⟨B un−1 , wnh − vnhk ⟩V ′ ×V + ⟨C θn , vnhk − vn ⟩V ′ ×V + ⟨C θn−1 , wnh − vnhk ⟩V ′ ×V
+ Rnhk + ψ(wnh ) − ψ(vn ) ≥ ⟨fn , wnh − vn ⟩V ′ ×V ,
∀ wnh ∈ V h ,
where for n = 1, . . . , N : tn n−1 hk hk Rn = B(tn − s) ε(u(s)) ds − k B(tn − tm ) ε(um ), −ε(en ) 0
H
m=0
n−1 hk + k B(tn − tm ) ε(um ), ε(wnh ) − ε(vn ) . m=0
H
Writing then hk hk vnhk − vn−1 vnhk − vn−1 + ; v˙ n = v˙ n − k k hk hk A vn = A vn − A vn + A vn , we obtain e − e n n−1 , en + ⟨A vn − A vnhk , en ⟩V ′ ×V H k vhk − vhk vn − vn−1 n n−1 ≤ ⟨A vnhk , wnh − vn ⟩V ′ ×V + v˙ n − , −en + , wnh − vn H H k k
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
9
hk hk + ⟨B un − B un−1 , −en ⟩V ′ ×V + ⟨B un−1 , wnh − vn ⟩V ′ ×V hk hk + ⟨C θn − C θn−1 , −en ⟩V ′ ×V + ⟨C θn−1 , wnh − vn ⟩V ′ ×V
+ Rnhk + ψ(wnh ) − ψ(vn ) − ⟨fn , wnh − vn ⟩V ′ ×V . Consider the quantity for n = 1, . . . , N , e − e n n−1 + ⟨A vn − A vnhk , en ⟩V ′ ×V . , en En = H k We have En ≥
1 ∥en ∥2H − ∥en−1 ∥2H + m A ∥en ∥2V . 2k
To continue, we denote in the sequel by Me := max ∥en ∥ H , 1≤n≤N
B0 := max ∥vn − wnh ∥ H , 1≤n≤N
B1 :=
B2 :=
N 2 v j − v j−1 − v˙ j , H k j=1 N
∥v j − whj ∥2V ;
Bˆ 2 :=
j=1
B3 :=
N
∥v j − whj ∥V ,
j=1
N −1
∥(v j − whj ) − (v j+1 − whj+1 )∥ H ,
j=1
I :=
N −1 j=1
B4 :=
N
tj
v−k
0
j 2 vi , V
i=1
∥u j − u j−1 ∥2V ,
A4 :=
j=1
B5 :=
N
N
∥θ j − θ j−1 ∥2F ,
j=1
|⟨B u j−1 , v j − whj ⟩V ′ ×V |,
A5 :=
j=1
F :=
N
N
|⟨C θ j−1 , v j − whj ⟩V ′ ×V |,
j=1
|⟨A v j , v j − whj ⟩V ′ ×V + ψ(v j ) − ψ(whj ) − ⟨f j , v j − whj ⟩V ′ ×V |.
j=1
Then we sum E j from j = 1 to j = n. After some algebraic manipulations, we obtain for any small ϵ0 > 0, for some constant c > 0, for n = 1, . . . , N , ∥en ∥2H + k
n
∥e j ∥2V ≤ c ∥e0 ∥2H + c ∥u0 − u0h ∥2V
j=1
+ c B02 + c k (B1 + B2 + B4 + B5 + I ) + c B3 Me + c k (A4 + A5 ) + c k F n n 2 h +ck (α hk ) + c k β hk j j ∥v j − w j ∥V j=1
+ ϵ0 k
n−1 j=0
j=1
∥ε j ∥2F + c k
j n−1 k ∥ei ∥2V , j=1
i=1
10
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
where for n = 1, . . . , N : tn n−1 hk hk B(tn − tm ) ε(um ) B(tn − s) ε(u(s)) ds − k αn := 0
βnhk
H
m=0
n−1 hk B(tn − tm ) ε(um ) . = k H
m=0
From the regularities (3.4)–(3.7), after some upper bound computations, we obtain for n = 1, . . . , N : k
n
2 2 (α hk j ) ≤ ck +ck
j=1
n−1
hk ∥um − um ∥2V
m=0
≤ c k 2 + c ∥u0 − u0h ∥2V + c k I + c k
j n−1 k ∥ei ∥2V ; j=1
i=1
and k
n
h β hk j ∥v j − w j ∥V ≤ c k
j=1
n−1
hk ∥um − um ∥2V + c k (B2 + Bˆ 2 )
m=0
≤
c ∥u0 − u0h ∥2V
j n−1 2 ∥ei ∥V + c k (B2 + Bˆ 2 ); +ck I +ck k j=1
i=1
and k I ≤ c k2; k A1 ≤ c k 2 ; k B1 ≤ c k 2 ; k A4 ≤ c k 2 ; k B4 ≤ c k 2 ; A5 + B5 + F ≤ c Bˆ 2 . Using then (3.14), we deduce for some constant c > 0, for n = 1, . . . , N , ∥en ∥2H + k
n
∥e j ∥2V ≤ c ∥e0 ∥2H + c ∥u0 − u0h ∥2V + c ∥ε0 ∥2F + c B02 + c A20
j=1
+ c k (B2 + Bˆ 2 ) + c k A2 + c k 2 + c B3 Me + c A3 Mε + c k
j n−1 2 k ∥ei ∥V . j=1
i=1
Using then Gronwall’s inequality (see for example the book [9] or the thesis [3] p.115 Lemma III.3.6), and again the estimation (3.14), we conclude that for some constant c > 0, and for n = 1, . . . , N , n n max ∥en ∥2H + k ∥e j ∥2V ; ∥εn ∥2F + k ∥ε j ∥2E j=1
≤
c ∥e0 ∥2H
j=1
+ c ∥u0 − u0h ∥2V
+ c ∥ε0 ∥2F
+ c B02 + c A20 + c k (B2 + Bˆ 2 ) + c k A2 + c k 2 + c B32 + c A23 .
This gives the estimation (3.13) stated in Theorem 3.2. The Gronwall inequality has been used in the following form (recall that k = T /N ). Let c, g0 , p1 , . . . , p N be non negative real numbers verifying pn ≤ cg0 + c
n−1
kp j ,
n = 1, . . . , N
j=1
(with the convention that for n = 1, we put
n−1 j=1
kp j = 0).
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
11
Then pn ≤ c(1 + cT ecT )g0 ,
n = 1, . . . , N .
The inequality (3.13) is a basis for error estimates for particular choice of the finite-dimensional subspace V h and under additional data and solution regularities. Let us consider Ω ⊂ Rd , d ∈ N∗ , a polygonal domain and V := H 1 (Ω )d ,
H := L 2 (Ω )d ,
E := H 1 (Ω ),
F := L 2 (Ω ).
Let T h be a regular finite element partition of Ω . Let V h ⊂ V and E h ⊂ E be the finite element space consisting of piecewise polynomials of degree ≤m −1, with m ≥ 2, according to the partition T h . Denote by ΠVh : H m (Ω )d → V h and Π Eh : H m (Ω ) → E h the corresponding finite element interpolation operator. Recall (see e.g. [4]) that: ∥w − ΠVh w∥ H l (Ω )d ≤ c h m−l |w| H m (Ω )d , ∀ w ∈ H m (Ω )d ; ∥η − Π Eh η∥ H l (Ω ) ≤ c h m−l |η| H m (Ω ) ,
∀ η ∈ H m (Ω ),
where l = 0 (for which H 0 = L 2 ) or l = 1. In the following we assume the additional data and solution regularities α+1 (Ω )d ; u0 ∈ H v ∈ C([0, T ]; H 2α+1 (Ω )d ), v˙ ∈ L 1 (0, T ; H α (Ω )d ); α+1 θ ∈ C([0, T ]; H (Ω )), θ˙ ∈ L 1 (0, T ; H α (Ω )).
(3.15)
Here α = m − 1 ≥ 1. We remark that for α = 1, the hypothesis (3.15) holds with the assumptions of Theorem 2.1 and v ∈ C([0, T ]; H 3 (Ω )d )
and
θ ∈ C([0, T ]; H 2 (Ω )).
Then we choose the elements u0h = ΠVh u0 ,
v0h = ΠVh v0 ,
whj = ΠVh v j ,
ηhj = Π Eh θ j ,
θ0h = Π Eh θ0 ,
and j = 1 · · · N.
From the assumptions (3.15), we have: ∥u0 − u0h ∥V ≤ c h α ,
∥e0 ∥ H ≤ c h 2α+1 ,
A0 ≤ c h α+1 ,
B0 ≤ c h 2α+1 ;
A3 ≤ c h α+1 ,
B3 ≤ c h 2α+1 ;
k A2 ≤ c h 2α ,
k B2 ≤ c h 4α ,
∥ε0 ∥ F ≤ c h α+1 ;
k Bˆ 2 ≤ c h 2α .
We have the following error estimate result. Theorem 3.3. We keep the assumptions of Theorem 3.2. Under the additional assumption (3.15), we obtain the error estimate for the corresponding discrete solution vnhk , θnhk , n = 1, . . . , N . 1/2 N max ∥vn − vnhk ∥ H + k ∥vn − vnhk ∥2V + max ∥θn − θnhk ∥ F 0≤n≤N
n=0
1/2 N hk 2 + k ∥θn − θn ∥ E ≤ c (h 2α + k). n=0
0≤n≤N
12
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
In particular, for α = 1, we have 1/2 1/2 N N hk hk hk 2 hk 2 + max ∥θn − θn ∥ F + k ≤ c (h 2 + k). max ∥vn − vn ∥ H + k ∥vn − vn ∥V ∥θn − θn ∥ E 0≤n≤N
0≤n≤N
n=0
n=0
4. Numerical computations Here we consider a typical example of thermal bilateral contact problem with Tresca’s friction law (see [1]). We provide numerical simulations derived from the discrete schemes in Section 3. The contact condition on Γ3 is bilateral, and satisfies (see e.g. [6,19]): u ν = 0, |σ τ | ≤ g, |σ τ | < g =⇒ u˙ τ = 0, on Γ3 × (0, T ). |σ τ | = g =⇒ u˙ τ = −λσ τ , for some λ ≥ 0, Here g represents the friction bound, i.e., the magnitude of the limiting friction traction at which slip begins, with g ∈ L ∞ (Γ3 ), g ≥ 0 a.e. on Γ3 . We deduce the admissible displacement space: V := {w ∈ H1 ; with w = 0 on Γ1 ; wν = 0 on Γ3 }, and the sub-differential contact function: ϕ(x, y) = g(x)|yτ (x) |
∀x ∈ Γ3 , y ∈ Rd ,
where yτ (x) := y − yν(x) ν(x), yν(x) := y · ν(x), with ν(x) the unit normal at x ∈ Γ3 . We have then ψ(v) := g |vτ | da, ∀v ∈ V Γ3
is well defined on V , and ψ is Lipschitz continuous on V (see [1]). For computations, we consider a rectangular open set, linear elastic and long memory viscoelastic operators, with a very tiny Γ1 . Ω = (0, L 1 ) × (0, L 2 );
Γ1 = {0} × [0, ϵ];
Γ2 = ({0}×]ϵ, L 2 ]) ∪ ([0, L 1 ] × {L 2 }) ∪ ({L 1 } × [0, L 2 ]); Γ3 =]0, L 1 [×{0}; EY κ EY (G τ )i j = (τ11 + τ22 ) δi j + τi j , 1 ≤ i, j ≤ 2, τ ∈ S2 ; 1+κ 1 − κ2 (A τ )i j = µ (τ11 + τ22 ) δi j + η τi j , 1 ≤ i, j ≤ 2, τ ∈ S2 ; (B(t) τ )i j = B1 (t) (τ11 + τ22 ) δi j + B2 (t) τi j ,
1 ≤ i, j ≤ 2, τ ∈ S2 , t ∈ [0, T ].
Here E Y is the Young modulus, κ the Poisson ratio of the material and δi j denotes the Kronecker symbol, and µ and η are viscosity constants. We use spaces of continuous piecewise affine functions V h ⊂ V and E h ⊂ E as families of approximating subspaces. For computations we considered the following data (IS unity): L 1 = L 2 = 1,
T =1
µ = 10, η = 10, f0 (x, t) = (0, −t), ci j = ki j = ke = 1, B1 (t) = B2 (t) = 10 u0 = (0, 0),
E Y = 2, κ = 0.1 f2 (x, t) = (1, 0), ∀t ∈ [0, T ] 1 ≤ i, j ≤ 2, q = 1 e ,
−2 −t
v0 = (0, 0),
∀t ∈ [0, T ] θ0 = 0.
In Fig. 1 is represented the initial configuration. Then we show in Fig. 2 the deformed configurations at final time, in the case of visco-elasticity of long memory, where the relaxation coefficients are positive decreasing, for two different types of Tresca’s friction bounds. For small friction bound, where g(x, 0) = x2 , 0 ≤ x ≤ 1, we observe on the contact
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
13
Fig. 1. Initial configuration.
Fig. 2. Long memory deformed configurations at final time, θ R (t) = 0, 0 ≤ t ≤ 1.
Fig. 3. Von Mises’ norm in long memory deformed configurations, θ R (t) = 0, 0 ≤ t ≤ 1.
surface a slip phenomena in the direction of the surface fraction on Γ2 : which means that the friction bound has been obtained, in the zone of the values of x near to 1. Whereas for large friction bound, e.g. for g(x, 0) = 20 x, 0 ≤ x ≤ 1, then slip in the direction of the traction may be harder to be realized and frictions appear more important. In Fig. 3, we compute the Von Mises norm. We observe that the norm is important in the neighborhood of the point (0, 0) for small friction bound, and in the neighborhood of the point (1, 0) for large friction bound. In Fig. 4, we show the influence of the different temperatures of the foundation on the temperature field of the body. We observe larger deformations of the body, for greater temperature of the foundation.
14
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
Fig. 4. Temperature field at final time in long memory deformed configurations, g(x, 0) = x2 , 0 ≤ x ≤ 1.
References [1] S. Adly, O. Chau, M. Rochdi, Solvability of a class of dynamical thermal contact problems with subdifferential conditions, Numer. Algebra Control Optim. 2 (1) (2012) 89–101. [2] B. Awbi, O. Chau, Quasistatic thermovisoelastic frictional contact problem with damped response, Appl. Anal. 83 (6) (2004) 635–648. [3] O. Chau, Analyse variationnelle et num´erique de quelques probl`emes aux limites en m´ecanique du contact (Ph.D. thesis), Perpignan, France, 2000. [4] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978. [5] L. De Lorenzis, I. Temizer, P. Wriggers, G. Zavarise, A large deformation frictional contact formulation using NURBS based isogeometric analysis, Internat. J. Numer. Methods Engrg. 87 (13) (2011) 1278–1300. [6] G. Duvaut, J.L. Lions, Les In´equations en M´ecanique et en Physique, Dunod, Paris, 1972. [7] C. Eck, Existence of solutions to a thermo-viscoelastic contact problem with Coulomb friction, Math. Models Methods Appl. Sci. 12 (10) (2002) 1491–1511. [8] D. Goeleven, D. Motreanu, Y. Dumont, M. Rochdi, Variational and Hemivariational Inequalities, Theory Methods and Applications, Volume I: Unilateral Analysis and Unilateral Mechanics, Kluwer Academic Publishers, 2003. [9] W. Han, M. Sofonea, Quasistatic contact problems in viscoelasticity and viscoplasticity, in: Studies in Advanced Mathematics, vol. 30, American Mathematical Society, Providence RI, International Press, Somerville MA, 2002. [10] S. Hueber, B.I. Wohlmuth, Thermo-mechanical contact problems on nonmatching meshes, Comput. Methods Appl. Mech. Engrg. 198 (15) (2009) 1338–1350. [11] R. Jeantet, T. Croguennec, P. Schuck, G. Brul´e, Science des Aliments, Lavoisier, Paris, 2006. [12] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity, SIAM, Philadelphia, 1988. [13] J.L. Lions, Quelques M´ethodes de R´esolution des Probl`emes aux Limites Non Lin´eaires, Dunod et Gauthier-Villars, 1969. [14] A. Matei, M. Sofonea, Variational inequalities with applications: A study of antiplane frictional contact problems, in: Adv. Mech. Math., vol. 18, Springer, 2009. [15] S. Migorski, Evolution hemivariational inequality for a class of dynamic viscoelastic non-monotone frictional contact problems, Comput. Math. Appl. 52 (5) (2006) 677–698. [16] S. Migorski, A. Ochal, M. Sofonea, Nonlinear inclusions and hemivariational inequalities. Models and analysis of contact problems, in: Advances in Mechanics and Mathematics, vol. 26, Springer, New York, 2013. [17] J. Neˇcas, I. Hlav´acˇ ek, Mathematical Theory of Elastic and Elasticoplastic Bodies: An Introduction, Elsevier, Amsterdam, 1981. [18] J. Nedoma, J. Stehlik, Mathematical and computational methods and algorithms in biomechanics, in: Human Skeletal Systems, vol. 12, John Wiley and Sons, 2011. [19] P.D. Panagiotopoulos, Inequality Problems in Mechanics and Applications, Birkh¨auser, Basel, 1985. [20] P.D. Panagiotopoulos, Hemivariational Inequalities, Applications in Mechanics and Engineering, Springer-Verlag, 1993. [21] M. Shillor, M. Sofonea, J.J. Telega, Models and analysis of quasistatic contact, in: Lecture Notes in Physis, vol. 655, Springer, Berlin, 2004. [22] M. Sofonea, W. Han, M. Shillor, Analysis and approximation of contact problems with adhesion or damage, in: Pure Appl. Math. (Boca Raton), vol. 276, Chapman & Hall/CRC, Boca Raton, 2006.
ScienceDirect Mathematics and Computers in Simulation (
)
– www.elsevier.com/locate/matcom
Original articles
Numerical treatment of a class of thermal contact problems Oanh Chau, Rachid Oujja ∗ Universit´e de La R´eunion, Laboratoire PIMENT EA4518, 97715 Saint-Denis Messag cedex 9 La R´eunion, France Received 13 March 2014; received in revised form 9 December 2014; accepted 22 December 2014
Abstract A mathematical problem modeling dynamical behavior of viscoelastic materials is studied. The problem is related to highly nonlinear and non-smooth phenomena like contact, friction and a class of dynamic contact problems with friction and thermal effects. The applied approach leads to a system involving variational inequalities and a fully discrete scheme for numerical approximations and the analysis of error estimate order is provided. Numerical computations are given so as to illustrate the mathematical model. c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Viscoelastic materials; Sub-differential contact condition; Heat transfer; Friction; Numerical simulations
1. Introduction Due to their growing complexity in the engineering sciences, contact problems for viscoelastic materials still remain a serious challenge. For example, contact phenomena are omnipresent in human skeletal systems, see for example [18]. Otherwise, many food materials used in process engineering are viscoelastic, see e.g. [11], and consequently, mathematical models can be very helpful in understanding various problems related to the product development, packing, transport, shelf life testing, thermal effects and heat transfer. It is thus important to study mathematical models that can be used to describe the dynamical behavior of a given viscoelastic material subjected to various highly nonlinear and even non-smooth phenomena like contact, friction and thermal effects. The earlier studies relating to contact problems in the framework of variational formulation and infinitesimal deformations could be found in the reference works [6,12,19]. Thermoviscoelastic contact problems by taking into account the evolution of the temperature parameter could be found in [2,7,10]. Since considerable progress and development have been made in contact mechanics, as widely attested the pioneering works in [9,14,22,21], including mathematical modeling, mathematical analysis, with specific studies on the error estimates analysis, numerical approximation and numerical simulations. In these works the authors show the cross-fertilization between numerous frictional new models arising in contact mechanics and various types of abstract variational inequalities. Further extensions to non convex contact conditions with non-monotone and possible multi-valued constitutive laws led to the active domain of non-smooth mechanics within the framework of the so-called hemivariational ∗ Corresponding author.
E-mail addresses: [email protected] (O. Chau), [email protected] (R. Oujja). http://dx.doi.org/10.1016/j.matcom.2014.12.007 c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.
2
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
inequalities. For a mathematical as well as mechanical treatment we refer to [8,15,16,20]. A large deformation formulation for frictional contact problem has been treated in [5]. This work constitutes in some sense a continuation paper of the results obtained in [1]. The work in [1] has been devoted to a mathematical study of a class of dynamical long memory viscoelastic thermal problems where the contact is governed by a general sub-differential condition. Qualitative results like existence and uniqueness result of weak solutions on displacement and temperature fields have been proved but no numerical approximations have been performed. Here we follow the latter work and propose a numerical scheme for the approximation of the solution fields so as to elaborate a general numerical analysis of error estimates. The theoretical results are then illustrated by different numerical simulations. The paper is organized as follows. In Section 2 we give a short description of the mathematical model and recall the main existence and uniqueness result. In Section 3, we introduce a fully discrete approximation scheme, derive and prove an optimal order error estimate under certain solution regularity assumptions. Finally, in Section 4, we present a set of numerical simulations showing the evolution of the displacement field, the temperature field as well as the Von Mises stress norm. 2. The contact problems
A viscoelastic body occupies the domain Ω with surface Γ that is partitioned into three disjoint measurable parts Γ1 , Γ2 and Γ3 , such that meas (Γ1 ) > 0. Let [0, T ] be the time interval of interest, where T > 0. The body is clamped on Γ1 × (0, T ) and therefore the displacement field vanishes there. We also assume that a volume force of density f0 acts in Ω × (0, T ) and that surface traction of density f2 acts on Γ2 × (0, T ). The body may come in contact with an obstacle, the foundation, over the potential contact surface Γ3 . The model of the contact is specified by a general sub-differential boundary condition, where thermal effects may occur in the frictional contact with the basis. We are interested in the dynamic evolution of the body. We denote by Sd the space of second order symmetric tensors on Rd (d = 1, 2, 3) in practical case, while “ · ” and |·| will represent the inner product and the Euclidean norm on Sd and Rd . Let Ω ⊂ Rd be a bounded domain with a Lipschitz boundary Γ and let ν denote the unit outer normal on Γ . Everywhere in the sequel the indexes i and j run from 1 to d, summation over repeated indices is implied and the index that follows a comma represents the partial derivative with respect to the corresponding component of the independent variable. We also use the following notation: d H = L 2 (Ω ) , H = {σ = (σi j )|σi j = σ ji ∈ L 2 (Ω ), 1 ≤ i, j ≤ d}, H1 = {u ∈ H |ε(u) ∈ H},
H1 = {σ ∈ H|Div σ ∈ H }.
Here ε : H1 −→ H and Div : H1 −→ H are the deformation and the divergence operators, respectively, defined by: 1 (u i, j + u j,i ), Div σ = (σi j, j ). 2 The spaces H , H, H1 and H1 are real Hilbert spaces endowed with the canonical inner products given by: (u, v) H = u i vi d x, (σ , τ )H = σi j τi j d x, ε(u) = (εi j (u)),
Ω
εi j (u) =
(u, v) H1 = (u, v) H + (ε(u), ε(v))H ,
Ω
(σ , τ )H1 = (σ , τ )H + (Div σ , Div τ ) H .
We recall that C denotes the class of continuous functions and C m , m ∈ N∗ is the set of m times differentiable functions. Finally D(Ω ) denotes the set of infinitely differentiable real functions with compact support in Ω ; and W m, p ,
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
3
m ∈ N, 1 ≤ p ≤ +∞ for the classical Sobolev spaces; and H0m (Ω ) := {w ∈ W m,2 (Ω ), w = 0 on Γ1 },
m ≥ 1.
The contact mechanical problem is now formulated. Problem. Q: Find a displacement field u : Ω × [0, T ] −→ Rd and a stress field σ : Ω × [0, T ] −→ Sd and a temperature field θ : Ω × [0, T ] −→ R+ such that for a.e. t ∈ (0, T ): t σ (t) = Aε(u(t)) ˙ + Gε(u(t)) + B(t − s) ε(u(s)) ds − θ (t) Ce in Ω (2.1) 0
u(t) ¨ = Div σ (t) + f0 (t) u(t) = 0
in Ω
(2.2)
on Γ1
σ (t)ν = f2 (t)
(2.3)
on Γ2
(2.4)
ϕ(w) − ϕ(u(t)) ˙ ≥ −σ (t)ν · (w − u(t)) ˙ ∂ u˙ i θ˙ (t) − div(K c ∇θ (t)) = −ci j (t) + q(t) in Ω ∂ xj
u(t) ∈ U,
−ki j
∂θ (t) n i = ke (θ (t) − θ R ) ∂ xj
θ (t) = 0 θ (0) = θ0 u(0) = u0 ,
∀w ∈ U on Γ3
on Γ3
(2.5) (2.6) (2.7)
on Γ1 ∪ Γ2
(2.8)
in Ω
(2.9)
u(0) ˙ = v0
in Ω .
(2.10)
Here, (2.1) is the long memory thermo-visco-elastic constitutive law of the body, σ the stress tensor, A is the viscosity operator, G for the elastic operator, Ce := (ci j ) represents the thermal expansion tensor, and B is the so called tensor of relaxation which defines the long memory of the material, as an important particular case, when B ≡ 0, we find again the usual visco-elasticity of short memory. Eq. (2.2) is the equation of motion where the mass density ϱ ≡ 1. The equation in (2.3) is the clamped condition and in (2.4) is the traction condition. On the contact surface, the general relation (2.5) is a sub-differential boundary condition, verifying D(Ω )d ⊂ U ⊂ H1 , where U represents the set of contact admissible test functions and D(Ω )d is the distribution space. σ ν denotes the Cauchy stress vector on the contact boundary and ϕ : Γ3 × Rd −→ R is a given function. Various situations may be modeled by such a condition, and a concrete example will be given below. The differential equation (2.6) describes the evolution of the temperature field, where K c := (ki j ) represents the thermal conductivity tensor, q(t) the density of volume heat sources. The associated temperature boundary condition is given by (2.7)–(2.8), where θ R is the temperature of the foundation, and ke is the heat exchange coefficient between the body and the obstacle. Finally in (2.9)–(2.10), u0 , v0 , θ0 represents the initial displacement, velocity and temperature, respectively. To derive the variational formulation of the mechanical problems (2.1)–(2.10) we need additional notations. Thus, let V denote the closed subspace of H1 defined by D(Ω )d ⊂ V = {v ∈ H1 | v = 0 on Γ1 } ∩ U ; E = {η ∈ H 1 (Ω ), η = 0 on Γ1 ∪ Γ2 },
F = L 2 (Ω ).
Since meas Γ1 > 0, Korn’s inequality holds: there exists C K > 0 which depends only on Ω and Γ1 such that ∥ε(v)∥H ≥ C K ∥v∥ H1
∀v ∈ V.
A proof of Korn’s inequality may be found in [17, p. 79]. On V we consider the inner product given by (u, v)V = (ε(u), ε(v))H
∀ u, v ∈ V,
4
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
and let ∥ · ∥V be the associated norm, i.e. ∥v∥V = ∥ε(v)∥H
∀ v ∈ V.
It follows that ∥ · ∥ H1 and ∥ · ∥V are equivalent norms on V and therefore (V, ∥ · ∥V ) is a real Hilbert space. Moreover, by Sobolev’s trace theorem, we have a constant C0 > 0 depending only on Ω , Γ1 , and Γ3 such that ∥v∥ L 2 (Γ3 ) ≤ C0 ∥v∥V
∀ v ∈ V.
In the study of the mechanical problem (2.1)–(2.10), we assume the following assumptions (see e.g. [1,13]). The viscosity operator A : Ω × Sd −→ Sd , (x, τ ) −→ (ai jkh (x) τkh ) is linear on the second variable and satisfies the usual properties of ellipticity and symmetry, i.e. (i) ai jkh ∈ W 1,∞ (Ω ); (ii) Aσ · τ = σ · Aτ ∀σ , τ ∈ Sd , a.e. in Ω ; (iii) there exists m A > 0 such that Aτ · τ ≥ m A |τ |2 ∀τ ∈ Sd , a.e. in Ω .
(2.11)
The elasticity operator G : Ω × Sd −→ Sd , (x, τ ) −→ (bi jkh (x) τkh ) is linear on the second variable, with bi jkh ∈ L ∞ (Ω ).
(2.12)
The relaxation tensor B : [0, T ] × Ω × Sd −→ Sd , (t, x, τ ) −→ (Bi jkh (t, x) τkh ) satisfies (i) Bi jkh ∈ W 1,∞ (0, T ; L ∞ (Ω )); (ii) B(t)σ · τ = σ · B(t)τ ∀σ , τ ∈ Sd , a.e. t ∈ (0, T ), a.e. in Ω .
(2.13)
We suppose the body forces and surface tractions satisfy f0 ∈ W 1,2 (0, T ; H ),
f2 ∈ W 1,2 (0, T ; L 2 (Γ2 )d ).
(2.14)
For the thermal tensors and the heat sources density, we suppose that Ce = (ci j ),
ci j = c ji ∈ L ∞ (Ω ),
q ∈ W 1,2 (0, T ; L 2 (Ω )).
(2.15)
The boundary thermic data satisfy ke ∈ L ∞ (Ω ; R+ ),
θ R ∈ W 1,2 (0, T ; L 2 (Γ3 )).
(2.16)
The thermal conductivity tensor verifies the usual symmetry and ellipticity. For some m K > 0 and for all (ξi ) ∈ Rd , K c = (ki j ),
ki j = k ji ∈ L ∞ (Ω ),
ki j ξi ξ j ≥ m K ξi ξi .
(2.17)
We assume that the initial data satisfy the conditions v0 ∈ V ∩ H02 (Ω )d ,
u0 ∈ V,
θ0 ∈ E ∩ H02 (Ω ).
On the contact surface, the following frictional contact function ψ(w) := ϕ(w) da Γ3
(2.18)
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
verifies (i) ψ : V −→ R is well defined, continuous and convex; (ii) there exists a sequence of differentiable convex functions (ψ ) : V −→ R such that ∀w ∈ L 2 (0, T ; V ), Tn T ψn (w(t)) dt −→ ψ(w(t)) dt, for n −→ +∞; 0
0
(iii) for all sequences (wn ) and w in W 1,2 (0, T ; V ) such that wn ⇀ w, w′n ⇀ w′ weakly in L 2 (0, T ; V ), T T then lim inf ψ (w (t)) dt ≥ ψ(w(t)) dt; n n n−→+∞ 0 0 (iv) ∀w ∈ V, (w = 0 on Γ3 =⇒ ∀n ∈ N, ψn′ (w) = 0V ′ ).
5
(2.19)
Here ψn′ (v) denotes the Fr´echet derivative of ψn at v. By using Green’s formula, we obtain the variational formulation of the mechanical problem Q in the abstract form as follows. Problem. QV : Find u : [0, T ] → V , θ : [0, T ] → E satisfying a.e. t ∈ (0, T ): ⟨u(t) ¨ + A u(t) ˙ + B u(t) + C θ (t), w − u(t)⟩ ˙ V ′ ×V t + ψ(w) − ψ(u(t)) ˙ + B(t − s) ε(u(s)) ds, ε(w) − ε(u(t)) ˙ 0
≥ ⟨f (t), w − u(t)⟩ ˙ ∀w ∈ V ; V ′ ×V ˙ + Q(t) in E ′ ; θ˙ (t) + K θ (t) = R u(t) u(0) = u0 , u(0) ˙ = v0 , θ (0) = θ0 .
H
Here the operators and functions A, B : V −→ V ′ , C : E −→ V ′ , K : E −→ E ′ , R : V −→ E ′ , f : [0, T ] −→ V ′ , and Q : [0, T ] −→ E ′ are defined by ∀v ∈ V , ∀w ∈ V , ∀τ ∈ E, ∀η ∈ E: ⟨A v, w⟩V ′ ×V = (A(εv), εw)H ; ⟨B v, w⟩V ′ ×V = (G(εv), εw)H ; ⟨Cτ, w⟩V ′ ×V = −(τ Ce , εw)H ; ⟨f (t), w⟩V ′ ×V = (f0 (t), w) H + (f2 (t), w)(L 2 (Γ2 ))d ; ′ ⟨Q(t), η⟩ E ×E = ke θ R (t) η da + q(t) η d x; Γ3
⟨K τ, η⟩
E ′ ×E
=
i, j=1
⟨R v, η⟩
Ω
d
∂τ ∂η ke τ · η da; dx + ki j ∂ x j ∂ xi Γ3 Ω
E ′ ×E
ci j
=− Ω
∂vi η d x. ∂x j
Let us recall now the following main mathematical result (see for details [1]): Theorem 2.1. Assume that (2.11)–(2.19) hold. Then there exists a unique solution {u, θ} to problem QV with the regularity: u ∈ W 2,2 (0, T ; V ) ∩ W 2,∞ (0, T ; H ) θ ∈ W 1,2 (0, T ; E) ∩ W 1,∞ (0, T ; F). 3. A numerical scheme Let {u, θ} be the unique solution of the problem QV , and introduce the velocity variable v(t) = u(t), ˙
∀t ∈ [0, T ].
6
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
Then t
u(t) = u0 +
v(s) ds,
∀t ∈ [0, T ].
0
From Theorem 2.1 we see that {v, θ} verify for all t ∈ [0, T ]: ⟨˙v(t) + A v(t) + B u(t) + C θ (t), w − v(t)⟩V ′ ×V t B(t − s) ε(u(s)) ds, ε(w) − ε(u(t)) ˙ + + ψ(w) − ψ(v(t)) 0 H ≥ ⟨f (t), w − v(t)⟩V ′ ×V , ∀ w ∈ V. ⟨θ˙ (t), η⟩ F + ⟨K θ (t), η⟩ E ′ ×E = ⟨Rv(t), η⟩ E ′ ×E + ⟨Q(t), η⟩ E ′ ×E , u(0) = u0 ,
v(0) = v0 ,
(3.1)
∀ η ∈ E.
(3.2)
θ (0) = θ0 ,
(3.3)
with the regularity:
v ∈ W 1,2 (0, T ; V ) ∩ W 1,∞ (0, T ; H ) θ ∈ W 1,2 (0, T ; E) ∩ W 1,∞ (0, T ; F).
(3.4)
We make the following additional assumptions on the solution and contact function: v ∈ W 2,2 (0, T ; H );
(3.5)
θ∈W
(3.6)
2,2
(0, T ; F);
ψ is Lipschitz continuous on V.
(3.7)
Before we continue, it is important to remark that proposition (3.7) is verified in the typical example of sub-differential condition we treated (see Section 4). Now let V h ⊂ V and E h ⊂ E be a family of finite dimensional subspaces, with h > 0 a discretization parameter. We divide the time interval [0, T ] into N equal parts: tn = n k, n = 0, 1, . . . , N , with the time step k = T /N . For a continuous function w ∈ C([0, T ]; X ) with values in a space X , we use the notation wn = w(tn ) ∈ X . Then from (3.1) and (3.2) we introduce the following fully discrete scheme. N N Problem. P hk . Find vhk = {vnhk }n=0 ⊂ V h , θ hk = {θnhk }n=0 ⊂ E h such that
v0hk = v0h ,
θ0hk = θ0h
(3.8)
and for n = 1, . . . , N , vhk − vhk n n−1
hk , wh − vnhk ⟩V ′ ×V + ⟨A vnhk , wh − vnhk ⟩V ′ ×V + ⟨B un−1 n−1 hk h hk hk h hk ′ h hk B(tn − tm ) ε(um ), ε(w ) − ε(vn ) + ⟨C θn−1 , w − vn ⟩V ×V + ψ(w ) − ψ(vn ) + k
k
, wh − vnhk
H
m=0
≥ ⟨fn , w
h
θ hk − θ hk n n−1 k
− vnhk ⟩V ′ ×V , , ηh
F
H
∀w ∈ V . h
h
+ ⟨K θnhk , ηh ⟩ E ′ ×E = ⟨R vnhk , ηh ⟩ E ′ ×E + ⟨Q n , ηh ⟩ E ′ ×E ,
(3.9) ∀ ηh ∈ E h
(3.10)
where hk unhk = un−1 + k vnhk ,
u0hk = u0h .
Here u0h ∈ V h , v0h ∈ V h , θ0h ∈ E h are suitable approximations of the initial values u0 , v0 , θ0 .
(3.11)
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
7
–
The main result of this section is the following. Theorem 3.2. We keep the assumptions of Theorem 2.1. Under the additional assumption (3.5)–(3.6), then for the unique solution vhk ⊂ V h , θ hk ⊂ E h of the discrete problem P hk , we have the following error estimate max
1≤n≤N
∥vn − vnhk ∥2H
1≤n≤N
n=1
c ∥u0 − u0h ∥2V
≤
N N hk 2 hk 2 hk 2 + k ∥vn − vn ∥V + max ∥θn − θn ∥ F + k ∥θn − θn ∥ E + c ∥v0 − v0hk ∥2H
n=1
+ c ∥θ0 − θ0h ∥2F
(3.12)
+ c max ∥vn − wnh ∥ H + c max ∥θn − ηnh ∥2F + c k 1≤n≤N
+c
N −1
1≤n≤N
∥(v j − whj ) − (v j+1 − whj+1 )∥ H
∥v j − whj ∥2V + c k
j=1
2 +c
j=1
+ c k2 + c k
N
N −1
N
∥θ j − ηhj ∥2E
j=1
∥θ j − ηhj − (θ j+1 − ηhj+1 )∥ F
2
j=1 N
∥v j − whj ∥V ,
(3.13)
j=1
where for j = 1, . . . , N , whj ∈ V h , ηhj ∈ E h are arbitrary. Proof. Denote by en = vn − vnhk ,
εn = θn − θnhk ,
0 ≤ n ≤ N,
the numerical solution errors. We begin by the estimate of (εn ). Let fix n = 1, . . . , N . From (3.2) and (3.10), we have hk θnhk − θn−1 θ˙n − , ηh + ⟨K θn − K θnhk , ηh ⟩ E ′ ×E = ⟨R vn − R vnhk , ηh ⟩ E ′ ×E , F k
Writing then θ˙n −
hk θnhk − θn−1
k
= θ˙n −
θn − θn−1 εn − εn−1 + , k k
and replacing ηh by ηnh − θn + εn we obtain ε − ε n n−1 , εn + ⟨K θn − K θnhk , εn ⟩ E ′ ×E F k hk θnhk − θn−1 = ⟨R vn − R vnhk , ηnh − θn + εn ⟩ E ′ ×E − θ˙n − , ηnh − θn F k θn − θn−1 hk h − ⟨K θn − K θn , ηn − θn ⟩ E ′ ×E − θ˙n − , εn . F k Consider the quantity for n = 1, . . . , N , ε − ε n n−1 Ξn := , εn + ⟨K θn − K θnhk , εn ⟩ E ′ ×E . F k We have Ξn ≥
1 ∥εn ∥2F − ∥εn−1 ∥2F + m K ∥εn ∥2E . 2k
∀ ηh ∈ E h .
8
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
Summing then Ξ j from j = 1 to j = n, and after some manipulation, we obtain n n n θ j − θ j−1 1 ˙ 2 ∥ε j ∥2E ≤ c ∥θ j − ηhj ∥2E ∥εn ∥2F − ∥ε0 ∥2F + θ j − +c F 2k k j=1 j=1 j=1 +c
n
∥e j ∥2V + c
j=1
n ε j − ε j−1 , θ j − ηhj . F k j=1
Denote now by Mε := max ∥εn ∥ F , 1≤n≤N
A0 := max ∥θn − ηnh ∥ F , 1≤n≤N
A1 :=
A2 :=
N θ j − θ j−1 ˙ 2 θ j − , F k j=1 N
∥θ j − ηhj ∥2E ,
j=1
A3 :=
N
∥θ j − ηhj − (θ j+1 − ηhj+1 )∥ F .
j=1
We deduce from the last inequality that for some constant c > 0 and for n = 1, . . . , N , ∥εn ∥2F + k
n
∥ε j ∥2E ≤ c ∥ε0 ∥2F + c A20 + c k (A1 + A2 ) + c A3 Mε + c k
j=1
n
∥e j ∥2V .
(3.14)
j=1
We now turn to the estimate of (en ). Let fix n = 1, . . . , N . Using (3.1) where we put t = tn and w = vnhk , and adding to (3.9) where wh = wnh , we have hk hk vn − vn−1 , wnh − vnhk + ⟨A vn , vnhk − vn ⟩V ′ ×V + ⟨A vnhk , wnh − vnhk ⟩V ′ ×V (˙vn , vnhk − vn ) H + k H hk hk + ⟨B un , vnhk − vn ⟩V ′ ×V + ⟨B un−1 , wnh − vnhk ⟩V ′ ×V + ⟨C θn , vnhk − vn ⟩V ′ ×V + ⟨C θn−1 , wnh − vnhk ⟩V ′ ×V
+ Rnhk + ψ(wnh ) − ψ(vn ) ≥ ⟨fn , wnh − vn ⟩V ′ ×V ,
∀ wnh ∈ V h ,
where for n = 1, . . . , N : tn n−1 hk hk Rn = B(tn − s) ε(u(s)) ds − k B(tn − tm ) ε(um ), −ε(en ) 0
H
m=0
n−1 hk + k B(tn − tm ) ε(um ), ε(wnh ) − ε(vn ) . m=0
H
Writing then hk hk vnhk − vn−1 vnhk − vn−1 + ; v˙ n = v˙ n − k k hk hk A vn = A vn − A vn + A vn , we obtain e − e n n−1 , en + ⟨A vn − A vnhk , en ⟩V ′ ×V H k vhk − vhk vn − vn−1 n n−1 ≤ ⟨A vnhk , wnh − vn ⟩V ′ ×V + v˙ n − , −en + , wnh − vn H H k k
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
9
hk hk + ⟨B un − B un−1 , −en ⟩V ′ ×V + ⟨B un−1 , wnh − vn ⟩V ′ ×V hk hk + ⟨C θn − C θn−1 , −en ⟩V ′ ×V + ⟨C θn−1 , wnh − vn ⟩V ′ ×V
+ Rnhk + ψ(wnh ) − ψ(vn ) − ⟨fn , wnh − vn ⟩V ′ ×V . Consider the quantity for n = 1, . . . , N , e − e n n−1 + ⟨A vn − A vnhk , en ⟩V ′ ×V . , en En = H k We have En ≥
1 ∥en ∥2H − ∥en−1 ∥2H + m A ∥en ∥2V . 2k
To continue, we denote in the sequel by Me := max ∥en ∥ H , 1≤n≤N
B0 := max ∥vn − wnh ∥ H , 1≤n≤N
B1 :=
B2 :=
N 2 v j − v j−1 − v˙ j , H k j=1 N
∥v j − whj ∥2V ;
Bˆ 2 :=
j=1
B3 :=
N
∥v j − whj ∥V ,
j=1
N −1
∥(v j − whj ) − (v j+1 − whj+1 )∥ H ,
j=1
I :=
N −1 j=1
B4 :=
N
tj
v−k
0
j 2 vi , V
i=1
∥u j − u j−1 ∥2V ,
A4 :=
j=1
B5 :=
N
N
∥θ j − θ j−1 ∥2F ,
j=1
|⟨B u j−1 , v j − whj ⟩V ′ ×V |,
A5 :=
j=1
F :=
N
N
|⟨C θ j−1 , v j − whj ⟩V ′ ×V |,
j=1
|⟨A v j , v j − whj ⟩V ′ ×V + ψ(v j ) − ψ(whj ) − ⟨f j , v j − whj ⟩V ′ ×V |.
j=1
Then we sum E j from j = 1 to j = n. After some algebraic manipulations, we obtain for any small ϵ0 > 0, for some constant c > 0, for n = 1, . . . , N , ∥en ∥2H + k
n
∥e j ∥2V ≤ c ∥e0 ∥2H + c ∥u0 − u0h ∥2V
j=1
+ c B02 + c k (B1 + B2 + B4 + B5 + I ) + c B3 Me + c k (A4 + A5 ) + c k F n n 2 h +ck (α hk ) + c k β hk j j ∥v j − w j ∥V j=1
+ ϵ0 k
n−1 j=0
j=1
∥ε j ∥2F + c k
j n−1 k ∥ei ∥2V , j=1
i=1
10
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
where for n = 1, . . . , N : tn n−1 hk hk B(tn − tm ) ε(um ) B(tn − s) ε(u(s)) ds − k αn := 0
βnhk
H
m=0
n−1 hk B(tn − tm ) ε(um ) . = k H
m=0
From the regularities (3.4)–(3.7), after some upper bound computations, we obtain for n = 1, . . . , N : k
n
2 2 (α hk j ) ≤ ck +ck
j=1
n−1
hk ∥um − um ∥2V
m=0
≤ c k 2 + c ∥u0 − u0h ∥2V + c k I + c k
j n−1 k ∥ei ∥2V ; j=1
i=1
and k
n
h β hk j ∥v j − w j ∥V ≤ c k
j=1
n−1
hk ∥um − um ∥2V + c k (B2 + Bˆ 2 )
m=0
≤
c ∥u0 − u0h ∥2V
j n−1 2 ∥ei ∥V + c k (B2 + Bˆ 2 ); +ck I +ck k j=1
i=1
and k I ≤ c k2; k A1 ≤ c k 2 ; k B1 ≤ c k 2 ; k A4 ≤ c k 2 ; k B4 ≤ c k 2 ; A5 + B5 + F ≤ c Bˆ 2 . Using then (3.14), we deduce for some constant c > 0, for n = 1, . . . , N , ∥en ∥2H + k
n
∥e j ∥2V ≤ c ∥e0 ∥2H + c ∥u0 − u0h ∥2V + c ∥ε0 ∥2F + c B02 + c A20
j=1
+ c k (B2 + Bˆ 2 ) + c k A2 + c k 2 + c B3 Me + c A3 Mε + c k
j n−1 2 k ∥ei ∥V . j=1
i=1
Using then Gronwall’s inequality (see for example the book [9] or the thesis [3] p.115 Lemma III.3.6), and again the estimation (3.14), we conclude that for some constant c > 0, and for n = 1, . . . , N , n n max ∥en ∥2H + k ∥e j ∥2V ; ∥εn ∥2F + k ∥ε j ∥2E j=1
≤
c ∥e0 ∥2H
j=1
+ c ∥u0 − u0h ∥2V
+ c ∥ε0 ∥2F
+ c B02 + c A20 + c k (B2 + Bˆ 2 ) + c k A2 + c k 2 + c B32 + c A23 .
This gives the estimation (3.13) stated in Theorem 3.2. The Gronwall inequality has been used in the following form (recall that k = T /N ). Let c, g0 , p1 , . . . , p N be non negative real numbers verifying pn ≤ cg0 + c
n−1
kp j ,
n = 1, . . . , N
j=1
(with the convention that for n = 1, we put
n−1 j=1
kp j = 0).
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
11
Then pn ≤ c(1 + cT ecT )g0 ,
n = 1, . . . , N .
The inequality (3.13) is a basis for error estimates for particular choice of the finite-dimensional subspace V h and under additional data and solution regularities. Let us consider Ω ⊂ Rd , d ∈ N∗ , a polygonal domain and V := H 1 (Ω )d ,
H := L 2 (Ω )d ,
E := H 1 (Ω ),
F := L 2 (Ω ).
Let T h be a regular finite element partition of Ω . Let V h ⊂ V and E h ⊂ E be the finite element space consisting of piecewise polynomials of degree ≤m −1, with m ≥ 2, according to the partition T h . Denote by ΠVh : H m (Ω )d → V h and Π Eh : H m (Ω ) → E h the corresponding finite element interpolation operator. Recall (see e.g. [4]) that: ∥w − ΠVh w∥ H l (Ω )d ≤ c h m−l |w| H m (Ω )d , ∀ w ∈ H m (Ω )d ; ∥η − Π Eh η∥ H l (Ω ) ≤ c h m−l |η| H m (Ω ) ,
∀ η ∈ H m (Ω ),
where l = 0 (for which H 0 = L 2 ) or l = 1. In the following we assume the additional data and solution regularities α+1 (Ω )d ; u0 ∈ H v ∈ C([0, T ]; H 2α+1 (Ω )d ), v˙ ∈ L 1 (0, T ; H α (Ω )d ); α+1 θ ∈ C([0, T ]; H (Ω )), θ˙ ∈ L 1 (0, T ; H α (Ω )).
(3.15)
Here α = m − 1 ≥ 1. We remark that for α = 1, the hypothesis (3.15) holds with the assumptions of Theorem 2.1 and v ∈ C([0, T ]; H 3 (Ω )d )
and
θ ∈ C([0, T ]; H 2 (Ω )).
Then we choose the elements u0h = ΠVh u0 ,
v0h = ΠVh v0 ,
whj = ΠVh v j ,
ηhj = Π Eh θ j ,
θ0h = Π Eh θ0 ,
and j = 1 · · · N.
From the assumptions (3.15), we have: ∥u0 − u0h ∥V ≤ c h α ,
∥e0 ∥ H ≤ c h 2α+1 ,
A0 ≤ c h α+1 ,
B0 ≤ c h 2α+1 ;
A3 ≤ c h α+1 ,
B3 ≤ c h 2α+1 ;
k A2 ≤ c h 2α ,
k B2 ≤ c h 4α ,
∥ε0 ∥ F ≤ c h α+1 ;
k Bˆ 2 ≤ c h 2α .
We have the following error estimate result. Theorem 3.3. We keep the assumptions of Theorem 3.2. Under the additional assumption (3.15), we obtain the error estimate for the corresponding discrete solution vnhk , θnhk , n = 1, . . . , N . 1/2 N max ∥vn − vnhk ∥ H + k ∥vn − vnhk ∥2V + max ∥θn − θnhk ∥ F 0≤n≤N
n=0
1/2 N hk 2 + k ∥θn − θn ∥ E ≤ c (h 2α + k). n=0
0≤n≤N
12
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
In particular, for α = 1, we have 1/2 1/2 N N hk hk hk 2 hk 2 + max ∥θn − θn ∥ F + k ≤ c (h 2 + k). max ∥vn − vn ∥ H + k ∥vn − vn ∥V ∥θn − θn ∥ E 0≤n≤N
0≤n≤N
n=0
n=0
4. Numerical computations Here we consider a typical example of thermal bilateral contact problem with Tresca’s friction law (see [1]). We provide numerical simulations derived from the discrete schemes in Section 3. The contact condition on Γ3 is bilateral, and satisfies (see e.g. [6,19]): u ν = 0, |σ τ | ≤ g, |σ τ | < g =⇒ u˙ τ = 0, on Γ3 × (0, T ). |σ τ | = g =⇒ u˙ τ = −λσ τ , for some λ ≥ 0, Here g represents the friction bound, i.e., the magnitude of the limiting friction traction at which slip begins, with g ∈ L ∞ (Γ3 ), g ≥ 0 a.e. on Γ3 . We deduce the admissible displacement space: V := {w ∈ H1 ; with w = 0 on Γ1 ; wν = 0 on Γ3 }, and the sub-differential contact function: ϕ(x, y) = g(x)|yτ (x) |
∀x ∈ Γ3 , y ∈ Rd ,
where yτ (x) := y − yν(x) ν(x), yν(x) := y · ν(x), with ν(x) the unit normal at x ∈ Γ3 . We have then ψ(v) := g |vτ | da, ∀v ∈ V Γ3
is well defined on V , and ψ is Lipschitz continuous on V (see [1]). For computations, we consider a rectangular open set, linear elastic and long memory viscoelastic operators, with a very tiny Γ1 . Ω = (0, L 1 ) × (0, L 2 );
Γ1 = {0} × [0, ϵ];
Γ2 = ({0}×]ϵ, L 2 ]) ∪ ([0, L 1 ] × {L 2 }) ∪ ({L 1 } × [0, L 2 ]); Γ3 =]0, L 1 [×{0}; EY κ EY (G τ )i j = (τ11 + τ22 ) δi j + τi j , 1 ≤ i, j ≤ 2, τ ∈ S2 ; 1+κ 1 − κ2 (A τ )i j = µ (τ11 + τ22 ) δi j + η τi j , 1 ≤ i, j ≤ 2, τ ∈ S2 ; (B(t) τ )i j = B1 (t) (τ11 + τ22 ) δi j + B2 (t) τi j ,
1 ≤ i, j ≤ 2, τ ∈ S2 , t ∈ [0, T ].
Here E Y is the Young modulus, κ the Poisson ratio of the material and δi j denotes the Kronecker symbol, and µ and η are viscosity constants. We use spaces of continuous piecewise affine functions V h ⊂ V and E h ⊂ E as families of approximating subspaces. For computations we considered the following data (IS unity): L 1 = L 2 = 1,
T =1
µ = 10, η = 10, f0 (x, t) = (0, −t), ci j = ki j = ke = 1, B1 (t) = B2 (t) = 10 u0 = (0, 0),
E Y = 2, κ = 0.1 f2 (x, t) = (1, 0), ∀t ∈ [0, T ] 1 ≤ i, j ≤ 2, q = 1 e ,
−2 −t
v0 = (0, 0),
∀t ∈ [0, T ] θ0 = 0.
In Fig. 1 is represented the initial configuration. Then we show in Fig. 2 the deformed configurations at final time, in the case of visco-elasticity of long memory, where the relaxation coefficients are positive decreasing, for two different types of Tresca’s friction bounds. For small friction bound, where g(x, 0) = x2 , 0 ≤ x ≤ 1, we observe on the contact
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
13
Fig. 1. Initial configuration.
Fig. 2. Long memory deformed configurations at final time, θ R (t) = 0, 0 ≤ t ≤ 1.
Fig. 3. Von Mises’ norm in long memory deformed configurations, θ R (t) = 0, 0 ≤ t ≤ 1.
surface a slip phenomena in the direction of the surface fraction on Γ2 : which means that the friction bound has been obtained, in the zone of the values of x near to 1. Whereas for large friction bound, e.g. for g(x, 0) = 20 x, 0 ≤ x ≤ 1, then slip in the direction of the traction may be harder to be realized and frictions appear more important. In Fig. 3, we compute the Von Mises norm. We observe that the norm is important in the neighborhood of the point (0, 0) for small friction bound, and in the neighborhood of the point (1, 0) for large friction bound. In Fig. 4, we show the influence of the different temperatures of the foundation on the temperature field of the body. We observe larger deformations of the body, for greater temperature of the foundation.
14
O. Chau, R. Oujja / Mathematics and Computers in Simulation (
)
–
Fig. 4. Temperature field at final time in long memory deformed configurations, g(x, 0) = x2 , 0 ≤ x ≤ 1.
References [1] S. Adly, O. Chau, M. Rochdi, Solvability of a class of dynamical thermal contact problems with subdifferential conditions, Numer. Algebra Control Optim. 2 (1) (2012) 89–101. [2] B. Awbi, O. Chau, Quasistatic thermovisoelastic frictional contact problem with damped response, Appl. Anal. 83 (6) (2004) 635–648. [3] O. Chau, Analyse variationnelle et num´erique de quelques probl`emes aux limites en m´ecanique du contact (Ph.D. thesis), Perpignan, France, 2000. [4] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978. [5] L. De Lorenzis, I. Temizer, P. Wriggers, G. Zavarise, A large deformation frictional contact formulation using NURBS based isogeometric analysis, Internat. J. Numer. Methods Engrg. 87 (13) (2011) 1278–1300. [6] G. Duvaut, J.L. Lions, Les In´equations en M´ecanique et en Physique, Dunod, Paris, 1972. [7] C. Eck, Existence of solutions to a thermo-viscoelastic contact problem with Coulomb friction, Math. Models Methods Appl. Sci. 12 (10) (2002) 1491–1511. [8] D. Goeleven, D. Motreanu, Y. Dumont, M. Rochdi, Variational and Hemivariational Inequalities, Theory Methods and Applications, Volume I: Unilateral Analysis and Unilateral Mechanics, Kluwer Academic Publishers, 2003. [9] W. Han, M. Sofonea, Quasistatic contact problems in viscoelasticity and viscoplasticity, in: Studies in Advanced Mathematics, vol. 30, American Mathematical Society, Providence RI, International Press, Somerville MA, 2002. [10] S. Hueber, B.I. Wohlmuth, Thermo-mechanical contact problems on nonmatching meshes, Comput. Methods Appl. Mech. Engrg. 198 (15) (2009) 1338–1350. [11] R. Jeantet, T. Croguennec, P. Schuck, G. Brul´e, Science des Aliments, Lavoisier, Paris, 2006. [12] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity, SIAM, Philadelphia, 1988. [13] J.L. Lions, Quelques M´ethodes de R´esolution des Probl`emes aux Limites Non Lin´eaires, Dunod et Gauthier-Villars, 1969. [14] A. Matei, M. Sofonea, Variational inequalities with applications: A study of antiplane frictional contact problems, in: Adv. Mech. Math., vol. 18, Springer, 2009. [15] S. Migorski, Evolution hemivariational inequality for a class of dynamic viscoelastic non-monotone frictional contact problems, Comput. Math. Appl. 52 (5) (2006) 677–698. [16] S. Migorski, A. Ochal, M. Sofonea, Nonlinear inclusions and hemivariational inequalities. Models and analysis of contact problems, in: Advances in Mechanics and Mathematics, vol. 26, Springer, New York, 2013. [17] J. Neˇcas, I. Hlav´acˇ ek, Mathematical Theory of Elastic and Elasticoplastic Bodies: An Introduction, Elsevier, Amsterdam, 1981. [18] J. Nedoma, J. Stehlik, Mathematical and computational methods and algorithms in biomechanics, in: Human Skeletal Systems, vol. 12, John Wiley and Sons, 2011. [19] P.D. Panagiotopoulos, Inequality Problems in Mechanics and Applications, Birkh¨auser, Basel, 1985. [20] P.D. Panagiotopoulos, Hemivariational Inequalities, Applications in Mechanics and Engineering, Springer-Verlag, 1993. [21] M. Shillor, M. Sofonea, J.J. Telega, Models and analysis of quasistatic contact, in: Lecture Notes in Physis, vol. 655, Springer, Berlin, 2004. [22] M. Sofonea, W. Han, M. Shillor, Analysis and approximation of contact problems with adhesion or damage, in: Pure Appl. Math. (Boca Raton), vol. 276, Chapman & Hall/CRC, Boca Raton, 2006.