A dual-phase-lag porous-thermoelastic problem with microtemperatures

: In this work, we consider a multi-dimensional dual-phase-lag problem arising in porous-thermoelasticity with microtemperatures. An existence and uniqueness result is proved by applying the semigroup of linear operators theory. Then, by using the finite element method and the Euler scheme, a fully discrete approximation is numerically studied, proving a discrete stability property and a priori error estimates. Finally, we perform some numerical simulations to demonstrate the accuracy of the approximation and the behavior of the solution in one-and two-dimensional problems.


Introduction
In the second part of the past century, there were developed different theories where the microstructure of elastic materials was taken into account. Two good books concerning these aspects were written by Eringen [1] and Ieşan [2]. In particular, Goodman and Cowin [3] established the theory of continuum for granular materials with intersticial voids. They proposed the bulk density of the materials as a product of the matrix density and the volume fraction. Later, Cowin and Nunziato [4][5][6] stated a theory of porous elastic solids modeling the deformations of materials with small voids distributed within them. Of course, this theory was extended to include the thermal effects [7]. Contributions in this theory became huge and we can recall a few of them [8][9][10][11][12][13][14][15][16][17][18][19]. In general, the microstructure components of elastic materials have deserved much attention. One of the possibilities could be the so-called "microtemperatures" [20,21].

Existence and uniqueness
In this section, we give an existence and uniqueness result for the problem determined by system (2.1)-(2.4), boundary conditions (2.5) and initial conditions (2.6).
To this end we will work in the Hilbert space: we define the inner product: R. It is worth noting that this inner product defines a norm which is equivalent to the usual one in the Hilbert space because of the previous assumptions.
Our problem can be written as where I is the identity operator and We note that the domain of our operator is the subset (u, v, ϕ, ψ, θ, ϑ, ξ, T, S, R) ∈ H such that It is clear that this is a dense subspace of our Hilbert space. On the other side, we have that for every U in the domain of the operator denoted by D(A).
To show that A generates a contractive semigroup, it is sufficient to prove that zero belongs to the resolvent of the operator.
To this end, we consider We need to prove that the equation has a solution in the domain of the operator. That is, we must solve the following system: We obtain the solution for v, ψ, R and ξ satisfying the regularity conditions and we also find that If we consider the last two equations we see that the right-hand side belongs to W −1,2 (B) × W −1,2 (B). On the other side, the bilinear form defines a bounded, and coercive, bilinear form in W 1,2 0 (B) × W 1,2 0 (B). In view of the Lax-Milgram lemma, we conclude the existence of (θ, T) ∈ W 1,2 0 (B) ×W 1,2 0 (B) satisfying the last two equations. Now, if we introduce this solution into the first two ones we obtain Au A similar argument shows the existence of (u, ϕ) ∈ W 1,2 0 (B) × W 1,2 0 (B). In fact, we can also prove an inequality of the type ∥U∥ ≤ C * ∥F∥ for a given positive constant C * > 0, and so, we have proved the following. Therefore, we can conclude the existence and uniqueness of solution to problem (2.1)-(2.6), that we state as follows.

