Analysis of a Poro-Thermo-Viscoelastic Model of Type III

: In this work, we numerically study a thermo-mechanical problem arising in poro-viscoelasticity with the type III thermal law. The thermomechanical model leads to a linear system of three coupled hyperbolic partial differential equations, and its weak formulation as three coupled parabolic linear variational equations. Then, using the ﬁnite element method and the implicit Euler scheme, for the spatial approximation and the discretization of the time derivatives, respectively, a fully discrete algorithm is introduced. A priori error estimates are proved, and the linear convergence is obtained under some suitable regularity conditions. Finally, some numerical results, involving one- and two-dimensional examples, are described, showing the accuracy of the algorithm and the dependence of the solution with respect to some constitutive parameters.


Introduction
Since the first works by Nunziato and Cowin [1,2], many research papers have been published involving the theory of elasticity of voids (see, for instance, [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and also the monographs [21,22]). The main idea of such a theory is the assumption that the mass at each material point is found as the product of the mass density by the volume fraction, adding a new variable in the constitutive equations. This theory has become very interesting because it has been useful in applications appearing in solid mechanics (rocks, woods or even bones).
In this paper, we also consider an improvement of the classical Fourier heat law to remove the well-known paradox of infinite wave speed. It is based on the work by Green and Nagdhi [23], where the thermal displacement was included as a new independent variable into the model (even of the elastic displacement), and it is usually called type III thermal law. Since then, many authors worked on this theory (see, e.g., [4,[24][25][26][27][28][29]). Here, our aim is to extend the results provided in [30], where a general thermoelastic law was studied from both mathematical (existence and uniqueness) and numerical (stability, a priori error estimates) points of view. We will consider the viscoelastic case and we will also include the porosity into the model. Therefore, using the well-known finite element method and the implicit Euler scheme, for the spatial approximation and the discretization of the time derivatives, we present a fully discrete algorithm, we perform an a priori error analysis, which leads to the linear convergence of the algorithm under some regularity conditions on the continuous solution, and we perform some numerical simulations, in one and two dimensions, to show the accuracy of the algorithm and the dependence of the solution on some constitutive parameters.
where u 0 , ϕ 0 and ψ 0 represent initial conditions for the displacement, the volume fraction and the thermal displacement. Therefore, the thermomechanical problem of a centrosymmetric, isotropic and homogeneous thermo-viscoelastic body within the type III thermal theory is written as follows (see [26,27,32]).
Here, v 0 , e 0 and θ 0 are given initial conditions for the velocity, the volume fraction speed and the temperature, respectively. As usual, the summation over repeated indexes is assumed, the partial derivative with respect to a variable is represented by a subscript preceded by a comma and one superposed dot denotes the first-order partial time derivative (two superposed dots represent the second-order partial time derivative).
In Equations (2)-(4), ρ is the mass density, J is the product of the mass density by the equilibrated inertia, λ and µ denote Lame's coefficients, λ * and µ * represent viscosity coefficients, a denotes the heat capacity, a 0 is the volume fraction diffusion, κ is the thermal diffusion and κ * is the viscous thermal diffusion. The remaining coefficients γ, β, ξ, m, d are model parameters. We also note that, for the sake of simplicity, we have neglected in Equations (2) and (4) the respective terms corresponding to the volume forces and the heat flux.
Problem VP. Find the velocity v : [0, T] → V, the volume fraction speed e : [0, T] → E and the temperature θ : [0, T] → E such that v(0) = v 0 , e(0) = e 0 , θ(0) = θ 0 and, for a.e. t ∈ (0, T), where the displacement, the volume fraction and the thermal displacement are then recovered from relations (1). Proceeding as in [27,31], we can state the existence of a unique solution to Problem VP. Details are omitted for the sake of reading. Theorem 1. Assume that the coefficients satisfy conditions (8) and the following regularity on the initial data: Then, Problem VP admits a unique solution with the regularity:

Fully Discrete Approximations: An a Priori Error Analysis
In this section, a finite element algorithm is shown for approximating solutions to Problem VP. This is done in two steps. First, to approximate the variational spaces V and E we define the finite element spaces V h and E h as follows, where we assume that Ω is a polyhedral domain and we denote by T h a regular triangulation of Ω. Moreover, the space of polynomials of global degree less or equal to 1 in element Tr is represented by P 1 (Tr). Here, h > 0 denotes the spatial discretization parameter. Secondly, in order to discretize the time derivatives we use a uniform partition of the time interval [0, T], that we denote by 0 = t 0 < t 1 < . . . < t N = T (k = T/N is the time step size). Moreover, for a continuous function f (t) 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 corresponding divided differences.
Using the classical implicit Euler scheme, the fully discrete approximation of Problem VP is the following.
Problem VP hk . Find the discrete velocity v hk = {v hk n } N n=0 ⊂ V h , the discrete volume fraction speed e hk = {e hk n } N n=0 ⊂ E h and the discrete temperature θ hk = {θ hk n } N n=0 ⊂ E h such that v hk 0 = v h 0 , e hk 0 = e h 0 , θ hk 0 = θ h 0 , and, for n = 1, . . . , N, where the discrete displacement, the discrete volume fraction and the discrete thermal displacement are then recovered from the relations: and the discrete initial conditions, denoted by u h 0 , v h 0 , ϕ h 0 , e h 0 , ψ h 0 and θ h 0 are given by Here, P 1h and P 2h are the respective projection operators over the finite element spaces V h and E h (see [33]).
It is easy to show that discrete problem VP hk admits a unique solution using well-known results on linear variational equations and conditions (8), so we omit the details.
The aim of this section is to obtain some a priori error estimates on the numerical errors u n − u hk n , v n − v hk n , ϕ n − ϕ hk n , e n − e hk n , ψ n − ψ hk n and θ n − θ hk n . First, we obtain the estimates on the velocity field. Subtracting Equation (9) at time t = t n for a test function w = w h ∈ V h ⊂ V and discrete Equation (15), we have, for all w h ∈ V h , and so, It is straightforward to show that where we have used that v hk n = δu hk n = (u hk n − u hk n−1 )/k and notations δv n = (v n − v n−1 )/k and δu n = (u n − u n−1 )/k. Therefore, we find that, for all Here and in what follows, C > 0 denotes a positive constant whose value may change for each expression, which depends on the continuous solution but it is independent of the discretization parameters h and k.
Secondly, we will find the error estimates on the volume fraction speed. Then, we subtract Equation (10) at time t = t n for a test function r = r h ∈ E h ⊂ E and discrete Equation (16) to obtain, for all r h ∈ E h , Thus, we find, for all r h ∈ E h , where we have used the fact that e hk n = δϕ hk n = (ϕ hk n − ϕ hk n−1 )/k and the notations δe n = (e n − e n−1 )/k and δϕ n = (ϕ n − ϕ n−1 )/k, it follows that, for all Finally, we get the estimates on the temperature field. Then, we subtract Equation (11) at time t = t n for a test function z = z h ∈ E h ⊂ E and discrete Equation (17) and so, it follows that, for all z h ∈ E h , Keeping in mind that where we have used the fact that θ hk n = δψ hk n = (ψ hk n − ψ hk n−1 )/k and the notations δθ n = (θ n − θ n−1 )/k and δψ n = (ψ n − ψ n−1 )/k, we obtain, for all z h ∈ E h , Combining now estimates (20)- (22) we find that Taking into account that where we have employed conditions (8), multiplying the previous estimates by k and summing up to n, it follows that, for all Keeping in mind again assumptions (8), we find that there exist two constants ζ 1 , ζ 2 > 0 such that γ/ξ < ζ 1 < (λ + µ)/γ and m/κ < ζ 2 < a 0 /m. Then, we have We will need the following discrete version of Gronwall's lemma (see, for instance, [34,35]). Lemma 1. Let {E n } N n=0 and {G n } N n=0 be two sequences of nonnegative real numbers satisfying, for a positive constant C 1 > 0 independent of G n and E n , where k is a positive constant. If we denote C = C 1 (1 + C 1 Te 2C 1 T ) and T = Nk, then we have applying Lemma 1 to estimates (23), we conclude the following a priori error estimates result.

Theorem 2.
Under the assumptions of Theorem 1, let us denote by (v, e, θ) the solution to Problem VP and by (v hk , e hk , θ hk ) the solution to Problem VP hk , then we obtain the following a priori error estimates, for all The previous error estimates can be used to derive the convergence order. As a particular case, under suitable additional regularity, the linear convergence is derived. Corollary 1. Let the assumptions of Theorem 2 still hold. If we assume that the solution to Problem VP has the additional regularity: and we use the finite element spaces V h and E h defined in (13) and (14), respectively, and the discrete initial conditions u h 0 , v h 0 , ϕ h 0 , e h 0 , ψ h 0 and θ h 0 given in (19), there exists a positive constant C > 0, independent of the discretization parameters h and k but depending on the continuous solution, such that

Numerical Simulations
In this final section, we will describe the algorithm and some numerical results involving oneand two-dimensional examples.

Numerical Algorithm
As a first step, given the solution u hk n−1 , v hk n−1 , ϕ hk n−1 , e hk n−1 , ψ hk n−1 and θ hk n−1 at time t n−1 , the velocity v hk n , the volume fraction speed e hk n and the temperature θ hk n are obtained from the coupled Equations (15)- (17). The corresponding linear system can be expressed in terms of a product variable x n and the resulting symmetric linear system is solved by using the well-known Cholesky method. Finally, the displacement field u hk n , the volume fraction ϕ hk n and the thermal displacement ψ hk n are updated using (18).
The numerical scheme was implemented on a Intel Core 3.2 GHz PC using MATLAB, and a typical one-dimensional run (h = k = 10 −3 ) took about 0.957 seconds of CPU time, meanwhile a two-dimensional run took about 2.74 s of CPU time.

A One-Dimensional Example: Numerical Convergence
To show the accuracy of the finite element approximations studied in the previous section, we will consider the following simpler one-dimensional problem: Problem P 1D .
where functions F i , i = 1, 2, 3, are given by, for all (x, t) ∈ (0, 1) × (0, 1), We point out that Problem P 1D is the one-dimensional version of Problem P using the following data: and all the initial conditions equal to x(x − 1).
We note that, in this case, the exact solution to Problem P 1D can be easily found and it has the following form: The numerical errors, for several values of the discretization parameters h and k, and given by are shown in Table 1. Moreover, the evolution of the error depending on the parameter h + k is depicted in Figure 1. The convergence of the numerical approximations is clearly observed, and its linear convergence, stated in Corollary 1, seems to be achieved.  If we assume now that there are not volume forces, and we use the following data: and the initial conditions, for all x ∈ (0, 1), taking the discretization parameters h = k = 0.001, the evolution in time of the discrete energy E hk n defined by is shown in Figure 2 using both natural and semi-log scales. We can see that the discrete energy converges to zero and an exponential decay seems to be achieved.

First Two-Dimensional Example: Application of a Surface Force
As a first two-dimensional example, we consider the square domain [0, 1] × [0, 1], which is assumed to be clamped on its left part {0} × [0, 1]. We also suppose that the thermal displacement and porosity vanish on the whole boundary, and we use the following expression for the mechanical surface force: f F (x, y, t) = (−ty 2 , 0) for x = 1, y ∈ (0, 1), which is assumed to be applied on the boundary Γ F = {1} × (0, 1), a part of the whole boundary ∂Ω. Even if this case was not studied in the previous sections, it is straightforward to extend the analysis to this more general case.
The following data have been employed in this simulation: and null initial conditions for all the variables. Taking the time discretization parameter k = 0.01, the deformation generated by this mechanical force is shown in Figure 3. In Figure 4 we plot the displacement in vector arrows (on the left-hand side) and its norm (on the right-hand side) at final time. As can be seen, there is a deformation along the horizontal axis. Moreover, the porosity (left) and the thermal displacement (right) are plotted in Figure 5 at final time. We note that both neglect on the boundary due to the null boundary conditions and they are generated by the deformation of the body.

Second Two-Dimensional Example: Dependence on the Type III Thermal Coefficient
In this second two-dimensional example, we will investigate the dependence on the type III thermal coefficient κ * .

Conclusions
In this paper, we numerically analyzed a dynamic poro-thermo-viscoelastic problem within the type III thermoelasticity. The weak form led to a linear system composed of parabolic variational equations written in terms of the velocity, the volume fraction speed and the temperature. Then, using the finite element method to approximate the spatial variable and the implicit Euler scheme to discretize the time derivatives, we introduced a fully discrete scheme. A priori error estimates were proved by using a discrete version of Gronwall's inequality. Finally, we presented a one-dimensional numerical simulation to show the convergence of the algorithm and the exponential decay of the discrete energy (Example 1), the effect of the application of a surface force (Example 2) and the dependence on the type III thermal coefficient (Example 3).