Finite Element Error Analysis of a Viscoelastic Timoshenko Beam with Thermodiffusion Effects

Abstract

In this paper, a thermomechanical problem involving a viscoelastic Timoshenko beam is analyzed from a numerical point of view. The so-called thermodiffusion effects are also included in the model. The problem is written as a linear system composed of two second-order-in-time partial differential equations for the transverse displacement and the rotational movement, and two first-order-in-time partial differential equations for the temperature and the chemical potential. The corresponding variational formulation leads to a coupled system of first-order linear variational equations written in terms of the transverse velocity, the rotation speed, the temperature and the chemical potential. The existence and uniqueness of solutions, as well as the energy decay property, are stated. Then, we focus on the numerical approximation of this weak problem by using the implicit Euler scheme to discretize the time derivatives and the classical finite element method to approximate the spatial variable. A discrete stability property and some a priori error estimates are shown, from which we can conclude the linear convergence of the approximations under suitable additional regularity conditions. Finally, some numerical simulations are performed to demonstrate the accuracy of the scheme, the behavior of the discrete energy decay and the dependence of the solution with respect to some parameters.

1. Introduction

The Timoshenko system, introduced by S. P. Timoshenko in 1921 [1], is a well-known mathematical framework that analyzes shear deformation and rotational bending effects on a beam. In this model, the two main unknowns are the transverse displacement and the neutral axis rotation caused by bending. Since this pioneering work, the asymptotic behavior of this type of problem has deserved much attention. In this way, in [2,3] the energy decay rates were analyzed when the past history was assumed; in [4], the so-called Kelvin–Voigt dissipation was considered; in [5], a new dynamic thermoviscoelastic contact problem was studied including a frictional contact; in [6], the uniform stability of a partially dissipative viscoelastic Timoshenko system was obtained; in [7], the energy decay was proved for a Timoshenko system with a memory type; the second sound effect was included in [8,9]. The nonexponential and polynomial energy decay was shown in [10,11] for a Kelvin–Voigt damping with heat conduction; the well-posedness and the energy decay (polynomial or exponential, depending on the speeds of wave propagation) was studied in [12] for a thermoelastic laminated Timoshenko beam; type III thermoelasticity was included in [13,14]; the long-time dynamics of a Timoshenko system modeling vibrations of beams with nonlinear localized damping mechanisms was derived in [15]; the global stability of Timoshenko systems was considered in [16]; and the well-posedness and regularity of the viscoelastic Timoshenko beam model, including a comparison with the widely used Euler–Bernoulli model, is provided in [17].

The second main aspect of the problem studied in this paper is the so-called thermodiffusion effect. It started being considered after the Second World War, and it refers to the coupling among the fields of strain, temperature and mass diffusion, leading to a random walk of a set of particles from regions of high concentrations to regions of lower concentrations. It is worth noting that there are many engineering applications, such as satellite problems, aircraft landing on water, the manufacturing of integrated circuits or oil extraction. This new diffusion theory was introduced by Nowacki (see [18]) in the 1970s and developed later by Sherief et al. [19]. Recently, this kind of problem has been considered by other authors. For instance, in [20], a thermodiffusive problem with voids was analyzed from the analytical point of view; a thermodiffusive model of elasticity was studied analytically and numerically in [21]; a contact problem with a thermodiffusive material (not beam) was considered in [22]; a thermodiffusive problem with two time delays was analyzed in [23]; and the exponential stability of a thermodiffusive Timoshenko beam was proved in [24].

In this paper, we continue the research started by Ramos et al. [25], where the asymptotic behavior of a viscoelastic Timoshenko system was considered, including some numerical simulations performed with a spectral discretization method. Here, we focus on the theoretical numerical analysis of this problem: a finite element approximation is introduced and an a priori error analysis is provided, including a discrete stability property. This leads us to obtain the convergence order or the scheme under some additional regularity assumptions. Some numerical simulations are also shown to demonstrate the linear convergence in an academical example with a known exact solution, and the decay of the discrete energy is compared for cases with different dissipative mechanisms in the viscoelasticity, following the examples depicted in [25]. To our knowledge, this is the first time where this type of thermodiffusive viscoelastic problem is considered from the numerical point of view.

This paper is organized in four sections. In Section 2, the mathematical model is described, including the conditions on the constitutive coefficients, and the analytical results provided in [25] are recalled and the variational formulation of the thermomechanical problem is derived. Then, its fully discrete approximation is introduced in Section 3 by discretizing the time derivatives using the implicit Euler scheme and approximating the spatial variable by the classical finite element method. The stability of the discrete solutions and a priori error estimates are proved, from which a convergence order can be shown under adequate regularity conditions. Finally, the numerical results are shown in Section 4.

2. The Mathematical Model: Analytical Results

In this section, we present the model and state the main analytical results, including the existence and uniqueness of solutions and an energy decay property (see [25] for further details).

Let us denote the spatial variable by x, which lies in the interval $[0,\ell ]$; and t is the time, which goes from 0 to T, where $T>0$ represents the final time. Additionally, the subscripts x and t indicate the spatial and time derivatives, respectively.

Denoting by $\phi$ the transverse displacement, $\psi$ the rotation of the neutral axis due to bending, P the chemical potential and $\theta$ the temperature field, the thermomechanical problem involving a viscoelastic beam with thermodiffusion effects is written as follows:

Find the transverse displacement $\phi :[0,\ell ]\times [0,T]\to \mathbb{R}$, the rotation field $\psi :[0,\ell ]\times [0,T]\to \mathbb{R}$, the chemical potential $P:[0,\ell ]\times [0,T]\to \mathbb{R}$ and the temperature field $\theta :[0,\ell ]\times [0,T]\to \mathbb{R}$, such that

$$\begin{array}{ccc}& & {\rho }_1{\phi }_{tt}-\kappa {({\phi }_x+\psi )}_x-{\mu }_1{({\phi }_x+\psi )}_{xt}=0\mathrm{in}(0,\ell )\times (0,T),\hfill \\ & & {\rho }_2{\psi }_{tt}-\alpha {\psi }_{xx}+\kappa ({\phi }_x+\psi )-{\gamma }_1{\theta }_x-{\gamma }_2{P}_x+{\mu }_1{({\phi }_x+\psi )}_t\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & -{\mu }_2{\psi }_{xxt}=0\mathrm{in}(0,\ell )\times (0,T),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & c{\theta }_t+d{P}_t-K{\theta }_{xx}-{\gamma }_1{\psi }_{xt}=0\mathrm{in}(0,\ell )\times (0,T),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & d{\theta }_t+r{P}_t-H{P}_{xx}-{\gamma }_2{\psi }_{xt}=0\mathrm{in}(0,\ell )\times (0,T),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & \phi (0,t)=\phi (\ell ,t)=\psi (0,t)=\psi (\ell ,t)=0\mathrm{for}\mathrm{a}.\mathrm{e}.t\in (0,T),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & \theta (0,t)=\theta (\ell ,t)=P(0,t)=P(\ell ,t)=0\mathrm{for}\mathrm{a}.\mathrm{e}.t\in (0,T),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & \phi (x,0)={\phi }_0\left(x\right),{\phi }_t(x,0)={\xi }_0\left(x\right),\psi (x,0)={\psi }_0\left(x\right)\mathrm{for}\mathrm{a}.\mathrm{e}.x\in (0,\ell ),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& & {\psi }_t(x,0)={\zeta }_0\left(x\right),P(x,0)=P_0\left(x\right),\theta (x,0)={\theta }_0\left(x\right)\mathrm{for}\mathrm{a}.\mathrm{e}.x\in (0,\ell ),\hfill \end{array}$$

(8)

where ${\rho }_1=\rho A$ and ${\rho }_2=\rho I$, with $\rho >0$ being the material density, A the cross-sectional area and I the second moment of the cross-sectional area; $\kappa$ is related to the modulus of rigidity and the transverse shear factor; $K>0$ is the heat coefficient; $H>0$ is the mass diffusion coefficient; $\alpha$ is related to the Young’s modulus; ${\mu }_1\ge 0$ and ${\mu }_2\ge 0$ are viscosity coefficients; and coefficients c, r and d satisfy the condition $cr-d^2>0$. Moreover, ${\phi }_0$, ${\xi }_0$, ${\psi }_0$, ${\zeta }_0$, $P_0$ and ${\theta }_0$ are given initial conditions.

In the rest of this section, the weak form of the thermomechanical problem (1)–(8) is obtained. Denote by $Y=L^2(0,\ell )$, and $(·,·)$ the space and its scalar product, with the usual norm ${\parallel ·\parallel }_Y$. Moreover, let us define the variational space $V=H_0^1(0,\ell )$ with norm ${\parallel ·\parallel }_V$.

Let us denote by $\xi ={\phi }_t$ the transverse velocity, and $\zeta ={\psi }_t$ the rotation speed. Therefore, integrating Equations (1)–(4) by parts leads to the following variational formulation of problem (1)–(8).

Find the transverse velocity $\xi :[0,T]\to V$, the rotation speed $\zeta :[0,T]\to V$, the chemical potential $P:[0,T]\to V$ and the temperature $\theta :[0,T]\to V$, such that $\xi \left(0\right)={\xi }_0$, $\zeta \left(0\right)={\zeta }_0$, $\theta \left(0\right)={\theta }_0$ and $P\left(0\right)=P_0$, and for a.e., $t\in (0,T)$ and $w,r,s,m\in V$:

$$\begin{array}{ccc}& & {\rho }_1({\xi }_t\left(t\right),w)+\kappa ({\phi }_x\left(t\right),w_x)-\kappa ({\psi }_x\left(t\right),w_{)}+{\mu }_1({\xi }_x\left(t\right),w_x)-{\mu }_1({\zeta }_x\left(t\right),w)=0,\hfill \\ & & {\rho }_2({\zeta }_t\left(t\right),r)+\alpha ({\psi }_x\left(t\right),r_x)+\kappa ({\phi }_x\left(t\right)+\psi \left(t\right),r)-({\gamma }_1{\theta }_x\left(t\right)+{\gamma }_2{P}_x\left(t\right),r)\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & +{\mu }_1({\xi }_x\left(t\right)+\zeta \left(t\right),r)+{\mu }_2({\zeta }_x\left(t\right),r_x)=0,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}& & c({\theta }_t\left(t\right),s)+d(P_t\left(t\right),s)+K({\theta }_x\left(t\right),s_x)-{\gamma }_1({\zeta }_x\left(t\right),s)=0,\hfill \end{array}$$

(11)

$$\begin{array}{ccc}& & d({\theta }_t\left(t\right),m)+r(P_t\left(t\right),m)+H(P_x\left(t\right),m_x)-{\gamma }_2({\zeta }_x\left(t\right),m)=0,\hfill \end{array}$$

(12)

where the transverse displacement and the rotation $\phi$ and $\psi$ are recovered from the equations:

$$\phi \left(t\right)={\int }_0^t\xi \left(s\right)ds+{\phi }_0,\psi \left(t\right)={\int }_0^t\zeta \left(s\right)ds+{\psi }_0.$$

(13)

Recently, in [25], the following existence and uniqueness result was proved.

Theorem 1.

If the initial conditions have the regularity

$$\begin{array}{c}{\phi }_0,{\psi }_0,{\xi }_0,{\zeta }_0,P_0,{\theta }_0\in H^2(0,\ell )\cap H_0^1(0,\ell ),\hfill \end{array}$$

then problem (9)–(13) admits a unique solution with the regularity

$$\phi ,\psi ,P,\theta \in C([0,T];H^2(0,\ell ))\cap C^1([0,T];V).$$

Moreover, if ${\mu }_1>0$, then the solutions to problem (9)–(13) decay exponentially, and if ${\mu }_1=0$ and ${\mu }_2>0$, the solutions decay polynomially.

3. An a Priori Error Analysis

In this section, we present the variational problem (9)–(13), introduced in the previous section, from the numerical point of view. We will carry this out in two steps. First, in order to obtain its spatial approximation, we consider a uniform partition of the interval $[0,\ell ]$ into M subintervals $a_0=0<a_1<\dots <a_M=\ell$, considering a uniform length $h=a_{i+1}-a_i=\ell /M$. Thus, the variational space V is approximated by the finite-dimensional space $V^h\subset V$, defined as

$$\begin{array}{c}{V}^h=\{w^h\in C\left([0,\ell ]\right);w_{{|}_{[a_i,a_{i+1}]}}^h\in P_1\left([a_i,a_{i+1}]\right)i=0,\dots ,M-1,w^h\left(0\right)=w^h(\ell )=0\}.\hfill \end{array}$$

Here, the space $P_1\left([a_i,a_{i+1}]\right)$ denotes the space of polynomials of a 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, and $h>0$ is the spatial discretization parameter. Moreover, we define the corresponding discrete initial conditions ${\phi }_0^h$, ${\xi }_0^h$, ${\psi }_0^h$, ${\zeta }_0^h$, $P_0^h$ and ${\theta }_0^h$ in the following form:

$$\begin{array}{c}{\phi }_0^h={\mathcal{P}}^h{\phi }_0,{\xi }_0^h={\mathcal{P}}^h{\xi }_0,{\psi }_0^h={\mathcal{P}}^h{\psi }_0,{\zeta }_0^h={\mathcal{P}}^h{\zeta }_0,P_0^h={\mathcal{P}}^h{P}_0,{\theta }_0^h={\mathcal{P}}^h{\theta }_0,\hfill \end{array}$$

where ${\mathcal{P}}^h$ is the classical finite element interpolation operator over $V^h$ (see [26]).

As a second step, we will discretize the time derivatives. Thus, the time interval $[0,T]$ is split into N subintervals, denoted by $0=t_0<t_1<\dots <t_N=T$, with a uniform time step size $k=T/N$ and nodes $t_n=nk$ for $n=0,1,\dots ,N$. As a notation, for a continuous function $z\left(t\right)$, let $z_n=z\left(t_n\right)$, and for a sequence ${\left\{z_n\right\}}_{n=0}^N$, let $\delta z_n=(z_n-z_{n-1})/k$ be its divided differences.

Therefore, using the well-known implicit Euler scheme, the fully discrete approximations of variational problem (9)–(13) are the following:

Find the discrete transverse velocity ${\left\{{\xi }_n^{hk}\right\}}_{n=0}^N\subset V^h$, the discrete rotation speed ${\left\{{\zeta }_n^{hk}\right\}}_{n=0}^N\subset V^h$, the discrete chemical potential ${\left\{P_n^{hk}\right\}}_{n=0}^N\subset V^h$ and the discrete temperature ${\left\{{\theta }_n^{hk}\right\}}_{n=0}^N\subset V^h$, such that ${\xi }_0^{hk}={\xi }_0^h$, ${\zeta }_0^{hk}={\zeta }_0^h$, $P_0^{hk}=P_0^h$ and ${\theta }_0^{hk}={\theta }_0^h$, and for $n=1,\dots ,N$ and $w^h,r^h,s^h,m^h\in V^h$,

$$\begin{array}{ccc}& & {\rho }_1(\delta {\xi }_n^{hk},w^h)+\kappa ({\phi }_{nx}^{hk},w_x^h)-\kappa ({\psi }_{nx}^{hk},w^h)+{\mu }_1({\xi }_{nx}^{hk},w_x^h)-{\mu }_1({\zeta }_{nx}^{hk},w^h)=0,\hfill \\ & & {\rho }_2(\delta {\zeta }_n^{hk},r^h)+\alpha ({\psi }_{nx}^{hk},r_x^h)+\kappa ({\phi }_{nx}^{hk}+{\psi }_n^{hk},r^h)-({\gamma }_1{\theta }_{nx}^{hk}+{\gamma }_2{P}_{nx}^{hk},r^h)\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & +{\mu }_1({\xi }_{nx}^{hk}+{\zeta }_n^{hk},r^h)+{\mu }_2({\zeta }_{nx}^{hk},r_x^h)=0,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & c(\delta {\theta }_n^{hk},s^h)+d(\delta P_n^{hk},s^h)+K({\theta }_{nx}^{hk},s_x^h)-{\gamma }_1({\zeta }_{nx}^{hk},s^h)=0,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & d(\delta {\theta }_n^{hk},m^h)+r(\delta P_n^{hk},m^h)+H(P_{nx}^{hk},m_x^h)-{\gamma }_2({\zeta }_{nx}^{hk},m^h)=0,\hfill \end{array}$$

(17)

where the discrete transverse displacement and the discrete rotation, ${\phi }_n^{hk}$ and ${\psi }_n^{hk}$, are recovered from the equations:

$${\phi }_n^{hk}=k\sum _{j=1}^n{\xi }_j^{hk}+{\phi }_0^h,{\psi }_n^{hk}=k\sum _{j=1}^n{\zeta }_j^{hk}+{\psi }_0^h.$$

(18)

By using the well-known Lax–Milgram lemma and the assumptions imposed on the constitutive parameters, it is easy to obtain that the fully discrete problem (14)–(18) has a unique solution.

We will now prove a discrete stability property for problem (14)–(18).

Lemma 1.

Let the assumptions of Theorem 1 hold. It follows that the sequences $\{{\phi }_n^{hk},{\xi }_n^{hk},{\psi }_n^{hk},{\zeta }_n^{hk},{\theta }^{hk},$$P_n^{hk}\}$, generated by discrete problem (14)–(18), satisfy the stability estimate

$${\displaystyle \parallel {\phi }_n^{hk}{\parallel }_V+\parallel {\psi }_n^{hk}{\parallel }_V+\parallel {\xi }_n^{hk}{\parallel }_Y+\parallel {\zeta }_n^{hk}{\parallel }_Y+\parallel P_n^{hk}{\parallel }_Y+{\parallel {\theta }_n^{hk}\parallel }_Y\le C,}$$

where C is a positive constant, which is independent of the discretization parameters h and k.

Proof. 

Here, we assume that ${\mu }_1>0$, for the sake of generality in the estimates obtained below. Otherwise, some of the calculations would be even easier.

Taking ${\xi }_n^{hk}$ as a test function in discrete variational Equation (14), we find that