Fully discrete approximations: a priori error analysis
In this section, we consider a fully discrete aproximation to a variational problem of the above thermomechanical problem by using the finite element method and the implicit Euler scheme.
Let Y = L 2 (B), H = [L 2 (B)] d and Q = [L 2 (B)] d×d , and denote by (·, ·) Y , (·, ·) H and (·, ·) Q the respective scalar products in these spaces, with corresponding norms ∥ · ∥ Y , ∥ · ∥ H and ∥ · ∥ Q . Moreover, let the variational spaces E and V be given as with respective scalar products (·, ·) E and (·, ·) V , and norms ∥ · ∥ E and ∥ · ∥ V . By using Green's formula and taking into account the above boundary conditions, we write the variational formulation of the thermomechanical problem in terms of variables v = (v i ) =u = (u i ) and ψ =φ, the temperature speed ϑ =θ, the temperature acceleration ξ =θ, the microtemperature speed where functions u, ϕ, ϑ, θ, S and T are then recovered from the relations: and T f > 0 denotes the final time of interest. Thus, we construct the finite element spaces V h ⊂ Vand E h ⊂ E given by where B is assumed to be a polyhedral domain, T h denotes a regular triangulation, in the sense of [29], of B, and P 1 (T r) represents the space of polynomials of global degree less or equal to 1 in element T r.
Here, h > 0 denotes the spatial discretization parameter. In order to discretize the time derivatives, we use a uniform partition of the time interval [0, and, for the sequence {z n } N n=0 , we denote by δz n = (z n − z n−1 )/k its corresponding divided differences.
Using the classical implicit Euler scheme, the fully discrete approximation of Problem VP is the following. Problem where discrete functions u hk n , ϕ hk n , ϑ hk n θ hk n , S hk n and T hk n are then recovered from the relations: and R h 0 are defined as follows: In the previous definitions, operators P 1h and P 2h are the projection operators over the finite element spaces V h and E h , respectively (see, for instance, [30]).
Using the previously given assumptions (2.7), applying the well-known Lax-Milgram lemma we can easily show that fully discrete problem V P hk has a unique solution.
In this section, the objective is to prove a main error estimates result regarding the approximation of the solution to Problem VP by the solution to Problem V P hk . So, first we have the following discrete stability property.
generated by Problem V P hk , 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 simplicity, in this proof we need to assume that τ 2 1 /2 = 1. We note that we can extend the analysis provided below to the general case doing some simple modifications.
First, we estimate the terms on the discrete function v hk n . Taking w h = v hk n as a test function in discrete variational equation (4.8) we obtain Now, we estimate the terms on the discrete function ψ hk n . Taking r h = ψ hk n as a test function in discrete variational equation (4.9) we have and, using the estimates Thirdly, we get the estimates on the discrete temperature acceleration ξ hk n . Thus, taking z h = ξ hk n as a test function in discrete variational equation (4.10) we obtain and therefore, taking into account that Finally, we get the estimates on the discrete microtemperatures acceleration R hk n . Then, taking η h = R hk n as a test function in discrete variational equation (4.11) it follows that b(δR hk n + τ 1 R hk n + S hk n , R hk n ) H + κ 6 (∇(T hk n + τ 2 S hk n ), ∇R hk n ) Q +(κ 4 + κ 5 )(div (T hk n + τ 2 S hk n ), div R hk n ) Y + κ 2 (T hk n + τ 2 S hk n , R hk n ) H = −κ 1 (∇(θ hk n + τ 1 ϑ hk n ), R hk n ) H − µ 2 (∇ψ hk n , R hk n ) H , and so, keeping in mind that Combining all these estimates, it leads to Multiplying the above estimates by k and summing up to n we get From the definition of θ hk n and T hk n we can easily show that Observing that and, using a discrete version of Gronwall's inequality (see [31]), we conclude the discrete stability property. □ Now, we will derive a main error estimates result on the numerical errors v n − v hk n , ψ n − ψ hk n , ξ n − ξ hk n and R n − R hk n . We have the following. Theorem 4.2. Under the assumptions of Lemma 4.1, if we denote by (v, ψ, ξ, R) the solution to Problem V P and by (v hk , ψ hk , ξ hk , R hk ) the solution to Problem V P hk , then we have the following a priori error estimates, for all where C > 0 is a positive constant which is be independent of the discretization parameters h and k, but depending on the continuous solution, (4.16) Proof. As we did in the proof of Lemma 4.1, we also assume that the time relaxation parameter τ 2 First, we obtain a priori error estimates on function v. Then, subtracting variational equation (4.1) at time t = t n for a test function w = w h ∈ V h ⊂ V and discrete variational equation (4.8), we have, for and so, we obtain Keeping in mind that using inequality (4.14) several times, after some easy calculations we get Secondly, we obtain the error estimates on function ψ. If we subtract variational equation (4.2) at time t = t n for a test function r = r h ∈ E h ⊂ E and discrete variational equation (4.9), we find that, for Then, we have, for all r h ∈ E h , and, taking into account that using again several times inequality (4.14) we get, for all r h ∈ E h , Thirdly, we derive the estimates on the temperature acceleration ξ. Therefore, subtracting variational equation (4.3) at time t = t n for a test function z = z h ∈ E h ⊂ E and discrete variational equation (4.10), for all z h ∈ E h we find that a(ξ n − δξ hk n + τ 1 (ξ n − ξ hk n ) + ϑ n − ϑ hk n , z h ) Y + κ(∇(θ n − θ hk n + τ 2 (ϑ n − ϑ hk n )), Thus, it follows that a(ξ n − δξ hk n + τ 1 (ξ n − ξ hk n ) + ϑ n − ϑ hk n , ξ n − ξ hk Using the following estimates applying several times inequality (4.14) we find that Finally, we get the error estimates for the microtemperature acceleration R. Then, if we subtract variational equation (4.4) at time t = t n for a test function η = η h ∈ V h ⊂ V and discrete variational equation (4.11), we have, for all ψ h ∈ V h , b(Ṙ n − δR hk n + τ 1 (R n − R hk n ) + S n − S hk n , η h ) H +(κ 4 + κ 5 )(div (T n − T hk n + τ 2 (S n − S hk n )), div η h ) Y +κ 6 (∇(T n − T hk n + τ 2 (S n − S hk n )), ∇η h ) Q + κ 2 (T n − T hk n + τ 2 (S n − S hk n ), η h ) H +κ 1 (∇(θ n − θ hk n + τ 1 (ϑ n − ϑ hk n )), η h ) H + µ 2 (∇(ψ n − ψ hk n ), η h ) H = 0, and therefore, we obtain b(Ṙ n − δR hk n + τ 1 (R n − R hk n ) + S n − S hk n , R n − R hk n ) H +κ 6 (∇(T n − T hk n + τ 2 (S n − S hk n )), ∇(R n − R hk n )) Q +(κ 4 + κ 5 )(div (T n − T hk n + τ 2 (S n − S hk n )), div (R n − R hk n )) Y +κ 2 (T n − T hk n + τ 2 (S n − S hk n ), Keeping in mind that using again inequality (4.14), after easy calculations it leads to Combining estimates (4.17)-(4.20) it follows that, for all w h , η h ∈ V h and r h , Multiplying the above estimates by k and summing up to n it follows that Taking into account that where similar estimates can be derived for the terms involving functions ψ, ξ and R, and using the estimates: where I n and J n are the integration errors defined in (4.15) and (4.16), respectively, applying a discrete version of Gronwall's inequality (see, again, [31]) we conclude the desired a priori error estimates. □ We note that the error estimates provided in Theorem 4.2 can be used to obtain the convergence order of the approximations provided by the fully discrete problem V P hk . As a particular case, if we assume the additional regularity:   .21), it follows the linear convergence of the approximations given by Problem VP hk ; that is, there exists a constant C > 0, independent of parameters h and k, such that

