Decay for strain gradient porous elastic waves

We study the one-dimensional problem for the linear strain gradient porous elasticity. Our aim is to analyze the behavior of the solutions with respect to the time variable when a dissipative structural mechanism is introduced in the system. We consider five different scenarios: hyperviscosity and viscosity for the displacement component and hyperviscoporosity, viscoporosity and weak viscoporosity for the porous component. We only apply one of these mechanisms at a time. We obtain the exponential decay of the solutions in the case of viscosity and a similar result for the viscoporosity. Nevertheless, in the hyperviscosity case (respectively hyperviscoporosity) the decay is slow and it can be controlled at least by t-1/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{-1/2}$$\end{document}. Slow decay is also expected for the weak viscoporosity in the generic case, although a particular combination of the constitutive parameters leads to the exponential decay. We want to emphasize the fact that the hyperviscosity (respectively hyperviscoporosity) is a stronger dissipative mechanism than the viscosity (respectively viscoporosity); however, in this situation, the second mechanism seems to be more “efficient” than the first one in order to pull along the solutions rapidly to zero. This is a striking fact that we have not seen previously at any other linear coupling system. Finally, we also present some numerical simulations by using the finite element method and the Newmark-β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} scheme to show the behavior of the energy decay of the solutions to the above problems, including a comparison between the hyperviscosity and the viscosity cases.


Introduction
It is known that the porous structure of one material can have a significant influence in the behavior of this material when it is exposed, for instance, to deformations. This is one of the reasons why porous elastic solids have been extensively studied. Nunziato and Cowin [26] put forth a nonlinear theory in which the skeletal or matrix material is elastic and the interstices are void of material. Later on, Cowin and Nunziato [4] derived the linear theory and Cowin [5] analyzed its viscoelastic behavior. For a thorough review of this theory, we refer the reader to the book of Ieşan [12].
On the other hand, some authors proposed the inclusion of higher-order gradients in the basic postulates of elasticity in order to obtain more detailed models for the configuration of the materials and their response to stimuli. As a matter or illustration, we cite the works of Green and Rivlin [14], Mindlin [22] and Toupin [31]. The theories including the second gradient of the displacement or the second gradient of the volume fraction field in the set of independent constitutive variables are now called strain gradient theories.
In this work, we study the one-dimensional problem for the linear strain gradient porous elasticity, theory that has been recently proposed by Ieşan [13]. Our main purpose is to analyze and to quantify the damping speed of the waves when we attach different types of dissipation in the system. To simplify, we distinguish only between exponential or slow decay of the solutions. The decay is said to be exponential if the energy of the system can be controlled by means of a negative exponential on the time variable. Otherwise, the decay is said to be slow, including the case in which the energy can be controlled by the inverse of a rational function. (This is usually called polynomial decay.) A porous elastic structure is determined by a macroscopic component (the elastic deformation) and a microscopic one (the porosity). Both components are coupled. It is interesting to know if the inclusion of one dissipation mechanism in one of the components is able to carry the entire structure to a state of quick decay or not. In fact, the time behavior of the solutions depends on three issues: the theory we work with, the dissipation mechanism, and the coupling between the macroscopic and the microscopic components. The first results in this line of research were obtained in 2003 by Quintanilla [30]. Since then, a lot of contributions can be found in the literature (see [18][19][20][21]23,25,27,28] among others), but without considering the new strain gradient assumption. In the generic case, the exponential decay can be guaranteed by choosing two well-combined dissipation mechanisms. Nevertheless, there are some singular cases in which a single mechanism is enough to get it, but it requires that the velocities of the elastic and of the porous waves coincide [1]. It has been also proved that the introduction of a suitable conservative heat conduction in the system leads to the exponential decay with a single dissipation mechanism [15,24]. Other results in the same direction have been shown depending on the kind of kernel considered when the dissipation depends on the history [7][8][9][10].
We want to point out that in 2015 Liu, Magaña and Quintanilla made a first approach to the strain gradient situation [17]. They considered second-order derivatives for the displacement in the constitutive equation for the hyperstress but only first-order derivatives in the gradient of the volume fraction. Applying the basic properties of thermomechanics, the coupling between both components was determined. They also showed the exponential decay with hyperviscosity, with viscosity and also with viscoporosity, and the slow decay in the presence of weak viscoporosity (when the dissipation depends on the variation of the volume fraction). These results were quite surprising in comparison with the known results for the classical situation.
In this paper, we consider fourth-order derivatives with respect to the spatial variable in both components of the system. The axioms of thermomechanics determine again the coupling, but, strikingly, the behavior of the solutions changes. We obtain polynomial (slow) decay with hyperviscosity and exponential decay with viscosity (respectively, hyperviscoporosity and viscoporosity), and slow decay with weak viscoporosity in the generic case, although there is a particular combination of the parameters of the system that leads this situation also to the exponential decay. We believe that these results are noteworthy because intuition says that the hyperviscosity is a stronger dissipation mechanism than the viscosity but, nevertheless, it seems that what really matters is how these mechanisms are coupled.
The structure of the paper is as follows. In Sect. 2, we state the basic equations we are going to work with. We only state the conservative structure, and we impose the boundary and initial conditions we use in all the systems of equations that we study later on. In Sect. 3, we introduce what a priori seems to be a very strong dissipation mechanism in the elasticity part. We call it hyperviscoelasticity and we prove that the solutions decay in a slow way. (In fact, we show that the decay can be controlled by t −1/2 .) In Sect. 4, we change the damping mechanism: we take now the first derivative of the displacement velocity with respect to the spatial variable. Surprisingly, we obtain now the exponential decay of the solutions. Section 5 is devoted to obtain similar results but for the porosity component. Three different dissipation mechanisms are analyzed there: hyperviscoporosity, viscoporosity and weak viscoporosity. We find slow decay for the first case, exponential decay for the second and slow decay again for the third, although we obtain a specific combination of the constitutive parameters that leads to the exponential decay. Finally, in Sect. 6, we describe some numerical simulations of the problems involving hyperviscosity, viscosity and weak viscoporosity. We show the evolution of the discrete energy in the three cases, including a comparison between the mechanical dissipation mechanisms.

