Analysis of two thermoelastic problems with the Green–Lindsay model

In this paper, we analyze, from the numerical point of view, two thermo-elastic problems involving the Green–Lindsay theory. The coupling term is different for each case, involving second order or first order spatial derivatives, respectively. The variational formulation leads to a linear coupled system which is written in terms of the velocity and temperature speed. An existence and uniqueness results and the exponential energy decay for the problem with the stronger coupling are recalled. The polynomial energy decay for the weaker coupling is then proved but using the theory of linear semigroups. Then, a fully discrete approximation is introduced using the finite element method and an implicit scheme. A discrete stability property and a main a priori error estimates result are shown, from which we can derive the linear convergence of the approximations. Finally, some numerical simulations are presented to demonstrate the accuracy of the algorithm, the discrete energy decay and the dependence on the relaxation parameter.


Introduction
In this work, we consider two thermoelastic models based on the Green-Lindsay theory (Green and Lindsay 1972). Therefore, it is worth recalling that these authors proposed a theory where the heat conduction was different than the parabolic version. In fact, since that parabolic version allows the propagation of instantaneous thermal waves (and so, it violates the so-called "causality principle"), several authors have proposed different theories trying to overcome this difficulty. Among the pioneering researchers we can cite the above commented proposal by Green and Lindsay (1972) that, even if it is based on Fourier law, we have an energy equation which brings us to a damped hyperbolic equation to describe the behavior of the temperature. There exists a huge number of contributions dealing with this theory [see the numerous references in the papers by Ignaczak (1999, 2000) and the books (Ignaczak and Ostoja-Starzewski 2009;Straughan 2011)]. It is well-known that, when we consider the usual theory of thermoelasticity of Green-Lindsay type, we conclude that the solutions decay through the time and, for the one-dimensional case, this decay is of exponential type. Moreover, if we consider the strain gradient elasticity (or the thermoelastic plate) with the heat conduction of parabolic type (Avalishvili et al. 2018;Aouadi et al. 2019;Aouadi and Moulahi 2015;Fernández-Sare et al. 2010;Liu and Renardy 1995;Liu and Quintanilla 2010;Lindsay and Straughan 1979;Liu et al. 2022;Shivay and Mukhopadhyay 2020), we find that the decay is also exponential in the one-dimensional case or in higher dimensions for the plate. It led recently to the question if it would also happen similarly for the heat conduction of the Green-Lindsay type. When we consider the systems of equations governing both problems we can appreciate a clear similarity, being its main difference the coupling term. In the case of the strain-gradient thermoelasticity, such coupling is weaker, obtaining that, for this case, the energy decay is slow (Quintanilla et al. 2023) and, even more, we prove that, in fact, it is polynomial (see the appendix provided in this work). For the other case, since the coupling is stronger, the energy decay is exponential (even for dimensions greater than one 1 ). In this paper, our aim is to develop a similar study from the numerical point of view, which constitutes one of the novelties of this work. We also think that it is relevant to understand the difference between the behavior of the solutions to these problems with different coupling mechanisms and so, we will study it numerically. Indeed, even if the discrete energy decay is always exponential, in this work we will show that it is much faster when the coupling mechanism is stronger and it leads to an exponential energy decay for the continuous problem.
The plan of this paper is as follows. The model equations for the two problems and the assumptions required on the constitutive parameters are presented in Sect. 2. The existence of a unique solution and the energy decay property, recently proved in Quintanilla et al. (2023), is recalled. Then, by using the classical finite element method and the implicit Euler scheme, a fully discrete approximation is introduced in Sect. 3. A discrete stability property and a priori error estimates are proved, and the linear convergence is derived under some additional regularity conditions on the continuous solution. Finally, some numerical simulations are presented in Sect. 4 to demonstrate the numerical convergence, the behavior of the discrete energy decay and the dependence on the coupling parameter.

