On existence and numerical approximation in phase-lag thermoelasticity with two temperatures

UPCommons Portal del coneixement obert de la UPC http://upcommons.upc.edu/e-prints This is a pre-copy-editing, author-produced PDF of an article accepted for publication in Discrete and continuous dynamical systems. Series B , following peer review. The definitive publisher-authenticated version: Campo, M. [et al.]. On existence and numerical approximation in phase-lag thermoelasticity with two temperatures. "Discrete and continuous dynamical systems. Series B", 26 Abril 2021. DOI 10.3934/dcdsb.2021130 is available online at: https://www.aimsciences.org/article/doi/10.3934/dcdsb.2021130 URL d'aquest document a UPCommons E-prints: https://upcommons.upc.edu/handle/2117/345198 Manuscript submitted to doi:10.3934/xx.xxxxxxx AIMS’ Journals Volume X, Number 0X, XX 200X pp. X–XX ON EXISTENCE AND NUMERICAL APPROXIMATION IN1 PHASE–LAG THERMOELASTICITY WITH TWO2 TEMPERATURES3 Marco Campo Departamento de Matem´aticas, ETS de Ingenieros de Caminos, Canales y Puertos Universidade da Coru˜na, Campus de Elvi˜na, 15071 A Coru˜na, Spain Maria I.M. Copetti Laborat´orio de An´alise Num´erica e Astrof´ısica, Departamento de Matem´atica Universidade Federal de Santa Maria, 97105-900, Santa Maria, RS, Brazil Jos´ e R. Fern´ andez∗ Departamento de Matem´atica Aplicada I, Universidade de Vigo Escola de Enxe˜ner´ıa de Telecomunicaci´on, Campus As Lagoas Marcosende s/n 36310 Vigo, Spain Ram´ on Quintanilla Departament de Matem`atiques, Universitat Polit`ecnica de Catalunya C. Colom 11, 08222 Terrassa, Barcelona, Spain (Communicated by the associate editor name) 2020 Mathematics Subject Classification. Primary: 35Q74, 74H10, 80A20, 65M15, 65M60; Secondary: 74F05. Key words and phrases. Phase-lag thermoelasticity, existence and uniqueness, finite elements, a priori estimates, numerical simulations. The work of M. Campo and J.R. Fern´andez has been partially supported by Ministerio de Ciencia, Innovaci´on y Universidades under the research project PGC2018-096696-B-I00 (FEDER, UE). The work of M.I.M. Copetti has been partially supported by the Brazilian institution CNPq (grant 304709/2017-4). The work of R. Quintanilla has been supported by Ministerio de Econom´ıa y Competitividad under the research project “An´alisis Matem´atico de Problemas de la Termomec´anica” (MTM2016-74934-P), (AEI/FEDER, UE), and Ministerio de Ciencia, Innovaci´on y Universidades under the research project “An´alisis matem´atico aplicado a la termomec´anica” (PID2019-105118GB-I00). The authors want to thank to the anonymous referees their useful comments which have allowed us to improve the paper. ∗Corresponding author: Jos´e R. Fern´andez. 1 2 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA Abstract. In this work we study from both variational and numerical points of view a thermoelastic problem which appears in the dual-phase-lag theory with two temperatures. An existence and uniqueness result is proved in the general case of different Taylor approximations for the heat flux and the inductive temperature. Then, in order to provide the numerical analysis, we restrict ourselves to the case of second-order approximations of the heat flux and first-order approximations for the inductive temperature. First, variational formulation of the corresponding problem is derived and an energy decay property is proved. Then, a fully discrete scheme is introduced by using the finite element method for the approximation of the spatial variable and the implicit Euler scheme for the discretization of the time derivatives. A discrete stability property is shown and a priori error estimates are provided, from which the linear convergence of the algorithm is derived under suitable additional regularity on the continuous solution. Finally, some numerical simulations are presented in two-dimensional numerical examples. 1. Introduction. It is widely accepted the Fourier formulation to describe the heat1 conduction. However, when we adjoin this relation with the usual energy equation:2 c˙ θ+ div q= 0, c > 0,(1) we arrive to the instantaneous propagation of heat. This is a drawback of the3 model because this fact is incompatible with basic axioms of the physics. In the4 above equation q= (qi) is the heat flux vector and θis the temperature. In order5 to overcome this difficulty, alternative proposals have been stated. In fact, in the6 late years different authors have developed new theories for the heat conduction.7 The most known is the hyperbolic damped equation [5] studied by Cattaneo and8 Maxwell which eliminates the drawback. Moreover, Green and Nagdhi [14,15]9 proposed three thermoelastic theories where, in each case, the heat conduction is10 described in alternative forms.11 In 1995, Tzou [33] introduced a theory in such a way that the heat flux and the12 gradient of the temperature have a delay in the constitutive equations. When delay13 parameters are taken into account, it is usual to speak about phase-lag theories.14 The constitutive equations proposed by Tzou are given by15 qi(x, t +τ1) = −Kθ,i(x, t +τ2), K > 0,(2) where τ1and τ2are the delay parameters which are assumed to be positive. As16 usual, the notation θ,i means the derivative of θwith respect to the variable xi, and17 repeated subscripts means summation. The derivative with respect to the time is18 denoted using a dot over the function. This equation suggests that the temperature19 gradient established across a material volume, at position xand time t+τ2, results20 in a heat flux to flow at a different time t+τ1. These time delays can be understood21 in terms of the microstructure of the material.22 Some time later, in 2007 Choudhuri [9] extended Tzou’s theory to propose that23 the heat flux is described using the following constitutive equation:24 qi(x, t +τ1) = −K1α,i(x, t +τ3)−K2θ,i(x, t +τ2),(3) where ˙α=θ. The variable αis called the thermal displacement, and the parameter25 τ3is another time delay parameter.26 These two aforementioned theories have several derivations when the heat flux27 and the gradients of the temperature and the thermal displacement are replaced28 by Taylor approximations. In fact, one can think that Choudhuri’s proposal tries29 to recover Green and Naghdi theories when different Taylor approximations are30 PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 3 considered. This new approach gives rise to different equations (depending on the1 selected Taylor approximation) to describe heat conduction that have been analyzed2 by many authors (see, for example, [1,3,16,20,23,27,28,29,30,31,32,35]).3 Unfortunately, the proposals of Tzou and Choudhuri lead to ill-posed problems in4 the sense of Hadamard. In fact, it can be shown that combining equation (2) (or (3))5 with the energy equation (1) leads to the existence of a sequence of elements in the6 point spectrum such that its real part tends to infinity [11]. At the same time, the7 Tzou’s theory is not compatible with the basic axioms of the thermomechanics [13].8 Therefore, we see that this theory cannot be accepted nor from the mathematical9 point of view neither from the thermodynamical point of view.10 In order to obtain a heat conduction theory with delays but without such an11 explosive behavior, Quintanilla [24,25] combined the delay parameters of Tzou12 and Choudhuri with the two-temperatures theory proposed by Chen and Gurtin13 [6,7,8,34]. The basic constitutive equation reads14 qi(x, t +τ1) = −K1β,i(x, t +τ3)−K2T,i(x, t +τ2),(4) where α=β−m∆β,θ=T−m∆Tand mis a positive constant. In this paper,15 we are going to consider the general case when the material is not isotropic, but16 we restrict our attention to the dual-phase-lag theory. We have the constitutive17 equation:18 qi(x, t +τ1) = −Kij T,j(x, t +τ2).(5) In fact, we are going to consider the Taylor approximation to the heat flux vector and the inductive temperature and we assume that q(x, t +τ1)≈a0q(x, t) + a1˙ q(x, t) + ...+anq(n)(x, t), T(x, t +τ2)≈b0T(x, t) + b1˙ T(x, t) + ...+blT(l)(x, t), where l≤n.1 19 This theory has also been extended to the thermoelasticity context [24,25]. To20 do so, one must assume the equation of motion:21 tji,j =ρ¨ui,(6) the energy equation:22 ˙η=−qi,i,(7) and the constitutive equations:23 tji =Cijkluk,l +βijθ, η=−βijui,j +cθ, (8) where tji represents the stress tensor, ηis the entropy, (ui) is the displacement24 vector, Cijkl and βij are constitutive tensors and ρand care the mass density and25 the thermal capacity, respectively.26 It is worth noting that these new thermomechanical theories have attracted much27 attention [2,12,26,18,19,21,22,35].28 Finally, we point out that, in this work, we restrict our attention to the ho-29 mogeneous case to make the calculations easier, but the extension to the case of30 nonhomogeneous materials is direct.31 The paper is outlined as follows. The thermomechanical problem with two tem-32 peratures and the general Taylor developments presented above is described in Sec-33 tion 2, with the assumptions on the constitutive data. Then, in Section 3it is34 1It is worth noting that we can recover the values of aiand bjin terms of τ1and τ2. 4 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA written as a Cauchy problem in a suitable Hilbert space, and an existence and1 uniqueness result is proved in Section 4. Next, a fully discrete approximation is2 introduced in Section 5, based on the finite element method to approximate the3 spatial domain and the backward Euler scheme to discretize the time derivatives.4 Under some conditions on the constitutive parameters, we prove that the discrete5 energy decays. Moreover, a priori error estimates are obtained, from which, under6 suitable additional regularity conditions, the linear convergence of the algorithm is7 deduced. Finally, some numerical simulations are presented in Section 6.8 2. Basic Equations and Assumptions. In this section, we recall the basic equa-9 tions and the assumptions under what we are going to work in this paper. The10 basic system of equations for the phase-lag thermoelasticity with two temperatures11 is given by the system:12 ρ¨ui= (Cijkluk,l +βij θ),j , cd dt(a0θ+a1˙ θ+... +anθ(n)) = (Kij(b0T,j +... +blT(l) ,j ),i) +βij d dt(a0ui,j +...anu(n) i,j ), (9) and the relation13 θ=T−m(KijT,i),j .(10) In this system of equations, Cijkl is the elasticity tensor, βij is the coupling term,14 Kij is the thermal conductivity, ρis the mass density, cis the thermal capacity,15 aiand bjare constants determined by the approximation we consider and mis a16 constant which is typical for the two temperatures theory. As usual, (ui) is the17 displacement vector and θand Tare the thermodynamical temperature and the18 inductive temperature, respectively.19 We are going to consider this system in a multi-dimensional domain Bin Rd 20 (d= 2,3) such that the boundary is smooth enough to apply the divergence the-21 orem. In order to define the problem, we need to assume the initial and bound-22 ary conditions and, to simplify the arguments, we assume homogeneous Dirichlet23 boundary conditions:24 ui(x, t) = T(x, t) = 0 for a.e. x∈∂B, t > 0.(11) We also impose the initial conditions:25 ui(x,0) = u0 i(x),˙ui(x,0) = v0 i(x) for a.e. x∈B, θ(x,0) = θ0(x),...,θ(n)(x,0) = θn(x)for a.e. x∈B. (12) As it is usual, we assume that the elasticity and the thermal conductivity tensors satisfy the symmetries Cijkl =Cklij , Kij =Kji. In this paper, we assume that the constitutive tensors are upper bounded and26 that27 •The mass density and the thermal conductivity are strictly positive; that is, ρ(x)≥ρ0>0, c(x)≥c0>0. •The elasticity tensor is positive definite. That is, there exists a positive constant Csuch that Cijkl ξij ξkl ≥Cξij ξij, for every tensor ξij .28 PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 5 •The thermal conductivity tensor is also positive definite. That is, there exists a positive constant Ksuch that Kijξiξj≥Kξiξi, for every vector ξi.1 •We also assume that parameter anis strictly positive.2 We now introduce the notation ˆ f=a0f+a1˙ f+...anf(n). System (9) can be written as3 ρ¨ ˆui= (Cijkl ˆuk,l +βij (a0θ+a1˙ θ+...+anθ(n))),j , cd dt a0θ+a1˙ θ+... +anθ(n)=Kij(b0T,j +...+blT(l) ,j ),i +βij d dt ˆui,j.(13) To make the notation easier we drop the hat in our system of equations. Therefore,4 we can write5 ρ¨ui= (Cijkluk,l +βij(a0θ+a1˙ θ+...+anθ(n))),j, cd dt a0θ+a1˙ θ+... +anθ(n)=Kij(b0T,j +...+blT(l) ,j ),i +βij ˙ui,j.(14) The aim of this paper is to study this system jointly with initial conditions (12)6 and boundary conditions (11).7 3. The Cauchy Problem. In this section, we will transform our initial-boundaryvalue problem into a Cauchy problem in a suitable Hilbert space. We will consider the space X=W1,2 0(B)×L2(B)×L2(B)×...×L2(B), where Wi,j 0and Liare the usual Sobolev spaces and the boldface means that this8 is the two or three times product.9 We denote by Ω = (ui, vi, θ{0}, θ{1},...,θ{n}) the elements in our Hilbert space and we define the inner product <Ω,Ω∗>=1 2ZB Wdx, where W=ρvivi∗+Cijklui,j u∗ k,l +α0θ{0}θ{0}∗ +...+αn−1θ{n−1}θ{n−1}∗ +c(a0θ{0}+...+anθ{n})(a0θ{0}∗ +...+anθ{n}∗), where αiare positive constants large enough to guarantee that Wdefines a positive definite bilinear form. We note that ||Ω||2 X=1 2ZB ρvivi+Cijklui,j uk,l +α0θ{0}θ{0}+...+αn−1θ{n−1}θ{n−1} +c(a0θ{0}+...+anθ{n})(a0θ{0}+...+anθ{n})dx. This norm is equivalent to the usual one in the Hilbert space.10 It is worth noting that θ(T) = T−m(KijT,j),i defines an isomorphism between L2(B) and W1,2 0(B)∩W2,2(B). We shall denote by Φ(θ) the inverse of this operator. We also note that ||θ||L2=ZB (T2+ 2mKij T,iT,j +m2[(KijT,i),j ]2)dx. Therefore, it is clear that the L2-norm of θis equivalent to the W2,2-norm of T.11 6 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA We shall define the following operators 2: Aiu=ρ−1(Cijkluk,l),j , Bi0θ{0}=ρ−1(βij a0θ{0}),j, Bi1θ{1}=ρ−1(βij a1θ{1}),j,...,Binθ{n}=ρ−1(βij anθ{n}),j, L0θ{0}= (can)−1(Kij b0Φ(θ{0}),j),i, L1θ{1}= (can)−1(Kij b1Φ(θ{1}),j),i −ca0θ{1}, ... Lnθ{n}= (can)−1(KijbnΦ(θ{n}),j),i −can−1θ{n}, Dv= (can)−1βij vi,j, A= (Ai),Bk= (Bik), and the matrix operator1 A=             0I0 0 0 ... 0 A0B0B1B2... Bn 0 0 0 I0... 0 . . . . . ... . . . . . . ... . . . . . . ... . 0 0 0 0 0 ... I 0D L0L1L2... Ln             ,(15) where Iis the identity operator. Our initial-boundary-value problem can be written as dω dt =Aω, ω0= (u0,v0, θ0, θ1, ..., , θn). It is worth noting that the domain of the operator Ais the set of ω∈ X such2 that Aω∈ X is a dense subspace of space X.3 4. The existence theorem. The aim of this section is to prove a theorem of exis-4 tence and uniqueness for the solutions to the Cauchy problem proposed previously.5 We will need a couple of lemmata.6 Lemma 4.1. There exists a positive constant Msuch that hAω, ωi ≤ M||ω||2 for every ωin the domain.7 Proof. If we take into account the evolution equations and the boundary conditions we see that hAω, ωi=1 2ZB D∗dx, where D∗=hCijkluk,l +βijn X k=0 akθ{k}i,j vi+Cijklvk,lui,j + n−1 X k=0 αkθ{k+1}θ{k} +hn−1 X k=0 akθ{k+1}+h n X k=0 KijbkΦ(θ{k}),ji,i − n−1 X k=0 akθ{k+1}+βijvi,j i n X k=0 akθ{k}. 2We recall that eventually vector bkcould be zero in case that l < k ≤n. PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 7 If we apply the divergence theorem we obtain D∗=α0θ{1}θ{0}+α1θ{2}θ{1}+...+αn−1θ{n}θ{n−1} + l X k=0 KijbkΦ(θ{k}),i,j n X s=0 asθ{s}. As l≤nand in view of the equivalence between L2-norm of θand the W2,2-norm of T, we find that there exists a positive constant M1such that ZB D∗dx≤M1ZB ((θ{0})2+...+ (θ{n})2)dx. As the inner product defines a bilinear form which is equivalent to the usual one in1 the Hilbert space, we conclude that the lemma is proved.2 Lemma 4.2. For λlarge enough the range of λI − A is the Hilbert space X.3 Proof. Let (u∗,v∗, θ{0}∗, θ{1}∗,...,θ{n}∗)∈ X, we have to show that for λlarge enough the system λu−v=u∗, λv−Au−B0θ{0}−...−Bnθ{n}=v∗, λθ{0}−θ{1}=θ{0}∗, λθ{1}−θ{2}=θ{1}∗, ... λθ{n−1}−θ{n}=θ{n−1}∗, λθ{n}−Dv−L0θ{0}−L1θ{1}−...−Lnθ{n}=θ{n}∗, has a solution.4 It follows that θ{n}=λθ{n−1}−θ{n−1}∗ =λ2θ{n−2}−λθ{n−2}∗ −θ{n−1}∗ =λnθ{0}−λn−1θ{0}∗−λn−2θ{1}∗−... −λθ{n−2}∗−θ{n−1}∗. In a similar way, we have θ{n−1}=λn−1θ{0}−λn−2θ{0}∗−... −λθ{n−3}∗−θ{n−2}∗, and, in general, we obtain θ{k}=λkθ{0}−λk−1θ{0}∗−... −λθ{k−2}∗−θ{k−1}∗. If we substitute these expressions in our system we see λ2u−Au−(B0+λB1+...+λnBn)θ{0}=F1, λn+1θ{0}−(L0+λL1+... +λnLn)θ{0}−λDu=F2, where F1=λu∗+v∗+P1(λ, θ{0}∗, θ{1}∗,...,θ{n−1}∗) and F2=θ{n}∗+Du∗+P2(λ, θ{0}∗, θ{1}∗,...,θ{n−1}∗). Here, P1and P2are polynomials in λ, but linear in the other components. It is clear that (F1, F2)∈W−1,2×L2. Therefore, to prove the lemma it will be sufficient to show that the bilinear form B[(u, θ{0}),(˜ u,˜ θ{0})] =<(λ2u−Au − n X i=0 λiBiθ{0},−λDu+λn+1θ{0}− n X j=0 λjLiθ{0})(˜ u,˜ θ{0})>L2×L2, 8 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA is a coercive and bounded bilinear form on W1,2 0×L2. In our case, we consider the weights given to multiply the first components by λand the second components by Q(λ) = n X i=0 λiai. We note that Q(λ) is positive for λlarge enough because anis strictly positive. It is clear that our product is bounded. On the other side, we have B[(u, θ{0}),(u, θ{0})] =ZB (Q(λ)λ2uiui+Q(λ)Cijkl ui,juk,l +λ[λn+1(θ{0})2− l X j=1 λjLjθ{0}θ{0}])dx. As Ljare bounded we can take λlarge enough to guarantee that this integral is1 equivalent to the inner product in the corresponding Sobolev space. Therefore, it2 leads to the coerciveness of the bilinear form and the lemma is proved.3 Theorem 4.3. The operator Adefined previously is the generator of a quasi-4 contractive semigroup.5 As a consequence, we have the following result.6 Theorem 4.4. Assume that conditions proposed previously are satisfied and that the7 initial conditions belong to the domain of the operator. Then, there exists a unique8 solution ω(t)which satisfies our system with the aforementioned initial conditions.9 Moreover, we know now that there is continuous dependence of the solutions10 with respect to the initial data.11 Since Ais the generator of a quasi-contractive semigroup, we can also obtain the12 existence and continuous dependence result when supply terms are imposed.13 5. Numerical approximation. In order to simplify the calculations, we assume14 that the material is homogeneous and isotropic, and we take n= 2 and l= 1.15 Hence, system (14) becomes16 ρ¨ui=µui,jj + (λ+µ)uj,ji +β(a0θ+a1˙ θ+a2¨ θ),i in B×(0, Tf), cd dt na0θ+a1˙ θ+a2¨ θo=K(b0∆T+b1∆˙ T) + βdiv vin B×(0, Tf),(16) where div represents the divergence operator and [0, Tf], Tf>0, is the time interval17 of interest.18 We also consider the following boundary and initial conditions:19 ui(x, t) = T(x, t) = 0 for a.e. (x, t)∈∂B ×(0, Tf),(17) ui(x,0) = u0 i(x),˙ui(x,0) = v0 i(x) for a.e. x∈∂B, (18) θ(x,0) = θ0(x),˙ θ(x,0) = θ1(x),¨ θ(x,0) = θ2(x) for a.e. x∈∂B. (19) In this section, we will assume the following conditions on the constitutive coef-20 ficients:21 ρ > 0, µ > 0, λ +µ > 0, a2>0, m > 0, K > 0, c > 0, a1>0, a0>0,b1>0, b2>0.(20) We note that conditions (20) are slightly stronger than those required in Section 222 but they are needed in the proof of the energy decay property for the variational23 solution and the a priori error estimates. Although we could weaken some of the24 conditions, we have imposed all for the sake of simplicity in the analysis.25 In order to simplify the writing and the calculations, in this section we will26 redefine constant mK as m, making an abuse of the notation.27 PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 15 Taking into account that (˙ vn−δvhk n,vn−vhk n)H≥(˙ vn−δvn,vn−vhk n)H +1 2kkvn−vhk nk2 H− kvn−1−vhk n−1k2 H, (div (un−uhk n),div (vn−vhk n))Y≥(div (un−uhk n),div ( ˙ un−δun))Y +1 2kkdiv (un−uhk n)k2 Y− kdiv (un−1−uhk n−1)k2 Y, (∇(un−uhk n),∇(vn−vhk n))Q≥(∇(un−uhk n),∇(˙ un−δun))Q +1 2kk∇(un−uhk n)k2 Q− k∇(un−1−uhk n−1)k2 Q, using Cauchy-Schwarz and Young’s inequalities it follows that, for all wh∈Vh,1 ρ 2kkvn−vhk nk2 H− kvn−1−vhk n−1k2 H+β(Rn−Rhk n,div (vn−vhk n))Y +λ+µ 2kkdiv (un−uhk n)k2 Y− kdiv (un−1−uhk n−1)k2 Y +µ 2kk∇(un−uhk n)k2 Q− k∇(un−1−uhk n−1)k2 Q ≤Ck˙ vn−δvnk2 H+kvn−whk2 V+k∇(un−uhk n)k2 Q+k˙ un−δunk2 V+kRn−Rhk nk2 Y +kdiv (un−uhk n)k2 Y+kvn−vhk nk2 H+ (δvn−δvhk n,vn−wh)H.(37) Now, we subtract variational equation (22) at time t=tnfor a test function η=ηh∈Eh⊂Eand discrete variational equation (32) to obtain, for all ηh∈Eh, c(˙ Rn−δRhk n, ηh)Y−β(div (vn−vhk n), ηh)Y−K(b0∆(Tn−Thk n) + b1∆(φn−φhk n), ηh)Y= 0, and so, we have, for all ηh∈Eh, c(˙ Rn−δRhk n, Rn−Rhk n)Y−K(b0∆(Tn−Thk n) + b1∆(φn−φhk n), Rn−Rhk n)Y −β(div (vn−vhk n), Rn−Rhk n)Y =c(˙ Rn−δRhk n, Rn−ηh)Y−K(b0∆(Tn−Thk n) + b1∆(φn−φhk n), Rn−ηh)Y −β(div (vn−vhk n), Rn−ηh)Y. Keeping in mind that (˙ Rn−δRhk n, Rn−Rhk n)Y≥(˙ Rn−δRn, Rn−Rhk n)Y+1 2knkRn−Rhk nk2 Y− kRn−1−Rhk n−1k2 Yo, (div (vn−vhk n), Rn−ηh)Y=−(vn−vhk n,∇(Rn−ηh))H, it follows that2 1 2knkRn−Rhk nk2 Y− kRn−1−Rhk n−1k2 Yo−K(b0∆(Tn−Thk n) + b1(φn−∆φhk n), Rn−Rhk n)Y −β(Rn−Rhk n,div (vn−vhk n))Y ≤Ck˙ Rn−δRnk2 Y+kRn−ηhk2 Y+k∇(Rn−ηh)k2 H+kRn−Rhk nk2 Y+kvn−vhk nk2 H +k∆(Tn−Thk n)k2 Y+k∆(φn−φhk n)k2 Y+ (δRn−δRhk n, Rn−ηh)Y.(38) 16 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA Combining estimates (37) and (38), multiplying the resulting estimates by kand1 summing up to n, we find that2 kvn−vhk nk2 H+kdiv (un−uhk n)k2 Y+k∇(un−uhk n)k2 Q+kRn−Rhk nk2 Y ≤Ck n X j=1 k˙ vj−δvjk2 H+kvj−wh jk2 V+k∇(uj−uhk j)k2 Q+k˙ uj−δujk2 V +kdiv (uj−uhk j)k2 Y+kvj−vhk jk2 H+ (δvj−δvhk j,vj−wh j)H+k˙ Rj−δRjk2 Y +kRj−ηh jk2 Y+k∇(Rj−ηh j)k2 H+k∆(Tj−Thk j)k2 Y+k∆(φj−φhk j)k2 Y +kRj−Rhk jk2 Y+ (δRj−δRhk j, Rj−ηh j)Y+Ckv0−v0hk2 H+kdiv (u0−u0h)k2 Y +k∇(u0−u0h)k2 Q+kR0−R0hk2 Y.(39) Finally, we subtract variational equation (22) at time t=tnand discrete variational equation (31) for a test function z=zh∈Wh⊂Wto obtain (Rn−Rhk n, zh)Y−(a2(ψn−ψhk n−m∆(ψn−ψhk n)), zh)Y −(a0(Tn−Thk n−m∆(Tn−Thk n)) + a1(φn−φhk n−m∆(φn−φhk n)), zh)Y= 0, and therefore, (Rn−Rhk n,∆(ψn−ψhk n))Y+a1(φn−φhk n−m∆(φn−φhk n)),∆(ψn−ψhk n))Y −(a0(Tn−Thk n−m∆(Tn−Thk n)) −(a2(ψn−ψhk n−m∆(ψn−ψhk n)),∆(ψn−ψhk n))Y = (Rn−Rhk n,∆(ψn−zh))Y+a1(φn−φhk n−m∆(φn−φhk n)),∆(ψn−zh))Y −(a0(Tn−Thk n−m∆(Tn−Thk n)) −(a2(ψn−ψhk n−m∆(ψn−ψhk n)),∆(ψn−zh))Y. Taking into account that −(a0(Tn−Thk n),∆(ψn−ψhk n))Y=a0(∇(Tn−Thk n),∇(ψn−ψhk n))H, −(a1(φn−φhk n),∆(ψn−ψhk n))Y=a1(∇(φn−φhk n),∇(ψn−ψhk n))H ≥a1(∇(φn−φhk n),∇(˙ φn−δφn))H+a1 2knk∇(φn−φhk n)k2 H− k∇(φn−1−φhk n−1)k2 Ho, −(a2(ψn−ψhk n),∆(ψn−ψhk n))Y=a2(∇(ψn−ψhk n),∇(ψn−ψhk n))H≥a2k∇(ψn−ψhk n)k2 H, (∆(φn−φhk n),∆(ψn−ψhk n))Y≥(∆(φn−φhk n),∆( ˙ φn−δφn))Y +1 2knk∆(φn−φhk n)k2 Y− k∆(φn−1−φhk n−1)k2 Yo, a2m(∆(ψn−ψhk n),∆(ψn−ψhk n))Y=a2mk∆(ψn−ψhk n)k2 Y, −(ψn−ψhk n,∆(ψn−zh))Y= (∇(ψn−ψhk n),∇(ψn−zh))H, −(φn−φhk n,∆(ψn−zh))Y= (∇(φn−φhk n),∇(ψn−zh))H, −(Tn−Thk n,∆(ψn−zh))Y= (∇(Tn−Thk n),∇(ψn−zh))H, we find that3 1 2knk∇(φn−φhk n)k2 H− k∇(φn−1−φhk n−1)k2 Ho+k∆(ψn−ψhk n)k2 Y+k∇(ψn−ψhk n)k2 H +1 2knk∆(φn−φhk n)k2 Y− k∆(φn−1−φhk n−1)k2 Yo ≤CkRn−Rhk nk2 Y+k∆(ψn−zh)k2 Y+k∇(˙ φn−δφn)k2 H+k∆( ˙ φn−δφn)k2 Y +k∇(φn−φhk n)k2 H+k∆(φn−φhk n)k2 Y+k∆(Tn−Thk n)k2 Y+k∇(ψn−zh)k2 H +k∇(Tn−Thk n)k2 H, PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 17 and therefore,1 k∇(φn−φhk n)k2 H+k∆(φn−φhk n)k2 Y+k n X j=1 nk∆(ψj−ψhk j)k2 Y+k∇(ψj−ψhk j)k2 Ho ≤Ck n X j=1 kRj−Rhk jk2 Y+k∆(ψj−zh j)k2 Y+k∇(˙ φj−δφj)k2 H+k∆( ˙ φj−δφj)k2 Y +k∇(φj−φhk j)k2 H+k∆(φj−φhk j)k2 Y+k∆(Tj−Thk j)k2 Y+k∇(ψj−zh j)k2 H +k∇(Tj−Thk j)k2 H+CkT1−T1hk2 H2(B).(40) Now, if we combine estimates (39) and (40) it follows that2 kvn−vhk nk2 H+kdiv (un−uhk n)k2 Y+k∇(un−uhk n)k2 Q+kRn−Rhk nk2 Y +k∇(φn−φhk n)k2 H+k∆(φn−φhk n)k2 Y+k n X j=1 nk∆(ψj−ψhk j)k2 Y+k∇(ψj−ψhk j)k2 Ho ≤Ck n X j=1 k˙ vj−δvjk2 H+kvj−wh jk2 V+k∇(uj−uhk j)k2 Q+k˙ uj−δujk2 V +kdiv (uj−uhk j)k2 Y+kRj−Rhk jk2 Y+k∇(Rj−ηh j)k2 H+kvj−vhk jk2 H +(δvj−δvhk j,vj−wh j)H+k˙ Rj−δRjk2 Y+kRj−ηh jk2 Y+k∆(Tj−Thk j)k2 Y +k∆(φj−φhk j)k2 Y+ (δRj−δRhk j, Rj−ηh j)Y+k∆(ψj−zh j)k2 Y+k∇(˙ φj−δφj)k2 H +k∆( ˙ φj−δφj)k2 Y+k∇(Tj−Thk j)k2 H+k∇(φj−φhk j)k2 H+Ckv0−v0hk2 H +kdiv (u0−u0h)k2 Y+k∇(u0−u0h)k2 Q+kR0−R0hk2 Y+kT1−T1hk2 H2(B). Now, we observe that3 k n X j=1 (δvj−δvhk j,vj−wh j)H= n X j=1 (vj−vhk j−(vj−1−vhk j−1),vj−wh j)H = (ρ(vn−vhk n),vn−wh n)H+ (ρ(v0h−v0),v1−wh 1)H + n−1 X j=1 (ρ(vj−vhk j),vj−wh j−(vj+1 −wh j+1))H, k n X j=1 (δRj−δRhk j, Rj−zh j)Y= n X j=1 (Rj−Rhk j−(Rj−1−Rhk j−1), Rj−zh j)Y = (Rn−Rhk n, Rn−zh n)Y+ (R0h−R0, R1−zh 1)Y + n−1 X j=1 (Rj−Rhk j, Rj−zh j−(Rj+1 −zh j+1))Y, k∇(Tn−Thk n)k2 H≤Ck∇(T0−T0h)k2 H+In+k n X j=1 k∇(φn−φhk n)k2 H, k∆(Tn−Thk n)k2 Y≤Ck∆(T0−T0h)k2 Y+Jn+k n X j=1 k∆(φn−φhk n)k2 Y, 18 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA where Inand Jnare the integration errors given by1 In=  Ztn 0 ∇φ(s)ds −k n X j=1 ∇φj   2 H, Jn=  Ztn 0 ∆φ(s)ds −k n X j=1 ∆φj   2 Y.(41) Using Poincar´e inequality for the inductive temperature and a discrete version2 of Gronwall’s inequality ([4]) we have the following.3 Theorem 5.4. Let the assumptions (20) hold. If we denote by (u,v, θ, e, ξ, T, φ, ψ)4 and (uhk,vhk, θhk, ehk, ξhk, T hk, φhk, ψhk)the respective solutions to problems V P5 and V Phk, then we have the following a priori error estimates, for all wh=6 {wh j}N j=0 ⊂Vhand ηh={ηh j}N j=0, zh={zh j}N j=0 ⊂Wh,7 max 0≤n≤Nnkvn−vhk nk2 H+kdiv (un−uhk n)k2 Y+k∇(un−uhk n)k2 Q+kRn−Rhk nk2 Y +k∇(φn−φhk n)k2 H+k∆(φn−φhk n)k2 Y+k∇(Tn−Thk n)k2 H+k∆(Tn−Thk n)k2 Yo +k N X j=1 nk∆(ψj−ψhk j)k2 Y+k∇(ψj−ψhk j)k2 Ho ≤Ck N X j=1 k˙ vj−δvjk2 H+kvj−wh jk2 V+k˙ uj−δujk2 V+k∇(Rj−ηh j)k2 H+k∆( ˙ φj−δφj)k2 Y 8 +k˙ Rj−δRjk2 Y+kRj−ηh jk2 Y+k∆(ψj−zh j)k2 Y+k∇(˙ φj−δφj)k2 H+Ij+Jj +Cmax 0≤n≤Nnkvn−wh nk2 H+kRn−zh nk2 Yo+C k N−1 X j=1 kRj−zh j−(Rj+1 −zh j+1)k2 Y +C k N−1 X j=1 kvj−wh j−(vj+1 −wh j+1)k2 H+Ckv0−v0hk2 H+ku0−u0hk2 V +kR0−R0hk2 Y+kT1−T1hk2 H2(B)+kT0−T0hk2 H2(B), where C > 0is a positive constant assumed to be independent of the discretization parameters hand kbut depending on the continuous solution, and the integration errors Ijand Jjare given by (41). We also recall the notations: R(t) = a0θ(t) + a1e(t) + a2ξ(t), Rhk n=a0θhk n+a1ehk n+a2ξhk n. We note that we can study the convergence order from the previous estimates.9 Therefore, as an example, we have the following result which states the linear con-10 vergence of the approximation under suitable additional regularity conditions (see11 [4] for details regarding the estimates of the non-usual finite element terms).12 Corollary 1. If we assume that the continuous solution to Problem V P has the regularity: u∈H3(0, T ;H)∩C1([0, T ]; [H2(B)]d)∩H2(0, T ;V), θ∈W2,∞(0, T ;H2(B)) ∩H3(0, T ;H1(B)), T ∈W2,∞([0, T ]; H3(B)) ∩H3(0, T ;H2(B)), PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 19 then the approximations provided by Problem V P hk are linearly convergent; i.e., there exists a positive constant C > 0such that max 0≤n≤Nnkvn−vhk nkH+kdiv (un−uhk n)kY+k∇(un−uhk n)kQ+kRn−Rhk nkY +k∇(φn−φhk n)kH+k∆(φn−φhk n)kY+k∇(Tn−Thk n)kH+k∆(Tn−Thk n)kYo≤C(h+k). Remark 2. If we assume the material to be viscoelastic, that is, if we replace equation (16) by ρ¨ui=µui,jj + (λ+µ)uj,ji +µ∗˙ui,jj + (λ∗+µ∗) ˙uj,ji +β(a0θ+a1˙ θ+a2¨ θ),i, where λ∗and µ∗are viscosity parameters, then the above error estimates can be1 improved.2 6. Numerical results. In order to verify the behavior of the numerical method3 described in the previous section, some numerical experiments have been performed4 in two-dimensional problems.5 6.1. Numerical scheme. First, we describe the numerical algorithm used to solve6 Problem V Phk. So, given vhk n−1,ξhk n−1and ψhk n−1, the discrete velocity, the discrete7 thermal acceleration and the discrete inductive thermal acceleration, vhk n,ξhk nand8 ψhk n, respectively, are the solution to the following coupled linear system:9 ρ(vhk n,wh)H+µk2(∇vhk n,∇wh)Q+k2(λ+µ)(div vhk n,div wh)Y +kβ((a0k2+a1k+a2)ξhk n,div wh)Y =ρ(vhk n−1,wh)H−µk(∇uhk n−1,∇wh)Q−k(λ+µ)(div uhk n−1,div wh)Y −βk(a0θhk n−1+ (a0k+a1)ehk n−1,div wh)Y, ((a0k2+a1k+a2)ξhk n−(a0k2+a1k+a2)ψhk n+mK(a0k2+a1k+a2)∆ψhk n, zh)Y = (−a0kθhk n−1−(a0k+a1)ehk n−1+a0kT hk n−1+ (a0k+a1)φhk n−1 −mK(a0k∆Thk n−1+a1∆φhk n−1), zh)Y, c((a0k2+a1k+a2)ξhk n, ηh)Y−K((b0k3+b1k2)∆ψhk n, ηh)Y−β(div vhk n, ηh)Y =c(a2ξhk n−1−a0kehk n−1, ηh)Y+kK(b0∆Thk n−1+ (b0k+b1)∆φhk n−1, ηh)Y, where the discrete displacement uhk n, thermal velocity ehk n, temperature θhk n, inductive thermal velocity φhk nand inductive temperature Thk nare then recovered from the relations: uhk n=kvhk n+uhk n−1, ehk n=kξhk n+ehk n−1, θhk n=kehk n+θhk n−1, φhk n=kψhk n+φhk n−1, T hk n=kφhk n+Thk n−1, We note that this discrete problem consists of three coupled symmetric linear10 equations, and so Cholesky’s method is used for the matrix factorization written in11 terms of a product variable.12 The numerical scheme was implemented using FreeFEM++ (see [17] for details)13 on a Intel Core i5-3337U @ 1.80GHz and a typical run (100 step times and 100014 nodes) took about 200 seconds of CPU time.15 20 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA 6.2. Numerical convergence. We consider the following academic problem:1 Problem Pex.Find the displacements u: [0,1] ×[0,1] ×[0,1] →R2, the temperature θ: [0,1] ×[0,1] ×[0,1] →Rand the inductive temperature T: [0,1] × [0,1] ×[0,1] →Rsuch that ¨ui−10ui,jj −20uj,ji −(θ+˙ θ+1 2¨ θ),i =Hiin (0,1) ×(0,1) ×(0,1), ˙ θ+¨ θ+1 2 ... θ−∆T−∆˙ T= div ˙ u+P, in (0,1) ×(0,1) ×(0,1), ui(x, y, t) = T(x, y, t) = 0 for i= 1,2and (x, y, t)∈∂([0,1] ×[0,1]) ×(0,1), ui(x, y, 0) = x2y2(1 −x)2(1 −y)2for i= 1,2and (x, y)∈[0,1] ×[0,1] ˙ui(x, y, 0) = x2y2(1 −x)2(1 −y)2for i= 1,2and (x, y)∈[0,1] ×[0,1], T(x, y, 0) = x2y2(1 −x)2(1 −y)2for (x, y)∈[0,1] ×[0,1], ˙ T(x, y, 0) = x2y2(1 −x)2(1 −y)2for (x, y)∈[0,1] ×[0,1], ¨ T(x, y, 0) = x2y2(1 −x)2(1 −y)2for (x, y)∈[0,1] ×[0,1], where the (artificial) body forces H= (H1, H2) and the heat supply Pare given2 by3 H1(x, y, t) = etx4y4−2x4y3−119 x4y2+ 120 x4y−20 x4−14 x3y4−292 x3y3+ 850 x3y2 −544 x3y+ 64 x3−341 x2y4+ 1162 x2y3−1397 x2y2+ 576 x2y−56 x2+ 426 x y4 −1012 x y3+ 738 x y2−152 x y + 12 x−96 y4+ 192 y3−96 y2, H2(x, y, t) = etx4y4−14 x4y3−341 x4y2+ 426 x4y−96 x4−2x3y4−292 x3y3+ 1162 x3y2 −1012 x3y+ 192 x3−119 x2y4+ 850 x2y3−1397 x2y2+ 738 x2y−96 x2+ 120 x y4 −544 x y3+ 576 x y2−152 x y −20 y4+ 64 y3−56 y2+ 12 y, P(x, y, t) = et5x2y2(x−1)2(y−1)2/2−9x2y2(y−1)2−9x2(x−1)2(y−1)2 −9y2(x−1)2(y−1)2−x2y2(2x−2)(y−1)2−x2y2(2y−2)(x−1)2−9x2y2(x−1)2 −18xy2(2x−2)(y−1)2−18x2y(2y−2)(x−1)2−2xy2(x−1)2(y−1)2 −2x2y(x−1)2(y−1)2. We note that Problem Pex corresponds to Problem (16)-(19) with the following4 data:5 B= (0,1) ×(0,1), Tf= 1, ρ = 1, λ =µ= 10, a0= 1, a1= 1, a2=1 2, β= 1, m = 1, b0= 1, b1= 1, K = 1, c = 1, u0 i(x, y) = v0 i(x, y) = x2y2(x−1)2(y−1)2for all (x, y)∈[0,1] ×[0,1], θ0(x, y) = θ1(x, y) = θ2(x, y) = x2y2(x−1)2(y−1)2for all (x, y)∈[0,1] ×[0,1]. The exact solution to Problem Pex has the following form:6 T(x, y, t) = etx2y2(x−1)2(y−1)2for (x, y, t)∈[0,1] ×[0,1] ×[0,1], ui(x, y, t) = etx2y2(x−1)2(y−1)2for i=1,2 and (x, y, t)∈[0,1] ×[0,1] ×[0,1]. The numerical errors given by max 0≤n≤Nnkvn−vhk nkH+kdiv (un−uhk n)kY+k∇(un−uhk n)kQ+kRn−Rhk nkY +k∇(φn−φhk n)kH+k∆(φn−φhk n)kY+k∇(Tn−Thk n)kH+k∆(Tn−Thk n)kYo, PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 21 nel ↓k→0.1 0.05 0.02 0.01 0.005 8 0.9251309 0.6235920 0.4613584 0.4168623 0.3989068 16 0.6973870 0.3659938 0.1719483 0.1109517 0.0833035 32 0.6726871 0.3381078 0.1391413 0.0734105 0.0411872 64 0.6725218 0.3349879 0.1340351 0.0677103 0.0347045 128 0.6732862 0.3355613 0.1338216 0.0668053 0.0334754 Table 1. Example 1: Numerical errors (×103) for some discretization parameters. and obtained for different discretization parameters nd and k, are depicted in Table1 1(being nd the number of subdivisions on each outer side of the square). Moreover,2 their evolution depending on the parameter h+kis plotted in Figure 1. We3 observe that the convergence of the numerical scheme is clearly obtained but the4 linear convergence, shown in Corollary 1, is not achieved.5 0 0.05 0.1 0.15 0.2 0.25 h+k 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Numerical errors 10-3 Asymptotic behavior Figure 1. Example 1: Asymptotic behavior of the numerical scheme. If we assume now that there are not volume forces nor heat supply, and we6 use as final time Tf= 2 s, (with the same data and mechanical initial conditions7 than in the previous example), taking the discretization parameters nd = 32 and8 k= 0.01,the evolution in time of the discrete energy Ehk n,defined by (34), is plotted9 in Figure 2in both usual and semilog scales. The energy converges to zero and10 an exponential decay seems to be achieved; however, we note that such behavior11 is not found in the continuous case (see, for instance, [19] in the analysis of the12 dual-phase-lag case).13 6.3. Dependence on the thermal coefficient m.In this second example, we14 study the dependence of the solution with respect to parameter m. In particular,15 we consider the domain B= (0,8) ×(0,1) and the final time Tf= 0.5. In these16 simulations we use the following data:17 ρ= 1, λ =µ= 10, a0= 1, a1= 1, a2=1 2, β = 1, b0= 1, K= 1, c = 1, b1= 1, 22 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 t 0 0.5 1 1.5 2 2.5 En hk 10-3 Discrete energy evolution 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 t -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 log10(En hk) Discrete energy evolution Figure 2. Example 1: Energy evolution in absolute and semilogarithmic scales. and the initial conditions, for all (x, y)∈[0,8] ×[0,1],1 u0 i(x, y) = v0 i(x, y) = 0, T0(x, y) = max{(x−3)(5 −x)y(1 −y),0}, T 1(x, y) = T2(x, y) = 0. We solve Problem V Phk with the time discretization parameter k= 0.001 and a2 fixed spatial finite element mesh. Thus, in Figure 3we plot the evolution in time3 of the temperature and inductive temperature at point x= (4,0.5) for different4 values of parameter m(varying between 0.5 and 0.005). As we can see, the inductive5 temperatures almost coincide and important differences appear for the temperature.6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 (4,0.5,t) m = 0.5 m = 0.2 m = 0.1 m = 0.05 m = 0.02 m = 0.01 m = 0.005 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t 0 0.05 0.1 0.15 0.2 0.25 T(4,0.5,t) m = 0.5 m = 0.2 m = 0.1 m = 0.05 m = 0.02 m = 0.01 m = 0.005 Figure 3. Example 2: Evolution in time of the temperature and inductive temperature at point x= (4,0.5) for different values of parameter m. Now, in Figure 4we plot the evolution in time of the temperature and inductive7 temperature at point x= (1,0.5) for the same values of parameter m. Again, the8 inductive temperatures are rather similar, although important differences are found9 for the temperatures.10 Finally, in Figure 5the evolution in time of the horizontal and vertical dis-11 placements at point x= (1,0.5) for the above values of parameter m. We can also12 observe the differences among the corresponding solutions.13 PHASE–LAG THERMOELASTICITY WITH TWO TEMPERATURES 23 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t -6 -4 -2 0 2 4 6 (1,0.5,t) 10-4 m = 0.5 m = 0.2 m = 0.1 m = 0.05 m = 0.02 m = 0.01 m = 0.005 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t -1 0 1 2 3 4 5 6 7 8 T(1,0.5,t) 10-5 m = 0.5 m = 0.2 m = 0.1 m = 0.05 m = 0.02 m = 0.01 m = 0.005 Figure 4. Example 2: Evolution in time of the temperature and inductive temperature at point x= (1,0.5) for different values of parameter m. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t -5 0 5 10 15 20 u1(1,0.5,t) 10-5 m = 0.5 m = 0.2 m = 0.1 m = 0.05 m = 0.02 m = 0.01 m = 0.005 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t -1 -0.5 0 0.5 1 1.5 u2(1,0.5,t) 10-5 m = 0.5 m = 0.2 m = 0.1 m = 0.05 m = 0.02 m = 0.01 m = 0.005 Figure 5. Example 2: Evolution in time of the horizontal and vertical displacements at point x= (1,0.5) for different values of parameter m. REFERENCES1 [1] I. A. Abdallah, Dual phase lag heat conduction and thermoelastic properties of a semi-infinite2 medium induced by ultrashort pulsed laser, Prog. Phys. 3 (2009) 60-63.3 [2] S. Banik, M. Kanoria, Effects of three-phase-lag on two temperatures generalized thermoelas-4 ticity for an infinite medium with a spherical cavity, Appl. Math. Mech. 33 (2012) 483-498.5 [3] K. Borgmeyer, R. Quintanilla, R. Racke, Phase-lag heat conduction: decay rates for limit6 problems and well-posedness, J. Evol. Equ. 14 (2014) 863-884.7 [4] M. Campo, J. R. Fern´andez, K. L. Kuttler, M. Shillor, J. M. Via˜no, Numerical analysis and8 simulations of a dynamic frictionless contact problem with damage, Comput. Methods Appl.9 Mech. Engrg. 196 (2006) 476–488.10 [5] C. Cattaneo, On a form of heat equation which eliminates the paradox of instantaneous11 propagation, C. R. Acad. Sci. Paris 247 (1958) 431–433.12 [6] P. J. Chen, M. E. Gurtin, On a theory of heat involving two temperatures, ZAMP-Z. Angew.13 Math. Phys. 19 (1968) 614–627.14 [7] P. J. Chen, M. E. Gurtin, W. O. Williams, A note on non-simple heat conduction, ZAMP-Z.15 Angew. Math. Phys. 19 (1968) 969–970.16 [8] P. J. Chen, M. E. Gurtin, W. O. Williams, On the thermodynamics of non-simple materials17 with two temperatures, ZAMP-Z. Angew. Math. Phys. 20 (1969) 107–112.18 24 CAMPO, COPETTI, FERN´ ANDEZ AND QUINTANILLA [9] S. K. R. Choudhuri, On a thermoelastic three-phase-lag model, J. Thermal Stresses 30 (2007)1 231-238.2 [10] P. G. Ciarlet, Basic error estimates for elliptic problems. In: Handbook of Numerical Analysis,3 P.G. Ciarlet and J.L. Lions eds., vol II (1993), 17-351.4 [11] M. Dreher, R. Quintanilla, R. Racke, Ill-posed problems in thermomechanics, Appl. Math.5 Lett. 22 (2009) 1374–1379.6 [12] M. A. Ezzat, A. S. El-Karamany, S. M. Ezzat, Two-temperature theory in magneto-7 thermoelasticity with fractional order dual-phase-lag heat transfer, Nuc. Eng. Des. 252 (2012)8 267–277.9 [13] M. Fabrizio, F. Franchi, Delayed thermal models: stability and thermodynamics, J. Thermal10 Stresses 37 (2014) 160–173.11 [14] A. E. Green, P. M. Naghdi, On undamped heat waves in an elastic solid, J. Thermal Stresses12 15 (1992) 253-264.13 [15] A. E. Green, P. M. Naghdi, Thermoelasticity without energy dissipation, J. Elasticity 3114 (1993) 189–208.15 [16] M. A. Hader, M. A. Al-Nimr, B. A. Abu Nabah, The Dual-Phase-Lag Heat Conduction Model16 in Thin Slabs Under a Fluctuating Volumetric Thermal Disturbance, Int. J. Thermophys. 2317 (2002) 1669–1680.18 [17] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (3-4) (2012) 251-265.19 [18] A. Maga˜na, A. Miranville, R. Quintanilla, On the stability in phase-lag heat conduction with20 two temperatures, J. Evol. Equations 18 (2018) 1697–1712.21 [19] A. Maga˜na, A. Miranville, R. Quintanilla, On the time decay in phase-lag thermoelasticity22 with two temperatures, Electronic Res. Arch. 27 (2019) 7–19.23 [20] A. Miranville, R. Quintanilla, A phase-field model based on a three-phase-lag heat conduction,24 Appl. Math. Optim. 63 (2011) 133-150.25 [21] S. Mukhopadhyay, R Prasad, R. Kumar, On the theory of Two-Temperature Thermoelaticity26 with Two Phase-Lags, J. Thermal Stresses 34 (2011) 352-365.27 [22] M. A. Othman, W. M. Hasona, E. M. Abd-Elaziz, Effect of rotation on micropolar generalized28 thermoelasticity with two temperatures using a dual-phase-lag model, Canadian J. Physics29 92 (2014) 149-158.30 [23] R. Quintanilla, Exponential stability in the dual-phase-lag heat conduction theory, J. Non-31 Equilib. Thermodyn. 27 (2002) 217–227.32 [24] R. Quintanilla, A Well-Posed Problem for the Dual-Phase-Lag Heat Conduction, J. Thermal33 Stresses 31 (2008) 260-269.34 [25] R. Quintanilla, A Well-Posed Problem for the Three-Dual-Phase-Lag Heat Conduction, J.35 Thermal Stresses 32 (2009) 1270-1278.36 [26] R. Quintanilla, P. M. Jordan, A note on the two-temperature theory with dual-phase-lag37 decay: some exact solutions, Mech. Research Comm. 36 (2009) 796–803.38 [27] R. Quintanilla, R. Racke, Qualitative aspects in dual-phase-lag thermoelasticity, SIAM J.39 Appl. Math. 66 (2006) 977-1001.40 [28] R. Quintanilla, R. Racke, A note on stability of dual-phase-lag heat conduction, Int. J. Heat41 Mass Transfer 49 (2006) 1209-1213.42 [29] R. Quintanilla, R. Racke, Qualitative aspects in dual-phase-lag heat conduction, Proc. Royal43 Society London A. 463 (2007) 659-674.44 [30] R. Quintanilla, R. Racke, A note on stability in three-phase-lag heat conduction, Int. J. Heat45 Mass Transfer 51 (2008) 24–29.46 [31] R. Quintanilla, R. Racke, Spatial behavior in phase-lag heat conduction, Differ. Integral Equ.47 28 (2015) 291-308.48 [32] S. A. Rukolaine, Unphysical effects of the dual-phase-lag model of heat conduction, Int. J.49 Heat Mass Transfer 78 (2014) 58–63.50 [33] D. Y. Tzou, A unified approach for heat conduction from macro to micro-scales, ASME J.51 Heat Transfer 117 (1995) 8–16.52 [34] W. E. Warren, P. J. Chen, Wave propagation in two temperatures theory of thermoelaticity,53 Acta Mech. 16 (1973) 83–117.54 [35] Y. Zhang, Generalized dual-phase lag bioheat equations based on nonequilibrium heat transfer55 in living biological tissues, Int. J. Heat Mass Transfer 52 (2009) 4829–4834.56 Received xxxx 20xx; revised xxxx 20xx.57