Quasistatic Porous-Thermoelastic Problems: An a Priori Error Analysis

: In this paper, we deal with the numerical approximation of some porous-thermoelastic problems. Since the inertial effects are assumed to be negligible, the resulting motion equations are quasistatic. Then, by using the ﬁnite element method and the implicit Euler scheme, a fully discrete approximation is introduced. We prove a discrete stability property and a main error estimates result, from which we conclude the linear convergence under appropriate regularity conditions on the continuous solution. Finally, several numerical simulations are shown to demonstrate the accuracy of the approximation, the behavior of the solution and the decay of the discrete energy.


Introduction
Since the first work by Cosserat and Cosserat [1], who introduced the theory of the so-called micropolar materials, there were several models to extend it to other type of elastic materials. For instance, assuming that the bulk density is equal to the product of the density matrix by the volume fraction, Cowin and Nunziato [2,3] introduced the theory of elastic solids with voids, where the pores are distributed within the material. Rocks, woods, biological materials (as bones), or soils are some examples of this kind of materials. Later, heat effects were also included and the number of papers dealing with these porous-thermoelastic problems is really huge (see, among others, [4][5][6][7][8][9]), and their numerical issues [10][11][12][13][14][15].
It is well-known that there is a paradox of instantaneous propagation of the thermal waves in the classical heat theory provided by the Fourier law. In order to overcome it, several proposals were made since 1990. One of them was proposed by Green and Naghdi [16], where a new variable, the thermal displacement, was considered. Here, we will study the type II theory, which is also called thermoelasticity without energy dissipation.
In this paper, we also assume that the process is quasi-static and so, we neglect the inertia term in the motion equation. This expression was intorduced by Mosconi [17] and we note that the system of equations that we obtain differs substantially from the original one, leading to a mixture of some parabolic and hyperbolic equations. One interesting question is, therefore, if the behavior of the solution with respect to the time variable remains unaltered in the quasi-static case. The number of contributions has increased over the last years (see, among others, [18][19][20][21]). In a recent article by Magaña and Quintanilla [22], they considered that this quasi-static hypothesis was also imposed in the volume fraction and in the temperature. So, it led to a set of different problems for which they proved the existence and uniqueness of solution and the slow or exponential energy decay. Here, we continue this research, aiming to provide an a priori error analysis of a fully discrete approximation and to perform some numerical simulations.
The paper is structured in the following way. In Section 2, we describe the thermomechanical problem following [22] and we recall an existence and uniqueness result. Then, in Section 3, we present the discrete problem, using the finite element method to approximate the spatial variable and the implicit Euler scheme to discretize the time derivatives. A discrete stability property and a priori error estimates are proved, from which the linear convergence of the algorithm is deduced under adequate additional regularity. Finally, in Section 4, some numerical results are presented to demonstrate the accuracy of the algorithm and the behavior of the solution.

The Thermo-Mechanical Problem and Its Variational Formulation
In this section, we describe briefly the problem and we recall an existence and uniqueness result (see [22] for further details).
Without loss of generality, we suppose that the spatial variable x lies in the interval [0, ] and that the time t goes from 0 to T, where T > 0 denotes the final time. Moreover, the subscript x indicates the spatial derivative and a dot over a variable represents the time derivative.
Defining the displacement field, the porosity field and the temperature field by u, ϕ and θ, respectively, the porous-thermo-mechanical problem is written as follows.
Find the displacement field u : where µ * , µ, b, β, J, δ, ξ, m, c and κ * are constitutive parameters satisfying some conditions which will be detailed later, and u 0 , v 0 , ϕ 0 , ψ 0 and θ 0 are given initial conditions. The influence of the dissipation constants µ * and κ * will be numerically studied in Section 4.2.
In what follows, we obtain the variational formulation of the thermo-mechanical problem (1)- (6). Let Y = L 2 (0, ), and denote by (·, ·) the scalar product in this space, with corresponding norm · Y . Moreover, let us define V = H 1 0 (0, ) with norm · V . Let us denote by v =u and ψ =φ the velocity and porosity speed. Therefore, integrating by parts it leads to the following variational formulation of problem (1)- (6).
Find the velocity field v : [0, T] → V, the porosity speed ψ : [0, T] → V and the temperature θ where the displacement and the porosity fields u and ϕ are recovered from the equations: The following existence and uniqueness result has been recently proved in [22].

Theorem 1.
Assume that the constitutive coefficients satisfy the following conditions: Then, there exists a unique solution to problem (1)- (6).
Regarding the energy decay, it was also proved in [22] that the energy of the system related to problem (1)-(6) decays in a slow way. However, if we replace the Fourier law by the type II thermal law (see Remark 1), then using the theory of linear semigroups it was shown in [22] that the energy decays in an exponential way.

