A MGT thermoelastic problem with two relaxation parameters

In this paper, we consider, from both analytical and numerical viewpoints, a thermoelastic problem. The so-called MGT model, with two different relaxation parameters, is used for both the displacements and the thermal displacement, leading to a linear coupled system made by two third-order in time partial differential equations. Then, using the theory of linear semi-groups the existence and uniqueness to this problem is proved. If we restrict ourselves to the one-dimensional case, the exponential decay of the energy is obtained assuming some conditions on the constitutive parameters. Then, using the classical finite element method and the implicit Euler scheme, we introduce a fully discrete approximation of a variational formulation of the thermomechanical problem. A main a priori error estimates result is shown, from which we conclude the linear convergence under suitable additional regularity conditions. Finally, we present some one-dimensional numerical simulations to demonstrate the convergence of the fully discrete approximation, the behavior of the discrete energy decay and the dependence on a coupling parameter.


Introduction
The study of the thermoviscoelastic theories has deserved much interest in the recent years.This is because, in many elastic materials, we can observe viscous mechanical aspects as well as the sensitivity to the thermal effects.We can see that the most general way to propose both aspects can be given by the Kelvin-Voigt dissipation mechanism for the viscosity and the Fourier heat conduction constitutive theory.Unfortunately, both aspects have a relevant drawback.The mechanical (thermal) waves for the theories based on the Kelvin-Voigt (Fourier) theory propagate instantaneously.That is, a mechanical (thermal) deformation imposed in a point in space is felt immediately at any other point in space.This fact violates in a relevant way the so-called causality principle.
This has been the reason why several alternatives theories have been proposed for the heat conduction.However, it is surprising that, although we have a similar effect when we consider the Kelvin-Voigt theory, this fact has not received particular attention.
The most known way to overcome the paradox of the infinite speed of propagation for the thermal waves is the Cattaneo-Maxwell proposition [6], which suggests the introduction of a relaxation parameter.This idea brings to Lord and Shulmann to propose their famous thermoelastic theory [19].After some time, many other theories were considered.In particular, we can recall the works of Green and Naghdi [12,13], who presented three new theories that they called type I, II and III, respectively.The difference among them comes from the choice of the independent constitutive variables.Linear version of type I recovers the classical theory of heat conduction, but types II and III became new theories which have deserved much attention over the last twenty-five years.The most general one is the type III theory, which contains the other two theories as limit cases.Unfortunately, type III theory has the same drawback as the

Basic equations
The aim of this section is to propose the basic equations governing the general theory of MGTthermoviscoelasticity.We are going to consider a bounded region B ⊂ R d , d = 1, 2, 3, with boundary smooth enough to apply the divergence theorem.In this sense, it is worth recalling that the evolution equations are Here, u i is the displacement, ρ is the mass density, t ij is the stress tensor, T 0 is the temperature (assumed to be uniform) in the reference configuration but, from now on, we will assume equal to one, η is the entropy and q i is the heat flux vector.
The constitutive equations are given by , where α is the thermal displacement, θ = α is the temperature, τ 1 and τ 2 are two positive constants, K ij is the thermal conductivity, C * ijrs is the elasticity tensor and C ijkl is the viscosity, c is the thermal capacity and K * ij is a tensor which is typical of the works related with the Green and Naghdi thermoelastic theories [12,13].
We impose the following symmetries on the previous tensors: After substitution of the constitutive equations into the evolution equations, we obtain the following system1 : We note that this coupling is new in the mathematical studies.We will consider this system of equations with homogeneous Dirichlet boundary conditions and imposing the initial conditions, for a.e.x ∈ B, (2.3) Sometimes, it is useful working with the homogeneous one-dimensional case.In this situation, our system of equations becomes: In order to make the calculations easier, in this case we assume that the boundary conditions are: where now the one-dimensional domain is taken as B = (0, π).We note that the last equation in the system of equations (2.1) can be written as (2.5) Therefore, the one-dimensional case becomes: We note that we can write the following equality: where In view of equality (2.6), it will be natural to assume that (i) ρ(x) ≥ ρ 0 > 0 and c(x) ≥ c 0 > 0. (ii) C * ijrs and K * ij are positive definite tensors.That is, there exist two positive constants C and K such that for every tensor ξ ij and vector η i .(iii) C ijrs and K ij are positive definite tensors.That is, there exist two positive constants C 1 and K 1 such that for every tensor ξ ij and vector η i .
We also assume that τ 1 ≥ τ 2 .We note that, for the one-dimensional homogeneous case, we can write the above assumptions as follows: The meaning of condition (i) is obvious.Condition (ii) can be interpreted in the context of the mathematical theory of thermoelastic stability.They are usually imposed.Condition (iii) guarantees that we will have mechanical and thermal dissipation.