$$\begin{array}{c}{\rho }_1(\delta {\xi }_n^{hk},{\xi }_n^{hk})+\kappa ({\phi }_{nx}^{hk},{\xi }_{nx}^{hk})-\kappa ({\psi }_{nx}^{hk},{\xi }_n^{hk})+{\mu }_1({\xi }_{nx}^{hk},{\xi }_{nx}^{hk})-{\mu }_1({\zeta }_{nx}^{hk},{\xi }_n^{hk})=0.\hfill \end{array}$$

Using Cauchy’s inequality,

$$ab\le ϵa^2+\frac{1}{4ϵ}{b}^2,a,b,ϵ\in \mathbb{R},ϵ>0,$$

(19)

keeping in mind that

$$\begin{array}{c}{\displaystyle {\rho }_1(\delta {\xi }_n^{hk},{\xi }_n^{hk})\ge \frac{{\rho }_1}{2k}\left\{\parallel {\xi }_n^{hk}{\parallel }_Y^2-{\parallel {\xi }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ {\displaystyle \kappa ({\phi }_{nx}^{hk},{\xi }_{nx}^{hk})\ge \frac{\kappa }{2k}\left\{\parallel {\phi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\phi }_{n-1x}^{hk}\parallel }_Y^2\right\},}\hfill \\ |{\mu }_1({\zeta }_{nx}^{hk},{\xi }_n^{hk})|=|{\mu }_1({\zeta }_n^{hk},{\xi }_{nx}^{hk})|\le C\parallel {\zeta }_n^{hk}{\parallel }_Y^2+\epsilon {\parallel {\xi }_{nx}^{hk}\parallel }_Y^2,\hfill \end{array}$$

where $\epsilon >0$ is assumed small enough, and the well-known Cauchy–Schwarz inequality, we obtain the following estimates for the transverse velocity:

$$\begin{array}{c}{\displaystyle \frac{{\rho }_1}{2k}\left\{\parallel {\xi }_n^{hk}{\parallel }_Y^2-{\parallel {\xi }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{\kappa }{2k}\left\{\parallel {\phi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\phi }_{n-1x}^{hk}\parallel }_Y^2\right\}+C{\parallel {\xi }_{nx}^{hk}\parallel }_Y^2}\hfill \\ \le C\left(\parallel {\zeta }_n^{hk}{\parallel }_Y^2+\parallel {\psi }_{nx}^{hk}{\parallel }_Y^2+{\parallel {\xi }_n^{hk}\parallel }_Y^2\right).\hfill \end{array}$$

(20)

Now, we obtain the estimates for the rotation speed. Thus, taking ${\zeta }_n^{hk}$ as a test function in discrete variational Equation (15), we have

$$\begin{array}{c}{\rho }_2(\delta {\zeta }_n^{hk},{\zeta }_n^{hk})+\alpha ({\psi }_{nx}^{hk},{\zeta }_{nx}^{hk})+\kappa ({\phi }_{nx}^{hk}+{\psi }_n^{hk},{\zeta }_n^{hk})-({\gamma }_1{\theta }_{nx}^{hk}+{\gamma }_2{P}_{nx}^{hk},{\zeta }_n^{hk})\hfill \\ +{\mu }_1({\xi }_{nx}^{hk}+{\zeta }_n^{hk},{\zeta }_n^{hk})+{\mu }_2({\zeta }_{nx}^{hk},{\zeta }_{nx}^{hk})=0.\hfill \end{array}$$

Taking into account that

$$\begin{array}{c}{\displaystyle {\rho }_2(\delta {\zeta }_n^{hk},{\zeta }_n^{hk})\ge \frac{{\rho }_2}{2k}\left\{\parallel {\zeta }_n^{hk}{\parallel }_Y^2-{\parallel {\zeta }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ {\displaystyle \alpha ({\psi }_{nx}^{hk},{\zeta }_{nx}^{hk})\ge \frac{\alpha }{2k}\left\{\parallel {\psi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1x}^{hk}\parallel }_Y^2\right\},}\hfill \\ {\displaystyle \kappa ({\psi }_n^{hk},{\zeta }_n^{hk})\ge \frac{\kappa }{2k}\left\{\parallel {\psi }_n^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ {\displaystyle {\mu }_2({\zeta }_{nx}^{hk},{\zeta }_{nx}^{hk})\ge 0,}\hfill \end{array}$$

and using again inequality (19) and the Cauchy–Schwarz inequality, after several simple calculations, it leads to the following estimates:

$$\begin{array}{c}{\displaystyle \frac{{\rho }_2}{2k}\left\{\parallel {\zeta }_n^{hk}{\parallel }_Y^2-{\parallel {\zeta }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{\alpha }{2k}\left\{\parallel {\psi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1x}^{hk}\parallel }_Y^2\right\}-({\gamma }_1{\theta }_{nx}^{hk}+{\gamma }_2{P}_{nx}^{hk},{\zeta }_n^{hk})}\hfill \\ {\displaystyle +\frac{\kappa }{2k}\left\{\parallel {\psi }_n^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1}^{hk}\parallel }_Y^2\right\}}\hfill \\ \le C\left(\parallel {\phi }_{nx}^{hk}{\parallel }_Y^2+ϵ\parallel {\xi }_{nx}^{hk}{\parallel }_Y^2+{\parallel {\zeta }_n^{hk}\parallel }_Y^2\right),\hfill \end{array}$$

(21)

where $\epsilon >0$ is assumed small enough.

Thirdly, we provide the estimates on the discrete temperature. Therefore, taking ${\theta }_n^{hk}$ as a test function in variational Equation (16), we obtain

$$\begin{array}{c}c(\delta {\theta }_n^{hk},{\theta }_n^{hk})+d(\delta P_n^{hk},{\theta }_n^{hk})+K({\theta }_{nx}^{hk},{\theta }_{nx}^{hk})-{\gamma }_1({\zeta }_{nx}^{hk},{\theta }_n^{hk})=0.\hfill \end{array}$$

By using inequality (19) several times, and keeping in mind that

$$\begin{array}{c}{\displaystyle c(\delta {\theta }_n^{hk},{\theta }_n^{hk})\ge \frac{c}{2k}\left\{\parallel {\theta }_n^{hk}{\parallel }_Y^2-\parallel {\theta }_{n-1}^{hk}{\parallel }_Y^2+{\parallel {\theta }_n^{hk}-{\theta }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ K({\theta }_{nx}^{hk},{\theta }_{nx}^{hk})\ge 0,\hfill \end{array}$$

we find the following estimates:

$$\begin{array}{c}{\displaystyle \frac{c}{2k}\left\{\parallel {\theta }_n^{hk}{\parallel }_Y^2-\parallel {\theta }_{n-1}^{hk}{\parallel }_Y^2+{\parallel {\theta }_n^{hk}-{\theta }_{n-1}^{hk}\parallel }_Y^2\right\}-{\gamma }_1({\zeta }_{nx}^{hk},{\theta }_n^{hk})+d(\delta P_n^{hk},{\theta }_n^{hk})\le 0.}\hfill \end{array}$$

(22)

Finally, we provide the estimates on the discrete chemical potential. Therefore, taking now as a test function $m^h=P_n^{hk}$ in discrete variational Equation (17), we have

$$d(\delta {\theta }_n^{hk},P_n^{hk})+r(\delta P_n^{hk},P_n^{hk})+H(P_{nx}^{hk},P_{nx}^{hk})-{\gamma }_2({\zeta }_{nx}^{hk},P_n^{hk})=0.$$

By again using inequality (19) several times, and keeping in mind that

$$\begin{array}{c}{\displaystyle r(\delta P_n^{hk},P_n^{hk})\ge \frac{r}{2k}\left\{\parallel P_n^{hk}{\parallel }_Y^2-\parallel P_{n-1}^{hk}{\parallel }_Y^2+{\parallel P_n^{hk}-P_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ H(P_{nx}^{hk},P_{nx}^{hk})\ge 0,\hfill \end{array}$$

we find the following estimates:

$$\begin{array}{c}{\displaystyle \frac{r}{2k}\left\{\parallel P_n^{hk}{\parallel }_Y^2-\parallel P_{n-1}^{hk}{\parallel }_Y^2+{\parallel P_n^{hk}-P_{n-1}^{hk}\parallel }_Y^2\right\}-{\gamma }_2({\zeta }_{nx}^{hk},P_n^{hk})+d(\delta {\theta }_n^{hk},P_n^{hk})\le 0.}\hfill \end{array}$$

(23)

Combining estimates (22) and (23), observing that

$$d(\delta {\theta }_n^{hk},P_n^{hk})+d(\delta P_n^{hk},{\theta }_n^{hk})=\frac{d}{k}\left\{({\theta }_n^{hk},P_n^{hk})-({\theta }_{n-1}^{hk},P_{n-1}^{hk})+({\theta }_n^{hk}-{\theta }_{n-1}^{hk},P_n^{hk}-P_{n-1}^{hk})\right\},$$

we find that

$$\begin{array}{c}{\displaystyle \frac{c}{2k}\left\{\parallel {\theta }_n^{hk}{\parallel }_Y^2-{\parallel {\theta }_{n-1}^{hk}\parallel }_Y^2\right\}-{\gamma }_1({\zeta }_{nx}^{hk},{\theta }_n^{hk})+\frac{r}{2k}\left\{\parallel P_n^{hk}{\parallel }_Y^2-{\parallel P_{n-1}^{hk}\parallel }_Y^2\right\}-{\gamma }_2({\zeta }_{nx}^{hk},P_n^{hk})}\hfill \\ {\displaystyle +\frac{d}{k}\left\{({\theta }_n^{hk},P_n^{hk})-({\theta }_{n-1}^{hk},P_{n-1}^{hk})\right\},}\hfill \end{array}$$

where we have used the constitutive condition $cr-d^2>0$.

Combining the previous estimates with (20) and (21), keeping in mind that

$$-{\gamma }_1({\zeta }_{nx}^{hk},{\theta }_n^{hk})-{\gamma }_2({\zeta }_{nx}^{hk},P_n^{hk})={\gamma }_1({\zeta }_n^{hk},{\theta }_{nx}^{hk})+{\gamma }_2({\zeta }_n^{hk},P_{nx}^{hk}),$$

we obtain

$$\begin{array}{c}{\displaystyle \frac{c}{2k}\left\{\parallel {\theta }_n^{hk}{\parallel }_Y^2-{\parallel {\theta }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{{\rho }_1}{2k}\left\{\parallel {\xi }_n^{hk}{\parallel }_Y^2-{\parallel {\xi }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{r}{2k}\left\{\parallel P_n^{hk}{\parallel }_Y^2-{\parallel P_{n-1}^{hk}\parallel }_Y^2\right\}}\hfill \\ {\displaystyle +\frac{\kappa }{2k}\left\{\parallel {\phi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\phi }_{n-1x}^{hk}\parallel }_Y^2\right\}+\frac{{\rho }_2}{2k}\left\{\parallel {\zeta }_n^{hk}{\parallel }_Y^2-{\parallel {\zeta }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{\alpha }{2k}\left\{\parallel {\psi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1x}^{hk}\parallel }_Y^2\right\}}\hfill \\ {\displaystyle +\frac{d}{k}\left\{({\theta }_n^{hk},P_n^{hk})-({\theta }_{n-1}^{hk},P_{n-1}^{hk})\right\}}\hfill \\ {\displaystyle \le C\left(\parallel {\psi }_n^{hk}{\parallel }_Y^2+\parallel {\phi }_{nx}^{hk}{\parallel }_Y^2+\parallel {\zeta }_n^{hk}{\parallel }_Y^2+\parallel {\psi }_{nx}^{hk}{\parallel }_Y^2+{\parallel {\xi }_n^{hk}\parallel }_Y^2\right).}\hfill \end{array}$$

Multiplying the above estimates by k and summing up to n, we have

$$\begin{array}{c}{\displaystyle c\parallel {\theta }_n^{hk}{\parallel }_Y^2+{\rho }_1\parallel {\xi }_n^{hk}{\parallel }_Y^2+r\parallel P_n^{hk}{\parallel }_Y^2+\kappa \parallel {\phi }_{nx}^{hk}{\parallel }_Y^2+{\rho }_2\parallel {\zeta }_n^{hk}{\parallel }_Y^2+\alpha {\parallel {\psi }_{nx}^{hk}\parallel }_Y^2+2d({\theta }_n^{hk},P_n^{hk})}\hfill \\ {\displaystyle \le Ck\sum _{j=1}^n\left(\parallel {\psi }_j^{hk}{\parallel }_Y^2+\parallel {\phi }_{jx}^{hk}{\parallel }_Y^2+\parallel {\zeta }_j^{hk}{\parallel }_Y^2+\parallel {\psi }_{jx}^{hk}{\parallel }_Y^2+{\parallel {\xi }_j^{hk}\parallel }_Y^2\right)+C(\parallel {\theta }_0^h{\parallel }_Y^2+{\parallel {\xi }_0^h\parallel }_Y^2}\hfill \\ +\parallel P_0^h{\parallel }_Y^2+\parallel {\phi }_0^h{\parallel }_V^2+\parallel {\zeta }_0^h{\parallel }_Y^2+{\parallel {\psi }_0^h\parallel }_V^2).\hfill \end{array}$$

Finally, since by using condition $cr-d^2>0$ again, it follows that

$$r\parallel P_n^{hk}{\parallel }_Y^2+c\parallel {\theta }_n^{hk}{\parallel }_Y^2+2d({\theta }_n^{hk},P_n^{hk})\ge C(\parallel P_n^{hk}{\parallel }_Y^2+\parallel {\theta }_n^{hk}{\parallel }_Y^2),$$

applying a discrete version of Gronwall’s inequality (see, e.g., [27]), we conclude the desired stability property. □

In the remainder of this section, we will derive some a priori error estimates on the numerical errors given by ${\xi }_n-{\xi }_n^{hk}$, ${\zeta }_n-{\zeta }_n^{hk}$, $P_n-P_n^{hk}$ and ${\theta }_n-{\theta }_n^{hk}$.

Theorem 2.

Let the assumptions of Theorem 1 still hold. If we denote by $(\phi ,\xi ,\psi ,\zeta ,P,\theta )$ the solution to problem (9)–(13) and by ${\{{\phi }_n^{hk},{\xi }_n^{hk},{\psi }_n^{hk},{\zeta }_n^{hk},P_n^{hk},{\theta }_n^{hk}\}}_{n=0}^N$ the solution to problem (14)–(18), then we have the following a priori error estimates, for all $w^h={\left\{w_j^h\right\}}_{j=0}^N,r^h={\left\{r_j^h\right\}}_{j=0}^N,$$s^h={\left\{s_j^h\right\}}_{j=0}^N,m^h={\left\{m_j^h\right\}}_{j=0}^N\subset V^h$,

$$\begin{array}{ccc}& & {\displaystyle \underset{0\le n\le N}{max}\{\parallel {\xi }_n-{\xi }_n^{hk}{\parallel }_Y^2+\parallel {\phi }_{nx}-{\phi }_{nx}^{hk}{\parallel }_Y^2+\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2+\parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{\parallel }_Y^2+{\parallel {\psi }_n-{\psi }_n^{hk}\parallel }_Y^2}\hfill \\ & & +\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2+{\parallel P_n-P_n^{hk}\parallel }_Y^2\}\hfill \\ & & \le Ck\sum _{j=1}^N(\parallel {\xi }_{tj}-\delta {\xi }_j{\parallel }_Y^2+\parallel {\phi }_{tj}-\delta {\phi }_j{\parallel }_V^2+\parallel {\xi }_j-w_j^h{\parallel }_V^2+{\parallel {\zeta }_{tj}-\delta {\zeta }_j\parallel }_Y^2\hfill \\ & & +\parallel {\psi }_{tj}-\delta {\psi }_j{\parallel }_V^2+\parallel {\zeta }_j-r_j^h{\parallel }_V^2+\parallel {\theta }_{tj}-\delta {\theta }_j{\parallel }_Y^2+{\parallel {\theta }_j-s_j^h\parallel }_V^2\hfill \\ & & +\parallel P_{tj}-\delta P_j{\parallel }_Y^2+\parallel P_j-m_j^h{\parallel }_V^2)+C\underset{0\le n\le N}{max}\parallel {\xi }_n-w_n^h{\parallel }_Y^2+C\underset{0\le n\le N}{max}{\parallel {\zeta }_n-r_n^h\parallel }_Y^2\hfill \end{array}$$

$$\begin{array}{ccc}& & {\displaystyle +C\underset{0\le n\le N}{max}\parallel P_n-m_n^h{\parallel }_Y^2+C\underset{0\le n\le N}{max}{\parallel {\theta }_n-s_n^h\parallel }_Y^2}\hfill \\ & & {\displaystyle +\frac{C}{k}\sum _{j=1}^{N-1}(\parallel {\xi }_j-w_j^h-({\xi }_{j+1}-w_{j+1}^h){\parallel }_Y^2+{\parallel {\zeta }_j-r_j^h-({\zeta }_{j+1}-r_{j+1}^h)\parallel }_Y^2}\hfill \\ & & +\parallel {\theta }_j-s_j^h-({\theta }_{j+1}-s_{j+1}^h){\parallel }_Y^2+\parallel P_j-m_j^h-(P_{j+1}-m_{j+1}^h){\parallel }_Y^2)+C({\parallel {\xi }_0-{\xi }_0^h\parallel }_Y^2\hfill \\ & & +\parallel {\phi }_0-{\phi }_0^h{\parallel }_V^2+\parallel {\zeta }_0-{\zeta }_0^h{\parallel }_Y^2+\parallel {\psi }_0-{\psi }_0^h{\parallel }_V^2+\parallel P_0-P_0^h{\parallel }_Y^2+{\parallel {\theta }_0-{\theta }_0^h\parallel }_Y^2),\hfill \end{array}$$

where C is again a positive constant, which does not depend on parameters h and k.

Proof. 

Again, for the sake of generality in the calculations presented below, we will assume that ${\mu }_1>0$.

As a first step, we obtain some error estimates for the transverse velocity. So, we subtract variational Equation (9) at time $t=t_n$ for a test function $w=w^h\in V^h\subset V$ and discrete variational Equation (14) to obtain, for all $w^h\in V^h$,

$$\begin{array}{c}{\rho }_1({\xi }_{tn}-\delta {\xi }_n^{hk},w^h)+\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},w_x^h)-\kappa ({\psi }_{nx}-{\psi }_{nx}^{hk},w^h)\hfill \\ +{\mu }_1({\xi }_{nx}-{\xi }_{nx}^{hk},w_x^h)-{\mu }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},w^h)=0.\hfill \end{array}$$

Therefore, it follows that for all $w^h\in V^h$,

$$\begin{array}{c}{\rho }_1({\xi }_{tn}-\delta {\xi }_n^{hk},{\xi }_n-{\xi }_n^{hk})+\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},{\xi }_{nx}-{\xi }_{nx}^{hk})-\kappa ({\psi }_{nx}-{\psi }_{nx}^{hk},{\xi }_n-{\xi }_n^{hk})\hfill \\ +{\mu }_1({\xi }_{nx}-{\xi }_{nx}^{hk},{\xi }_{nx}-{\xi }_{nx}^{hk})-{\mu }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\xi }_n-{\xi }_n^{hk})\hfill \\ ={\rho }_1({\xi }_{tn}-\delta {\xi }_n^{hk},{\xi }_n-w^h)+\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},{\xi }_{nx}-w_x^h)-\kappa ({\psi }_{nx}-{\psi }_{nx}^{hk},{\xi }_n-w^h)\hfill \\ +{\mu }_1({\xi }_{nx}-{\xi }_{nx}^{hk},{\xi }_{nx}-w_x^h)-{\mu }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\xi }_n-w^h).\hfill \end{array}$$

Keeping in mind that

$$\begin{array}{c}{\rho }_1({\xi }_{tn}-\delta {\xi }_n^{hk},{\xi }_n-{\xi }_n^{hk})={\rho }_1({\xi }_{tn}-\delta {\xi }_n,{\xi }_n-{\xi }_n^{hk})+{\rho }_1(\delta {\xi }_{tn}-\delta {\xi }_n^{hk},{\xi }_n-{\xi }_n^{hk}),\hfill \\ {\displaystyle {\rho }_1(\delta {\xi }_{tn}-\delta {\xi }_n^{hk},{\xi }_n-{\xi }_n^{hk})\ge \frac{{\rho }_1}{2k}\left\{\parallel {\xi }_n-{\xi }_n^{hk}{\parallel }_Y^2-{\parallel {\xi }_{n-1}-{\xi }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ \kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},{\xi }_{nx}-{\xi }_{nx}^{hk})=\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},{\phi }_{tnx}-\delta {\phi }_{nx})+\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},\delta {\phi }_{nx}-\delta {\phi }_{nx}^{hk}),\hfill \\ {\displaystyle \kappa ({\phi }_{nx}-{\phi }_{nx}^{hk},\delta {\phi }_{nx}-\delta {\phi }_{nx}^{hk})\ge \frac{\kappa }{2k}\left\{\parallel {\phi }_{nx}-{\phi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\phi }_{n-1x}-{\phi }_{n-1x}^{hk}\parallel }_Y^2\right\},}\hfill \\ |{\mu }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\xi }_n-{\xi }_n^{hk})|\le C\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2+\epsilon {\parallel {\xi }_{nx}-{\xi }_{nx}^{hk}\parallel }_Y^2,\hfill \end{array}$$