Basic equations
First of all, we recall the evolution and constitutive equations which govern the theory we are going to deal with. We follow the guidelines proposed by Ieşan [13]. As we consider several dissipation mechanisms, in this section we only state the conservative structure. Later, in each section, we write the constitutive equation (or equations) that we conveniently modify to introduce the dissipation.
Our analysis is focused on the one-dimensional problem, whose evolution equations are Here, u is the displacement, ϕ is the fraction of volume, τ is the stress, μ is the hyperstress, χ is the equilibrated stress vector, σ is the equilibrated hyperstress tensor and g is the equilibrated body force. As usual, ρ stands for the mass density and J for the product of the mass density by the equilibrated inertia, and both are assumed to be positive. The primary constitutive equations are given by The conditions for the constitutive coefficients a, b, β, k 1 , γ, α, d, k 2 and ξ will be stated in the following section. In fact, in view of the field equations, we will introduce some other notation to simplify the writing.
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 ∞.
The following set of boundary and initial conditions are imposed for all the different systems that we analyze: and, for a.e. x ∈ (0, π), We intend to introduce dissipative mechanisms in the system and determine how the solutions decay with respect to the time variable.

Hyperviscoelasticity
We introduce a first dissipative mechanism in the elasticity. We call it hyperviscoelasticity because it is the second derivative of the displacement velocity with respect to the spatial variable. To be precise, we assume that with k * 1 > 0, while the other constitutive equations remain unaltered. With the above assumption, if we substitute the constitutive equations into the evolution equations we obtain the system of field equations: where, to simplify the notation, we set η = γ − β and δ = α − 2d. The constitutive coefficients satisfy the following conditions: Basically, the above assumptions guarantee the elastic stability of the material. These conditions are assumed for all the systems we are going to study in this paper. The existence of solutions that do not decay is clear, but if the average of the initial condition for u 0 and ϕ 0 vanishes, then we avoid this possibility. Remark 3.1. It is not difficult to see that this system has undamped solutions. Take, for example, u = 0 and ϕ = e ωt cos(nx). Substituting it into (3.1), we obtain from the first equation b + n 2 η = 0, and, from the second one, Therefore, taking appropriate values of ω and a specific combination of b and η, the above expressions can be a solution to system (3.1).
Hence, throughout this section and in the rest of the paper, we will assume that b + n 2 η = 0 for all n ∈ N. In fact, we need to impose η = 0 (as we set in conditions (3.2)).
We transform our initial-boundary problem (3.1) into a more abstract problem in an appropriate Hilbert space with an adequate inner product.
Let us denote v =u and ψ =φ. We consider the Hilbert space (u, v, ϕ, ψ) and U * = (u * , v * , ϕ * , ψ * ) are two elements of H, we define its inner product as As usual, a superposed bar is used to denote the conjugate of a complex number. It is worth recalling that this product is equivalent to the usual product in the Hilbert space H.
Using D i to notate the i-derivative with respect to the spatial variable x, we can rewrite system (3.1) as follows:u With the above notation, our initial-boundary value problem can be written as where (u 0 , v 0 , ϕ 0 , ψ 0 ) are the initial conditions (2.2) and A is the following 4 × 4 matrix: In this matrix, I denotes the identity operator. The domain of the operator A, which will be denoted by D(A), is given by the elements We prove first the existence and uniqueness of solutions. To do so, we have to show that the operator is dissipative and that 0 belongs to the resolvent of A.
Using the inner product defined above and taking into account the assumed boundary conditions, we obtain which proves that operator A is dissipative. Let us remark that the boundary conditions play an important role to obtain this result (and similar results for other matrix operators that we use later). It can be proved that the general solutions to system (3.3) are given by the semigroup of contractions generated by the operator A. Proof. For any F = (f 1 , f 2 , f 3 , f 4 ) ∈ H we will find U ∈ H such that AU = F. Writing this condition term by term, we get: We will solve the above system using the expressions of f i as Fourier series. That means that we write f i = f i n sin nx for i = 1, 2 and f j = f j n cos nx for j = 3, 4, with On the other hand, the solutions we are looking for can be written also as Fourier series with unknown coefficients: If we substitute these expressions into the system, we get straightforwardly that v n = f 1 n and ψ n = f 3 n . Moreover, for each n a new system of equations is obtained: n . The solution of this system is given by where which is strictly positive for all n due to the hypotheses over the constitutive coefficients. The only doubt can be found in the term that goes with n 4 . However, we notice that Hence, it is not difficult to see that n 4 u 2 n < ∞ and n 4 ϕ 2 n < ∞. It remains to show that but easy calculations give that, for each n, where p 6 (n) and p 6 (n) are polynomials of degree six on n and a 4 is a 4 /n 2 . Therefore, it can be seen that both linear combinations belong to L 2 . Finally, taking into account the solutions obtained for u n , v n , ϕ n and ψ n , it can be shown that The fact that the operator A is dissipative, jointly with the above lemma and the Lumer-Phillips theorem, proves the existence and uniqueness of solutions. We write this result in the following theorem.
We prove now that the solutions to system (3.1) do not decay exponentially. To do so, we prove that there exists a solution of the form such that (ω) > − for all positive small enough. This fact implies that we can find a solution ω as near to the imaginary axis as we desire and, hence, it is impossible to have uniform exponential decay on the solutions to the problem determined by (3.1), with conditions (2.1) and (2.2).
Imposing that u and ϕ are of the form u = A 1 e ωt sin(nx) and ϕ = A 2 e ωt cos(nx), the following homogeneous system on the unknowns A 1 and A 2 is obtained: This linear system will have nontrivial solutions if, and only if, the determinant of the coefficients matrix is null. We denote by p(x) the fourth-degree polynomial obtained from the determinant of the coefficients matrix once ω is replaced by x, and by a i its coefficients for i = 0, 1, 2, 3, 4. These coefficients depend on the parameters of system (3.1) and on n. To be precise: To prove that there are roots of p(x) as near to the complex axis as desired is equivalent to show that, for any > 0, there are roots of p(x) located on the right-hand side of the vertical line (X) = − . If we make a translation, this fact is equivalent to show that polynomial p(x − ) has a root with positive real part.
We use the Routh-Hurwitz theorem (see Dieudonné [6]), which states that, if a 0 > 0, then all the roots of polynomial have negative real part if, and only if, a 4 and all the leading diagonal minors of matrix The coefficients of p(x − ) are given by The third leading minor of the Routh-Hurwitz matrix is a sixteenth-degree polynomial on n whose main coefficient is negative for n large enough. To be precise, if we denote by L i the leading minors of this matrix, where p 14 (n) is a fourteenth-degree polynomial on n. Therefore, it is clear that, for n large enough, L 3 will be negative provided that is sufficiently small. The above argument proves the slow decay of the solutions to system (3.1). We can be more specific and prove that, in fact, the solutions decay polynomially. We use the characterization given by Borichev and Tomilov [2], which we recall in the following theorem.
. Then, given α > 0, the following conditions are equivalent: Unfortunately, we cannot prove these conditions in a straightforward way. We have to decompose the Hilbert space H as the direct sum of two subspaces: Notice that K N is invariant under the semigroup. That means that the solutions starting at K N always belong to K N . A solution to system (3.1), U (t), can also be decomposed as the sum of two elements: As U 1 (t) belongs to a finite-dimensional subspace, if all the eigenvalues have negative real part, the exponential decay of U 1 (t) is guaranteed and, therefore, the polynomial decay is also satisfied. Proof. Imposing as above that u and ϕ are of the form u = A 1 e ωt sin(nx) and ϕ = A 2 e ωt cos(nx), we obtain the same linear homogeneous system of equations and the same polynomial p(x) with coefficients a i for i = 0, 1, 2, 3, 4. Straightforward computations show that the leading minors (denoted by M i below) corresponding to polynomial p(x) are all positive: By hypothesis, b + ηn 2 = 0 for all n ∈ N, and a and k 1 are positive.
We can show the two conditions of the Borichev and Tomilov characterization for the part of the solutions in K, the non-finite-dimensional subspace. We prove first that the imaginary axis is contained in the resolvent of A. Lemma 3.6. Let A be the matrix operator defined before. Therefore, iR ⊂ (A).
Proof. The proof has three steps. The first two ones refer to the operator A and are quite standard (for details see, for instance, [16], page 25). We concentrate in the third one, which is specific for each case.
Here, it reads as follows: suppose that the statement of this lemma is not true. Therefore, there exist a sequence of real numbers λ n with λ n → ∈ R, |λ n | < | | and a sequence of vectors U n = (u n , v n , ϕ n , ψ n ) in D(A) and with unit norm such that (iλ n I − A)U n → 0.
We prove now the second condition of Borichev and Tomilov's characterization taking α = 2.
Lemma 3.7. Let A be the above matrix operator. Then, Suppose that the statement of the lemma is not true. Then, there exist a sequence of real numbers, λ n , with |λ n | → ∞ and a sequence of unit norm vectors in the domain of A, U n = (u n , v n , ϕ n , ψ n ), such that λ 2 n (iλ n I − A)U n → 0. Writing this condition term by term, we get Selecting the real part of the product λ 2 n (iλ n I − A)U n , U n and taking into account (3.4), we get λ n D 2 v n → 0. Hence, it will also be λ n D 2 u n → 0.
We repeat the argument we did in the Proof of Lemma 4.5. First, we multiply convergence (3.11) by Dϕ n and notice that we can remove the λ 2 n because the expression inside the parentheses clearly tends to zero. We obtain again (3.9). In this case, it follows that Du n , D 4 ϕ n ≈ iJλ n k 2 Du n , ψ n = iJ k 2 λ n Du n , ψ n → 0, and so, we find that Dv n , D 4 ϕ n → 0. This argument shows that U n cannot be of unit norm, which finishes the proof of this lemma.
From the above results, we can state the following theorem.

