A generalized finite element method for linear thermoelasticity

# A generalized finite element method for linear thermoelasticity

###### Abstract.

We propose and analyze a generalized finite element method designed for linear quasistatic thermoelastic systems with spatial multiscale coefficients. The method is based on the local orthogonal decomposition technique introduced by Målqvist and Peterseim (Math. Comp., 83(290): 2583–2603, 2014). We prove convergence of optimal order, independent of the derivatives of the coefficients, in the spatial -norm. The theoretical results are confirmed by numerical examples.

11footnotetext: Department of Mathematical Sciences, Chalmers University of Technology and University of Gothenburg, SE-412 96 Göteborg, Sweden.22footnotetext: Supported by the Swedish Research Council and the Swedish Foundation for Strategic Research.

## 1. Introduction

In many applications the expansion and contraction of a material exposed to temperature changes are of great importance. To model this phenomenon a system consisting of an elasticity equation describing the displacement coupled with an equation for the temperature is used, see, e.g., [6]. The full system consists of a hyperbolic elasticity equation coupled with a parabolic equation for the temperature, see [8] for a comprehensive treatment of this formulation. If the inertia effects are negligible, the hyperbolic term in the elasticity equation can be removed. This leads to an elliptic-parabolic system, often referred to as quasistatic. This formulation is discussed in, for instance, [22, 25]. In some settings it is justified to also remove the parabolic term, which leads to an elliptic-elliptic system, see, e.g., [22, 25]. Since the thermoelastic problem is formally equivalent to the system describing poroelasticity, several papers on this equation are also relevant, see, e.g., [5, 24].

In this paper we study the quasistatic case. Existence and uniqueness of a solution to this system are discussed in [22] within the framework of linear degenerate evolution equations in Hilbert spaces. It is also shown that this system is essentially of parabolic type. Existence and uniqueness are also treated in [25] (only two-dimensional problems) and in [23, 21] some results on the thermoelastic contact problem are presented. The classical finite element method for the thermoelastic system is analyzed in [10, 25], where convergence rates of optimal order are derived for problems with solution in or higher.

When the elastic medium of interest is strongly heterogeneous, like composite materials, the coefficients are highly varying and oscillating. Commonly, such coefficients are said to have multiscale features. For these problems classical polynomial finite elements, as in [10, 25], fail to approximate the solution well unless the mesh width resolves the data variations. This is due to the fact that a priori bounds of the error depend on (at least) the spatial -norm of the solution. Since this norm depends on the derivative of the diffusion coefficient, it is of order if the coefficient oscillates with frequency . To overcome this difficulty, several numerical methods have been proposed, see for instance [4, 3, 15, 17, 14].

In this paper we suggest a generalized finite element method based on the techniques introduced in [17], often referred to as local orthogonal decomposition. This method builds on ideas from the variational multiscale method [14, 15], where the solution space is split into a coarse and a fine part. The coarse space is modified such that the basis functions contain information from the diffusion coefficient and have support on small patches. With this approach the basis functions have good approximation properties locally. In [17] the technique is applied to elliptic problems with an arbitrary positive and bounded diffusion coefficient. One of the main advantages is that no assumptions on scale separation or periodicity of the coefficient are needed. Recently, this technique has been applied to several other problems, for instance, semilinear elliptic equations [12], boundary value problems [11], eigenvalue problems [18], linear and semilinear parabolic equations [16], and the linear wave equation [1].

The method we propose in this paper uses generalized finite element spaces similar to those used [17] and [13], together with a correction building on the ideas in [11, 15]. We prove convergence of optimal order that does not depend on the derivatives of the coefficients. We emphasize that by avoiding these derivatives, the a priori bound does not contain any constant of order , although coefficients are highly varying.

In Section 2 we formulate the problem of interest, in Section 3 we first recall the classical finite element method for thermoelasticity and then we define the new generalized finite element method. In Section 4 we perform a localization of the basis functions and in Section 5 we analyze the error. Finally, in Section 6 we present some numerical results.

## 2. Problem formulation

