On the mixtures of MGT viscoelastic solids

: In this paper, we study, from both analytical and numerical points of view, a problem involving a mixture of two viscoelastic solids. An existence and uniqueness result is proved using the theory of linear semigroups. Exponential decay is shown for the one-dimensional case. Then, fully discrete approximations are introduced using the finite element method and the implicit Euler scheme. Some a priori error estimates are obtained and the linear convergence is derived under suitable regularity conditions. Finally, one- and two-dimensional numerical simulations are presented to demonstrate the convergence, the discrete energy decay and the behavior of the solution.


Introduction
It is well-known that to describe several kinds of materials we cannot use the usual theories of elasticity or viscoelasticity. In fact, in the second part of the past century, several generalized models were considered to describe the behavior of materials. If different components interact within a material, it is suitable to use the so-called theory of mixtures. The easiest example of mixtures could be metallic alloys, but many other applications can be found in the theory of composites [1]. We cite classical references where these kind of materials are described [2][3][4][5][6][7][8][9][10]. Another application could be the interaction in saturated soils [11][12][13][14][15] or the water-heat coupling process during freezing in a saline soil [16]. We also recall the first proposition with Lagrangian description, and where the independent variables were the gradient of each displacement and the relative displacement, was suggested by Bedford and Stern [5,6].
The theory of mixtures is currently under deep study in the scientific community. If we restrict our attention to the case of two interacting continua, the displacement of each component's particles depends on the material point and the time, and we will assume that the particles under consideration are at the same position at initial time.
We can find in the literature many theories about Kelvin-Voigt mixtures. These theories are known, but they address the instantaneous propagation of the mechanical waves, which is not compatible with the causality principle. Therefore, it should be useful to provide an alternative theory for the behavior of viscoelastic mixtures. When we introduce delay parameters, the problem becomes hyperbolic and the mechanical waves propagate at a finite speed, which is more consistent from the physical point of view. This process is similar to the one used in the heat conduction theory between the classical Fourier law and the theory of Maxwell-Cattaneo.
This mechanism has been observed previously in the case of a single viscoelastic material and it has been successfully noted that we should substitute the usual parabolic equation by a hyperbolic equation of the Moore-Gibson-Thompson type [17][18][19][20][21][22][23][24]. In fact, in a recent contribution [25], the authors proposed a system of equations describing a mixture of a viscous solid (of Moore-Gibson-Thompson type) and an elastic solid.
In the present paper, we want to continue the work of [25], studying the problem determined by the mixture of two viscoelastic solids. Furthermore, it is known that a material described by the Moore-Gibson-Thompson equation requires a "time delay parameter". Here, we will consider the general case where the "time delay parameter" of each constituent can be different. This is important because it will propose some new relevant difficulties in the study of the mathematical system.
In the next section, we propose the basic equations and the assumptions we are going to focus on. In Section 3 we give a functional description of the problem and we obtain a result about the existence and uniqueness of a quasi-contractive semigroup. In order to simplify the mathematical analysis, we restrict our attention to the one-dimensional case in Section 4 and we provide an exponential decay result (under suitable conditions). Later, in Section 5 we study numerically a variational formulation of this problem by using the classical finite element method to approximate the spatial variable and the implicit Euler scheme to discretize the time derivatives. Some a priori error estimates are proved, from which the linear convergence is obtained under suitable regularity conditions on the continuous solution. Finally, some one-and two-dimensional numerical simulations are presented in Section 6 to demonstrate the convergence of the approximation, the discrete energy decay and the behavior of the solution.

Basic equations
In this section, we propose the basic equations that we want to study in this article. Let us denote by Ω the three-dimensional domain occupied by the mixture, and let x and t be the spatial and time variables, respectively. As usual, we will assume that indexes i, j, k, l vary between 1 and 3, and we will use the repeated index summation.