Viscoelasticity
We introduce now a dissipative mechanism in the elasticity that, intuitively, is weaker than the previous one because we take only the first derivative of the displacement velocity with respect to x. Let us assume that τ = au x + bϕ + βϕ xx + a * u x , with a * > 0, while the other constitutive equations remain unaltered.
Substituting the constitutive equations into the evolution equations, we obtain a new system of field equations: Conditions (3.2) are assumed for the system coefficients. The boundary and initial conditions (2.1) and (2.2) are also assumed for the above system.
We want to highlight the fact that, in this system, the dissipation term is given by a second-order derivative with respect to x, while, in the previous section, it was given by a fourth-order derivative. Nevertheless, we will prove that the solutions to this system decay exponentially.
We still assume that b + n 2 η = 0 for all n ∈ N. The same Hilbert space is considered, with the same inner product. We can rewrite system (4.1) as follows: . We denote by B the matrix operator corresponding to this system: Therefore, system (4.1) can be written as where (u 0 , v 0 , ϕ 0 , ψ 0 ) are the initial conditions (2.2). The domain of B is given by the elements We prove first the existence and uniqueness of solutions. The operator B is dissipative and a direct calculation gives Proof. We proceed as in the proof of Lemma 3.2. For any F = (f 1 , f 2 , f 3 , f 4 ) ∈ H, we will find U ∈ H such that BU = F, or equivalently, we will find a solution to the system: v = f 1 , We write f i = f i n sin nx for i = 1, 2 and f j = f j n cos nx for j = 3, 4, with n 4 (f i n ) 2 < ∞ for i = 1, 3 and (f i n ) 2 < ∞ for i = 2, 4. We make an abuse of the notation, and we write again u = u n sin nx, v = v n sin nx, ϕ = ϕ n cos nx, ψ = ψ n cos nx for the solutions. It is clear that v n = f 1 n and ψ n = f 3 n . Simplifying, the following system of equations is obtained for each n: −an 2 u n − bnϕ n − k 1 n 4 u n − ηn 3 ϕ n = ρf 2 n + a * f 1 n n 2 , −ηn 3 u n − bnu n − δn 2 ϕ n − k 2 n 4 ϕ n − ξϕ n = Jf 4 n . The solution of this system is given by where a 4 is the independent term of the polynomial we have seen in the proof of Lemma 3.2. Hence, it is not difficult to see that n 4 u 2 n < ∞ and n 4 ϕ 2 n < ∞. It remains to show that k 1 D 4 u + ηD 3 ϕ ∈ L 2 and ηD 3 u − k 2 D 4 ϕ ∈ L 2 . Easy calculations give that, for each n, k 1 n 4 u n + ηn 3 ϕ n = −a * k 1 k 2 f 1 n n 8 + p 6 (n) a 4 and −ηn 3 u n − k 2 n 4 ϕ n = Jk 1 k 2 f 4 n n 6 + p 5 (n) a 4 , where p 6 (n) and p 5 (n) are polynomials of degree six and five on n, respectively, and a 4 is a 4 /n 2 . Therefore, it can be seen that both linear combinations belong to L 2 . Finally, taking into account the solutions obtained for u n , v n , ϕ n and ψ n , it can be shown that where K is a constant independent of U .
Therefore, the existence and uniqueness of solutions is clear. We write this result in the following theorem. To prove the exponential decay of the solutions we need to split again H in two subspaces and to decompose a solution to system (4.1) as the sum of two elements, U (t) = U 1 (t) + U 2 (t), as we did in Sect. 3.
Again, if U 1 (t) belongs to a finite-dimensional subspace and all the eigenvalues have negative real part, the exponential decay of U 1 (t) is guaranteed. Proof. Imposing u = A 1 e ωt sin(nx) and ϕ = A 2 e ωt cos(nx), the following homogeneous system on the unknowns A 1 and A 2 is obtained: k 1 n 4 + (a + ωa * )n 2 + ρω 2 ηn 3 + bn ηn 3 + bn k 2 n 4 + δn 2 + Jω 2 + ξ We make again an abuse of the notation and we denote by a i , for i = 0, 1, 2, 3, 4, the coefficients of the fourth-degree polynomial obtained from the determinant of the coefficients matrix once ω is replaced by x: a 0 = ρJ, a 1 = Ja * n 2 , a 2 = (Jk 1 + k 2 ρ) n 4 + (aJ + δρ)n 2 + ξρ, a 3 = a * k 2 n 6 + a * δn 4 + a * ξn 2 , As we obtained before, it is clear that a i > 0 for i = 0, 1, 2, 3, 4. A direct calculation shows that all the leading minors of the Routh-Hurwitz matrix are positive. To be precise: Proposition 4.3 shows the exponential decay of U 1 (t). We study now U 2 (t). To prove the exponential decay, we use the characterization given by Huang [11] or Prüss [29]. We recall it below. We split these conditions in two separate lemmata.

