Numerical analysis of a problem of elasticity with several dissipation mechanisms

In this work, we numerically study a problem including several dissipative mechanisms. A particular case involving the symmetry of the coupling matrix and three mechanisms is considered, leading to the exponential decay of the corresponding solutions. Then, a fully discrete approximation of the general case in two dimensions is introduced by using the finite element method and the implicit Euler scheme. A priori error estimates are obtained and the linear convergence is derived under some appropriate regularity conditions on the continuous solution. Finally, some numerical simulations are performed to illustrate the numerical convergence and the behavior of the discrete energy depending on the number of dissipative mechanisms.


Introduction
Time decay of solutions in thermomechanics has deserved much attention. One aspect corresponds to the thermoelasticity. Dafermos [4] started this study showing the asymptotic stability for the one-dimensional case and the existence of undamped isothermal solutions for dimension greater than one. Several authors proved the exponential decay later [16,22] as well as the polynomial decay (but not exponential) in higher dimensions [12,14,15]. Similar results were also obtained for other thermoelastic theories which are not based on the Fourier law [1,[18][19][20][21].
Recently, in [6] Fernández and Quintanilla studied the problem determined by the elasticity with n 2 dissipative mechanisms, where n is the dimension of the domain. They showed that, whenever the matrix of the coupling coefficients is of the rank n 2 , the decay of the solutions is exponential, but different issues need to be considered. One could be if the number n 2 of the dissipative mechanisms is optimal. Another one could be to clarify how the number of dissipative coupling mechanisms affects to The work of R. Quintanilla is supported by the Ministerio de Ciencia, Innovación y Universidades under the research project "Análisis matemático aplicado a la termomecánica" (Ref. PID2019-105118GB-I00). the decay of the solutions. This paper is addressed to deal with these two questions.
The contribution [6] applies to the linearized theory [5,9,10,13], where we note that few symmetries are satisfied by the constitutive coefficients. More symmetries are hold in the linear case. We will prove that, whenever the matrix of coupling coefficients satisfies the symmetries A l ij = A l ji (see system (4)), we can reduce the necessary number of coupling dissipation mechanisms to obtain the exponential decay. We point out that the linear theory satisfies these conditions. To be precise, we only need three dissipative mechanisms in dimension two and six in dimension three.
The mathematical aim of the contribution [6] was analytical. So, it will be suitable to provide a numerical study to the problem of the elasticity with several dissipative mechanisms. In this sense, we will introduce some fully discrete approximations of the general case involving four dissipative mechanisms (the two temperatures and the two mass diffusions), providing a priori error estimates which lead to the linear convergence of the algorithm. Later, we will carry out some simulations (always in dimension two) referring to the different problems when the coupling matrix has rank one, two, three and four, respectively. Since we work with simulations, we will be able to see the (numerical) exponential decay in all cases. This result should be interpreted as an approximation since the simulations correspond to numerical approximations, but we note that the introduction of an increasing number of dissipative mechanisms provides that the rate of the exponential decay (obtained from the numerical approach) also increases. This fact is not surprising, but this is the first mathematical contribution in this aspect. We show a numerical proof that the decay of the solution increases when we increase the number of coupling dissipative mechanisms.
The plan of this article is as follows: in the next section we recall the basic equations of the problem as well as the results obtained in [6] in this sense. Then, we focus on the case that we consider three dissipation mechanisms and the coupling coefficients satisfy a certain symmetry and we verify that, in this case, there is also an exponential decay. The numerical study of the problem is carried out in the fourth section, obtaining some a priori error estimates and deriving the linear convergence under adequate regularity conditions. In the fifth section, we carry out some simulations when we incorporate different dissipation mechanisms and we check how the decay rate increases as we introduce new mechanisms.

Basic equations and assumptions
In this section, we describe the model and the required assumptions on the data, and we recall some results regarding the existence and uniqueness or the energy decay (see [6] for further details).
Let us consider a two-dimensional bounded domain B with a sufficiently smooth boundary which is made of an elastic material coupled with (in general) four dissipative mechanisms. From a physical point of view, we can think that they are determined by two temperatures and two mass diffusion processes. For the sake of simplicity in the modeling, we assume that the material is homogeneous and isotropic in the mechanical, thermal and mass diffusion parts, but we do not impose this requirement on the coupling coefficients. Therefore, the general system is written as follows (see [6]): where i, j, k = 1, 2 and l, q, p = 1, … , 4 , u i denotes the displacement, 1 , 2 are the two temperatures, 3 , 4 are the other two dissipative variables which can be seen as two mass diffusions (see [13]), is the mass density, , are the Lamé constants, ij are related with the thermal and diffusive conductivity, m ij are related with the thermal and diffusive capacity, A l ij are the coupling terms and l kj corresponds to the coefficients associated to the relative temperatures (or concentrations).
We will impose the boundary conditions, for a.e.
where n j is the outward normal vector to the boundary B , for i = 1, 2 and k = 1, … , 4 , and some initial conditions, for a.e. x ∈ B, for i = 1, 2 and k = 1, … , 4 , where functions u 0 1 , u 0 2 , v 0 1 , v 0 2 , 0 1 , 0 2 , 0 3 and 0 4 represent the initial data of the problem. We note that conditions (2) imply that there is not displacement on the boundary and that the dissipative mechanisms are isolated with respect to the exterior.
We will impose the following conditions on the data: inite, that is, there exist three positive constants C 1 , C 2 and C 3 such that Let W 1,2 0 (B) and L 2 (B) be the well-known Sobolev spaces and Then, in [6] we proved the following.
Moreover, if the coupling matrix defined, from the coefficients given in system (1), as has rank 4, then the solutions to problem (1-3) decay in an exponential way.

Case A l ij = A l ji and three dissipative mechanisms
As we pointed out before, four coupling mechanisms are sufficient to guarantee the exponential decay for the elasticity. The natural question is if we can relax the number of couplings. A possible answer is that it will be positive whenever we assume that A l ij = A l ji for every i, j, l and we consider three dissipative mechanisms. The proof of this fact is very easy from the analysis proposed in [6]. For this reason, we are going to sketch the proofs and to follow the similar steps proposed there.

The problem
In this case, we consider the system (see Ieşan [11]): where i, j, k = 1, 2 and l, q, p = 1, … , 3 , which is similar to system (1). It is relevant to say that, in contrast to the reference [6], we assume that l, p, q can be 1, 2 or 3 and 2 12 = − 1 12 = l but the other combinations vanish. From a physical point of view, we can think that the dissipative mechanisms are determined by the two temperatures and one mass dissipation mechanism.
We will assume again conditions (i)-(iii) but now i, j vary between 1 and 3, without the condition on l i since it now leads to a unique positive value l.
We combine system (4) with the initial conditions (3) and the boundary conditions (2), but in the case that k = 1, 2, 3 for a.e. x ∈ B.
We can study the problem (4), (3) and (2) in the Hilbert space: We can define the inner product associated with the norm: where e ij = 1 2 (u i,j + u j,i ). As usual, the bar denotes the conjugated complex.
If U = (u i , v i , d ) we can write our problem in the following abstract from: and n dl m lp = dp , with dp being the Kronecker symbol.
At is was pointed out in [6] the domain of the operator A is made by the elements of the Hilbert space H such that v i ∈ W 1,2 0 (B) , u i + ( + )u j,ji ∈ L 2 (B) and i ∈ W 2,2 (B) ∩ L 2 * (B) . Therefore, the domain of the operator is dense.
At the same time, we have that Following the arguments used in the proof of Lemma 2 in [6] we can develop (word by word) a proof in the sense that zero belongs to the resolvent of the operator A in our case. Therefore, we can obtain the following.

Theorem 2
The operator A generates a C 0 -semigroup of contractions.

Exponential decay
In this subsection, we prove the exponential decay of the solutions to our problem. It is relevant to say that we assume that A l ij = A l ji , for every i, j = 1, 2 and l = 1, 2, 3 , and that the matrix has rank 3.
Then, we get the exponential decay of the solutions that we state as follows.

Theorem 3
The solutions to problem (2-4) decay in an exponential way.
Proof The proof follows the same steps of the arguments used in the proof of Theorem 2 in [6]. That is, we need to prove that the imaginary axis is contained in the resolvent of the operator A and that the condition In order to show the first condition, we use a contradiction argument. Let us assume that this condition does not hold and we will arrive to a contradiction. Assuming that there exists ≠ 0 such that i belongs to the spectrum, there will exist a sequence of real numbers n → and a sequence of elements U n of unit norm, at the domain of the operator A, such that where From the energy inequality, we see that in → 0 in W 1,2 (B) for i = 1, 2, 3 . Then, we obtain where We multiply by A j pq e pqn the jth-convergence and we note that as far as the H 2 −norm of −1 n u in is bounded. It then follows that If we use the Korn inequality we also obtain that we deduce that v pn → 0 in L 2 (B) and we arrive to a contradiction.
We also note that we can use the same argument to prove the asymptotic condition and the theorem is

Remark 1
We have developed the analysis in the two-dimensional case. For the three-dimensional one, we would need that the corresponding coupling matrix has rank six.

Remark 2
Two recent papers [7,8] have extended the decay results obtained in [6] to the type III and the hyperbolic cases. It is worth saying that we can also extend Theorem 3 to these alternative theories. That is, it is possible to obtain the exponential decay in two dimensions for the type III or the hyperbolic cases in the case that A l ij = A l ji and the matrix (5) has rank three.

Numerical analysis of a fully discrete scheme
In the rest of this work, we will simplify the system (1) considerably in order to better demonstrate our results. So, we will take particular values of tensors A l ik , m lp , lp and pq to obtain the following system: In this case, we will impose the following conditions on the constitutive coefficients in order to guarantee assumptions (i)-(iv) 1 : In this section, we present a fully discrete approximation of the simplified problem given by the system (6) with boundary conditions (2) and initial conditions (3), and we perform an a priori error analysis. As usual, we will study the deformation of the body over a finite time interval [0, T], for a given final time First, we provide the variational formulation of problem (6), (2) and (3). So, let us denote the two components of the velocity by v 1 =u 1 and v 2 =u 2 , respectively. Let Y = L 2 (B) and H = [L 2 (B)] 2 and let (⋅, ⋅) Y and ‖ ⋅ ‖ Y (resp. (⋅, ⋅) H and ‖ ⋅ ‖ H ) be the usual scalar product and the norm defined in Y (resp. H). Moreover, let V be the variational space given by V = H 1 0 (B) and denote by ‖ ⋅ ‖ V its usual norm. Multiplying the equations of system (6) by adequate test functions, using the Green's formula and the boundary conditions (2), we obtain the following weak form.
Find the first component of the velocity 1 We note that we assume that l 2 = 0 , but the analysis could be extended without difficulties to the case l 2 ≠ 0 whenever where u(t) = (u 1 (t), u 2 (t)) represents the displacement field and its components are recovered from the relations: Now, we provide the fully discrete approximation of the variational problem (8)(9)(10)(11)(12)(13)(14). This is done in two steps. First, in order to obtain the spatial approximation we assume that the domain B is polyhedral and let us denote by T h a regular finite element triangulation in the sense of [3]. Thus, we construct the finite dimensional space V h ⊂ V given by In the above definition, P 1 (Tr) represents the space of polynomials of degree less or equal to one in the element Tr, i.e. the finite element space V h is made of continuous and piecewise affine functions, and h > 0 is the spatial discretization parameter. Moreover, we assume that the discrete initial conditions, denoted by u 0h i , v 0h i and 0h j (for i = 1, 2 and j = 1, … , 4 ), are given by where P h is the classical finite element interpolation operator over V h (see, e.g., [3]).
In order to obtain the discretization of the time derivatives, we use a uniform partition of the time interval [0, T], denoted by 0 = t 0 < t 1 < ⋯ < t N = T , with step size k = T∕N and nodes t n = n k for n = 0, 1, … , N . For a continuous function z(t), we use the notation z n = z(t n ) and, for the sequence {z n } N n=0 , we denote by z n = (z n − z n−1 )∕k its corresponding divided differences.
where u hk,n = (u hk,n 1 , u hk,n 2 ) represents the discrete displacement field and its components are recovered from the relations: It is worth noting that, under conditions (7), by using Lax-Milgram lemma it is easy to show that the above fully discrete problem has a unique solution.
In the rest of this section, our aim will be to obtain some a priori error estimates for the numerical errors (̇v n where v n 1 = (v n 1 − v n−1 1 )∕k and u n 1 = (u n 1 − u n−1 1 )∕k , after some algebraic manipulations and using several times Cauchy-Schwarz inequality and Cauchy's inequality ab ≤ a 2 + 1 4 b 2 , we find that, for all w h ∈ V h , Proceeding in a similar form, we obtain the error estimates for the second component of the velocity field, for all r h ∈ V h , where v n 2 = (v n 2 − v n−1 2 )∕k and u n 2 = (u n 2 − u n−1 2 )∕k. Combining the previous two estimates and taking into account that it leads, for all w h , r h ∈ V h , Now, we obtain the error estimates for the first temperature. If we subtract variational Eq. (10) at time t = t n for a test function s = s h ∈ V , and discrete variational Eq. (19), we obtain, for all s h ∈ V h , and therefore, we find that, for all s h ∈ V h , Since using again several times Cauchy-Schwarz inequality and the above commented Cauchy's inequality, we obtain, for all s h ∈ V h , where n 1 = ( n 1 − n−1 1 )∕k. In a similar way, we obtain the error estimates for the other temperature 2 and the two mass diffusions 3 and 4 (we skip the details for the sake of clarity): Combining all the estimates of the two components of the velocity, the two temperatures and the two mass diffusions, it follows that, for all Multiplying the above estimates by k and summing up to n, we find that, for all Finally, taking into account that where similar estimates can be found for v 2 , 1 , 2 , 3 and 4 , applying a discrete version of Gronwall's inequality (see [2]) we obtain the following a priori error estimates result.

where C is a positive constant which does not depend on parameters h and k.
It is worth noting that, from the above estimates we can derive the convergence order of the approximations provided by the fully discrete problem (17)(18)(19)(20)(21)(22)(23). For instance, under the additional regularity, for i = 1, 2 and l = 1, 2, 3, 4, the algorithm is linearly convergent, that is, there exists a positive constant C, independent of the discretization operators h and k, such that

Numerical results
In this final section, we present a numerical example to show the accuracy of the approximations and the dependence on the number of dissipation mechanisms. The numerical scheme obtained in Problem VP hk was implemented on a 3.2 GHz PC using MATLAB, and typical run (using parameters h = √ 2∕8 and k = 0.01 ) took about 0.12 s of CPU time.
First, we will demonstrate the numerical convergence of the algorithm. As a simple example, in order to show the accuracy of the approximations, we consider the two-dimensional domain B = (0, 1) × (0, 1) and we will divide it into 2nd 2 triangles.
We have employed the following data in these simulations: By using the following initial conditions, for all (x, y) ∈ B = (0, 1) × (0, 1), considering homogeneous Dirichlet boundary conditions in the boundaries and the (artificial) supply terms, for all (x, y, t) ∈ B × (0, 1), the exact solution to the above two-dimensional problem can be found and it has the form, for (x, y, t) ∈ B × [0, 1]: Thus, the approximation errors estimated by are shown in Table 1 (multiplied by 10) for some values of the discretization parameters h and k. Moreover, by using the diagonal of this table the evolution of the error depending on the parameter h + k is plotted in Fig. 1. As can be seen, the convergence of the algorithm is achieved, although the linear convergence obtained in the previous section under some additional regularity conditions, does not seem to be found. Now, our aim is to compare the effect of the number of dissipation mechanisms in the evolution of the discrete energy. We assume now that there are not supply terms, and we use the final time T = 10 , the following data (common for the four problems described later): and the initial conditions, for all (x, y) ∈ B, We define the following four problems: -Problem 1 the problem with a unique dissipation mechanism. In this case, we only consider the temperature 1 and the data 1 = 1. -Problem 2 the problem with two dissipation mechanisms. Now, we consider the two temperatures 1 and 2 and the data 1 = 2 = l 1 = 1. -Problem 3 the problem with three dissipation mechanisms. The case with the two temperatures 1 and 2 , and one mass diffusion 3 , where we use the data 1 = 2 = l 1 = 3 = 1.
Taking the discretization parameter k = 0.001 and a very fine finite element mesh, the evolution in time of the discrete energy of the above four problems is plotted in Fig. 2 (in both natural and semi-log scales). As can be seen, every energy converges to zero and an exponential decay seems to be achieved in the four cases. Moreover, it is clear that, when the number of dissipation mechanisms increases, the energy decay seems to be faster.

Conclusions
In this paper, we have analyzed a problem of elasticity with several dissipative mechanisms from the analytical and numerical points of view. From the analytical point of view, we have seen that, if we assume that the coupling coefficients satisfy the symmetry A l ij = A l ji , we only need three dissipative mechanisms (in dimension two) to guarantee the exponential energy decay; however, to know if three mechanisms are sufficient when this symmetry does not hold is an open question yet.
From the numerical point of view, we have provided a fully discrete approximation of a weak form of the above thermomechanical problem, and we have proved a main a priori error estimates result. In a numerical example, we have also shown the discrete energy decay of the solutions and we have seen that, when the number of dissipative mechanisms increases, the rate of the decay of the solutions also increases.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Conflicts of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.