On the mixtures of MGT viscoelastic solids

1.   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.

2.   Basic equations

In this section, we propose the basic equations that we want to study in this article. Let us denote by $\Omega$ the three-dimensional domain occupied by the mixture, and let $\mathit{\boldsymbol{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.

2.1.   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$, $\rho_1, \, \rho_2$ are the mass density of every constituent, $\sigma_{ij}, \, \tau_{ij}$ 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 Ⅲ 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 $\tau_1$ and $\tau_2$, which can be different for the sake of generality, assumed to satisfy the condition $\tau_1\geq \tau_2$. The case $\tau_2\geq \tau_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. $\mathit{\boldsymbol{x}}\in \Omega$,

We will consider homogeneous Dirichlet boundary conditions. Therefore, we impose that

To study problems (2.4)–(2.6), it will be useful to manipulate the second equation of the system (2.4). We note that

and

and so, we can write the second equation of the system (2.4) as

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 $\Omega = (0, \ell)$, with $\ell > 0$.

2.2.   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

where

Here, $\rho_2^* = \rho_2\tau_2\tau_1^{-1}$, $\overline{A}_{ijkl} = A_{ijkl}^*-\tau_1 A_{ijkl}$, $\overline{B}_{ijkl} = B_{ijkl}^*-\tau_1 B_{ijkl}$ and $\overline{C}_{ijkl} = C_{ijkl}^*-\tau_1 C_{ijkl}$.

In view of these relations it will be natural to assume

${\rm{i}})$ $\rho_1(\mathit{\boldsymbol{x}})\geq \rho_1^0 > 0$ and $\rho_2(\mathit{\boldsymbol{x}})\geq \rho_2^0 > 0$.

${\rm{ii}})$ There exists a positive constant $C_1$ such that

for every tensors $\xi_{ij}$ and $\eta_{kl}$.

${\rm{iii}})$ There exists a positive constant $C_2$ such that

for every tensors $\xi_{ij}$ and $\eta_{ij}$.

${\rm{iv}})$ There exists a positive constant $C_3$ such that

for every vector $\xi_{i}$.

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 ⅳ) 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

ⅰ$^*$) $\rho_1 > 0$, $\rho_2 > 0$ and $a > 0$.

ⅱ$^*$) The matrices

are positive definite. That is, $AC > B^2$, $A > 0$, $\overline{A}\overline{C} > \overline{B}^2$ and $\overline{A} > 0$.

Remark 1. When $B = B^* = 0$ we can propose an alternative approach. We assume that, in this case, $\tau_2 > \tau_1$. We note that

and

Therefore, we have

Then, we can obtain the equality

where

In the previous definitions, we have used the notation:

3.   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 ⅰ)–ⅳ).

We first consider the Hilbert space

where $\mathit{\boldsymbol{W}}_0^{1, 2}(\Omega) = [W_0^{1, 2}(\Omega)]^3$ and $\mathit{\boldsymbol{L}}^2(\Omega) = [L^2(\Omega)]^3$, which is equipped by an inner product defined as

where $U = (u_i, v_i, z_i, w_i, y_i, r_i)$ and $U^* = (u_i^*, v_i^*, z_i^*, w_i^*, y_i^*, r_i^*)$, 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 $\mathcal{H}$.

We can write our problem in the following form:

where the matrix differential operator $\mathcal{A}$ is defined as

and its terms are given by

We note that the domain of the operator $\mathcal{A}$ is the subspace of $U\in\mathcal{H}$ such that $\mathcal{A}U\in \mathcal{H}$. It is clear that this is a dense subspace.

Lemma 1. There exists a positive constant $K_1$ such that

for every $U\in Dom(\mathcal{A})$.

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.

Lemma 2. There exists a positive constant $K_2$ such that $K_2 \mathit{\boldsymbol{I}}-\mathcal{A}$ is exhaustive.