Lemma 4.5.
Let B be the matrix operator defined above. Therefore, iR ⊂ (B).
Proof. We suppose then that there exist a sequence of real numbers λ n with λ n → , |λ n | < | | and a sequence of vectors U n = (u n , v n , ϕ n , ψ n ) in D(B) and, with unit norm, such that (iλ n I − B)U n → 0. If we write the above expression term by term, we obtain the following conditions: Selecting the real part of the product (iλ n I − B)U n , U n and taking into account (4.3), it is clear that Dv n → 0 and, hence, it follows that λ n Du n → 0. Let us multiply expression (4.4) by u n : iλ n ρ v n , u n − a D 2 u n , u n − b Dϕ n , u n + η D 3 ϕ n , u n + k 1 D 4 u n , u n − a * D 2 v n , u n → 0. Now, we remove the first term (that tends to zero), integrate by parts on the other terms and it follows that a Du n , Du n + b ϕ n , Du n − η D 2 ϕ n , Du n + k 1 D 2 u n , D 2 u n + a * Dv n , Du n → 0, which yields D 2 u n → 0. Notice that D 2 ϕ n is bounded because ϕ ∈ H 2 . Notice also that D 2 v n → 0. We multiply again expression (4.4) but now by Dϕ n . We remove the terms with D 2 u n and D 2 v n , and we obtain: From (4.4), iλ n ϕ n ≈ ψ n . Hence, integrating by parts, we get iλ n ρ v n , Dϕ n = −ρ v n , iλ n Dϕ n ≈ ρ Dv n , ψ n → 0.
Since the second and third terms tend to zero, we now see that − Jψ n , iλ n ϕ n ≈ −J ψ n , ψ n → 0, which finishes the proof because this shows that vector U n cannot be of unit norm. Proof. Notice that throughout the proof of the previous lemma we only make use of the fact that λ n does not tend to zero, but it does not depend on λ n tending to a finite number or to infinity.
As a consequence of the above lemmata, we have the following result.

