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.
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 ellipticparabolic 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 ellipticelliptic 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 twodimensional 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
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
Under the assumption that the displacement gradients are small, the (linearized) strain tensor is given by
Assuming further that the material is isotropic, Hooke’s law gives the (total) stress tensor, see e.g. [21] and the references therein,
where is the dimensional identity matrix, is the thermal expansion coefficient, and and are the so called Lamé coefficients given by
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.2)  
(2.3)  
(2.4)  
(2.5)  
(2.6)  
(2.7) 
where is the heat conductivity parameter, which is assumed to be rapidly varying in space.
Remark 2.1.
Assumptions.
We make the following assumptions on the data

, symmetric,

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

, , , 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)  
(2.9) 
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) 
where (resp. ) depends on (resp. and ). Similarly, there are constants (resp. ) depending on the bound (resp. ) such that
(2.11) 
Furthermore, we use the following notation for the energy norms induced by the bilinear forms
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.
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
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 (P1P1 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 timedependent problem locking can occur if is neglected in (2.2) and P1P1 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 P1P1 discretization in this paper.
The classical finite element method with a backward Euler scheme in time reads; for find and , such that
(3.1)  
(3.2) 
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
and we thus define to be the solution to
(3.3) 
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 (P2P1 elements). It is, however, noted that the problem is still stable using P1P1 elements. In this paper we use P1P1 elements and derive error bounds in the norm, of optimal order, for both the displacement and the temperature.
Theorem 3.2.
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.
Proof.
From (3.1)(3.2) and the initial data (3.3) we deduce that the following relation must hold for
(3.7)  
(3.8) 
By choosing and and adding the resulting equations we have
(3.9) 
Note that the coupling terms cancel. By using CauchySchwarz and Young’s inequality we can bound
Multiplying (3.9) by , summing over , and using (2.10) gives
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)  
(3.11) 
Now choose and and add the resulting equations to get
Multiplying by and using CauchySchwarz and Young’s inequality gives
Summing over and using (2.10) now gives
Here we use summation by parts to get
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 ,
The observation that completes the bound (3.5).
Now assume and and note that the following holds for ,
Choosing , and adding the resulting equations gives
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 CauchySchwarz and Young’s inequality
Note that , where we use that if . Summing over now gives
To bound the last sum we choose , in (3.10)(3.11), now with and . Adding the resulting equations gives
Multiplying by and gives after using CauchySchwarz inequality
Summing over and using (2.10) thus gives
To bound the last sum in this estimate we choose , in (3.7)(3.8) and multiply by to get
Summing over and using (2.11) gives
(3.12) 
It remains to bound , , and . For this purpose we recall that and use (3.12) for to get
Finally, we have that
and thus (3.6) follows. ∎
3.2. Generalized finite element
In this section we shall derive a generalized finite element method. First we define