where $\epsilon >0$ is assumed small enough, using several time inequalities (19), we obtain the following estimates for the transverse velocity field, for all $w^h\in V^h$,

$$\begin{array}{ccc}& & \frac{{\rho }_1}{2k}\left\{\parallel {\xi }_n-{\xi }_n^{hk}{\parallel }_Y^2-{\parallel {\xi }_{n-1}-{\xi }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{\kappa }{2k}\left\{\parallel {\phi }_{nx}-{\phi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\phi }_{n-1x}-{\phi }_{n-1x}^{hk}\parallel }_Y^2\right\}\hfill \\ & & \le C(\parallel {\xi }_{tn}-\delta {\xi }_n{\parallel }_Y^2+\parallel {\phi }_{tn}-\delta {\phi }_n{\parallel }_V^2+\parallel {\xi }_n-w^h{\parallel }_V^2+\parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{)\parallel }_Y^2+{\parallel {\xi }_n-{\xi }_n^{hk}\parallel }_Y^2\hfill \\ & & +\parallel {\phi }_{nx}-{\phi }_{nx}^{hk}{\parallel }_Y^2+(\delta {\xi }_n-\delta {\xi }_n^{hk},{\xi }_n-w^h)).\hfill \end{array}$$

(24)

Now, we obtain the error estimates for the rotation speed. Thus, subtracting variational Equation (10) at time $t=t_n$ for a test function $r=r^h\in V^h\subset V$ and discrete variational Equation (15), we obtain, for all $r^h\in V^h$,