An a Priori Error Analysis
In this section, we study a fully discrete approximation of the porous-thermo-mechanical problem (1)-(6) described in the previous section. First, in order to obtain the spatial approximation, we assume that the interval [0, ] is divided into M subintervals a 0 = 0 < a 1 < . . . < a M = of length h = a i+1 − a i = /M. Therefore, to approximate the variational space V, we define the finite dimensional space V h ⊂ V given by where P 1 ([a i , a i+1 ]) represents the space of polynomials of degree less or equal to one in the subinterval [a i , a i+1 ]; i.e., the finite element space V h is made of continuous and piecewise affine functions. Here, h > 0 denotes the spatial discretization parameter. Furthermore, let the discrete initial conditions u h 0 , v h 0 , ϕ h 0 , ψ h 0 and θ h 0 be defined as where P h is the classical finite element interpolation operator over V h (see [23]). Secondly, we consider a uniform partition of the time interval [0, T], denoted by 0 = t 0 < t 1 < . . . < t N = T, with step size k = T/N and nodes t n = n k for n = 0, 1, . . . , N.
Therefore, using the well-known implicit Euler scheme, the fully discrete approximations of the above variational problem are the following.
Find the discrete velocity v hk = {v hk n } N n=0 ⊂ V h , the discrete porosity speed ψ hk = {ψ hk n } N n=0 ⊂ V h and the discrete temperature and, for all w h , r h , s h ∈ V h and n = 1, . . . , N, where the discrete displacement and the discrete porosity u hk n and ϕ hk n are now recovered from the equations: We note that, using the well-known Lax Milgram lemma (see, for instance, [23]) and the required assumptions on the constitutive parameters, it is straightforward to obtain that discrete problem (11)-(14) has a unique solution.

Remark 1.
Proceeding as in the derivation of variational problems (7)-(10) and (11)- (14) we could study the rest of the problems considered in [22]. In particular, we could study the following problem which involves the type II thermal law, which is written as: where α represents the thermal displacement which satisfiesα = θ, l is a new constitutive parameter, and u 0 , v 0 , ϕ 0 , ψ 0 , α 0 and θ 0 are given initial conditions. Therefore, the resulting variational problem is written as follows.
Find the velocity field v : where the displacement, the porosity and the thermal displacement u, ϕ and α are recovered from the equations: So this problem is approximated by the discrete problem defined by the following. Find the discrete velocity v hk = {v hk n } N n=0 ⊂ V h , the discrete porosity speed ψ hk = {ψ hk n } N n=0 ⊂ V h and the discrete temperature and, for all w h , r h , s h ∈ V h and n = 1, . . . , N, where the discrete displacement u hk n , the discrete porosity ϕ hk n and the discrete thermal displacement α hk n are now recovered from the equations: We will prove now a discrete stability property for problem (11)- (14).

Lemma 1.
Let the assumptions of Theorem 1 hold. It follows that the sequences {u hk , v hk , ϕ hk , ψ hk , θ hk }, generated by discrete problem (11)- (14), satisfy the stability estimate: where C is a positive constant which is independent of the discretization parameters h and k.

Proof.
Taking v hk n as a test function in variational Equation (11) we find that and, using the Cauchy's inequality we get the following estimates for the velocity field: Now, we deal with the porosity speed. Thus, taking ψ hk n s a test function in variational Equation (12) we have taking into account that and using again inequality (15) after several simple calculations it leads to the following estimates Finally, we obtain the estimates on the discrete temperature. Therefore, taking θ hk n as a test function in variational Equation (13) we get and using inequality (15) several times and keeping in mind that we find the following estimates: Combining estimates (17)-(18) we obtain Multiplying the above estimates by k and summing up to n we have Combining estimates (16) and (19) it follows that Now, using the definition of the discrete displacements we easily find that and applying a discrete version of Gronwall's inequality (see, e.g., [24]) we conclude the desired stability property.
In the rest of the section, we will prove some a priori error estimates on the numerical errors v n − v hk n , u n − u hk n , ψ n − ψ hk n , ϕ n − ϕ hk n and θ n − θ hk n .

Theorem 2.
Let the assumptions of Theorem 1 still hold. If we denote by (u, v, ϕ, ψ, θ) the solution to problem (7)-(10) and by (u hk , v hk , ϕ hk , ψ hk , θ hk ) the solution to problem (11)- (14), then we have the following a priori error estimates, for all where C is again a positive constant which does not depend on parameters h and k, and I n is the integration error given by Proof. As a first step, we get some error estimates for the velocity field. So, we subtract variational Equation (7) at time t = t n for a test function w = w h ∈ V h ⊂ V and discrete variational Equation (11) to obtain, for all w h ∈ V h , Therefore, it follows that, for all w h ∈ V h , . Using several times inequality (15) we get the following estimates for the velocity field, for all w h ∈ V h , Now, we obtain the error estimates for the porosity speed. Thus, subtracting variational Equation (8) at time t = t n for a test function r = r h ∈ V h ⊂ V and discrete variational Equation (12) we get, for all r h ∈ V h , and so, for all r h ∈ V h we have Keeping in mind that using inequality (15) we obtain the following estimates, for all r h ∈ V h , Finally, we obtain the error estimates for the temperature field. Then, we subtract variational Equation (9) at time t = t n for a test function s = s h ∈ V h ⊂ V and discrete variational Equation (13) we get, for all s h ∈ V h , and so, we find that, for all s h ∈ V h , Taking into account that Now, we combine estimates (22) and (23) to obtain Multiplying the above estimates by k and summing up to n it follows that and, together with estimates (21), it leads to Finally, keeping in mind that where I n is the integration error defined in (20), using again a discrete version of Gronwall's inequality (see [24]) we obtain the desired a priori error estimates.
The estimates provided in the above theorem can be used to obtain the convergence order of the approximations given by discrete problem (11)- (14). Hence, as an example, if we assume the additional regularity: we obtain the linear convergence of the algorithm applying some results on the approximation by finite elements (see [23]) and previous estimates already derived in [24], that is, we can prove that there exists a positive constant C > 0, independent of the discretization parameters h and k, such that max 0≤n≤N v n − v hk n V + u n − u hk n V + ψ n − ψ hk n Y + ϕ n − ϕ hk n V + θ n − θ hk n Y ≤ C(h + k).