Dissipation in the porosity
In this section, we want to introduce different dissipation mechanisms in the porosity component. In fact, we will study three of them. We will develop an analysis quite similar to the ones used in Sects. 3 and 4. Some parts of the analysis are, mutatis mutandis, equal to the previous ones, and, for this reason, we only write the main results. In the following subsections, we denote by A the matrix operator corresponding to each system. We think that it will not cause any misunderstanding. Moreover, in each subsection, to obtain the system of field equations, we write only the changes we need to impose in one (or more) of the constitutive equations to introduce the dissipation. The other constitutive equations remain unaltered.

Hyperviscoporosity
We assume that with k * 2 > 0. The system of field equations is given by (5.1) Proof. The proof follows the same scheme as we did in Sect. 3. It is worth noting that, if A denotes the matrix operator corresponding to system (5.1), therefore Using the Routh-Hurwitz theorem for the polynomial p(x − ) corresponding to system (5.1), we obtain that the third leading minor is which gives the slow energy decay. Using the Borichev and Tomilov characterization, we obtain the polynomial rate of decay.

Viscoporosity
We change now the dissipation mechanism. In the constitutive equations, we consider xx , which give rise to the following system of field equations: where δ * = α * − 2d * > 0. Proof. The proof follows the same scheme used in Section 4. In this case, we apply Huang's characterization. If A is the matrix operator obtained from system (5.2), then which is used in the proof of the equivalent here to Lemmata 4.5 and 4.6.
It is worth noting that systems (5.1) and (5.2) are quite "symmetric" to systems (3.1) and (4.1), respectively, and hence it is not surprising at all that the solutions behave analogously.