$$\begin{array}{c}{\rho }_2({\zeta }_{tn}-\delta {\zeta }_n^{hk},r^h)+\alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},r_x^h)-({\gamma }_1({\theta }_{nx}-{\theta }_{nx}^{hk})+{\gamma }_2(P_{nx}-P_{nx}^{hk}),r^h)\hfill \\ +\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk}+{\psi }_n-{\psi }_n^{hk},r^h)+{\mu }_1({\xi }_{nx}-{\xi }_{nx}^{hk}+{\zeta }_n-{\zeta }_n^{hk},r^h)+{\mu }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},r_x^h)=0,\hfill \end{array}$$

and so, for all $r^h\in V^h$, we have

$$\begin{array}{c}{\rho }_2({\zeta }_{tn}-\delta {\zeta }_n^{hk},{\zeta }_n-{\zeta }_n^{hk})+\alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},{\zeta }_{nx}-{\zeta }_{nx}^{hk})+\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk}+{\psi }_n-{\psi }_n^{hk},{\zeta }_n-{\zeta }_n^{hk})\hfill \\ -({\gamma }_1({\theta }_{nx}-{\theta }_{nx}^{hk})+{\gamma }_2(P_{nx}-P_{nx}^{hk}),{\zeta }_n-{\zeta }_n^{hk})\hfill \\ +{\mu }_1({\xi }_{nx}-{\xi }_{nx}^{hk}+{\zeta }_n-{\zeta }_n^{hk},{\zeta }_n-{\zeta }_n^{hk})+{\mu }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\zeta }_{nx}-{\zeta }_{nx}^{hk})\hfill \\ ={\rho }_2({\zeta }_{tn}-\delta {\zeta }_n^{hk},{\zeta }_n-r^h)+\alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},{\zeta }_{nx}-r_x^h)+\kappa ({\phi }_{nx}-{\phi }_{nx}^{hk}+{\psi }_n-{\psi }_n^{hk},{\zeta }_n-r^h)\hfill \\ -({\gamma }_1({\theta }_{nx}-{\theta }_{nx}^{hk})+{\gamma }_2(P_{nx}-P_{nx}^{hk}),{\zeta }_n-r^h)\hfill \\ +{\mu }_1({\xi }_{nx}-{\xi }_{nx}^{hk}+{\zeta }_n-{\zeta }_n^{hk},{\zeta }_n-r^h)+{\mu }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\zeta }_{nx}-r_x^h).\hfill \end{array}$$

Keeping in mind that

$$\begin{array}{c}{\rho }_2({\zeta }_{tn}-\delta {\zeta }_n^{hk},{\zeta }_n-{\zeta }_n^{hk})={\rho }_2({\zeta }_{tn}-\delta {\zeta }_n,{\zeta }_n-{\zeta }_n^{hk})+{\rho }_2(\delta {\zeta }_{tn}-\delta {\zeta }_n^{hk},{\zeta }_n-{\zeta }_n^{hk}),\hfill \\ {\displaystyle {\rho }_2(\delta {\zeta }_{tn}-\delta {\zeta }_n^{hk},{\zeta }_n-{\zeta }_n^{hk})\ge \frac{{\rho }_2}{2k}\left\{\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2-{\parallel {\zeta }_{n-1}-{\zeta }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \\ \alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},{\zeta }_{nx}-{\zeta }_{nx}^{hk})=\alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},{\psi }_{tnx}-\delta {\psi }_{nx})+\alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},\delta {\psi }_{nx}-\delta {\psi }_{nx}^{hk}),\hfill \\ {\displaystyle \alpha ({\psi }_{nx}-{\psi }_{nx}^{hk},\delta {\psi }_{nx}-\delta {\psi }_{nx}^{hk})\ge \frac{\alpha }{2k}\left\{\parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1x}-{\psi }_{n-1x}^{hk}\parallel }_Y^2\right\},}\hfill \\ {\displaystyle \kappa ({\psi }_n-{\psi }_n^{hk},\delta {\psi }_n-\delta {\psi }_n^{hk})\ge \frac{\kappa }{2k}\left\{\parallel {\psi }_n-{\psi }_n^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1}-{\psi }_{n-1}^{hk}\parallel }_Y^2\right\},}\hfill \end{array}$$

using inequality (19) and the Cauchy–Schwarz inequality again, we obtain the following estimates for all $r^h\in V^h$:

$$\begin{array}{ccc}& & \frac{{\rho }_2}{2k}\left\{\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2-{\parallel {\zeta }_{n-1}-{\zeta }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{\alpha }{2k}\left\{\parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1x}-{\psi }_{n-1x}^{hk}\parallel }_Y^2\right\}\hfill \\ & & +\frac{\kappa }{2k}\left\{\parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{\parallel }_Y^2-{\parallel {\psi }_{n-1x}-{\psi }_{n-1x}^{hk}\parallel }_Y^2\right\}\hfill \\ & & -({\gamma }_1({\theta }_{nx}-{\theta }_{nx}^{hk})+{\gamma }_2(P_{nx}-P_{nx}^{hk}),{\zeta }_n-{\zeta }_n^{hk})\hfill \\ & & \le C(\parallel {\zeta }_{tn}-\delta {\zeta }_n{\parallel }_Y^2+\parallel {\psi }_{tn}-\delta {\psi }_n{\parallel }_V^2+\parallel {\zeta }_n-r^h{\parallel }_V^2+\parallel {\phi }_{nx}-{\phi }_{nx}^{hk}{\parallel }_Y^2+{\parallel {\zeta }_n-{\zeta }_n^{hk}\parallel }_Y^2\hfill \\ & & +\parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{\parallel }_Y^2+{\parallel {\psi }_n-{\psi }_n^{hk}\parallel }_Y^2+(\delta {\zeta }_n-\delta {\zeta }_n^{hk},{\zeta }_n-r^h)).\hfill \end{array}$$

(25)

Thirdly, we derive some estimates for the temperature field. Then, we subtract variational Equation (11) at time $t=t_n$ for a test function $s=s^h\in V^h\subset V$ and discrete variational Equation (16), and we obtain, for all $s^h\in V^h$,

$$\begin{array}{c}c({\theta }_{tn}-\delta {\theta }_n^{hk},s^h)+d(P_{tn}-\delta P_n^{hk},s^h)+K({\theta }_{nx}-{\theta }_{nx}^{hk},s_x^h)-{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},s^h)=0,\hfill \end{array}$$

and so, we find that for all $s^h\in V^h$,

$$\begin{array}{c}c({\theta }_{tn}-\delta {\theta }_n^{hk},{\theta }_n-{\theta }_n^{hk})+d(P_{tn}-\delta P_n^{hk},{\theta }_n-{\theta }_n^{hk})+K({\theta }_{nx}-{\theta }_{nx}^{hk},{\theta }_{nx}-{\theta }_{nx}^{hk})\hfill \\ -{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\theta }_n-{\theta }_n^{hk})\hfill \\ =c({\theta }_{tn}-\delta {\theta }_n^{hk},{\theta }_n-s^h)+d(P_{tn}-\delta P_n^{hk},{\theta }_n-s^h)+K({\theta }_{nx}-{\theta }_{nx}^{hk},{\theta }_{nx}-s_x^h)\hfill \\ -{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\theta }_n-s^h).\hfill \end{array}$$

Taking into account that

$$\begin{array}{c}c({\theta }_{tn}-\delta {\theta }_n^{hk},{\theta }_n-{\theta }_n^{hk})=c({\theta }_{tn}-\delta {\theta }_n,{\theta }_n-{\theta }_n^{hk})+c(\delta {\theta }_n-\delta {\theta }_n^{hk},{\theta }_n-{\theta }_n^{hk}),\hfill \\ {\displaystyle c(\delta {\theta }_n-\delta {\theta }_n^{hk},{\theta }_n-{\theta }_n^{hk})\ge \frac{{\rho }_2}{2k}\{\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2-{\parallel {\theta }_{n-1}-{\theta }_{n-1}^{hk}\parallel }_Y^2}\hfill \\ +\parallel {\theta }_n-{\theta }_n^{hk}-({\theta }_{n-1}-{\theta }_{n-1}^{hk}){\parallel }_Y^2\},\hfill \\ -{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\theta }_n-s^h)={\gamma }_1({\zeta }_n-{\zeta }_n^{hk},{\theta }_{nx}-s_x^h),\hfill \end{array}$$

we have, for all $s^h\in V^h$,

$$\begin{array}{ccc}& & \frac{c}{2k}\left\{\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2-\parallel {\theta }_{n-1}-{\theta }_{n-1}^{hk}{\parallel }_Y^2+{\parallel {\theta }_n-{\theta }_n^{hk}-({\theta }_{n-1}-{\theta }_{n-1}^{hk})\parallel }_Y^2\right\}\hfill \\ & & +d(\delta P_n-\delta P_n^{hk},{\theta }_n-{\theta }_n^{hk})-{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\theta }_n-{\theta }_n^{hk})\hfill \\ & & \le C(\parallel {\theta }_{tn}-\delta {\theta }_n{\parallel }_Y^2+\parallel {\theta }_n-s^h{\parallel }_V^2+\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2+\parallel P_{tn}-\delta P_n{\parallel }_Y^2+{\parallel {\theta }_n-{\theta }_n^{hk}\parallel }_Y^2\hfill \\ & & +(\delta {\theta }_n-\delta {\theta }_n^{hk},{\theta }_n-s^h)+(\delta P_n-\delta P_n^{hk},{\theta }_n-s^h)).\hfill \end{array}$$

(26)

Finally, we obtain the estimates for the chemical potential. Thus, subtracting variational Equation (12) at time $t=t_n$ for a test function $m=m^h\in V^h\subset V$ and discrete variational Equation (17), we find that for all $m^h\in V^h$,

$$\begin{array}{c}d({\theta }_{tn}-\delta {\theta }_n^{hk},m^h)+r(P_{tn}-\delta P_n^{hk},m^h)+H(P_{nx}-P_{nx}^{hk},m_x^h)-{\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},m^h)=0,\hfill \end{array}$$

and so, we have, for all $m^h\in V^h$,

$$\begin{array}{c}d({\theta }_{tn}-\delta {\theta }_n^{hk},P_n-P_n^{hk})+r(P_{tn}-\delta P_n^{hk},P_n-P_n^{hk})+H(P_{nx}-P_{nx}^{hk},P_{nx}-P_{nx}^{hk})\hfill \\ -{\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},P_n-P_n^{hk})\hfill \\ =d({\theta }_{tn}-\delta {\theta }_n^{hk},P_n-m^h)+r(P_{tn}-\delta P_n^{hk},P_n-m^h)+H(P_{nx}-P_{nx}^{hk},P_{nx}-m_x^h)\hfill \\ -{\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},P_n-m^h).\hfill \end{array}$$

Taking into account that

$$\begin{array}{c}r(P_{tn}-\delta P_n^{hk},P_n-P_n^{hk})=r(P_{tn}-\delta P_n,P_n-P_n^{hk})+r(\delta P_n-\delta P_n^{hk},P_n-P_n^{hk}),\hfill \\ {\displaystyle r(\delta P_n-\delta P_n^{hk},P_n-P_n^{hk})\ge \frac{r}{2k}\{\parallel P_n-P_n^{hk}{\parallel }_Y^2-{\parallel P_{n-1}-P_{n-1}^{hk}\parallel }_Y^2}\hfill \\ +\parallel P_n-P_n^{hk}-(P_{n-1}-P_{n-1}^{hk}){\parallel }_Y^2\},\hfill \\ {\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},P_n-m^h)={\gamma }_2({\zeta }_n-{\zeta }_n^{hk},P_{nx}-m_x^h),\hfill \end{array}$$

we have, for all $m^h\in V^h$,

$$\begin{array}{ccc}& & \frac{r}{2k}\left\{\parallel P_n-P_n^{hk}{\parallel }_Y^2-\parallel P_{n-1}-P_{n-1}^{hk}{\parallel }_Y^2+{\parallel P_n-P_n^{hk}-(P_{n-1}-P_{n-1}^{hk})\parallel }_Y^2\right\}\hfill \\ & & +d(\delta {\theta }_n-\delta {\theta }_n^{hk},P_n-P_n^{hk})-{\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},P_n-P_n^{hk})\hfill \\ & & \le C(\parallel P_{tn}-\delta P_n{\parallel }_Y^2+\parallel P_n-m^h{\parallel }_V^2+\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2+{\parallel {\theta }_{tn}-\delta {\theta }_n\parallel }_Y^2\hfill \\ & & +\parallel P_n-P_n^{hk}{\parallel }_Y^2+(\delta {\theta }_n-\delta {\theta }_n^{hk},P_n-m^h)+(\delta P_n-\delta P_n^{hk},P_n-m^h)).\hfill \end{array}$$

(27)

Observing that

$$\begin{array}{c}{\displaystyle d(\delta {\theta }_n-\delta {\theta }_n^{hk},P_n-P_n^{hk})+d(\delta P_n-\delta P_n^{hk},{\theta }_n-{\theta }_n^{hk})}\hfill \\ {\displaystyle =\frac{d}{k}\{({\theta }_n-{\theta }_n^{hk},P_n-P_n^{hk})-({\theta }_{n-1}-{\theta }_{n-1}^{hk},P_{n-1}-P_{n-1}^{hk})}\hfill \\ +({\theta }_n-{\theta }_n^{hk}-({\theta }_{n-1}-{\theta }_{n-1}^{hk}),P_n-P_n^{hk}-(P_{n-1}-P_{n-1}^{hk}))\},\hfill \\ c\parallel {\theta }_n-{\theta }_n^{hk}-({\theta }_{n-1}-{\theta }_{n-1}^{hk}){\parallel }_Y^2+r{\parallel P_n-P_n^{hk}-(P_{n-1}-P_{n-1}^{hk})\parallel }_Y^2\hfill \\ +2d({\theta }_n-{\theta }_n^{hk}-({\theta }_{n-1}-{\theta }_{n-1}^{hk}),P_n-P_n^{hk}-(P_{n-1}-P_{n-1}^{hk}))\ge 0,\hfill \end{array}$$