The model
We first recall those corresponding to a mixture of materials. The evolution equations are: where a dot represents the first-order time derivatives, two dots the second-order and three dots the third-order, a subscript j after a comma denotes the spatial derivative with respect to variable x j , ρ 1 , ρ 2 are the mass density of every constituent, σ i j , τ i j are the partial stress tensors, p i is the internal body force and u i and w i are the displacements of each constituent. In the case that we consider a mixture of Kelvin-Voigt viscoelastic solids, the constitutive equations are (usually) assumed: In this work, we assume the following symmetries on these tensors: If we introduce the constitutive equations into the evolution equations, we obtain a system of equations such that it allows the propagation of mechanical waves instantaneously. This fact violates the causality principle and it is very similar to what happens in the case of the heat conduction of Fourier (or type III Green-Naghdi) theory. A way to overcome this difficulty is the introduction of delay parameters as the Maxwell-Cattaneo heat equation (or the Moore-Gibson-Thompson one). In this work, we want to study what happens when we follow a similar proposition in our case and we will assume that That is, we introduce two relaxation parameters τ 1 and τ 2 , which can be different for the sake of generality, assumed to satisfy the condition τ 1 ≥ τ 2 . The case τ 2 ≥ τ 1 would follow a similar procedure.
If we introduce the new constitutive equations (2.2) 3 and (2.3) into the evolution equations, taking into account that we obtain the following system: In order to work with this system, we need to impose initial and boundary conditions. We will assume that, for a.e. x ∈ Ω, (2.5) We will consider homogeneous Dirichlet boundary conditions. Therefore, we impose that u i (x, t) = w i (x, t) = 0 for a.e. x ∈ ∂Ω and t ≥ 0. (2.6) To study problems (2.4)-(2.6), it will be useful to manipulate the second equation of the system (2.4). We note that ρ 2 (τ 2 ... and so, we can write the second equation of the system (2.4) as (2.7) In several parts of this paper, we restrict our attention to the one-dimensional homogeneous case in order to clarify better the analysis. The system of equations is written in this case as where the domain is now Ω = (0, ℓ), with ℓ > 0.

The assumptions
To understand in a clear way the assumptions that we need to impose in the constitutive terms, we will need to consider the energy equation. We have that In view of these relations it will be natural to assume i) ρ 1 (x) ≥ ρ 0 1 > 0 and ρ 2 (x) ≥ ρ 0 2 > 0.
ii) There exists a positive constant C 1 such that for every tensors ξ i j and η kl . iii) There exists a positive constant C 2 such that for every tensors ξ i j and η i j . iv) There exists a positive constant C 3 such that These assumptions are natural in the context of the MGT problems. The meaning of condition i) is clear. Assumptions ii) and iii) are needed to guarantee the positivity of the internal energy, and assumption iv) is a usual condition in the study of mixtures.
It is worth noting that, in the case of one-dimensional homogeneous materials, these conditions can be written as and A > 0.
Remark 1. When B = B * = 0 we can propose an alternative approach. We assume that, in this case, τ 2 > τ 1 . We note that Therefore, we have Then, we can obtain the equality In the previous definitions, we have used the notation:

Existence of solutions
The aim of this section is to obtain an existence and uniqueness result for the problem determined by the system (2.4) and conditions (2.5) and (2.6) under the assumptions i)-iv).
We first consider the Hilbert space where W 1,2 0 (Ω) = [W 1,2 0 (Ω)] 3 and L 2 (Ω) = [L 2 (Ω)] 3 , which is equipped by an inner product defined as , and the bar over the elements in the Hilbert space denotes the complex conjugate. It is easy to show that this product is equivalent to the usual one in the Hilbert space H.
We can write our problem in the following form: where the matrix differential operator A is defined as and its terms are given by We note that the domain of the operator A is the subspace of U ∈ H such that AU ∈ H. It is clear that this is a dense subspace. Lemma 1. There exists a positive constant K 1 such that Proof. In view of the energy equality we see that where K is a positive constant, and, because of the definition of the inner product, we obtain the required inequality.
There exists a positive constant K 2 such that K 2 I − A is exhaustive.
Proof. Let us consider F = ( f 1 , . . . , f 6 ) ∈ H. We have to see that the system admits a solution in the domain of the operator. We note that we can write where the elements F 1 are F 2 are given by and they are obtained after a direct substitution of one equation into the other ones. As a consequence of the previous regularities, we find that (F 1 , On the other side, if we define the bilinear form it is clear that, for K 2 large enough, B is a coercive and bounded bilinear form on W 1,2 0 (Ω) × W 1,2 0 (Ω). In view of the Lax-Milgram theorem, we obtain the existence of (u, w) satisfying the system (3.1). Then, we also obtain the existence of (v, z, y, r) and the lemma is proved.
Remark 2. We note that Lemma 2 also holds when K 2 = 0. Therefore, zero belongs to the resolvent of the operator.
If we use now the Lumer-Phillips corollary to the Hille-Yosida theorem we find the following. Theorem 1. The operator A generates a quasi-contractive semigroup.
Therefore, we also have. Theorem 2. The problems (2.4)-(2.6) admits a unique solution. In fact, for every U 0 ∈ Dom(A) there exists a unique solution with the following regularity: Remark 3. It is clear that we could also prove the existence of a quasi-contractive semigroup determined by the scalar product associated to the energy E 2 (t).

The one-dimensional case: energy decay
In this section, we are interested on the decay of solutions to our problem. In order to make the calculations easier, we restrict our attention to the homogeneous one-dimensional case taking the domain Ω = (0, ℓ), ℓ > 0, and we deal with system (2.8). We also consider the one-dimensional version of the initial condition (2.5) and the boundary condition (2.6). Therefore, the energy equation is satisfied, where If we impose now that assumptions i * ) and ii * ) hold, then the matrices associated with coefficients A, B and C, and A, B and C, are positive definite and so, the resulting semigroup associated to operator A is quasi-contractive. In fact, we could prove that where K is a positive constant, and that the zero belongs to the resolvent of the operator.
If we denote then we can consider the bilinear form defined by the matrix: where we recall that τ 1 − τ 2 ≥ 0. We note that after a simple calculation we find that It is worth noting that, if matrix M is positive definite, then Re⟨AU, U⟩ ≤ 0 for every U ∈ Dom (A) and, therefore, the semigroup associated to operator A is contractive. Remark 4. M is positive definite if and only the following conditions hold: We see that whenever τ 1 − τ 2 is small enough in comparison with mλ * 1 , and assuming that a and τ 1 are not very large, we can guarantee that M is positive definite. Of course, to give the explicit set where these conditions hold is a very difficult task, since we are involving nonlinear algebraic equations with several parameters. It seems much easier to check, for each particular case, if the conditions are satisfied. However, we note that Finally, we have the following. Theorem 3. If the matrix (4.1) is positive definite, the solutions to the one-dimensional problem decay in an exponential way. That is, there exist two positive constants N and ω such that Proof. To prove this result we will use the well-known characterization of the exponentially stable semigroups. We know that it is sufficient to prove that the imaginary axis is contained in the resolvent of the operator and that the asymptotic condition holds (see [26]). A well-known argument (see, again, [26]) shows that, if we assume the existence of λ ∈ R such that iλ belongs to the spectrum, then there exist a sequence of real numbers λ n → λ 0, and a sequence of unit norm vectors U n = (u n , v n , z n , w n , y n , r n ) at the domain of the operator A such that ∥(iλ n I − A)U n ∥ → 0.
In view of the dissipation inequality we see that v n , y n → 0 in W 1,2 (0, ℓ) and r n → 0 in L 2 (0, ℓ). Then, we also obtain that λ n w n , λ n u n → 0 in W 1,2 (0, ℓ). Now, if we multiply the third convergence by z n we also obtain that z n → 0 in L 2 (0, ℓ). We have found a contradiction to suppose that the thesis was not correct. Then, we have proved that the imaginary axis is contained in the resolvent of the operator. Now, we want to prove that the asymptotic condition (4.2) holds. Let us assume that this condition is not fulfilled. Then, there exist a sequence of real numbers λ n → ∞ and a sequence of unit norm vectors at the domain of the operator A such that the condition holds. From this point, we can follow the same argument we used before because the only condition is λ 0. In short, the theorem has been proved.

Remark 5.
It is suitable to consider the case τ 1 = τ 2 . In this case, we cannot assume that the matrix M is positive definite. Nevertheless, we can obtain the same conclusion. Although we cannot obtain that r n → 0 from the dissipation inequality, we can see that y n , v n → 0 in W 1,2 0 (0, ℓ) and then, λ n u n , λ n w n → 0 in W 1,2 0 (0, ℓ). From the convergences third and sixth we can obtain that z n , r n → 0 in L 2 (0, ℓ) after multiplication by v n and y n , respectively. It is worth noting that this argument can be used to prove that the imaginary axis is contained in the resolvent and that the asymptotic condition (4.2) holds. Therefore, the solutions to the problem determined by the system (2.8), when τ 1 = τ 2 , and the one-dimensional version of the initial condition (2.5) and the boundary condition (2.6) decay in an exponential way whenever conditions i * ) and ii * ) hold.