Numerical results
In this section, we describe some numerical results obtained solving one-and two-dimensional problems. First, the numerical convergence and the discrete energy decay are shown in the onedimensional example. Secondly, the numerical dependence on the coupling parameter µ 2 is studied in a two-dimensional example.

First example: numerical convergence in a one-dimensional example
As an academical example, in order to show the accuracy of the approximations the following onedimensional problem is considered: The following data have been used: and the artificial supply terms F i , i = 1, 2, 3, 4, are given by, for a.e. (x, t) ∈ (0, 1) × (0, 1), We note that the numerical scheme was implemented on a 3.2 Ghz PC using MATLAB, and a typical run (h = k = 0.001) took about 0.13 seconds of CPU time.
In this case, the exact solution to this one-dimensional problem can be easily calculated and it has the following form, for (x, t) ∈ [0, 1] × [0, 1]: Thus, the approximation errors estimated by max 0≤n≤N ∥v n − v hk n ∥ Y + ∥u n − u hk n ∥ E + ∥ϕ n − ϕ hk n ∥ E + ∥ψ n − ψ hk n ∥ Y + ∥θ n − θ hk n ∥ E +∥ϑ n − ϑ hk n ∥ E + ∥ξ n − ξ hk n ∥ Y + ∥T n − T hk n ∥ E + ∥S n − S hk n ∥ E + ∥R n − R hk n ∥ Y  are presented in Table 1 for different values of h and k. Moreover, the evolution of the error depending on the parameter h + k is plotted in Figure 1. We notice that the convergence of the algorithm is clearly observed, and the linear convergence, stated in Theorem 4.3, is achieved. If we assume now that there are not supply terms and we use the final time T f = 30, the same data than in the previous example and the initial conditions: u 0 = v 0 = ϕ 0 = ψ 0 = θ 0 = ϑ 0 = ξ 0 = 0, T 0 (x) = R 0 (x) = S 0 (x) = x(x − 1) for a.e. x ∈ (0, 1), taking the discretization parameters h = k = 0.001, the evolution in time of the discrete energy E hk n = ||v hk n || 2 Y + 2||u hk n || 2 E + ∥ψ hk n ∥ 2 Y + ∥ϕ hk n ∥ 2 E + ∥ξ hk n ∥ 2 Y + 3∥ϑ hk n ∥ 2 E + 3∥θ hk n ∥ 2 E +2∥R hk n ∥ 2 Y + 2∥S hk n ∥ 2 E + 2∥T hk n ∥ 2 E is plotted in Figure 2 (in both natural and semi-log scales). As can be seen, it converges to zero and an exponential decay seems to be achieved. We consider now the two-dimensional domain B = (0, 1) × (0, 1). In this case, our aim is to study the dependence on the coupling parameter µ 2 .
Using the time discretization parameter k = 0.01 and the uniform finite element mesh obtained dividing the unit square into 2048 triangles, in Figure 3 we plot the norm of function u and the microtemperatures, at middle point x = (0.5, 0.5), for some values of parameter µ 2 . It seems that, for the highest value µ 2 = 20, an oscillation is produced.  Finally, in Figure 4 we plot function ϕ and the temperature θ at final time for the value of parameter µ 2 = 20. We can note that they have been generated by the resulting deformation.