thanks to condition $cr>d^2$, and combining estimates (26) and (27), we obtain

$$\begin{array}{ccc}& & \frac{c}{2k}\left\{\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2-{\parallel {\theta }_{n-1}-{\theta }_{n-1}^{hk}\parallel }_Y^2\right\}+\frac{r}{2k}\left\{\parallel P_n-P_n^{hk}{\parallel }_Y^2-{\parallel P_{n-1}-P_{n-1}^{hk}\parallel }_Y^2\right\}\hfill \\ & & +\frac{d}{k}\left\{(P_n-P_n^{hk},{\theta }_n-{\theta }_n^{hk})-(P_{n-1}-P_{n-1}^{hk},{\theta }_{n-1}-{\theta }_{n-1}^{hk})\right\}\hfill \\ & & -{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\theta }_n-{\theta }_n^{hk})-{\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},P_n-P_n^{hk})\hfill \\ & & \le C(\parallel {\theta }_{tn}-\delta {\theta }_n{\parallel }_Y^2+\parallel {\theta }_n-s^h{\parallel }_V^2+\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2+\parallel P_{tn}-\delta P_n{\parallel }_Y^2+{\parallel {\theta }_n-{\theta }_n^{hk}\parallel }_Y^2\hfill \\ & & +\parallel P_n-m^h{\parallel }_V^2+(\delta {\theta }_n-\delta {\theta }_n^{hk},{\theta }_n-s^h)+(\delta P_n-\delta P_n^{hk},{\theta }_n-s^h)+{\parallel P_n-P_n^{hk}\parallel }_Y^2\hfill \\ & & +(\delta P_n-\delta P_n^{hk},P_n-m^h)+(\delta {\theta }_n-\delta {\theta }_n^{hk},P_n-m^h)).\hfill \end{array}$$

Combining the above estimates with (24) and (25), keeping in mind that

$$\begin{array}{c}-{\gamma }_1({\zeta }_{nx}-{\zeta }_{nx}^{hk},{\theta }_n-{\theta }_n^{hk})-{\gamma }_2({\zeta }_{nx}-{\zeta }_{nx}^{hk},P_n-P_n^{hk})\hfill \\ ={\gamma }_1({\zeta }_n-{\zeta }_n^{hk},{\theta }_{nx}-{\theta }_{nx}^{hk})+{\gamma }_2({\zeta }_n-{\zeta }_n^{hk},P_{nx}-P_{nx}^{hk}),\hfill \end{array}$$

multiplying the resulting estimates by k and summing up to n, it follows that, for all ${\left\{w_j^h\right\}}_{j=0}^n,{\left\{r_j^h\right\}}_{j=0}^n,{\left\{s_j^h\right\}}_{j=0}^n,{\left\{m_j^h\right\}}_{j=0}^n\subset V^h$,

$$\begin{array}{c}{\displaystyle {\rho }_1\parallel {\xi }_n-{\xi }_n^{hk}{\parallel }_Y^2+\kappa \parallel {\phi }_{nx}-{\phi }_{nx}^{hk}{\parallel }_Y^2+{\rho }_2\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y^2+\alpha \parallel {\psi }_{nx}-{\psi }_{nx}^{hk}{\parallel }_Y^2+\kappa {\parallel {\psi }_n-{\psi }_n^{hk}\parallel }_Y^2}\hfill \\ +c\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2+r{\parallel P_n-P_n^{hk}\parallel }_Y^2+2d(P_n-P_n^{hk},{\theta }_n-{\theta }_n^{hk})\hfill \\ {\displaystyle \le Ck\sum _{j=1}^n(\parallel {\xi }_{tj}-\delta {\xi }_j{\parallel }_Y^2+\parallel {\phi }_{tj}-\delta {\phi }_j{\parallel }_V^2+\parallel {\xi }_j-w_j^h{\parallel }_V^2+\parallel {\psi }_{jx}-{\psi }_{jx}^{hk}{\parallel }_Y^2+{\parallel {\theta }_j-{\theta }_j^{hk}\parallel }_Y^2}\hfill \\ +\parallel {\xi }_j-{\xi }_j^{hk}{\parallel }_Y^2+\parallel {\phi }_{jx}-{\phi }_{jx}^{hk}{\parallel }_Y^2+(\delta {\xi }_j-\delta {\xi }_j^{hk},{\xi }_j-w_j^h)+\parallel {\zeta }_{tj}-\delta {\zeta }_j{\parallel }_Y^2+{\parallel P_j-P_j^{hk}\parallel }_Y^2\hfill \\ +\parallel {\psi }_{tj}-\delta {\psi }_j{\parallel }_V^2+\parallel {\zeta }_j-r_j^h{\parallel }_V^2+(\delta {\zeta }_j-\delta {\zeta }_j^{hk},{\zeta }_j-r_j^h)+\parallel {\theta }_{tj}-\delta {\theta }_j{\parallel }_Y^2+{\parallel {\theta }_j-s_j^h\parallel }_V^2\hfill \\ +\parallel P_{tj}-\delta P_j{\parallel }_Y^2+{\parallel P_j-m_j^h\parallel }_V^2+(\delta {\theta }_j-\delta {\theta }_j^{hk},{\theta }_j-s_j^h)+(\delta P_j-\delta P_j^{hk},{\theta }_j-s_j^h)\hfill \\ +(\delta {\theta }_j-\delta {\theta }_j^{hk},P_j-m_j^h)+(\delta P_j-\delta P_j^{hk},P_j-m_j^h))+C(\parallel {\xi }_0-{\xi }_0^h{\parallel }_Y^2+{\parallel {\phi }_0-{\phi }_0^h\parallel }_V^2\hfill \\ +\parallel {\zeta }_0-{\zeta }_0^h{\parallel }_Y^2+\parallel {\psi }_0-{\psi }_0^h{\parallel }_V^2+\parallel P_0-P_0^h{\parallel }_Y^2+{\parallel {\theta }_0-{\theta }_0^h\parallel }_Y^2).\hfill \end{array}$$

Now, again using condition $cr>d^2$, we easily find that

$$c\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2+r\parallel P_n-P_n^{hk}{\parallel }_Y^2+2d(P_n-P_n^{hk},{\theta }_n-{\theta }_n^{hk})\ge C(\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y^2+\parallel P_n-P_n^{hk}{\parallel }_Y^2).$$

Finally, keeping in mind that

$$\begin{array}{c}{\displaystyle k\sum _{j=1}^n(\delta {\xi }_j-\delta {\xi }_j^{hk},{\xi }_j-w_j^h)=({\xi }_n-{\xi }_n^{hk},{\xi }_n-w_n^h)+({\xi }_0^h-{\xi }_0,{\xi }_1-w_1^h)}\hfill \\ {\displaystyle +\sum _{j=1}^{n-1}({\xi }_j-{\xi }_j^{hk},{\xi }_j-w_j^h-({\xi }_{j+1}-w_{j+1}^h)),}\hfill \\ {\displaystyle k\sum _{j=1}^n(\delta {\zeta }_j-\delta {\zeta }_j^{hk},{\zeta }_j-r_j^h)=({\zeta }_n-{\zeta }_n^{hk},{\zeta }_n-r_n^h)+({\zeta }_0^h-{\zeta }_0,{\zeta }_1-r_1^h)}\hfill \\ {\displaystyle +\sum _{j=1}^{n-1}({\zeta }_j-{\zeta }_j^{hk},{\zeta }_j-r_j^h-({\zeta }_{j+1}-r_{j+1}^h)),}\hfill \\ {\displaystyle k\sum _{j=1}^n(\delta {\theta }_j-\delta {\theta }_j^{hk},{\theta }_j-s_j^h)=({\theta }_n-{\theta }_n^{hk},{\theta }_n-s_n^h)+({\theta }_0^h-{\theta }_0,{\theta }_1-s_1^h)}\hfill \\ {\displaystyle +\sum _{j=1}^{n-1}({\theta }_j-{\theta }_j^{hk},{\theta }_j-s_j^h-({\theta }_{j+1}-s_{j+1}^h)),}\hfill \end{array}$$

$$\begin{array}{c}{\displaystyle k\sum _{j=1}^n(\delta P_j-\delta P_j^{hk},P_j-m_j^h)=(P_n-P_n^{hk},P_n-m_n^h)+(P_0^h-P_0,P_1-m_1^h)}\hfill \\ {\displaystyle +\sum _{j=1}^{n-1}(P_j-P_j^{hk},P_j-m_j^h-(P_{j+1}-m_{j+1}^h)),}\hfill \\ {\displaystyle k\sum _{j=1}^n(\delta P_j-\delta P_j^{hk},{\theta }_j-s_j^h)=(P_n-P_n^{hk},{\theta }_n-s_n^h)+(P_0^h-P_0,{\theta }_1-s_1^h)}\hfill \\ {\displaystyle +\sum _{j=1}^{n-1}(P_j-P_j^{hk},{\theta }_j-s_j^h-({\theta }_{j+1}-s_{j+1}^h)),}\hfill \\ {\displaystyle k\sum _{j=1}^n(\delta {\theta }_j-\delta {\theta }_j^{hk},P_j-m_j^h)=({\theta }_n-{\theta }_n^{hk},P_n-m_n^h)+({\theta }_0^h-{\theta }_0,P_1-m_1^h)}\hfill \\ {\displaystyle +\sum _{j=1}^{n-1}({\theta }_j-{\theta }_j^{hk},P_j-m_j^h-(P_{j+1}-m_{j+1}^h)),}\hfill \end{array}$$