Weak viscoporosity
Finally, we consider an even weaker dissipation mechanism. We take Therefore, the following system of field equations is obtained: where ξ * > 0. Proof. The existence and uniqueness part can be proved as in the previous sections. The slow decay can also be showed following the same methods: The third leading minor of the Routh-Hurwitz technique is which is negative for n large enough when Jk 1 = ρk 2 provided that is sufficiently small. We concentrate now in the exponential decay, which is the difficult part because it is quite different from the ones we have done previously. We suppose that Jk 1 = ρk 2 .
If A denotes the matrix operator obtained from system (5.3), a direct calculation gives We prove both conditions of Huang's characterization. First of all, as in Lemma 3.6, we consider a sequence of unit norm vectors in the domain of the operator and we write term by term the convergences: iλ n Jψ n − ηD 3 u n + bDu n − δD 2 ϕ n + k 2 D 4 ϕ n + ξϕ n + ξ * ψ n → 0, in L 2 . (5.7) From (5), we get ψ n → 0 and, hence, λ n ϕ n → 0 in L 2 . We denote by m n the function such that D 2 m n = ϕ n and Dm n is zero at the boundary. Notice that, in particular, there exists a real number C such that m n ≤ C ϕ n and λ n ϕ n → 0 implies λ n m n → 0 and, moreover, λ n Dm n → 0.
We remove from (5.7) the terms that tend to zero, and we multiply the remaining part by m n : iJλ n ψ n , m n − η D 3 u n , m n + b Du n , m n − δ D 2 ϕ n , m n + k 2 D 4 ϕ n , m n → 0.
Integrating by parts and taking into account that Du n is bounded, the above expression becomes iJψ n , λ n m n + η D 2 u n , Dm n − δ ϕ n , ϕ n − k 2 Dϕ n , Dϕ n → 0.
The three first terms tend to zero and, therefore, it is clear that Dϕ n → 0. We multiply (5.7) by ϕ n , and we obtain: Integrating again by parts, using that Dϕ n → 0 and that D 2 u n is bounded we get D 2 ϕ n → 0. We remove from (5.5) and (5.7) the terms which tend to zero, and we multiply the remaining parts by Dϕ n and by Du n , respectively. We get: iλ n ρv n , Dϕ n − a D 2 u n , Dϕ n + η D 3 ϕ n , Dϕ n + k 1 D 4 u n , Dϕ n → 0 and iλ n Jψ n , Du n − η D 3 u n , Du n + b Du n , Du n − δ D 2 ϕ n , Du n + k 2 D 4 ϕ n , Du n → 0.
Using the previous results, the first expression reduces to On the other hand, from convergence (5.6) we find that ψ n ∼ iλ n ϕ n , and convergence (5.9) becomes Finally, applying the hypothesis Jk 1 = ρk 2 , we obtain η D 2 u n , D 2 u n + b Du n , Du n → 0, which implies that Du n → 0 and D 2 u n → 0 for n large enough. Notice that we do not distinguish between λ n being bounded or not because this does not matter in the proof. The only relevant point is that λ n does not tend to zero.