The thermoelastic model
In this section, we describe the thermomechanical problems that we will study in the paper. We refer the reader to the recently published work (Quintanilla et al. 2023) for further details. Here, we restrict, for the sake of simplicity in the numerical simulations, to the onedimensional case, although the extension of the numerical analysis, presented in the next section, to the multi-dimensional setting is not difficult. Therefore, we will consider the following linear systems: We note that the unique difference is that, in the coupling term of the first equation of both systems, there is a stronger coupling in the first system. Here, u and θ denote the displacement and the temperature deviation in a rod occupying the interval (0, ), > 0, ρ is the mass density, κ is the thermal conductivity, d is the thermal capacity, μ is the elastic modulus in system (1) but the hyperelastic modulus in system (2) and b is the elastic modulus in system (2). Moreover, α is a relaxation parameter usually used in the Green-Lindsay theory which satisfies the condition αd − h > 0, and we also assume that all parameters ρ, μ, |a|, α, h, d, κ and b are positive. We will study the deformation of the rod over a finite time interval (0, T ), with T > 0 being the final time.
To obtain a well-posed problem, we need to impose boundary and initial conditions to the above systems. Therefore, we assume the following Dirichlet-type boundary conditions for the displacements and Neumann-type boundary conditions for the temperature, that is, and the initial conditions, for a.e. x ∈ (0, ), In the above conditions, functions u 0 , v 0 , θ 0 and ϑ 0 represent the initial data of the problems. Moreover, to guarantee the uniqueness of solutions to the problems we need to assume that functions θ 0 and ϑ 0 have null average.
Hence, we will consider two similar problems: the first one defined by system (1) with boundary conditions (3) and initial conditions (4), and the second one defined by system (2) with boundary conditions (3) and initial conditions (4). Since no doubt the numerical analysis of both problems is rather similar, we will focus on the first problem since it has some additional difficulties due to the stronger coupling between the equations of the system.
The existence and uniqueness, as well as the exponential decay, were recently studied in Quintanilla et al. (2023) and they are summarized as follows.
Theorem 1 Under the condition αd − h > 0, then problem (1), (3) and (4) admits a unique solution which is exponentially stable. In a similar way, we also conclude that problem (2)-(4) admits a unique solution with the same regularity but, in this case, the energy decay is slow (not exponential at least).
The proof of the above theorem is shown in the recent paper (Quintanilla et al. 2023) and it is based on the use of the theory of linear semigroups, adequate energy functionals and some a priori estimates. However, following the work done in Bazarra et al. (2022) for another MGT problem, in fact we could prove that this energy decay is polynomial (see the Appendix A for details).