Numerical Results
In this section, we will describe show some numerical solutions to the previously stated problems. We start by showing the numerical convergence of the algorithm to follow with some examples that check some theoretical statements made in [22]. In all the examples presented, we consider a domain of length = 1.

Numerical Convergence
To show the numerical convergence of the algorithm we propose an analytical solution to the problem such that: To ensure this solution is achieved, we introduce to the original problem artificial supply terms, defined as follows: For this verification of the convergence we choose a final time as T = 1 and the following parameters: After solving the problem for different discretization parameters h and k, we compute the numerical error at the final time step using the known solution: The errors obtained for the discretization parameters that were tested are summarized in Table 1. The main diagonal of this table is plotted in Figure 1. These results confirm the linear convergence of the algorithm employed.

Parametric Study of the Dissipation Parameters
In this problem there are two dissipation mechanisms present: one associated with µ * in Equation (1) and another associated with κ * in Equation (3). We will study the influence of these parameters on the decay velocity of the energy. Unless stated otherwise the values for the remaining parameters used for these cases are: µ = 10, b = 1, β = 1, J = 1, δ = 10, ξ = 1, m = 1, c = 1, µ * = 1, κ * = 1.
We recall from [22] the expression for the discrete energy of the system: We start by studying the influence of µ * . In Figure 2, we show the evolution of the energy with time for several values of this parameter in both natural (left) and semi-log (right) scales. It can be seen that, after some time, the influence is minimal, reaching a significantly small energy at around time t = 70. On the semi-log graph, the slopes for all the cases are very similar, indicating they have close exponents on the energy decay rate. We repeat the study for the second dissipation parameter κ * . We study now more values, since the variability is higher than in the previous case. We show our findings in Figure 3, again in natural and semi-log scales. We can see that for all the values of κ * , except for κ * = 1, the decay is slower than before. In addition, there is a maximum decay rate around the value of κ * = 1 for this selection of parameters, since both for higher and lower values than this one the decay rates are slower. This behavior is best shown in the semi-log scale, where the slope of the line associated with κ * = 1 is much higher than for the others.

Velocity of the Type II Model
In Theorem 3.7 of [22], the authors proved the exponential decay of the energy of the type II thermal law when coupled with the quasi-static displacement and the porosity.
They differentiate several cases, depending on the values of the parameters l and β. In particular, when l = 0, the decay is slow in general. However, for the particular case where J κ = c δ the authors suggest that the exponential decay is achieved again. In the discrete problem, it is not possible to see a non-exponential energy decay. However, we expect to see a difference in the decay rate when the previous condition is fulfilled.
In order to test the two different cases we choose different values for c. With c = 1 and c = 3, the condition is not verified, while for c = 2 it is. The results obtained are shown in Figure 4. In natural scale (left), it is not clear that the case where c = 2 is faster, but in the zoom around time 180 it can be seen that the dashed line in red overtakes the green line. The expected behavior is clearer in the semi-log figure, where the slope of the dashed red line is higher than the other two.  Figure 4. Evolution of the energy in natural (left) and semi-log (right) scales. In the case where c = 2 the decay velocity is expected to be higher because J κ = c δ.

Conclusions
In this paper, we have studied, from the numerical point of view, several porousthermoelastic problems assuming that the inertia effects are negligible (and so, the problems became quasi-static). As an example, we have focused on the study of the quasi-static porous-thermoelastic problem where the thermal law was modeled by the classical Fourier law. By using the implicit Euler scheme and the finite element method to discretize both the time derivatives and spatial variable, fully discrete approximations have been introduced. A discrete stability property and a priori error estimates have been proved, and the linear convergence has been deduced under some adequate additional regularity. Finally, we have provided some numerical simulations to show the numerical convergence of the algorithm in Section 4.1, where the linear convergence seems to be achieved. Then, in Section 4.2 a parametric study has been performed with respect to dissipative parameters µ * and κ * . Finally, the problem involving the so-called type II thermal law has been considered. A detailed analysis on some constitutive coefficients, as suggested in [22], has been shown.
At this time, the lack of experimental results does not allow to show which model is closer to reality. Both models behave similarly from a numerical point of view. In addition, the performance and complexity is similar in both cases. Therefore, the selection of the most appropriate model will depend on the specific application.