Numerical behavior
In this section, we study a fully discrete approximation of a variational version of the above mechanical problems. So, we introduce its variational formulation. Let Y = L 2 (0, π), and denote by (·, ·) the scalar product in this space, with corresponding norm · .
We replace boundary conditions (2.1) by the following ones: Therefore, integrating by parts we derive the following variational formulation for the problems studied in the previous sections.
We note that we omit the analysis of the problems involving hyperviscoporosity and viscoporosity cases because they are similar to the hyperviscosity and viscosity ones, respectively. Now, we provide the fully discrete approximation of the previous weak problem. This is done in two steps. First, 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 , and so, to approximate the variational space represents the space of polynomials of degree less or equal to three in the subinterval [a i , a i+1 ]; i.e., the finite element space V h is made of C 1 and piecewise cubic functions. Here, h > 0 denotes the spatial discretization parameter. Furthermore, let the discrete initial conditions u h 0 , v h 0 , φ h 0 and ψ h 0 be defined as where P h is the classical finite element interpolation operator over V h (see [3]).
Therefore, using the well-known Newmark-β scheme, the fully discrete approximations of the above variational problem are the following.
Find the discrete displacement u hk = {u hk n } N n=0 ⊂ V h and the discrete porosity function where the discrete velocity, the discrete porosity speed, the discrete acceleration and the discrete porosity accelerationu hk n ,φ hk n ,ü hk n andφ hk n are now recovered from the relations: We note that the first time iteration is done using the implicit Euler scheme, and so, the accelerations at time t 1 are obtained asü hk It is straightforward to obtain that this fully discrete problem has a unique solution applying the well-known Lax-Milgram lemma and the required assumptions on the constitutive parameters.
In all the numerical simulations described below, we have used the following data: T = 7000, ρ = 1, k 1 = 1, a = 1, b = 1, η = 1, J = 1, k 2 = 1, δ = 2, ξ = 2, and the initial conditions, for all x ∈ (0, 1): We note that, for the sake of simplicity in the numerical implementation, we have assumed that the length of the beam is 1 (instead of π). Moreover, we have chosen the discretization parameters h = 0.025 and k = 10 −3 and the Newmark−β coefficients α = 0.25 and β = 0.5. In the first example, we solve the discrete problem assuming that k * 1 = ξ * = 0 and varying parameter a * between 0.01 and 100 (which corresponds to the numerical resolution of system (4.1) with boundary conditions (6.1)). In Fig. 1, we plot the evolution in time of the discrete energy given by in both normal and semi-log scales. ZAMP As can be seen, the theoretical asymptotic exponential behavior of the energy can be clearly seen for all the coefficients, although we can also appreciate that, when parameter a * increases (higher than 50), the energy decay seems to reduce. A possible explanation for this finding could be the fact that the dissipation mechanism of the beam becomes too rigid, and so, the dissipation is strongly affected (see also the zoom part shown on the left-hand side).
Secondly, we consider the dependence of the solution with respect to parameter k * 1 assuming now that ξ * = a * = 0 (i.e., it corresponds to the numerical resolution of system (3.1) with boundary conditions (6.1)). Therefore, the evolution in time of the discrete energy given above is shown in Fig. 2 for some values of parameter k * 1 (k * 1 = 0.1, 1, 10). As can be clearly seen, an asymptotic exponential behavior is again observed for the discrete energy. Although we have proved theoretically that it should decay as t −1/2 , we note that it cannot be found in the numerical simulations because, in this case, the variational space has a finite dimension and so, all the eigenvalues of the corresponding operator (the eigenvalues of the matrix system) have real part. Therefore, the energy decay is always exponential. Comparison between the solution obtained with the dissipation mechanism a * = 1 and ξ * = k * 1 = 0 (second order) and the dissipation mechanism k * 1 = 1 and ξ * = a * = 0 (fourth order) One interesting issue in this experiment is the dependence on parameter k * 1 because, when it increases, the dissipation mechanism becomes rigid (as in the previous example) and the energy decay is slower; however, it remains to be understood what happens when this parameter becomes smaller. In Fig. 3, we plot the energy decay for a large number of solutions with parameter k * 1 varying between 0.01 and 0.3. As can be seen on the left-hand side, the energy curve decreases when the parameter k * 1 increases until value k * 1 = 0.04, and then, it starts to increase again. In order to analyze easily this behavior, on the right-hand side we plot the values of the energy at time t = 50, and we can clearly appreciate how this minimum is achieved. Now, the aim is to compare both dissipation mechanisms with the same values for the equivalent constitutive parameters (so, we have used values a * = 1 and k * 1 = 1 for each case). The comparison of the energy decay is shown in Fig. 4. As can be seen, the energy decay is clearly faster when the second-order dissipation mechanism is considered (case a * = 1 and ξ * = k * 1 = 0). However, the asymptotic energy decay for the fourth-order case (k * 1 = 1 and ξ * = a * = 0) is also exponential. Finally, we analyze the dependence of the solution with respect to parameter ξ * assuming now that k * 1 = a * = 0 (i.e., it corresponds to the numerical resolution of system (5.3) with boundary conditions (6.1)). Thus, the evolution in time of the discrete energy given above is shown in Fig. 5 for some values of parameter ξ * (ξ * = 0.01, 0.1, 1, 10, 100).
As can be seen, an asymptotic exponential behavior is again found for the discrete energy. Although in Theorem 5.3 it is shown that this behavior depends on the condition Jk 1 = ρk 2 , as in the previous example, we can conclude that this is not required for a finite-dimensional setting. Therefore, the energy decay is always exponential as in the previous cases involving the mechanical dissipation. Now, we focus on the dependence on parameter ξ * . As in the case of a fourth-order dissipation mechanism, when it increases the mechanism becomes rigid and the energy decay is slower; however, it remains to be explained what happens when this parameter becomes high. In Fig. 6, we plot the energy decay for a large number of solutions with parameter ξ * varying between 10 and 6000. On the left-hand side, we can appreciate that the smaller curve is found for the value ξ * = 100. Moreover, on the right-hand side this discrete energy is shown at time t = 10 for these values where the minimum is clearly seen.