Fully discrete approximations and an a priori error analysis
In this section, we present an a priori error analysis of a fully discrete scheme for solving problem (1), (3) and (4). As we commented above, the second problem described in the previous section can be analyzed analogously.
First, we need to provide a variational formulation of this problem. Therefore, let us denote Y = L 2 (0, ) and denote by (·, ·) and · the inner product and the norm in this space, respectively. Moreover, let E = H 1 (0, ) and V = H 2 0 (0, ). Integrating by parts system (1) and using boundary conditions (3) we obtain the following weak form of the thermomechanical problem.
Find the velocity v : [0, T ] → V and the temperature speed ϑ : where the displacements u(t) and the temperature θ(t) are recovered from the relations: Remark 1 We note that the variational formulation of the problem defined by system (2) is rather similar. The unique difference is that we must replace variational equations (5) and (6) by the following ones: Now, we will obtain the approximation of problem (5)-(7). We will proceed in two steps. First, we approximate the problem in space. Thus, let us assume that the interval [0, ] is divided into M subintervals of length h = /M, with nodes a 0 = 0 < · · · < a M = , and construct the finite element spaces: where P r ([a i , a i+1 ]) represents the space of polynomials of degree r in the subinterval The discrete initial conditions u 0h , v 0h , θ 0h and ϑ 0h are approximations of the respective initial conditions u 0 , v 0 , θ 0 and ϑ 0 defined as Here, we denote by P h 1 and P h 2 the interpolation operators over the finite element spaces V h and E h , respectively (see Ciarlet 1991 for details).
Second, to provide the time discretization, we consider a uniform partition of the time Therefore, using the well-known implicit Euler scheme we can introduce the following fully discrete problem.
Find the discrete velocity v hk = {v hk n } N n=0 ⊂ V h and the discrete temperature speed , and for all n = 1, . . . , N and where the discrete displacements u hk n and the discrete temperature θ hk n are updated from the relations: Using the assumptions on the constitutive coefficients and applying Lax-Milgram lemma, it is easy to prove that the above discrete problem has a unique solution.
Remark 2 Again, the numerical approximation of the problem obtained by using system (2) is similar to the previous one. The unique difference is that we must replace the discrete variational equations (12) and (13) by the following ones (as we also did in (8), (9)): Now, we will provide a discrete stability property that we summarize in the following.
Lemma 1 Under the assumptions of Theorem 1, the sequences {u hk , v hk , θ hk , ϑ hk }, generated by discrete problem (12)-(14), satisfy the stability estimate: where C is a positive constant assumed to be independent of the discretization parameters h and k.
Proof For the sake of clarity in the calculations, we assume that α = 1. It is straightforward to extend this result to the general case.
Taking as a test function w h = v hk n in discrete variational equation (12) we obtain Using the estimates Now, taking the discrete variational equation (13) for a test function ξ h = ϑ hk n we have (hδϑ hk n + dϑ hk n , ϑ hk Keeping in mind that Combining estimates (17) and (18) Multiplying the above estimates by k and summing up to n it leads v hk Finally, keeping in mind that applying several times Cauchy-Schwarz inequality and Cauchy's inequality ab ≤ a 2 + 1 4 b 2 , a, b, ∈ R with > 0, and using a discrete version of Gronwall's inequality (see Campo et al. (2006)), we obtain the discrete stability.
In the rest of the section, we will obtain some a priori error estimates on the numerical errors v n − v hk n and ϑ n − ϑ hk n , which we state in the following result.
Theorem 2 Let the assumptions of Theorem 1 hold. If we denote by (u, v, θ, ϑ) the solution to problem (5)-(7)s and by (u hk , v hk , θ hk , ϑ hk ) the solution to problem (12)- (14), then we have the following a priori error estimates, for all where C is a positive constant which does not depend on parameters h and k, I j is the integration error given by and, for a Hilbert space X , let · X denote the norm in X .
Proof In this proof, to simplify the calculations, we will consider again that α = 1. We note that we can modify the arguments used below to the general case with some minor changes. First, we obtain the error estimates on the velocity v n − v hk n . Thus, we subtract variational equation (5) for a test function w = w h ∈ V h ⊂ V , at time t = t n , and discrete variational equation (12) where we have used the notations δv n = (v n − v n−1 )/k and δu n = (u n − u n−1 )/k, applying several times Cauchy-Schwarz inequality and the previously commented Cauchy's inequality, we obtain the following error estimates for the velocity field, for all where, here and in what follows, C will represent a positive constant which depends on the constitutive coefficients, but it does not depend on the discretization parameters h and k, and whose value may change even within the same line. Secondly, we derive the a priori error estimates on the temperature speed ϑ n − ϑ hk n . Subtracting variational equation (6), for a test function ξ = ξ h ∈ E h ⊂ E at time t = t n , and discrete variational equation (13), we have, for all ξ h ∈ E h , and so we obtain, for all ξ h ∈ E h ,

Keeping in mind that
Combining estimates (21) and (22) we find that, for all Multiplying the above estimates by k and summing up to n we have, for all where I j is the integration error given in (20), using again a discrete version of Gronwall's inequality (see Campo et al. (2006)) we conclude the desired a priori error estimates.
We note that we can use the above a priori error estimates to derive the convergence order of the approximations under additional regularity conditions on the continuous solution. Therefore, if we assume that we obtain the following.
Corollary 1 Under the additional regularity conditions (23) and the assumptions of Theorem 2, we find that the approximations obtained by problem (12)-(14) are linearly convergent; that is, there exists a positive constant C, independent of the discretization parameters h and k, such that

Remark 3
We note that we could analyze in a similar way variational problem (8), (6) and (7), approximated by the fully discrete problem (15), (13) and (14). Proceeding as in the proof of Lemma 1, we could derive the same dicrete stability property, and, following Theorem 2, we could prove a priori error estimates (19), which would lead to the linear convergence of the approximations under additional regularity conditions (23). However, we omit the details for the sake of clarity in the presentation of this work.

The numerical scheme
In this subsection, we describe the algorithm for solving problem (1), (3) and (4). Thus, given the solution u hk n−1 , v hk n−1 , θ hk n−1 and ϑ hk n−1 at time t n−1 , the discrete velocity v hk n and the dicrete temperature speed ϑ hk n are the solution to the following linear system: Using the well-known commercial code MATLAB, this algorithm was implemented on a 3.2 Ghz PC, and we note that a typical run (h = k = 0.001) took about 1.92 seconds of CPU time.

Numerical convergence in a simple problem
As a first example, we solve problem (1), (3) and (4) to show the accuracy of the approximations. Then, we have used the following data in the simulations: Defining the following initial conditions, for all x ∈ (0, 1), considering the homogeneous Dirichlet boundary conditions (3), if we add the (artificial) supply terms for each equation of system (1), for all (x, t) ∈ (0, 1) × (0, 1), the exact solution to this slightly modified version of problem (1), (3) and (4) can be easily calculated and it has the form, for (x, t) ∈ [0, 1] × [0, 1]: We note that the analysis of this problem is similar to the one shown in the previous section.
Therefore, in Table 1 we show the approximation errors estimated as for several values of the discretization parameters h and k and, in Fig. 1, we plot the evolution of the error depending on the parameter h + k. As can be clearly seen, the convergence of the algorithm is observed, and the linear convergence, stated in Corollary 1, is achieved. If we assume that there are not supply terms, and we use the data: T = 60, ρ = 5, μ = 2, a = 2, α = 0.5, h = 1, d = 10, κ = 1. and the initial conditions, for all x ∈ (0, 1), taking the discretization parameter h = 0.001 and k = 0.001, the evolution in time of the discrete energy (see Quintanilla et al. 2023) defined as: is plotted in Fig. 2 (in both natural and semi-log scales). We note that we have represented their evolution (in natural scale), on the left-hand side, until time t = 10 to see better the shape of the curves. Even if these decays seem to be exponential, on the right-hand side (semi-log scale) we can observe that, after time t = 30, the decay is faster for the secondorder coupling corresponding to Problem 1. We recall that, for this problem, it was proved that the energy decay was exponential meanwhile it was polynomial for Problem 2.

Dependence of the solution on the coupling parameter a
In this example, we will investigate the dependence on parameter a for the solution to the thermomechanical problem defined by system (1) with initial conditions (4) and boundary conditions (3), which we name as Problem 1, and for the solution to the thermomechanical problem defined by system (2) with initial conditions (4) and boundary conditions (3), which we name as Problem 2.
In these simulations, we have used the following data: T = 1, ρ = 10, μ = 2, α = 0.5, h = 1, d = 10, κ = 3, and the initial conditions, for all x ∈ (0, 1), Taking the discretization parameters h = k = 0.001, the solution to Problem 1 is plotted in Figs. 3 and 4 for some values of parameter a. Regarding the displacements and the velocities (see Fig. 3), we can see that both displacements and velocities have a quadratic form, being larger when parameter a increases. If we consider the temperature and the temperature speed (see Fig. 4), we observe that the shape for this largest value of the coupling parameter a is drastically different. By using the same data that in the previous simulations, the solution to Problem 2 is shown in Figs. 5 and 6 for those values of parameter a. We can appreciate that, for values of parameter a less than 1, the displacements and velocities (see Fig. 5) the solutions are very similar. However, for the largest value of this parameter, it seems that there is a delay on the solution. The temperature and the temperature speed are plotted in Fig. 6 and we can see that both are quadratic, being rather similar.