Existence of solutions
In this section, we propose an existence and uniqueness result for the solutions to the problem defined by system (2.1) with the modified equation (2.5), the boundary conditions (2.2) and the initial conditions (2.3).
We will work on the Hilbert space: Here, W 1,2 0 (B) and L 2 (B) are the usual Sobolev spaces, and In this space H, we will consider the inner product defined as ) and the bar over an element of the Hilbert space represents its complex conjugated.
In this situation, we can write our problem defined by system (2.1) with the modified equation (2.5), the boundary conditions (2.2) and the initial conditions (2.3) as the following Cauchy problem: where the matrix operator A is defined as In this matrix operator, the elements are given by The domain of the operator A is made by the elements U ∈ H satisfying AU ∈ H.In fact, we note that it is subspace of elements U ∈ H such that It is clear that it is a dense subspace of the Hilbert space H.At the same time, we can obtain the following properties.

Lemma 3.1. There exists a positive constant L such that
for every U ∈ Dom(A).
Proof.If we recall the energy equality (2.6), it is straightforward to see that the relation (3.2) holds.

Lemma 3.2. Zero belongs to the resolvent of operator A.
Proof.
) ∈ H, we need to prove that the equation has a solution.We obtain the system: from which, after easy algebraic manipulations, we conclude that we must solve the simpler system: It is straightforward to see that the right-hand side of the previous system is in In view of the assumptions on C * ijkl and K * ij , we obtain that this system admits a solution in the space From Lemmas 3.1 and 3.2, we can prove that A generates a quasi-contractive semigroup by using the Lumer-Phillips corollary applied to Hille-Yosida theorem.Therefore, problem (3.1) has a unique solution.This is obtained in the following result.
Theorem 3.3.The operator A generates a C 0 -semigroup of contractions in the space H.Moreover, for any initial data U (0) in the domain of the operator A, we conclude that there exists at least one solution to Cauchy problem (3.1) with the regularity: Remark 3.4.It is possible to obtain an existence result under more general assumptions.Even if we could extend this comment to the non-homogeneous case, we now assume that the material is homogeneous.Indeed, we also have the equality: and In view of this equality, under the assumption proposed previously and that l ij is a positive definite tensor, we can define the inner product in H: dv.
In this case, we can see that there exists a positive constant but we do not need the assumptions required above.Again, the domain is dense and we can also prove that zero belongs to the resolvent of this operator A. Therefore, we can obtain the existence of solutions under the assumptions (i)-(iii), but changing the condition on K ij by l ij , which is a weaker condition.

One-dimensional case
In this section, we restrict our analysis to the one-dimensional homogeneous case.
If we define the inner product Therefore, if we assume that the bilinear form given by the matrix is positive definite we conclude that Re AU, U ≤ 0 for every U ∈ Dom(A).
Remark 4.1.Matrix (4.1) is positive definite if and only if the following conditions hold: We see that whenever τ 1 − τ 2 is small enough, the previous conditions hold.
It is worth noting that, under the assumptions proposed in this section, we can prove that zero belongs to the resolvent of the operator.
Therefore, we have shown the following.
Theorem 4.2.Under the previous assumptions (i * ) and that the matrix (4.1) is positive definite, the solutions decay in an exponential way.
Proof.Since the semigroup is contractive, we can use the semigroup theory of linear operators as well as the characterization of the exponentially stable semigroups obtained by Pruss (and other authors), which can be recalled in the book of Liu and Zheng [18].Therefore, in order to prove the theorem it is sufficient to show that the imaginary axis is contained in the resolvent of the operator A and that the asymptotic condition holds.
We will prove the first condition.In the case that this is not fulfilled, there exist a sequence of real numbers λ n → λ = 0 and a sequence of unit norm vectors In view of the assumption on the dissipation D(t), we have that and therefore, If we multiply the above third convergence by v n , we also obtain that a n → 0 in L 2 (B).This is a contradiction because we assumed that U n were unit norm vectors.
To prove the asymptotic condition, we can use a similar argument.