Let , , be a polygonal/polyhedral domain describing the reference configuration of an elastic body. For a given time we let denote the displacement field and the temperature. To impose Dirichlet and Neumann boundary conditions, we let and denote two disjoint segments of the boundary such that . The segments and are defined similarly.

We use to denote the inner product in and for the corresponding norm. Let denote the classical Sobolev space with norm and let denote the dual space to . Furthermore, we adopt the notation for the Bochner space with the norm

 ∥v∥Lp([0,T];X) =(∫T0∥v∥pXdt)1/p,1≤p<∞, ∥v∥L∞([0,T];X) =esssup0≤t≤T∥v∥X,

where is a Banach space equipped with the norm . The notation is used to denote . The dependence on the interval and the domain is frequently suppressed and we write, for instance, for . We also define the following subspaces of

 V1:={v∈(H1(Ω))d:v=0 on ΓuD},V2:={v∈H1(Ω):v=0 on ΓθD}.

Under the assumption that the displacement gradients are small, the (linearized) strain tensor is given by

 ε(u)=12(∇u+∇u⊺).

Assuming further that the material is isotropic, Hooke’s law gives the (total) stress tensor, see e.g. [21] and the references therein,

 ¯σ=2με(u)+λ(∇⋅u)I−αθI,

where is the -dimensional identity matrix, is the thermal expansion coefficient, and and are the so called Lamé coefficients given by

 μ=E2(1+ν),λ=Eν(1+ν)(1−2ν),

where denotes Young’s elastic modulus and denotes Poisson’s ratio. The materials of interest are strongly heterogeneous which implies that , , and are rapidly varying in space.

The linear quasistatic thermoelastic problem takes the form

 (2.1) −∇⋅(2με(u)+λ∇⋅uI−αθI) =f, in (0,T]×Ω, (2.2) ˙θ−∇⋅κ∇θ+α∇⋅˙u =g, in (0,T]×Ω, (2.3) u =0, in (0,T]×ΓuD, (2.4) ¯σ⋅n =0, in (0,T]×ΓuN. (2.5) θ =0, on (0,T]×ΓθD, (2.6) ∇θ⋅n =0, on (0,T]×ΓθN. (2.7) θ(0) =θ0, in Ω,

where is the heat conductivity parameter, which is assumed to be rapidly varying in space.

###### Remark 2.1.

For simplicity we have assumed homogeneous boundary data (2.3)-(2.6). However, using techniques similar to the ones used in [11, 13] the analysis in this paper can be extended to non-homogeneous situations.

###### Assumptions.

We make the following assumptions on the data

1. , symmetric,

 0<κ1:=essinfx∈Ωinfv∈Rd∖{0}κ(x)v⋅vv⋅v,∞>κ2:=esssupx∈Ωsupv∈Rd∖{0}κ(x)v⋅vv⋅v,
2. , and

 0<μ1:=essinfx∈Ωμ(x)≤esssupx∈Ωμ(x)=:μ2<∞.

Similarly, the constants , and are used to denote the corresponding upper and lower bounds for and .

3. , , , and .

To pose a variational form we multiply the equations (2.1) and (2.2) with test functions from and and using Green’s formula together with the boundary conditions (2.3)-(2.6) we arrive at the following weak formulation [10]. Find and , such that,

 (2.8) (σ(u):ε(v1))−(αθ,∇⋅v1) =(f,v1), ∀v1∈V1, (2.9) (˙θ,v2)+(κ∇θ,∇v2)+(α∇⋅˙u,v2) =(g,v2), ∀v2∈V2,

and the initial value is satisfied. Here we use to denote the effective stress tensor and we use to denote the Frobenius inner product of matrices. Using Korn’s inequality we have the following bounds, see, e.g., [7],

 (2.10) cσ∥v1∥2H1≤(σ(v1):ϵ(v1))≤Cσ∥v1∥2H1,∀v1∈V1

where (resp. ) depends on (resp. and ). Similarly, there are constants (resp. ) depending on the bound (resp. ) such that

 (2.11) cκ∥v2∥2H1≤(κ∇v2,∇v2)≤Cκ∥v2∥2H1,∀v2∈V2.