Fully discrete approximations: a priori error estimates
In this section, we introduce a fully discrete approximation of the multidimensional problem studied in the Sections 2 and 3, and we perform an a priori error analysis. However, we consider now the deformation of the body Ω over a finite time interval [0, T ], for a given final time T > 0.

Fully discrete approximations
First, we provide the variational formulation of problems (2.4)-(2.6). Let us introduce the velocities of the first and second constituents, denoted by v =u and y =ẇ, and the accelerations of the first and second constituents, denoted by z =ü =v and r =ẅ =ẏ, respectively. Moreover, let V be the variational space given by V = [H 1 0 (Ω)] 3 . In order to simplify the writing, and for the sake of clarity, we define the following operators, for By using boundary condition (2.6), we derive the variational formulation of problem (2.4), with the modified Eq (2.7), initial condition (2.5) and boundary condition (2.6).
Find the acceleration of the first constituent z : [0, T ] → V and the acceleration of the second constituent r : [0, T ] → V such that z(0) = z 0 = (z 0 i ), r(0) = r 0 = (r 0 i ) and, for a.e. t ∈ [0, T ] and for all ξ, η ∈ V, where v 0 = (v 0 i ), y 0 = (y 0 i ), u 0 = (u 0 i ) and w 0 = (w 0 i ). Now, we introduce a fully discrete approximation of problems (5.1)-(5.4) and we provide an a priori error analysis. In order to obtain the spatial approximation, we assume that the domain Ω has a finite element triangulation T h assumed to be regular in the sense of [27] and so, we define the finite element space: In the above definition, P 1 (T r) is the space of polynomials of degree less or equal to one in the element T r, i.e., the finite element space V h is made of continuous and piecewise affine functions. Here, h > 0 represents the spatial discretization parameter. Moreover, we assume that the discrete initial conditions, denoted by u 0h , v 0h , z 0h , w 0h , y 0h and r 0h , are given by where P h is the classical finite element interpolation operator over V h (see, again, [27]). In order to discretize the time derivatives, we consider a partition of the time interval [0, T ], denoted by 0 = t 0 < t 1 < · · · < t N = T , and we use a uniform partition with step size k = T/N and nodes t n = n k for n = 0, 1, . . . , N. For a continuous function z(t), we use the notation z n = z(t n ) and, for the sequence {z n } N n=0 , we denote by δz n = (z n − z n−1 )/k its corresponding divided differences. Therefore, using the implicit Euler scheme, the fully discrete approximations of problems (5.1)-(5.4) are considered as follows.
Find the discrete acceleration of the first constituent z hk = {z hk n } N n=0 ⊂ V h and the discrete acceleration of the second constituent r hk = {r hk n } N n=0 ⊂ V h such that z hk 0 = z 0h , r hk 0 = r 0h and, for all n = 1, . . . , N and ξ h , η h ∈ V h , where the discrete velocities and the discrete displacements of the first and second constituents are recovered from the relations: v hk n = k n j=1 z hk j + v 0h , y hk n = k n j=1 r hk j + y 0h , (5.9) u hk n = k n j=1 v hk j + u 0h , w hk n = k n j=1 y hk j + w 0h . (5.10) Using assumptions i)-iv) and the well-known Lax-Milgram lemma, it is easy to prove that the above fully discrete problem admits a unique discrete solution.