A fully discrete scheme: a priori error estimates
In this section, we will introduce a fully discrete approximation of the problem defined by system (2.1), boundary conditions (2.2) and initial conditions (2.3) over a finite time interval [0, T ], T > 0. First, we need to write this problem in its variational form.Therefore, let us denote Moreover, for a Hilbert space X, let (•, •) X and • X be the inner product and norm in X, respectively.
Therefore, multiplying the equations of system (2.1) by adequate test functions belonging to the spaces V and E, respectively, and using boundary conditions (2.2), we obtain the following variational formulation of this thermo-mechanical problem, written in terms of the acceleration a(t) = (a i (t)), the thermal acceleration φ(t), the velocity v(t) = (v i (t)) and the temperature θ(t).
Find the acceleration a : [0, T ] → V and the thermal acceleration φ : [0, T ] → E such that a(0) = a 0 , φ(0) = φ 0 ,and for a.e.t ∈ (0, T ) and where the operators C, C * , K and K * are given by and the velocity v(t), the temperature θ(t), the displacements u(t) and the thermal displacements α(t) are recovered from the relations: Now, we introduce a fully discrete approximation of problem (5.1)-(5.2).We will proceed in two steps.First, we approximate the problem in space.Thus, let us assume that B is a polyhedral domain and construct the finite element spaces: where T h is a regular finite element triangulation of the domain B (in the sense of [7]), and P 1 (T r) represents the space of affine functions in T r.Moreover, as usual, parameter h > 0 denotes the mesh size.The discrete initial conditions u 0h , v 0h , a 0h , α 0h , θ 0h and φ 0h are approximations of the respective initial conditions u 0 , v 0 , a 0 , α 0 , θ 0 and φ 0 defined as (5.4) ZAMP 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, again, [7]).
Secondly, in order to provide the time discretization, we consider a uniform partition of the time interval [0, T ], denoted by 0 = t 0 < t 1 < • • • < t N = T , where k = T /N is the time step size.If f is a continuous function, we denote f n = f (t n ) and, for the sequence {z n } N n=0 , let δz n = (z n − z n−1 )/k be its divided differences.
Therefore, using the well-known implicit Euler scheme we can introduce the following fully discrete problem.
Find the discrete acceleration a hk = {a hk n } N n=0 ⊂ V h and the discrete thermal acceleration φ hk = {φ hk n } N n=0 ⊂ E h such that a hk 0 = a 0h , φ hk 0 = φ 0h , and for all n = 1, . . ., N and where the discrete velocity v hk n , the discrete temperature θ hk n , the discrete displacements u hk n and the discrete thermal displacements α hk n are updated from the relations: (5.6) Using assumptions (i)-(iii) and applying Lax-Milgram lemma, it is easy to prove that the above discrete problem has a unique solution.
In the rest of the section, we will obtain some a priori error estimates on the numerical errors a n − a hk n and φ n − φ hk n , which we state in the following result.
where C is a positive constant which does not depend on parameters h and k, and the integration errors I 1 j and I 2 j are defined as (5.7) Proof.In this proof, in order to simplify the calculations, we will consider that τ 1 = τ 2 = 1.We note that we can modify the arguments used below to the general case with some minor changes.First, we will obtain the error estimates on the acceleration term a n − a hk n .Thus, we subtract variational equation (5.1) 1 for a test function w = w h ∈ V h ⊂ V , at time t = t n , and discrete variational equation (5.5) 1 to find that Therefore, it follows that, for all ) Y , where we have used assumptions (i)-(iii), applying several times Cauchy-Schwarz inequality and Cauchy's inequality ab ≤ a 2 + 1 4 b 2 , a, b, ∈ R with > 0, we obtain the following error estimates for the acceleration terms: where, here and in what follows, C will represent a positive constant which depends on the constitutive tensors and coefficients, but it does not depend on the discretization parameters h and k, and whose value may change even within the same line.Now, we obtain the error estimates on the thermal acceleration term φ n − φ hk n .Subtracting variational equation (5.1) 2 , for a test function ξ = ξ h ∈ E h ⊂ E at time t = t n , and discrete variational equation (5.5) 2 it follows that Therefore, we obtain that, for all Combining estimates (5.8) and (5.9) and taking into account that Multiplying the above estimates by k and summing up to n, we have, for all We note that it is easy to prove that where I 1 j and I 2 j are the integration errors defined in (5.7).Thus, using a discrete version of Gronwall's inequality (see [5]) 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 (5.10) we obtain the following.
Corollary 5.2.Under the additional regularity conditions (5.10) and the assumptions of Theorem 5.1, we find that the approximations obtained by problem (5.5)-(5.6)are linearly convergent; that is, there exists a positive constant C, independent of the discretization parameters h and k, such that