Furthermore, we use the following notation for the energy norms induced by the bilinear forms

 ∥v1∥2σ:=(σ(v1):ε(v1)), v1∈V1,∥v2∥2κ:=(κ∇v2∇v2), v2∈V2

Existence and uniqueness of a solution to (2.8)-(2.9) have been proved in [22, 25]. There are also some papers on the solution to contact problems, see [2, 23].

###### Theorem 2.2.

Assume that 1-3 hold and that is sufficiently smooth. Then there exist and such that , , , and satisfying (2.8)-(2.9) and the initial condition .

###### Remark 2.3.

We remark that the equations (2.1)-(2.7) also describe a poroelastic system. In this case denotes the fluid pressure, the permeability and viscosity of the fluid.

## 3. Numerical approximation

In this section is we first recall some properties of the classical finite element method for (2.8)-(2.9). In subsection 3.2 we propose a new numerical method built on the ideas from [17]. The localization of this method is treated in Section 4.

### 3.1. Classical finite element

First, we need to define appropriate finite element spaces. For this purpose we let { be a family of shape regular triangulations of with the mesh size , for . Furthermore, we denote the largest diameter in the triangulation by . We now define the classical piecewise affine finite element spaces

 V1h ={v∈(C(¯Ω))d:v=0 on ΓuD,v|K is a polynomial of degree≤1,∀K∈Th}, V2h ={v∈C(¯Ω):v=0 on ΓθD,v|K is a polynomial of degree≤1,∀K∈Th}.

For the discretization in time we consider, for simplicity, a uniform time step such that for and .

###### Remark 3.1.

The classical linear elasticity equation can in some cases suffer from locking effects when using continuous piecewise linear polynomials in both spaces (P1-P1 elements). These typically occur if is close to (Poisson locking) or if the thickness of the domain is very small (shear locking). In the coupled time-dependent problem locking can occur if is neglected in (2.2) and P1-P1 elements are used. The locking produces artificial oscillations in the numerical approximation of the temperature (or pressure) for early time steps. However, it shall be noted that in the case when is not neglected, this locking effect does not occur, see [20]. Thus, we consider a P1-P1 discretization in this paper.

The classical finite element method with a backward Euler scheme in time reads; for find and , such that

 (3.1) (σ(unh):ε(v1))−(αθnh,∇⋅v1) =(fn,v1), ∀v1∈V1h, (3.2) (¯∂tθnh,v2)+(κ∇θnh,∇v2)+(α∇⋅¯∂tunh,v2) =(gn,v2), ∀v2∈V2h,

where and similarly for . The right hand sides are evaluated at time , that is, and . Given initial data and the system (3.1)-(3.2) is well posed [10]. We assume that is a suitable approximation of . For we note that is uniquely determined by (2.8) at , that is, fulfills the equation

 (σ(u(0)):ε(v1))−(αθ0,∇⋅v1)=(f0,v1),∀v1∈V1,

and we thus define to be the solution to

 (3.3) (σ(u0h):ε(v1))−(αθ0h,∇⋅v1)=(f0,v1),∀v1∈V1h.

The following theorem is a consequence of [10, Theorem 3.1]. The convergence rate is optimal for the two first norms. However, it is not optimal for the -norm . In [10] this is avoided by using second order continuous piecewise polynomials for the displacement (P2-P1 elements). It is, however, noted that the problem is still stable using P1-P1 elements. In this paper we use P1-P1 elements and derive error bounds in the -norm, of optimal order, for both the displacement and the temperature.

###### Theorem 3.2.

Let be the solution to (2.8)-(2.9) and be the solution to (3.1)-(3.2). Then for

 ∥un−unh∥H1+(n∑m=1τ∥θm−θmh∥2H1)1/2+∥θn−θnh∥≤Cϵ−1(h+τ),

where is of order if the material varies on a scale of size .

Note that the constant involved in this error bound contains derivatives of the coefficients. Hence, convergence only takes place when the mesh size is sufficiently small (). Throughout this paper, it is assumed that is small enough and and are referred to as reference spaces for the solution. Similarly, and are referred to as reference solutions. In Section 5 this solution is compared with the generalized finite element solution. We emphasize that the generalized finite element solution is computed in spaces of lower dimension and hence not as computationally expensive.

In the following theorem we prove some regularity results for the finite element solution.

###### Theorem 3.3.

Let and be the solution to (3.1)-(3.2). Then the following bound holds

 (3.4) (n∑j=1τ∥¯∂tujh∥2H1)1/2 +(n∑j=1τ∥¯∂tθjh∥2)1/2+∥θnh∥H1 ≤C(∥g∥L∞(L2)+∥˙f∥L∞(H−1)+∥θ0h∥H1)

If , then for

 (3.5) ∥¯∂tunh∥H1 +∥¯∂tθnh∥+(n∑j=1τ∥¯∂tθjh∥2H1)1/2

If and , then for

 (3.6) ∥¯∂tunh∥H1+∥¯∂tθnh∥+t1/2n∥¯∂tθnh∥H1≤Ct−1/2n∥θ0h∥H1.
###### Proof.

From (3.1)-(3.2) and the initial data (3.3) we deduce that the following relation must hold for

 (3.7) (σ(¯∂tunh):ε(v1))−(α¯∂tθnh,∇⋅v1) =(¯∂tfn,v1), ∀v1∈V1h, (3.8) (¯∂tθnh,v2)+(κ∇θnh,∇v2)+(α∇⋅¯∂tunh,v2) =(gn,v2), ∀v2∈V2h.

By choosing and and adding the resulting equations we have

 (3.9) cσ2∥¯∂tunh∥2H1+12∥¯∂tθnh∥2+(κ∇θnh,∇¯∂tθnh)≤C(∥gn∥2+∥¯∂tfn∥2H−1).

Note that the coupling terms cancel. By using Cauchy-Schwarz and Young’s inequality we can bound

 τ(κ∇θnh,∇¯∂tθnh)=∥κ1/2∇θnh∥2−(κ∇θnh,∇θn−1h)≥12∥θnh∥2κ−12∥θn−1h∥2κ.

Multiplying (3.9) by , summing over , and using (2.10) gives

 n∑j=1τ∥¯∂tujh∥2H1+n∑j=1τ∥¯∂tθjh∥2+∥θnh∥2H1 ≤Cn∑j=1τ(∥gj∥2+∥¯∂tfj∥2H−1) +C∥θ0h∥H1,

which is bounded by the right hand side in (3.4).

For the bound (3.5) we note that the following relation must hold for

 (3.10) (σ(¯∂tunh):ε(v1))−(α¯∂tθnh,∇⋅v1) =(¯∂tfn,v1), ∀v1∈V1h, (3.11) (¯∂2tθnh,v2)+(κ∇¯∂tθnh,∇v2)+(α∇⋅¯∂2tunh,v2) =(¯∂tgn,v2), ∀v2∈V2h.

Now choose and and add the resulting equations to get

 (σ(¯∂tunh):ε(¯∂2tunh))+(¯∂2tθnh,¯∂tθnh) +(κ∇¯∂tθnh,∇¯∂tθnh) =(¯∂tfn,¯∂2tunh)+(¯∂tgn,¯∂tθnh).

Multiplying by and using Cauchy-Schwarz and Young’s inequality gives

 12∥¯∂tunh∥2σ+12∥¯∂tθnh∥2+Cτ∥¯∂tθnh∥2H1 ≤12∥¯∂tθn−1h∥2+12∥¯∂tun−1h∥2σ +τ(¯∂tfn,¯∂2tunh)+C∥¯∂tgn∥2H−1.

Summing over and using (2.10) now gives

 ∥¯∂tunh∥2H1+∥¯∂tθnh∥2+n∑j=2τ∥¯∂tθjh∥2H1 ≤C(∥¯∂tu1h∥2H1+∥¯∂tθ1h∥2 +n∑j=2τ((¯∂tfj,¯∂2tujh)+∥¯∂tgj∥2H−1)).

Here we use summation by parts to get

 n∑j=2τ(¯∂tfj,¯∂2tujh) =(¯∂tfn,¯∂tunh)−(¯∂tf1,¯∂tu1h)−n∑j=2τ(¯∂2tfj,¯∂tuj−1h) ≤C(max1≤j≤n∥¯∂tfj∥H−1+n∑j=2τ∥¯∂2tfj∥H−1)max1≤j≤n∥¯∂tujh∥H1,

and can now be kicked to the left hand side.

To estimate and we choose and in (3.7)-(3.8) for . We thus have, since ,

 ∥¯∂tu1h∥2H1+∥¯∂tθ1h∥2+1τ∥θ1h∥2H1≤C(∥¯∂tf1∥2H−1+∥g1∥2).

The observation that completes the bound (3.5).

Now assume and and note that the following holds for ,

 (σ(¯∂2tunh):ε(v1))−(α¯∂2tθnh,∇⋅v1) =0, ∀v1∈V1h, (¯∂2tθnh,v2)+(κ∇¯∂tθnh,∇v2)+(α∇⋅¯∂2tunh,v2) =0, ∀v2∈V2h.

Choosing , and adding the resulting equations gives

 (σ(¯∂2tunh):ε(¯∂2tunh))+(¯∂2tθnh,¯∂2tθnh)+(κ∇¯∂tθnh,∇¯∂2tθnh)=0,

where, again, the coupling terms cancel. The two first terms on the left hand side are positive and can thus be ignored. Multiplying by and gives after using Cauchy-Schwarz and Young’s inequality

 t2n∥¯∂tθnh∥2κ≤t2n−1∥¯∂tθn−1h∥2κ+(t2n−t2n−1)∥¯∂2tθn−1h∥2κ.

Note that , where we use that if . Summing over now gives

 t2n∥¯∂tθnh∥2κ≤t21∥¯∂tθ1h∥2κ+3n∑j=2τtj−1∥¯∂tθj−1h∥2κ.

To bound the last sum we choose , in (3.10)-(3.11), now with and . Adding the resulting equations gives

 (¯∂2tθnh,¯∂tθnh)+(κ∇¯∂tθnh,∇¯∂tθnh)+(σ(¯∂tunh):ε(¯∂2tunh))=0,

Multiplying by and gives after using Cauchy-Schwarz inequality

 tn2 ∥¯∂tunh∥2σ+tn2∥¯∂tθnh∥2+cκτtn∥¯∂tθnh∥2H1 ≤tn−12∥¯∂tun−1h∥2σ+tn−12∥¯∂tθn−1h∥2+τ2∥¯∂tun−1h∥2σ+τ2∥¯∂tθn−1h∥2.

Summing over and using (2.10) thus gives

 cσtn2∥¯∂tunh∥2H1+tn2∥¯∂tθnh∥2+n∑j=2τtj∥¯∂tθjh∥2H1 ≤Cσt12∥¯∂tu1h∥2H1+t12∥¯∂tθ1h∥2+Cn∑j=2τ(∥¯∂tuj−1h∥2H1+∥¯∂tθj−1h∥2).

To bound the last sum in this estimate we choose , in (3.7)-(3.8) and multiply by to get

 cστ∥¯∂tunh∥2H1+τ∥¯∂tθnh∥2+12∥θnh∥2κ≤12∥θn−1h∥2κ.

Summing over and using (2.11) gives

 (3.12) Cn∑j=1τ(∥¯∂tθjh∥2+∥¯∂tujh∥2H1)+cκ2∥θnh∥2H1≤Cκ2∥θ0h∥2H1.

It remains to bound , , and . For this purpose we recall that and use (3.12) for to get

 t1∥¯∂tu1h∥H1 +t1∥¯∂tθ1h∥2+t21∥¯∂tθ1h∥2H1 ≤C(τ(∥¯∂tu1h∥2H1+∥¯∂tθ1h∥2)+∥θ1h∥2H1+∥θ0h∥2H1)≤C∥θ0∥2H1.

Finally, we have that

 tn∥¯∂tunh∥2H1+tn∥¯∂tθnh∥2≤C∥θ0∥2H1,t2n∥¯∂tθnh∥2H1≤C∥θ0∥2H1,

and thus (3.6) follows. ∎

### 3.2. Generalized finite element

In this section we shall derive a generalized finite element method. First we define