A priori error estimates
Now, we aim to obtain some a priori error estimates on the numerical errors z n − z hk n and r n − r hk n . We have the following main a priori error estimates result.
where M > 0 is a positive constant which is be independent of the discretization parameters h and k, but depending on the continuous solution, δz j = (z j − z j−1 )/k, δv j = (v j − v j−1 )/k, δr j = (r j − r j−1 )/k and δy j = (y j − y j−1 )/k, and I 1 j and I 2 j are the integration errors given by Proof. First, we obtain some estimates for the acceleration of the first constituent z n . Then, we subtract variational equation (5.1) at time t = t n for a test function ξ = ξ h ∈ V h ⊂ V and discrete variational equation (5.7) to obtain, for all ξ h ∈ V h , and so, we have, for all ξ h ∈ V h , Taking into account that after easy algebraic manipulations, using several times Cauchy-Schwarz inequality and Cauchy's inequality it follows that, for all ξ h ∈ V h , (Ω)] 3 + A(u n − u hk n , δv n − δv hk n ) + B(w n − w hk n , δv n − δv hk n ) +A * (v n − v hk n , δv n − δv hk n ) + B * (y n − y hk n , δv n − δv hk n ) Proceeding in a similar form, we obtain the following error estimates for the acceleration of the second constituent, for all η h ∈ V h , ∥r n − r hk n ∥ 2 [L 2 (Ω)] 3 − ∥r n−1 − r hk n−1 ∥ 2 [L 2 (Ω)] 3 + B(u n − u hk n , δy n − δy hk n ) + C(w n − w hk n , δy n − δy hk n ) +B * (v n − v hk n , δy n − δy hk n ) + C * (y n − y hk n , δy n − δy hk n ) +∥v n − v hk n ∥ 2 V + ∥y n − y hk n ∥ 2 V + (δr n − δr hk n , r n − η h ) [L 2 (Ω)] 3 .
Multiplying the above estimates by k and summing up to n, we find that, for all Now, keeping in mind that (thanks again to property iii)), and, therefore, we have, for all Finally, taking into account that where I 1 j and I 2 j are the integration errors defined in (5.11) and (5.12), respectively, applying a discrete version of Gronwall's inequality (see, again, [28]) we conclude the desired error estimates.
The above main error estimates can be used to derive the convergence order of the approximations given by problems (5.7)-(5.10). Therefore, as an example, if we assume the following additional regularity: using the classical results on the approximation by finite elements and the regularities of the initial conditions (see, for instance, [27]), we have the following. ∥z n − z hk n ∥ [L 2 (Ω)] 3 + ∥v n − v hk n ∥ V + ∥u n − u hk n ∥ V + ∥r n − r hk n ∥ [L 2 (Ω)] 3 + ∥y n − y hk n ∥ V +∥w n − w hk n ∥ V ≤ M(h + k).

Numerical results
In this final section, we present some numerical results, obtained after an implementation of the fully discrete problem in MATLAB, involving examples in one and two dimensions. Our aim here is to show the numerical convergence of the approximations, the discrete energy decay and the behavior of the solution.
First, we note that, given the solution u hk n−1 , v hk n−1 , z hk n−1 , w hk n−1 , y hk n−1 and r hk n−1 at time t n−1 , the discrete accelerations at time t = t n , z hk n and r hk n , are obtained by solving the discrete linear system, for all τ 2 τ 1 (τ 1 r hk n + kr hk n , η h ) [L 2 (Ω)] 3 + k 3 C(r hk n , η h ) + k 2 C * (r hk n , η h ) + k 3 a(r hk n , η h ) + τ 1 k 2 a(r hk n , η h ) We note that the above numerical scheme was implemented on a 3.2 Ghz PC using MATLAB, and a typical one-dimensional run (h = k = 0.001) took about 0.23 seconds of CPU time.
6.1. First example: numerical convergence and energy decay in a 1D problem As an academical example, in order to show the accuracy of the approximations the following onedimensional problem is considered.
Regarding the boundary conditions, we assume that the boundary Γ is decomposed into three disjoint parts Γ D , Γ N and Γ F , being Γ D = {0} × [0, 1], Γ N = {1} × [0, 1] and Γ F = Γ − Γ D − Γ N . We note that the analysis shown in the previous section can be easily extended to this slightly more general case. Therefore, on the part Γ D we assume that the body is clamped and so, u(0, y, t) = w(0, y, t) = 0 for y ∈ [0, 1], t ∈ [0, 1], meanwhile on the Neumann part Γ N we assume that a surface force G is applied in the form: G(1, y, t) = (−y, 0) for t ∈ [0, 1], y ∈ [0, 1], and boundary Γ F is assumed to be traction-free.
Using the time discretization parameter k = 10 −2 , the norm of the displacements of the first constituent is plotted at final time T = 1 in Figure 3(left) over the deformed configuration. We can see that the body bends due to the application of the surface force and, due to the clamping conditions, the more displaced areas are located on the top right part, since the force is linearly increasing with the value of the Y-coordinate. Moreover, the norm of these displacements is also shown by using arrows in Figure 3(right). We can clearly appreciate the orientation due the boundary conditions. Finally, in Figure 4 we plot the displacements and the velocity of the second constituent at final time over the deformed configuration. We note that they are generated by the resulting deformation and the coupling with the displacement and velocity of the first constituent.