again using a discrete version of Gronwall’s inequality (see [27]), we obtain the desired a priori error estimates. □

The estimates presented in the theorem above can be utilized to determine the convergence order of the approximations provided by the discrete problem (14)–(18). Hence, provided we assume the following regularity:

$$\begin{array}{c}\phi ,\psi \in H^3(0,T;Y)\cap H^2(0,T;V)\cap C^1([0,T];H^2(0,\ell )),\hfill \\ \theta ,P\in H^2(0,T;Y)\cap H^1(0,T;V)\cap C([0,T];H^2(0,\ell )),\hfill \end{array}$$

(28)

by applying some results on the approximation by finite elements (refer to [26]) and utilizing previous estimates derived in [27], we can demonstrate the linear convergence of the algorithm. Specifically, we can establish the existence of a positive constant $C>0$, which remains independent of the discretization parameters h and k, such that

$$\begin{array}{c}{\displaystyle \underset{0\le n\le N}{max}\{\parallel {\xi }_n-{\xi }_n^{hk}{\parallel }_Y+\parallel {\phi }_n-{\phi }_n^{hk}{\parallel }_V+\parallel {\zeta }_n-{\zeta }_n^{hk}{\parallel }_Y+{\parallel {\psi }_n-{\psi }_n^{hk}\parallel }_V}\hfill \\ +\parallel {\theta }_n-{\theta }_n^{hk}{\parallel }_Y+{\parallel P_n-P_n^{hk}\parallel }_Y\}\le C(h+k).\hfill \end{array}$$

4. Numerical Results

In this section, we study some problems solved with the numerical approximation presented previously. We start by checking numerically that the convergence rate of the algorithm matches the one obtained theoretically. Then, in order to analyze the energy decay, we solve three variants of the model: the fully viscoelastic problem, as well as the partially viscoelastic I and partially viscoelastic II problems.

We note that the problem solved in [25] is slightly different, since the authors impose Neumann boundary conditions for $\psi$. With our definition of the model, their problem is not compatible with our boundary conditions. However, by taking into account their boundary conditions, we were able to replicate their results.

Finally, it is also worth noting that the numerical resolution of problem (14)–(18) was implemented in a 4-core 3.40 GHz computer with 16 Gb of RAM by using the well-known software MATLAB.

4.1. Error Convergence

In order to assess the convergence rate of the implementation (and also to check that it is correct), we solve a problem with a known analytical solution. We select the following parameters, which will be kept constant throughout all the numerical experiments unless specified otherwise:

$$\begin{array}{c}\ell =1;T=1;{\rho }_1=0.1;{\rho }_2=1;{\mu }_1=2;{\mu }_2=3;{\gamma }_1=2.6;{\gamma }_2=7.5;\hfill \\ \kappa =5;K=9;\alpha =7;c=10;d=2;r=2;H=8.\hfill \end{array}$$

To be able to obtain the analytical solution, we manufacture some artificial supply terms:

$$\begin{array}{c}{f}_{\phi }=\frac{-1}{10}x{\mathrm{e}}^t\left(x^5-993x^4+828x^3+1319x^2-1485x+330\right),\hfill \\ f_{\psi }=\frac{-3}{2}x{\mathrm{e}}^t\left(13x^5-3121x^4+7144x^3-4977x^2+821x+120\right),\hfill \\ f_{\theta }=\frac{3}{5}x{\mathrm{e}}^t\left(-500x^5+1578x^4+2805x^3-8344x^2+5361x-900\right),\hfill \\ f_P=\frac{5}{2}x{\mathrm{e}}^t\left(-88x^5+318x^4+9201x^3-19004x^2+11493x-1920\right).\hfill \end{array}$$

With these supply terms, the solution for the problem variables are, for all $x\in (0,1)$ and $t\in (0,1)$,

$$\begin{array}{c}\phi (x,t)=x^3{\mathrm{e}}^t{\left(1-x\right)}^3,\psi (x,t)=3x^3{\mathrm{e}}^t{\left(1-x\right)}^3,\hfill \\ \theta (x,t)=10x^3{\mathrm{e}}^t{\left(1-x\right)}^3,P(x,t)=100x^3{\mathrm{e}}^t{\left(1-x\right)}^3.\hfill \end{array}$$

The boundary conditions are homogeneous Dirichlet conditions for all the variables, and the initial conditions are obtained by substituting $t=0$ in the previous expressions.

Under these conditions, we compute the numerical errors for different mesh sizes and timesteps. The results are compiled in

Table 1. Numerical errors for some values of h and k.

$\mathit{h}\downarrow \mathit{k}\to$	$5\times {10}^{-2}$	$1\times {10}^{-2}$	$5\times {10}^{-3}$	$1\times {10}^{-3}$	$5\times {10}^{-4}$	$1\times {10}^{-4}$	$5\times {10}^{-5}$	$1\times {10}^{-5}$
$1\times {10}^{-1}$	0.156721	0.157573	0.157723	0.157851	0.157867	0.157887	0.157882	0.157883
$5\times {10}^{-2}$	0.057486	0.056802	0.056903	0.057011	0.057026	0.057038	0.057040	0.057041
$2\times {10}^{-2}$	0.021111	0.018167	0.017962	0.017963	0.017974	0.017985	0.017987	0.017988
$1\times {10}^{-2}$	0.012934	0.008800	0.008431	0.008225	0.008225	0.008233	0.008234	0.008235
$5\times {10}^{-3}$	0.009643	0.004798	0.004291	0.003966	0.003938	0.003931	0.003931	0.003932
$2\times {10}^{-3}$	0.008115	0.002616	0.002032	0.001604	0.001559	0.001531	0.001529	0.001529
$1\times {10}^{-3}$	0.007801	0.001953	0.001311	0.000853	0.000801	0.000763	0.000759	0.000757
$2\times {10}^{-4}$	0.007684	0.001578	0.000822	0.000262	0.000203	0.000160	0.000155	0.000151
$1\times {10}^{-4}$	0.007680	0.001560	0.000790	0.000195	0.000131	$8.55\times {10}^{-5}$	$8.02\times {10}^{-5}$	$7.60\times {10}^{-5}$

, and the main diagonal of this table is plotted in Figure 1. These results agree with Theorem 2, and they show the linear convergence of the algorithm.

4.2. Energy Decay

To numerically check the decay rates for the energy, we performed some experiments with the parameters, boundary conditions and initial conditions listed in the previous section. However, in this case, we did not have any supply term. In all the cases, we refined the mesh and the timestep until the solution stopped changing. This resulted in 51 elements ($h=0.02$) and a timestep of $k={10}^{-4}$ (except in the partially viscoelastic I case; see below).

The results for the fully viscoelastic case (${\mu }_1=2,{\mu }_2=3$) are plotted in Figure 2, both in natural (left) and semilogarithmic (right) axes. The straight line in the semilogarithmic plot (after a brief initial transient) shows that the energy decay is exponential.

When ${\mu }_1=0$ and ${\mu }_2=3$, the problem becomes partially viscoelastic (of type I). In this case, much finer timesteps are required to obtain the solution. This effect is clearly seen in Figure 3, where we represent the solution obtained for different timesteps (again in natural and semilogarithmic axes). A timestep of $k={10}^{-6}$ is required to obtain a converged solution. Also, as the theoretical results predict, the decay rate of the energy is not exponential (it is polynomial).

Finally, in the partially viscoelastic II problem (${\mu }_2=0$ and ${\mu }_1=2$) shown in Figure 4, we can see that we recover the exponential decay; however, since the dissipation mechanism is weaker, some oscillations remain the whole time.

5. Conclusions

In this work, a viscoelastic Timoshenko beam problem is analyzed from the numerical point of view. The thermodiffusion effects were included into the model, continuing the recently published paper [25], where the existence and uniqueness of solution and the energy decay were studied from the analytical point of view, and some numerical results were described by using an spectral method.

A finite element approximation was also derived by using the implicit Euler scheme to discretize the time derivatives, and a discrete stability property and an a priori error analysis were proved. From these results, the linear convergence of the approximation was concluded under some additional regularity conditions.

Finally, some numerical simulations were performed to demonstrate the predicted theoretical convergence in an academical example, and the discrete energy behavior depending on the dissipative mechanisms (the viscosity terms) was studied. As a conclusion, it was shown numerically that if the model was fully viscoelastic or partially viscoelastic of type II, then the decay was exponential, but if it was partially viscoelastic of type I (${\mu }_1=0$ and ${\mu }_2\ne 0$), then the decay was polynomial. These results are compared directly with those shown in [25].

As a possible future work, we could consider the extension to the multidimensional setting. From our opinion, the theoretical analysis would be rather similar, but the numerical implementation could lead to some problems of CPU time, because the case corresponding to partial type I viscoelasticity needs very refined timesteps.