Numerical results
In this final section, we present the numerical scheme implemented in MATLAB for solving problem (5.5)-(5.6),and we show some numerical examples to demonstrate the accuracy of the approximations, the behavior of the discrete energy decay and the dependence on the coupling coefficient β.

Numerical scheme for the one-dimensional problem
As a first step, given the solution This numerical scheme was implemented on a 3.2 GHz PC using MATLAB, and a typical run (using parameters h = k = 0.001) took about 0.17 s of CPU time

First example: numerical convergence
As a first simpler example, in order to show the accuracy of the approximations the following problem is considered over the domain B = (0, 1).
with the following data: By using the following initial conditions, for all x ∈ B = (0, 1), considering homogeneous Dirichlet boundary conditions and the (artificial) supply terms, for all (x, t) ∈ (0, 1) × (0, 1), the exact solution to the above one-dimensional problem can be easily calculated and it has the form, for (x, t) ∈ [0, 1] × [0, 1]: u(x, t) = α(x, t) = e t x(x − 1).are presented in Table 1 for several values of the discretization parameters h and k.Moreover, the evolution of the error depending on the parameter h + k is plotted in Fig. 1.We notice that the convergence of the algorithm is clearly observed, and the linear convergence, stated in Corollary 5.2, is achieved.If we assume that there are not supply terms, and we use the final time T = 500, the data taking the discretization parameters h = k = 0.001, the evolution in time of the discrete energy given by is plotted in Fig. 2 (in both natural and semi-log scales).As it is demonstrated in this figure, it converges to zero and an exponential decay seems to be achieved.
Then, taking the discretization parameters h = 0.001 and k = 0.001, we plot the solution to problem (5.5)-(5.6) in Figs. 3 and 4 for some values of parameter β.As can be seen in Fig. 3, the displacements and velocities are rather similar for all the values but, if we focus on the accelerations, some oscillations appear for the largest value of the parameter.In Fig. 4, as expected we can see a big dependence on the value of the parameter.The reason is that these thermal displacements are produced by the deformation and so, they increase when the value of the coupling coefficient is greater.

Author contributions
The three authors have contributed equally to this manuscript.

Fig. 4 .
Fig. 4. Example 2: Thermal displacement, temperature and thermal acceleration for different values of β 1 and φ hk n−1 at time t n−1 , variables a hk n and φ hk n are obtained by solving the discrete linear system, for all w

Table 1 .
Example 1: Numerical errors for some values of Example 1: Asymptotic constant error Fig. 2. Example 1: Evolution in time of the discrete energy (natural and semi-log scales)