Proof. Let us consider $\mathcal{F} = (\mathit{\boldsymbol{f}}_1, \ldots, \mathit{\boldsymbol{f}}_6)\in \mathcal{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 $(\mathit{\boldsymbol{F}}_1, \mathit{\boldsymbol{F}}_2)\in \mathit{\boldsymbol{W}}^{-1, 2}(\Omega)\times \mathit{\boldsymbol{W}}^{-1, 2}(\Omega).$

On the other side, if we define the bilinear form

where

it is clear that, for $K_2$ large enough, $\mathcal{B}$ is a coercive and bounded bilinear form on $\mathit{\boldsymbol{W}}_0^{1, 2}(\Omega)\; \times \; \mathit{\boldsymbol{W}}_0^{1, 2}(\Omega)$. In view of the Lax-Milgram theorem, we obtain the existence of $(\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{w}})$ satisfying the system (3.1). Then, we also obtain the existence of $(\mathit{\boldsymbol{v}}, \mathit{\boldsymbol{z}}, \mathit{\boldsymbol{y}}, \mathit{\boldsymbol{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 $\mathcal{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\in Dom(\mathcal{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)$.

4.   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 $\Omega = (0, \ell)$, $\ell > 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 $\overline{A}$, $\overline{B}$ and $\overline{C}$, are positive definite and so, the resulting semigroup associated to operator $\mathcal{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

and

then we can consider the bilinear form defined by the matrix:

where we recall that $\tau_1-\tau_2\geq 0$.

We note that after a simple calculation we find that

and

It is worth noting that, if matrix $\mathcal{M}$ is positive definite, then $Re\langle \mathcal{A}U, U\rangle\leq 0$ for every $U\in \hbox{Dom }(\mathcal{A})$ and, therefore, the semigroup associated to operator $\mathcal{A}$ is contractive.

Remark 4. ${\mathcal {M}}$ is positive definite if and only the following conditions hold:

We see that whenever $\tau_1-\tau_2$ is small enough in comparison with $m\lambda_1^*$, and assuming that $a$ and $\tau_1$ are not very large, we can guarantee that ${\mathcal {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

is a sufficient condition to guarantee that ${\mathcal {M}}$ is positive definite.

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 $\omega$ such that

for every $U\in \hbox{Dom }(\mathcal{A})$.

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 $\lambda\in \mathbb{R}$ such that $i\lambda$ belongs to the spectrum, then there exist a sequence of real numbers $\lambda_n\rightarrow \lambda\neq 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 $\mathcal{A}$ such that

That is, we have the following convergences:

In view of the dissipation inequality we see that $v_n, y_n \rightarrow 0$ in $W^{1, 2}(0, \ell)$ and $r_n\rightarrow 0$ in $L^2(0, \ell)$. Then, we also obtain that $\lambda_nw_n, \lambda_n u_n \rightarrow 0$ in $W^{1, 2}(0, \ell)$. Now, if we multiply the third convergence by $z_n$ we also obtain that $z_n\rightarrow 0$ in $L^2(0, \ell)$. 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 $\lambda_n\rightarrow \infty$ and a sequence of unit norm vectors at the domain of the operator $\mathcal{A}$ such that the condition holds. From this point, we can follow the same argument we used before because the only condition is $\lambda\neq 0$. In short, the theorem has been proved.

Remark 5. It is suitable to consider the case $\tau_1 = \tau_2$. In this case, we cannot assume that the matrix $\mathcal{M}$ is positive definite. Nevertheless, we can obtain the same conclusion. Although we cannot obtain that $r_n\rightarrow 0$ from the dissipation inequality, we can see that $y_n, v_n\rightarrow 0$ in $W^{1, 2}_0(0, \ell)$ and then, $\lambda_n u_n, \lambda_n w_n\rightarrow 0$ in $W^{1, 2}_0(0, \ell)$. From the convergences third and sixth we can obtain that $z_n, r_n\rightarrow 0$ in $L^2(0, \ell)$ 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 $\tau_1 = \tau_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 ⅰ$^*$) and ⅱ$^*$) hold.

5.   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 $\Omega$ over a finite time interval $[0, T]$, for a given final time $T > 0$.

5.1.   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 $\mathit{\boldsymbol{v}} = \dot{ \mathit{\boldsymbol{u}}}$ and $\mathit{\boldsymbol{y}} = \dot{ \mathit{\boldsymbol{w}}}$, and the accelerations of the first and second constituents, denoted by $\mathit{\boldsymbol{z}} = \ddot{ \mathit{\boldsymbol{u}}} = \dot{ \mathit{\boldsymbol{v}}}$ and $\mathit{\boldsymbol{r}} = \ddot{ \mathit{\boldsymbol{w}}} = \dot{ \mathit{\boldsymbol{y}}}$, respectively. Moreover, let $V$ be the variational space given by $V = [H^1_0(\Omega)]^3$.

In order to simplify the writing, and for the sake of clarity, we define the following operators, for all $\mathit{\boldsymbol{u}} = (u_i), \, \mathit{\boldsymbol{v}} = (v_i)\in V,$

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 $\mathit{\boldsymbol{z}}:[0, T]\rightarrow V$ and the acceleration of the second constituent $\mathit{\boldsymbol{r}}:[0, T]\rightarrow V$ such that $\mathit{\boldsymbol{z}}(0) = \mathit{\boldsymbol{z}}^0 = (z_i^0)$, $\mathit{\boldsymbol{r}}(0) = \mathit{\boldsymbol{r}}^0 = (r_i^0)$ and, for a.e. $t\in [0, T]$ and for all $\mathit{\boldsymbol{\xi}}, \mathit{\boldsymbol{\eta}}\in V$,

Here, the velocities and displacements of the first and second constituents are recovered from the relations:

where $\mathit{\boldsymbol{v}}^0 = (v_i^0)$, $\mathit{\boldsymbol{y}}^0 = (y_i^0)$, $\mathit{\boldsymbol{u}}^0 = (u_i^0)$ and $\mathit{\boldsymbol{w}}^0 = (w_i^0)$.

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 $\overline{\Omega}$ has a finite element triangulation $\mathcal{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(Tr)$ is the space of polynomials of degree less or equal to one in the element $Tr$, 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 $\mathit{\boldsymbol{u}}^{0h}$, $\mathit{\boldsymbol{v}}^{0h}$, $\mathit{\boldsymbol{z}}^{0h}$, $\mathit{\boldsymbol{w}}^{0h}$, $\mathit{\boldsymbol{y}}^{0h}$ and $\mathit{\boldsymbol{r}}^{0h}$, are given by

where $\mathcal{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 < \cdots < 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, \dots, N$. For a continuous function $z(t)$, we use the notation $z_n = z(t_n)$ and, for the sequence $\{z_n\}_{n = 0}^N$, we denote by $\delta 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 $\mathit{\boldsymbol{z}}^{hk} = \{ \mathit{\boldsymbol{z}}_n^{hk}\}_{n = 0}^N \subset V^h$ and the discrete acceleration of the second constituent $\mathit{\boldsymbol{r}}^{hk} = \{ \mathit{\boldsymbol{r}}_n^{hk}\}_{n = 0}^N \subset V^h$ such that $\mathit{\boldsymbol{z}}_0^{hk} = \mathit{\boldsymbol{z}}^{0h}$, $\mathit{\boldsymbol{r}}_0^{hk} = \mathit{\boldsymbol{r}}^{0h}$ and, for all $n = 1, \ldots, N$ and $\mathit{\boldsymbol{\xi}}^h, \mathit{\boldsymbol{\eta}}^h\in V^h$,

where the discrete velocities and the discrete displacements of the first and second constituents are recovered from the relations:

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.

5.2.   A priori error estimates

Now, we aim to obtain some a priori error estimates on the numerical errors $\mathit{\boldsymbol{z}}_n- \mathit{\boldsymbol{z}}_n^{hk}$ and $\mathit{\boldsymbol{r}}_n- \mathit{\boldsymbol{r}}_n^{hk}$. We have the following main a priori error estimates result.

Theorem 4. Under the assumptions ⅰ)–ⅳ), if we denote by $(\mathit{\boldsymbol{z}}, \mathit{\boldsymbol{r}})$ the solution to problems (5.1)–(5.4) and by $(\mathit{\boldsymbol{z}}^{hk}, \mathit{\boldsymbol{r}}^{hk}$ the solution to problems (5.7)–(5.10), then we have the following a priori error estimates, for all $\{ \mathit{\boldsymbol{\xi}}_n^{h}\}_{n = 0}^N, \, \{ \mathit{\boldsymbol{\eta}}_n^{h}\}_{n = 0}^N\subset V^h$,

where $M > 0$ is a positive constant which is be independent of the discretization parameters $h$ and $k$, but depending on the continuous solution, $\delta \mathit{\boldsymbol{z}}_j = (\mathit{\boldsymbol{z}}_j- \mathit{\boldsymbol{z}}_{j-1})/k$, $\delta \mathit{\boldsymbol{v}}_j = (\mathit{\boldsymbol{v}}_j- \mathit{\boldsymbol{v}}_{j-1})/k$, $\delta \mathit{\boldsymbol{r}}_j = (\mathit{\boldsymbol{r}}_j- \mathit{\boldsymbol{r}}_{j-1})/k$ and $\delta \mathit{\boldsymbol{y}}_j = (\mathit{\boldsymbol{y}}_j- \mathit{\boldsymbol{y}}_{j-1})/k$, and $I_j^1$ and $I_j^2$ are the integration errors given by

Proof. First, we obtain some estimates for the acceleration of the first constituent $\mathit{\boldsymbol{z}}_n$. Then, we subtract variational equation (5.1) at time $t = t_n$ for a test function $\mathit{\boldsymbol{\xi}} = \mathit{\boldsymbol{\xi}}^h\in V^h\subset V$ and discrete variational equation (5.7) to obtain, for all $\mathit{\boldsymbol{\xi}}^h\in V^h$,

and so, we have, for all $\mathit{\boldsymbol{\xi}}^h\in 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 $\mathit{\boldsymbol{\xi}}^h\in V^h$,

Proceeding in a similar form, we obtain the following error estimates for the acceleration of the second constituent, for all $\mathit{\boldsymbol{\eta}}^h\in V^h$,

Combining the previous two estimates, we get, for all $\mathit{\boldsymbol{\xi}}^h, \, \mathit{\boldsymbol{\eta}}^h\in V^h$,

Now, we observe that

and, thanks to assumptions ⅲ), it follows that

Multiplying the above estimates by $k$ and summing up to $n$, we find that, for all $\{ \mathit{\boldsymbol{\xi}}_j^h\}_{j = 0}^n, \{ \mathit{\boldsymbol{\eta}}_j^h\}_{j = 0}^n\subset V^h$,

Now, keeping in mind that (thanks again to property iii)),

and, therefore, we have, for all $\{ \mathit{\boldsymbol{\xi}}_j^h\}_{j = 0}^n, \{ \mathit{\boldsymbol{\eta}}_j^h\}_{j = 0}^n\subset V^h$,

Finally, taking into account that

where $I_j^1$ and $I_j^2$ 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.

Corollary 1. Let the assumptions of Theorem 4 hold. Under the additional regularity (5.13) it follows that the approximations obtained by problems (5.7)–(5.10) are linearly convergent; that is, there exists a positive constant $M$, independent of the discretization parameters $h$ and $k$, such that

6.   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 $\mathit{\boldsymbol{u}}_{n-1}^{hk}$, $\mathit{\boldsymbol{v}}_{n-1}^{hk}$, $\mathit{\boldsymbol{z}}_{n-1}^{hk}$, $\mathit{\boldsymbol{w}}_{n-1}^{hk}$, $\mathit{\boldsymbol{y}}_{n-1}^{hk}$ and $\mathit{\boldsymbol{r}}_{n-1}^{hk}$ at time $t_{n-1}$, the discrete accelerations at time $t = t_n$, $\mathit{\boldsymbol{z}}_n^{hk}$ and $\mathit{\boldsymbol{r}}_n^{hk}$, are obtained by solving the discrete linear system, for all $\mathit{\boldsymbol{\xi}}^h, \mathit{\boldsymbol{\eta}}^h\in V^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 one-dimensional problem is considered.

with the following data:

By using the following initial conditions, for all $x\in(0, 1)$,

considering homogeneous Dirichlet boundary conditions in the boundaries and the (artificial) supply terms, for all $(x, t)\in(0, 1)\times (0, 1)$,

the exact solution to the above one-dimensional problem can be easily calculated and it has the form, for $(x, t)\in[0, 1]\times[0, 1]$:

Thus, the approximation errors estimated by

are shown in Table 1 for several values of the discretization parameters $h$ and $k$. The evolution of the error, depending on the parameter $h+k$, is plotted in Figure 1. The convergence of the algorithm is clearly observed, and the linear convergence, stated in Corollary 1, is achieved.

If we assume now that there are not supply terms, and we use the final time $T = 10$, the data

and the initial conditions, for all $x\in (0, 1)$,

taking the discretization parameters $h = k = 0.0001$, the evolution in time of the discrete energy

is plotted in Figure 2 (in both natural and semi-log scales). As demonstrated in this figure, it converges to zero and an exponential decay seems to be found.

6.2.   Second example: a two-dimensional problem

Now, we consider that the body occupies the two-dimensional domain $\Omega = (0, 1)\times(0, 1)$, and we use the following data:

with the initial conditions, for all $(x, y)\in [0, 1]\times [0, 1]$,

Regarding the boundary conditions, we assume that the boundary $\Gamma$ is decomposed into three disjoint parts $\Gamma_D$, $\Gamma_N$ and $\Gamma_F$, being $\Gamma_D = \{0\}\times [0, 1]$, $\Gamma_N = \{1\}\times [0, 1]$ and $\Gamma_F = \Gamma-\Gamma_D-\Gamma_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 $\Gamma_D$ we assume that the body is clamped and so,

meanwhile on the Neumann part $\Gamma_N$ we assume that a surface force $G$ is applied in the form:

and boundary $\Gamma_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.

Acknowledgments

The authors would like to express their gratitude to the anonymous referees, whose comments improved the quality of the paper.

This paper is part of the project PID2019-105118GB-I00, funded by the Spanish Ministry of Science, Innovation and Universities and FEDER "A way to make Europe".

Conflict of interest

The authors declare there is no conflict of interest.