Conclusions
We have analyzed the time decay for the solutions to the system of partial differential equations that models the behavior of porous elastic materials when fourth-order derivatives with respect to the spatial variable are considered in both components, the displacement and the porosity, and one dissipation mechanism is present in the system. Let us briefly summarize what we have found for each case: (1) Hyperviscoelasticity: polynomial (slow) decay, controlled by t −1/2 .
(3) Hyperviscoporosity: polynomial (slow) decay, controlled by t −1/2 . (4) Viscoporosity: exponential decay. (5) Weak viscoporosity: slow decay in the generic case and exponential decay in a specific situation. These behaviors differ from the ones known for the classical theory, where, generically, two dissipation mechanisms are needed (one at the macrostructure level and another one at the microstructure) to guarantee the exponential decay. They differ also from the ones obtained for the strain gradient situation with only high-order derivatives in the elastic component of the structure.
Finally, we have performed some numerical simulations to analyze this theoretical behavior. Therefore, using the finite element method and the Newmark-β scheme we have implemented a numerical algorithm in MATLAB for the solution of the hyperviscoelastic, the viscoelastic and the weak viscoporosity cases (the remaining two cases are similar to the previous ones). We have found that, for every problem, the discrete energy decay is always exponential, but it is significantly faster for the viscoelastic case (second order). We have also seen that the energy decay has a minimum value depending on the constitutive parameter which is different for each case.