Analysis of two thermoelastic problems with the Green–Lindsay model

Abstract

In this paper, we analyze, from the numerical point of view, two thermo-elastic problems involving the Green–Lindsay theory. The coupling term is different for each case, involving second order or first order spatial derivatives, respectively. The variational formulation leads to a linear coupled system which is written in terms of the velocity and temperature speed. An existence and uniqueness results and the exponential energy decay for the problem with the stronger coupling are recalled. The polynomial energy decay for the weaker coupling is then proved but using the theory of linear semigroups. Then, a fully discrete approximation is introduced using the finite element method and an implicit scheme. A discrete stability property and a main a priori error estimates result are shown, from which we can derive the linear convergence of the approximations. Finally, some numerical simulations are presented to demonstrate the accuracy of the algorithm, the discrete energy decay and the dependence on the relaxation parameter.

1 Introduction

In this work, we consider two thermoelastic models based on the Green–Lindsay theory (Green and Lindsay 1972). Therefore, it is worth recalling that these authors proposed a theory where the heat conduction was different than the parabolic version. In fact, since that parabolic version allows the propagation of instantaneous thermal waves (and so, it violates the so-called “causality principle”), several authors have proposed different theories trying to overcome this difficulty. Among the pioneering researchers we can cite the above commented proposal by Green and Lindsay (1972) that, even if it is based on Fourier law, we have an energy equation which brings us to a damped hyperbolic equation to describe the behavior of the temperature. There exists a huge number of contributions dealing with this theory [see the numerous references in the papers by Hetnarski and Ignaczak (1999, 2000) and the books (Ignaczak and Ostoja-Starzewski 2009; Straughan 2011)]. It is well-known that, when we consider the usual theory of thermoelasticity of Green–Lindsay type, we conclude that the solutions decay through the time and, for the one-dimensional case, this decay is of exponential type. Moreover, if we consider the strain gradient elasticity (or the thermoelastic plate) with the heat conduction of parabolic type (Avalishvili et al. 2018; Aouadi et al. 2019; Aouadi and Moulahi 2015; Fernández-Sare et al. 2010; Liu and Renardy 1995; Liu and Quintanilla 2010; Lindsay and Straughan 1979; Liu et al. 2022; Shivay and Mukhopadhyay 2020), we find that the decay is also exponential in the one-dimensional case or in higher dimensions for the plate. It led recently to the question if it would also happen similarly for the heat conduction of the Green–Lindsay type. When we consider the systems of equations governing both problems we can appreciate a clear similarity, being its main difference the coupling term. In the case of the strain-gradient thermoelasticity, such coupling is weaker, obtaining that, for this case, the energy decay is slow (Quintanilla et al. 2023) and, even more, we prove that, in fact, it is polynomial (see the appendix provided in this work). For the other case, since the coupling is stronger, the energy decay is exponential (even for dimensions greater than one¹). In this paper, our aim is to develop a similar study from the numerical point of view, which constitutes one of the novelties of this work. We also think that it is relevant to understand the difference between the behavior of the solutions to these problems with different coupling mechanisms and so, we will study it numerically. Indeed, even if the discrete energy decay is always exponential, in this work we will show that it is much faster when the coupling mechanism is stronger and it leads to an exponential energy decay for the continuous problem.

The plan of this paper is as follows. The model equations for the two problems and the assumptions required on the constitutive parameters are presented in Sect. 2. The existence of a unique solution and the energy decay property, recently proved in Quintanilla et al. (2023), is recalled. Then, by using the classical finite element method and the implicit Euler scheme, a fully discrete approximation is introduced in Sect. 3. A discrete stability property and a priori error estimates are proved, and the linear convergence is derived under some additional regularity conditions on the continuous solution. Finally, some numerical simulations are presented in Sect. 4 to demonstrate the numerical convergence, the behavior of the discrete energy decay and the dependence on the coupling parameter.

2 The thermoelastic model

In this section, we describe the thermomechanical problems that we will study in the paper. We refer the reader to the recently published work (Quintanilla et al. 2023) for further details. Here, we restrict, for the sake of simplicity in the numerical simulations, to the one-dimensional case, although the extension of the numerical analysis, presented in the next section, to the multi-dimensional setting is not difficult. Therefore, we will consider the following linear systems:

$$\begin{aligned}{} & {} \left. \begin{array}{l} \rho u_{tt}+\mu u_{xxxx}-a(\theta _{xx}+\alpha \theta _{txx})=0,\\ \overline{h}\theta _{tt}+d\theta _t-\kappa \theta _{xx}+a u_{txx}=0, \end{array}\right\} \end{aligned}$$

(1)

$$\begin{aligned}{} & {} \left. \begin{array}{l} \rho u_{tt}+\mu u_{xxxx}-bu_{xx}-a(\theta _{x}+\alpha \theta _{tx}),\\ \overline{h}\theta _{tt}+d\theta _t-\kappa \theta _{xx}-a u_{tx}=0. \end{array}\right\} \end{aligned}$$

(2)

We note that the unique difference is that, in the coupling term of the first equation of both systems, there is a stronger coupling in the first system. Here, u and $\theta$ denote the displacement and the temperature deviation in a rod occupying the interval $(0,\ell )$, $\ell >0$, $\rho$ is the mass density, $\kappa$ is the thermal conductivity, d is the thermal capacity, $\mu$ is the elastic modulus in system (1) but the hyperelastic modulus in system (2) and b is the elastic modulus in system (2). Moreover, $\alpha$ is a relaxation parameter usually used in the Green–Lindsay theory which satisfies the condition $\alpha d-\overline{h}>0$, and we also assume that all parameters $\rho ,\, \mu ,\, |a|,\, \alpha ,\, \overline{h},\,d,\,\kappa$ and b are positive.

We will study the deformation of the rod over a finite time interval (0, T), with $T>0$ being the final time.

To obtain a well-posed problem, we need to impose boundary and initial conditions to the above systems. Therefore, we assume the following Dirichlet-type boundary conditions for the displacements and Neumann-type boundary conditions for the temperature, that is,

$$\begin{aligned} \begin{array}{l} u(0,t)=u(\ell ,t)=u_x(0,t)=u_x(\ell ,t)=0 \quad \hbox {for a.e. } t\in (0,T),\\ \theta _x(0,t)=\theta _x(\ell ,t)=0 \quad \hbox {for a.e. } t\in (0,T), \end{array} \end{aligned}$$

(3)

and the initial conditions, for a.e. $x\in (0,\ell )$,

$$\begin{aligned} u(x,0)=u^0(x),\quad u_t(x,0)=v^0(x),\quad \theta (x,0)=\theta ^0(x),\quad \theta _t(x,0)=\vartheta ^0(x). \end{aligned}$$

(4)

In the above conditions, functions $u^0$, $v^0$, $\theta ^0$ and $\vartheta ^0$ represent the initial data of the problems. Moreover, to guarantee the uniqueness of solutions to the problems we need to assume that functions $\theta ^0$ and $\vartheta ^0$ have null average.

Hence, we will consider two similar problems: the first one defined by system (1) with boundary conditions (3) and initial conditions (4), and the second one defined by system (2) with boundary conditions (3) and initial conditions (4). Since no doubt the numerical analysis of both problems is rather similar, we will focus on the first problem since it has some additional difficulties due to the stronger coupling between the equations of the system.

The existence and uniqueness, as well as the exponential decay, were recently studied in Quintanilla et al. (2023) and they are summarized as follows.

Theorem 1

Under the condition $\alpha d-\overline{h}>0$, then problem (1), (3) and (4) admits a unique solution

$$\begin{aligned} \begin{array}{l} u\in C^1([0,T];H^2_0(0,\ell ))\cap C^2([0,T];L^2(0,\ell )),\\ \theta \in { C([0,T];H^2_0(0,\ell ))\cap }C^1([0,T];H^1(0,\ell ))\cap C^2([0,T];L^2(0,\ell )) \end{array} \end{aligned}$$

which is exponentially stable. In a similar way, we also conclude that problem (2)–(4) admits a unique solution with the same regularity but, in this case, the energy decay is slow (not exponential at least).

The proof of the above theorem is shown in the recent paper (Quintanilla et al. 2023) and it is based on the use of the theory of linear semigroups, adequate energy functionals and some a priori estimates. However, following the work done in Bazarra et al. (2022) for another MGT problem, in fact we could prove that this energy decay is polynomial (see the Appendix A for details).

3 Fully discrete approximations and an a priori error analysis

In this section, we present an a priori error analysis of a fully discrete scheme for solving problem (1), (3) and (4). As we commented above, the second problem described in the previous section can be analyzed analogously.

First, we need to provide a variational formulation of this problem. Therefore, let us denote $Y=L^2(0,\ell )$ and denote by $(\cdot ,\cdot )$ and $\Vert \cdot \Vert$ the inner product and the norm in this space, respectively. Moreover, let $E=H^1(0,\ell )$ and $V=H^2_0(0,\ell )$.

Integrating by parts system (1) and using boundary conditions (3) we obtain the following weak form of the thermomechanical problem.

Find the velocity $v:[0,T]\rightarrow V$ and the temperature speed $\vartheta :[0,T]\rightarrow E$ such that $v(0)=v^0$, $\vartheta (0)=\vartheta ^0$, and for a.e. $t\in (0,T)$ and $w\in V,\, \xi \in E$,

$$\begin{aligned}{} & {} \displaystyle \rho (v_t(t),w)+\mu (u_{xx}(t),w_{xx})-a(\theta (t)+\alpha \theta _t(t),w_{xx})=0, \end{aligned}$$

(5)

$$\begin{aligned}{} & {} { (\overline{h}\vartheta _t(t)+d\vartheta (t),\xi )}+\kappa (\theta _x(t),\xi _x)+a(v_{xx}(t),\xi )=0, \end{aligned}$$

(6)

where the displacements u(t) and the temperature $\theta (t)$ are recovered from the relations:

$$\begin{aligned} \displaystyle \displaystyle u(t)=\int _0^t v(s)\,\textrm{d}s+u^0,\quad \theta (t)=\int _0^t\vartheta (s)\,\textrm{d}s+\theta ^0. \end{aligned}$$

(7)

Remark 1

We note that the variational formulation of the problem defined by system (2) is rather similar. The unique difference is that we must replace variational equations (5) and (6) by the following ones:

$$\begin{aligned}{} & {} \rho (v_t(t),w)+\mu (u_{xx}(t),w_{xx}){+b(u_{x}(t),w_x)}+a(\theta (t)+\alpha \theta _t(t),w_{x})=0, \end{aligned}$$

(8)

$$\begin{aligned}{} & {} {(\overline{h}\vartheta _t(t)+d\vartheta (t),\xi )}+\kappa (\theta _x(t),\xi _x)-a(v_{x}(t),\xi )=0. \end{aligned}$$

(9)

Now, we will obtain the approximation of problem (5)–(7). We will proceed in two steps. First, we approximate the problem in space. Thus, let us assume that the interval $[0,\ell ]$ is divided into M subintervals of length $h=\ell /M$, with nodes $a_0=0<\cdots <a_M=\ell$, and construct the finite element spaces:

(10)

where $P_r([a_i,a_{i+1}])$ represents the space of polynomials of degree r in the subinterval $[a_i,a_{i+1}]$.

The discrete initial conditions $u^{0h}$, $v^{0h}$, $\theta ^{0h}$ and $\vartheta ^{0h}$ are approximations of the respective initial conditions $u^{0}$, $v^{0}$, $\theta ^{0}$ and $\vartheta ^{0}$ defined as

$$\begin{aligned} u^{0h}=\mathcal {P}_1^{h} u^0,\quad v^{0h}=\mathcal {P}_1^{h} v^0, \quad \theta ^{0h}=\mathcal {P}_2^{h} \theta ^0,\quad \vartheta ^{0h}=\mathcal {P}_2^{h} \vartheta ^0. \end{aligned}$$

(11)

Here, we denote by $\mathcal {P}_1^{h}$ and $\mathcal {P}_2^{h}$ the interpolation operators over the finite element spaces $V^h$ and $E^h$, respectively (see Ciarlet 1991 for details).

Second, to provide the time discretization, we consider a uniform partition of the time interval [0, T], denoted by $0=t_0<t_1<\ldots <t_N=T$, where $k=T/N$ is the time step size. If f is a continuous function, we denote $f_n=f(t_n)$ and, for the sequence $\{z_n\}_{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 we can introduce the following fully discrete problem.

Find the discrete velocity $v^{hk}=\{v_n^{hk}\}_{n=0}^N\subset V^h$ and the discrete temperature speed $\vartheta ^{hk}=\{\vartheta _n^{hk}\}_{n=0}^N\subset E^h$ such that $v_0^{hk}=v^{0h}$, $\vartheta _0^{hk}=\vartheta ^{0h}$, and for all $n=1,\ldots ,N$ and $w^h\in V^h,\, \xi ^h\in E^h$,

$$\begin{aligned}{} & {} \displaystyle \rho (\delta {v}_n^{hk},w^h)+\mu ((u_n^{hk})_{xx},w^h_{xx})-a(\theta _n^{hk}+\alpha \vartheta _n^{hk},w^h_{xx})=0, \end{aligned}$$

(12)

$$\begin{aligned}{} & {} (\overline{h}\delta \vartheta _n^{hk}+d\vartheta _n^{hk},\xi ^h)+\kappa ((\theta _n^{hk})_x,\xi ^h_x)+a((v_n^{hk})_{xx},\xi ^h)=0, \end{aligned}$$

(13)

where the discrete displacements $u_n^{hk}$ and the discrete temperature $\theta _n^{hk}$ are updated from the relations:

$$\begin{aligned} \displaystyle u_n^{hk}=k\sum _{j=1}^n v_j^{hk}+u^{0h},\quad \theta _n^{hk}=k\sum _{j=1}^n \vartheta _j^{hk}+\theta ^{0h}. \end{aligned}$$

(14)

Using the assumptions on the constitutive coefficients and applying Lax-Milgram lemma, it is easy to prove that the above discrete problem has a unique solution.

Remark 2

Again, the numerical approximation of the problem obtained by using system (2) is similar to the previous one. The unique difference is that we must replace the discrete variational equations (12) and (13) by the following ones (as we also did in (8), (9)):

$$\begin{aligned}{} & {} \rho (\delta {v}_n^{hk},w^h)+\mu ((u_n^{hk})_{xx},w^h_{xx})+{ b((u_n^{hk})_x,w^h_x)}+a(\theta _n^{hk}+\alpha \vartheta _n^{hk},w^h_{x})=0, \end{aligned}$$

(15)

$$\begin{aligned}{} & {} (\overline{h}\delta \vartheta _n^{hk}+d\vartheta _n^{hk},\xi ^h)+\kappa ((\theta _n^{hk})_x,\xi ^h_x)-a((v_n^{hk})_{x},\xi ^h)=0. \end{aligned}$$

(16)

Now, we will provide a discrete stability property that we summarize in the following.

Lemma 1

Under the assumptions of Theorem 1, the sequences $\{u^{hk}, v^{hk},\theta ^{hk},\vartheta ^{hk}\}$, generated by discrete problem (12)–(14), satisfy the stability estimate:

$$\begin{aligned} \begin{array}{l} \displaystyle \Vert v_n^{hk}\Vert ^2+ \Vert (u_n^{hk})_{xx}\Vert ^2+ \Vert \vartheta _n^{hk}\Vert ^2+ \Vert (\theta _n^{hk})_x\Vert ^2\le C,\end{array} \end{aligned}$$

where C is a positive constant assumed to be independent of the discretization parameters h and k.

Proof

For the sake of clarity in the calculations, we assume that $\alpha =1$. It is straightforward to extend this result to the general case.

Taking as a test function $w^h=v_n^{hk}$ in discrete variational equation (12) we obtain

$$\begin{aligned} \rho (\delta {v}_n^{hk},v_n^{hk})+\mu ((u_n^{hk})_{xx},(v_n^{hk})_{xx})-a(\theta _n^{hk}+\vartheta _n^{hk},(v_n^{hk})_{xx})=0. \end{aligned}$$

Using the estimates

$$\begin{aligned} \begin{array}{l} \displaystyle \rho (\delta {v_n^{hk}},v_n^{hk})\ge \frac{\rho }{2k}\left\{ \Vert v_n^{hk}\Vert ^2-\Vert v_{n-1}^{hk}\Vert \right\} ,\\[0.2cm] \displaystyle \mu ((u_n^{hk})_{xx}, (v_n^{hk})_{xx})\ge \frac{\mu }{2k}\left\{ \Vert (u_n^{hk})_{xx}\Vert ^2-\Vert (u_{n-1}^{hk})_{xx}\Vert ^2\right\} , \end{array} \end{aligned}$$

we find that

$$\begin{aligned}{} & {} \displaystyle \frac{\rho }{2k}\left\{ \Vert v_n^{hk}\Vert ^2-\Vert v_{n-1}^{hk}\Vert \right\} +\frac{\mu }{2k}\left\{ \Vert (u_n^{hk})_{xx}\Vert ^2-\Vert (u_{n-1}^{hk})_{xx}\Vert ^2\right\} \nonumber \\{} & {} \quad -a(\theta _n^{hk}+\vartheta _n^{hk},(v_n^{hk})_{xx})\le 0. \end{aligned}$$

(17)

Now, taking the discrete variational equation (13) for a test function $\xi ^h=\vartheta _n^{hk}$ we have

$$\begin{aligned} (\overline{h}\delta \vartheta _n^{hk}+d\vartheta _n^{hk},\vartheta _n^{hk})+\kappa ((\theta _n^{hk})_x,(\vartheta _n^{hk})_x)+a((v_n^{hk})_{xx},\vartheta _n^{hk})=0. \end{aligned}$$

Keeping in mind that

$$\begin{aligned} \begin{array}{l} \displaystyle (\overline{h}\delta {\vartheta _n^{hk}},\vartheta _n^{hk})\ge \frac{\overline{h}}{2k}\left\{ \Vert \vartheta _n^{hk}\Vert ^2-\Vert \vartheta _{n-1}^{hk}\Vert \right\} ,\\ \displaystyle \kappa ((\theta _n^{hk})_{x}, (\vartheta _n^{hk})_{x})\ge \frac{\kappa }{2k}\left\{ \Vert (\theta _n^{hk})_{x}\Vert ^2-\Vert (\theta _{n-1}^{hk})_{x}\Vert ^2\right\} , \end{array} \end{aligned}$$

it follows that

$$\begin{aligned}{} & {} \displaystyle \frac{\overline{h}}{2k}\left\{ \Vert \vartheta _n^{hk}\Vert ^2-\Vert \vartheta _{n-1}^{hk}\Vert \right\} +\frac{\kappa }{2k}\left\{ \Vert (\theta _n^{hk})_{x}\Vert ^2-\Vert (\theta _{n-1}^{hk})_{x}\Vert ^2\right\} \nonumber \\{} & {} \quad +a((v_n^{hk})_{xx},\vartheta _n^{hk})\le 0. \end{aligned}$$

(18)

Combining estimates (17) and (18) we have

$$\begin{aligned}{} & {} \displaystyle \frac{\rho }{2k}\left\{ \Vert v_n^{hk}\Vert ^2-\Vert v_{n-1}^{hk}\Vert \right\} +\frac{\mu }{2k}\left\{ \Vert (u_n^{hk})_{xx}\Vert ^2-\Vert (u_{n-1}^{hk})_{xx}\Vert ^2\right\} -a(\theta _n^{hk},(v_n^{hk})_{xx})\\{} & {} \quad + \displaystyle \frac{\overline{h}}{2k}\left\{ \Vert \vartheta _n^{hk}\Vert ^2-\Vert \vartheta _{n-1}^{hk}\Vert \right\} +\frac{\kappa }{2k}\left\{ \Vert (\theta _n^{hk})_{x}\Vert ^2-\Vert (\theta _{n-1}^{hk})_{x}\Vert ^2\right\} \le 0. \end{aligned}$$

Multiplying the above estimates by k and summing up to n it leads

$$\begin{aligned} \begin{array}{l} \displaystyle \Vert v_n^{hk}\Vert ^2+ \Vert (u_n^{hk})_{xx}\Vert ^2+ \Vert \vartheta _n^{hk}\Vert ^2+ \Vert (\theta _n^{hk})_{x}\Vert ^2-ak\sum _{j=1}^n(\theta _j^{hk},(v_j^{hk})_{xx})\\ \qquad \le C \left( \Vert v^{0h}\Vert ^2+\Vert (u^{0h})_{xx}\Vert ^2+\Vert \vartheta ^{0h}\Vert ^2+\Vert (\theta ^{0h})_x\Vert ^2\right) . \end{array} \end{aligned}$$

Finally, keeping in mind that

$$\begin{aligned} \begin{array}{l} \displaystyle k\sum _{j=1}^n(\theta _j^{hk},(v_j^{hk})_{xx})=\sum _{j=1}^n(\theta _j^{hk},(u_j^{hk})_{xx}-(u_{j-1}^{hk})_{xx})\\ \qquad \displaystyle = \sum _{j=1}^{n-1}(\theta _j^{hk}-\theta _{j+1}^{hk},(u_j^{hk})_{xx})+(\theta _n^{hk},(u_n^{hk})_{xx}) -(\theta _1^{hk},(u^{0h})_{xx})\\ \qquad \displaystyle = -k\sum _{j=1}^{n-1}(\vartheta _{j}^{hk},(u_j^{hk})_{xx})+(\theta _n^{hk},(u_n^{hk})_{xx}) -(\theta _1^{hk},(u^{0h})_{xx}),\\ \displaystyle \Vert \theta _n^{hk}\Vert ^2\le C\left( \Vert \theta ^{0h}\Vert ^2+k\sum _{j=1}^n \Vert \vartheta ^{hk}_j\Vert ^2\right) , \end{array} \end{aligned}$$

applying several times Cauchy–Schwarz inequality and Cauchy’s inequality $ab\le \epsilon a^2+\frac{1}{4\epsilon }b^2$, $a,b,\epsilon \in \mathbb {R}$ with $\epsilon >0$, and using a discrete version of Gronwall’s inequality (see Campo et al. (2006)), we obtain the discrete stability. $\square$

In the rest of the section, we will obtain some a priori error estimates on the numerical errors $v_n-v_n^{hk}$ and $\vartheta _n-\vartheta _n^{hk}$, which we state in the following result.

Theorem 2

Let the assumptions of Theorem 1 hold. If we denote by $(u,v,\theta ,\vartheta )$ the solution to problem (5)–(7)s and by $(u^{hk},v^{hk},\theta ^{hk},\vartheta ^{hk})$ the solution to problem (12)-(14), then we have the following a priori error estimates, for all $\{w_n^{h}\}_{n=0}^N\subset V^h,\, \{\xi _n^{h}\}_{n=0}^N\subset E^h$,

$$\begin{aligned}{} & {} \displaystyle \max _{0\le n\le N}\Big \{\Vert v_n-v_n^{hk}\Vert ^2+\Vert (u_n-u_n^{hk})_{xx}\Vert ^2 +\Vert (\theta _n-\theta _n^{hk})_x\Vert ^2+\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2\Big \}\nonumber \\{} & {} \le Ck\sum _{j=1}^N\Big [\Vert v_{t}(t_j)-\delta v_j\Vert ^2+\Vert u_{t}(t_j)-\delta u_j\Vert ^2_V+\Vert v_j-w_j^{h}\Vert _{V}^2+\Vert \vartheta _{t}(t_j)-\delta \vartheta _j\Vert ^2\nonumber \\{} & {} \quad +\Vert \theta _{t}(t_j)-\delta \theta _j\Vert _{E}^2+\Vert \vartheta _j-\xi _j^{h}\Vert _{E}^2\Big ]\nonumber \\{} & {} \quad +\frac{C}{k}\sum _{j=1}^{N-1}\Big [\Vert v_j-w_j^h-(v_{j+1}-w_{j+1}^h)\Vert ^2+ \Vert \vartheta _j-\xi _j^h-(\vartheta _{j+1}-\xi _{j+1}^h)\Vert ^2+I_j\Big ]\nonumber \\{} & {} \quad +C\max _{0\le n\le N}\Vert v_n-w_n^h\Vert ^2+C\max _{0\le n\le N}\Vert \vartheta _n-\xi _n^h\Vert ^2\nonumber \\{} & {} \quad +C\Big (\Vert v^0-v^{0h}\Vert ^2+\Vert u^0-u^{0h}\Vert _{V}^2+\Vert \theta ^0-\theta ^{0h}\Vert _{E}^2 +\Vert \vartheta ^0-\vartheta ^{0h}\Vert ^2\Big ), \end{aligned}$$

(19)

where C is a positive constant which does not depend on parameters h and k, $I_j$ is the integration error given by

$$\begin{aligned} I_j=\left\| \int _0^{t_j}\vartheta (s)\,\textrm{d}s-k\sum _{j=1}^l\vartheta _l\right\| ^2 \end{aligned}$$

(20)

and, for a Hilbert space X, let $\Vert \cdot \Vert _X$ denote the norm in X.

Proof

In this proof, to simplify the calculations, we will consider again that $\alpha =1$. We note that we can modify the arguments used below to the general case with some minor changes.

First, we obtain the error estimates on the velocity $v_n-v_n^{hk}$. Thus, we subtract variational equation (5) for a test function $w=w^h\in V^h\subset V$, at time $t=t_n$, and discrete variational equation (12) to find that, for all $w^h\in V^h$,

$$\begin{aligned} \displaystyle \rho (v_{t}(t_n)-\delta {v}_n^{hk},w^h)+\mu ((u_n-u_n^{hk})_{xx},w^h_{xx})-a(\theta _n-\theta _n^{hk}+\vartheta _n-\vartheta _n^{hk},w^h_{xx})=0. \end{aligned}$$

Then, we find that, for all $w^h\in V^h$,

$$\begin{aligned}{} & {} \displaystyle \rho (v_{t}(t_n)-\delta {v}_n^{hk},v_n-v_n^{hk})+\mu ((u_n-u_n^{hk})_{xx},(v_n-v_n^{hk})_{xx})\\{} & {} \qquad -a(\theta _n-\theta _n^{hk}+\vartheta _n-\vartheta _n^{hk},(v_n-v_n^{hk})_{xx})\\{} & {} \quad \displaystyle =\rho (v_{t}(t_n)-\delta {v}_n^{hk},v_n-w^h)+\mu ((u_n-u_n^{hk})_{xx},(v_n-w^h)_{xx})\\{} & {} \qquad -a(\theta _n-\theta _n^{hk}+\vartheta _n-\vartheta _n^{hk},(v_n-w^h)_{xx}). \end{aligned}$$

Taking into account that

$$\begin{aligned}{} & {} \displaystyle \left( \delta v_n-\delta v_n^{hk},v_n-v_n^{hk}\right) \ge \frac{1}{2k}\Big \{\Vert v_n-v_n^{hk}\Vert ^2-\Vert v_{n-1}-v_{n-1}^{hk}\Vert ^2\Big \},\\{} & {} \displaystyle ((u_n-u_n^{hk})_{xx},(v_n-v_n^{hk})_{xx})= ((u_n-u_n^{hk})_{xx},(u_{t}(t_n)-\delta u_n)_{xx})\\{} & {} \quad + ((u_n-u_n^{hk})_{xx},(\delta u_n-\delta u_n^{hk})_{xx}),\\{} & {} \displaystyle \left( (u_n-u_n^{hk})_{xx},(\delta u_n-u_n^{hk})_{xx}\right) \ge \frac{1}{2k}\Big \{\Vert (u_n-u_n^{hk})_{xx}\Vert ^2-\Vert (u_{n-1}-u_{n-1}^{hk})_{xx}\Vert ^2\Big \}, \end{aligned}$$

where we have used the notations $\delta v_n=(v_n-v_{n-1})/k$ and $\delta u_n=(u_n-u_{n-1})/k$, applying several times Cauchy–Schwarz inequality and the previously commented Cauchy’s inequality, we obtain the following error estimates for the velocity field, for all $w^h\in V^h$,

$$\begin{aligned}{} & {} \displaystyle \frac{\rho }{2k}\Big \{\Vert v_n-v_n^{hk}\Vert ^2-\Vert v_{n-1}-v_{n-1}^{hk}\Vert ^2\Big \}\nonumber \\{} & {} \qquad -a(\theta _n-\theta _n^{hk}+\vartheta _n-\vartheta _n^{hk},(v_n-v_n^{hk})_{xx}) \nonumber \\{} & {} \qquad +\frac{\mu }{2k}\Big \{\Vert (u_n-u_n^{hk})_{xx}\Vert ^2-\Vert (u_{n-1}-u_{n-1}^{hk})_{xx}\Vert ^2\Big \}\nonumber \\{} & {} \quad \displaystyle \le C\Big ( \Vert v_{t}(t_n)-\delta v_n\Vert ^2+\Vert v_n-v_n^{hk}\Vert ^2+\Vert v_n-w^h\Vert _V^2+\Vert (u_n-u_n^{hk})_{xx}\Vert ^2\nonumber \\{} & {} \qquad +\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2+\Vert \theta _n-\theta _n^{hk}\Vert ^2+\Vert u_{t}(t_n)-\delta u_n\Vert _V^2\nonumber \\{} & {} \qquad + (\delta {v}_n-\delta {v^{hk}_n},v_n-w^h)\Big ), \end{aligned}$$

(21)

where, here and in what follows, C will represent a positive constant which depends on the constitutive coefficients, but it does not depend on the discretization parameters h and k, and whose value may change even within the same line.

Secondly, we derive the a priori error estimates on the temperature speed $\vartheta _n-\vartheta _n^{hk}$. Subtracting variational equation (6), for a test function $\xi =\xi ^h\in E^h\subset E$ at time $t=t_n$, and discrete variational equation (13), we have, for all $\xi ^h\in E^h$,

$$\begin{aligned} \begin{array}{l} \displaystyle (\overline{h}\vartheta _{t}(t_n)-\delta \vartheta _n^{hk}+d(\vartheta _n-\vartheta _n^{hk}),\xi ^h)+\kappa ((\theta _n-\theta _n^{hk})_x,\xi ^h_x)+a((v_n-v_n^{hk})_{xx},\xi ^h)=0, \end{array} \end{aligned}$$

and so we obtain, for all $\xi ^h\in E^h$,

$$\begin{aligned}{} & {} \displaystyle (\overline{h}\vartheta _{t}(t_n)-\delta \vartheta _n^{hk}+d(\vartheta _n-\vartheta _n^{hk}),\vartheta _n-\vartheta _n^{hk})+ \kappa ((\theta _n-\theta _n^{hk})_x,(\vartheta _n-\vartheta _n^{hk})_x)\\{} & {} \qquad +a((v_n-v_n^{hk})_{xx},\vartheta _n-\vartheta _n^{hk})\\{} & {} \quad \displaystyle = (\overline{h}\vartheta _{t}(t_n)-\delta \vartheta _n^{hk}+d(\vartheta _n-\vartheta _n^{hk}),\vartheta _n-\xi ^h)+ \kappa ((\theta _n-\theta _n^{hk})_x,(\vartheta _n-\xi ^{h})_x)\\{} & {} \qquad +a((v_n-v_n^{hk})_{xx},\vartheta _n-\xi ^{h}). \end{aligned}$$

Keeping in mind that

$$\begin{aligned}{} & {} \displaystyle \left( \overline{h}\delta \vartheta _n-\delta \vartheta _n^{hk},\vartheta _n-\vartheta _n^{hk}\right) _\ge \frac{\overline{h}}{2k}\Big \{\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2-\Vert \vartheta _{n-1}-\vartheta _{n-1}^{hk}\Vert ^2\Big \},\\{} & {} \kappa ((\theta _n-\theta _n^{hk})_x,(\vartheta _n-\vartheta _n^{hk})_x)= \kappa ((\theta _n-\theta _n^{hk})_x,(\theta _{t}(t_n)-\delta \theta _n)_x)\\{} & {} \quad +\kappa ((\theta _n-\theta _n^{hk})_x,(\delta {\theta }_n-\delta \theta _n^{hk})_x),\\{} & {} \displaystyle \kappa ((\theta _n-\theta _n^{hk})_x,(\delta {\theta }_n-\delta \theta _n^{hk})_x)\ge \frac{\kappa }{2k}\Big \{\Vert (\theta _n-\theta _n^{hk})_x\Vert ^2-\Vert (\theta _{n-1}-\theta _{n-1}^{hk})_x\Vert ^2\Big \}, \end{aligned}$$

we have, for all $\xi ^h\in E^h,$

$$\begin{aligned}{} & {} \displaystyle \frac{\overline{h}}{2k}\Big \{\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2-\Vert \vartheta _{n-1}-\vartheta _{n-1}^{hk}\Vert ^2\Big \}\nonumber \\{} & {} \displaystyle \qquad +\frac{\kappa }{2k}\Big \{\Vert (\theta _n-\theta _n^{hk})_x\Vert ^2-\Vert (\theta _{n-1}-\theta _{n-1}^{hk})_x\Vert ^2\Big \}\nonumber \\{} & {} \qquad +a((v_n-v_n^{hk})_{xx},\vartheta _n-\vartheta _n^{hk})-a((v_n-v_n^{hk})_{xx},\vartheta _n-\xi ^{h})\nonumber \\{} & {} \quad \displaystyle \le C\Big ( \Vert \vartheta _{t}(t_n)-\delta \vartheta _n\Vert ^2+\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2 +\Vert \vartheta _n-\xi ^h\Vert _E^2+\Vert (\theta _n-\theta _n^{hk})_x\Vert ^2\nonumber \\{} & {} \qquad +\Vert \theta _{t}(t_n)-\delta \theta _n\Vert ^2_E + (\delta {\vartheta }_n-\delta {\vartheta ^{hk}_n},\vartheta _n-\xi ^h)\Big ). \end{aligned}$$

(22)

Combining estimates (21) and (22) we find that, for all $w^h\in V^h,\, \xi ^h\in E^h$,

$$\begin{aligned}{} & {} \displaystyle \frac{\rho }{2k}\Big \{\Vert v_n-v_n^{hk}\Vert ^2-\Vert v_{n-1}-v_{n-1}^{hk}\Vert ^2\Big \} -a(\theta _n-\theta _n^{hk},(\delta u_n-\delta u_n^{hk})_{xx})\\{} & {} +\frac{\mu }{2k}\Big \{\Vert (u_n-u_n^{hk})_{xx}\Vert ^2-\Vert (u_{n-1}-u_{n-1}^{hk})_{xx}\Vert ^2\Big \} -a((\delta u_n-\delta u_n^{hk})_{xx},\vartheta _n-\xi ^{h})\\{} & {} +\frac{\overline{h}}{2k}\Big \{\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2-\Vert \vartheta _{n-1}-\vartheta _{n-1}^{hk}\Vert ^2\Big \} +\frac{\kappa }{2k}\Big \{\Vert (\theta _n-\theta _n^{hk})_x\Vert ^2-\Vert (\theta _{n-1}-\theta _{n-1}^{hk})_x\Vert ^2\Big \} \\{} & {} \displaystyle \le C\Big ( \Vert v_{t}(t_n)-\delta v_n\Vert ^2+\Vert v_n-v_n^{hk}\Vert ^2+\Vert v_n-w^h\Vert _V^2+\Vert (u_n-u_n^{hk})_{xx}\Vert ^2\\{} & {} +\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2+\Vert \theta _n-\theta _n^{hk}\Vert ^2+\Vert u_{t}(t_n)-\delta u_n\Vert _V^2 + (\delta {v}_n-\delta {v^{hk}_n},v_n-w^h)\\{} & {} +\Vert \vartheta _{t}(t_n)-\delta \vartheta _n\Vert ^2+\Vert \vartheta _n-\xi ^h\Vert _E^2+\Vert \theta _{t}(t_n)-\delta \theta _n\Vert ^2_E + (\delta {\vartheta }_n-\delta {\vartheta ^{hk}_n},\vartheta _n-\xi ^h) \Big ), \end{aligned}$$

Multiplying the above estimates by k and summing up to n we have, for all $\{w_j^h\}_{j=1}^n\subset V^h, \, \{\xi _j^h\}_{j=1}^n\subset E^h$,

$$\begin{aligned}{} & {} \displaystyle \Vert v_n-v_n^{hk}\Vert ^2 -k\sum _{j=1}^n(\theta _j-\theta _j^{hk},(\delta u_j-\delta u_j^{hk})_{xx})+\Vert (u_n-u_n^{hk})_{xx}\Vert ^2\\{} & {} +\Vert \vartheta _n-\vartheta _n^{hk}\Vert ^2- k\sum _{j=1}^n ((\delta u_j-\delta u_j^{hk})_{xx},\vartheta _j-\xi _j^{h})+\Vert (\theta _n-\theta _n^{hk})_x\Vert ^2\\{} & {} \displaystyle \le Ck\sum _{j=1}^n\Big ( \Vert v_{t}(t_j)-\delta v_j\Vert ^2_V+\Vert v_j-v_j^{hk}\Vert ^2+\Vert v_j-w_j^h\Vert _V^2+\Vert (u_j-u_j^{hk})_{xx}\Vert ^2\\{} & {} +\Vert \vartheta _j-\vartheta _j^{hk}\Vert ^2+\Vert \theta _j-\theta _j^{hk}\Vert ^2+\Vert u_{t}(t_j)-\delta u_j\Vert _V^2 + (\delta {v}_j-\delta {v^{hk}_j},v_j-w^h_j)\\{} & {} +\Vert \vartheta _{t}(t_j)-\delta \vartheta _j\Vert ^2+\Vert \vartheta _j-\xi _j^h\Vert _E^2+\Vert \theta _{t}(t_j)-\delta \theta _j\Vert ^2_E + (\delta {\vartheta }_j-\delta {\vartheta ^{hk}_j},\vartheta _j-\xi _j^h)\Big ). \end{aligned}$$

Since

$$\begin{aligned}{} & {} \displaystyle k\sum _{j=1}^n (\delta v_j-\delta v_j^{hk},v_j-w_j^h)=\sum _{j=1}^n (v_{j}-v_{j}^{hk}-(v_{j-1}-v_{j-1}^{hk}),v_j-v_j^h)\\{} & {} \quad =(v_n-v_n^{hk},v_n-w_n^h)+(v^{0h}-v^0,v_1-w_1^h)\\{} & {} \qquad \displaystyle + \sum _{j=1}^{n-1} (v_j-v_j^{hk},v_j-w_j^h-(v_{j+1}-w_{j+1}^h)),\\{} & {} \displaystyle k\sum _{j=1}^n (\delta \vartheta _j-\delta \vartheta _j^{hk},\vartheta _j-\xi _j^h)=\sum _{j=1}^n (\vartheta _{j}-\vartheta _{j}^{hk}-(\vartheta _{j-1}-\vartheta _{j-1}^{hk}),\vartheta _j-\xi _j^h)\\{} & {} \displaystyle \quad =(\vartheta _n-\vartheta _n^{hk},\vartheta _n-\xi _n^h)+(\vartheta ^{0h}-\vartheta ^0,\vartheta _1-\xi _1^h)\\{} & {} \qquad \displaystyle + \sum _{j=1}^{n-1} (\vartheta _j-\vartheta _j^{hk},\vartheta _j-\xi _j^h-(\vartheta _{j+1}-\xi _{j+1}^h)),\\{} & {} \displaystyle k\sum _{j=1}^n(\theta _j-\theta _j^{hk},(\delta u_j-\delta u_j^{hk})_{xx})= \sum _{j=1}^n(\theta _j-\theta _j^{hk},(u_j-u_j^{hk}-(u_{j-1}-u_{j-1}^{hk}))_{xx})\\{} & {} \quad \displaystyle =(\theta _n-\theta _n^{hk},(u_n-u_n^{hk})_{xx})+(\theta _1^{hk}-\theta _1,(u^0-u^{0h})_{xx})\\{} & {} \qquad \displaystyle + k\sum _{j=1}^{n-1}\Big ( (\vartheta _j-\vartheta _{j}^{hk},(u_j-u_j^{hk})_{xx})+ (\delta \theta _j-\theta _{t}(t_j),(u_j-u_j^{hk})_{xx})\Big ),\\{} & {} \displaystyle k\sum _{j=1}^n ((\delta u_j-\delta u_j^{hk})_{xx},\vartheta _j-\xi _j^{h})= \sum _{j=1}^n (( u_j-u_j^{hk}-(u_{j-1}-u_{j-1}^{hk}))_{xx},\vartheta _j-\xi _j^{h})\\{} & {} \quad \displaystyle =((u_n-u_n^{hk})_{xx},\vartheta _n-\xi _n^h)+((u^{0h}-u^0)_{xx},\vartheta _1-\xi _1^h)\\{} & {} \qquad \displaystyle + \sum _{j=1}^{n-1} ((u_j-u_j^{hk})_{xx},\vartheta _j-\xi _j^h-(\vartheta _{j+1}-\xi _{j+1}^h)),\\{} & {} \displaystyle \Vert \theta _n-\theta _n^{hk}\Vert ^2\le C\left( \Vert \theta ^0-\theta ^{0h}\Vert ^2+I_j+k\sum _{j=1}^n\Vert \vartheta _j-\vartheta _j^{hk}\Vert ^2\right) , \end{aligned}$$

where $I_j$ is the integration error given in (20), using again a discrete version of Gronwall’s inequality (see Campo et al. (2006)) we conclude the desired a priori error estimates.

We note that we can use the above a priori error estimates to derive the convergence order of the approximations under additional regularity conditions on the continuous solution. Therefore, if we assume that

$$\begin{aligned} \begin{array}{l} u\in H^3(0,T;Y)\cap W^{1,\infty }(0,T;H^3(0,\ell ))\cap H^2(0,T;V),\\ \theta \in H^3(0,T;Y)\cap W^{1,\infty }(0,T;H^2(0,\ell ))\cap H^2(0,T;E), \end{array} \end{aligned}$$

(23)

we obtain the following.

Corollary 1

Under the additional regularity conditions (23) and the assumptions of Theorem 2, we find that the approximations obtained by problem (12)–(14) are linearly convergent; that is, there exists a positive constant C, independent of the discretization parameters h and k, such that

$$\begin{aligned} \displaystyle \max _{0\le n\le N}\Big \{\Vert v_n-v_n^{hk}\Vert +\Vert u_n-u_n^{hk}\Vert _V +\Vert \theta _n-\theta _n^{hk}\Vert _E+\Vert \vartheta _n-\vartheta _n^{hk}\Vert \Big \}\le C(h+k). \end{aligned}$$

Remark 3

We note that we could analyze in a similar way variational problem (8), (6) and (7), approximated by the fully discrete problem (15), (13) and (14). Proceeding as in the proof of Lemma 1, we could derive the same dicrete stability property, and, following Theorem 2, we could prove a priori error estimates (19), which would lead to the linear convergence of the approximations under additional regularity conditions (23). However, we omit the details for the sake of clarity in the presentation of this work.

4 Numerical results

Now, we present the numerical algorithm that we have used for solving problem (12)–(14) and we show some numerical examples.

4.1 The numerical scheme

In this subsection, we describe the algorithm for solving problem (1), (3) and (4). Thus, given the solution $u_{n-1}^{hk}$, $v_{n-1}^{hk}$, $\theta _{n-1}^{hk}$ and $\vartheta _{n-1}^{hk}$ at time $t_{n-1}$, the discrete velocity $v_n^{hk}$ and the dicrete temperature speed $\vartheta _n^{hk}$ are the solution to the following linear system:

$$\begin{aligned}{} & {} \displaystyle \rho ({v}_n^{hk},w^h)+\mu k^2((v_n^{hk})_{xx},w^h_{xx})=\rho ({v}_{n-1}^{hk},w^h) -\mu k((u_{n-1}^{hk})_{xx},w^h_{xx})\\{} & {} \quad +a k(\theta _n^{hk}+\alpha \vartheta _n^{hk},w^h_{xx})\quad \forall w^h\in V^h,\\{} & {} \displaystyle (\overline{h}\vartheta _n^{hk}+d k \vartheta _n^{hk},\xi ^h) +\kappa k^2((\vartheta _n^{hk})_x,\xi ^h_x)=(\overline{h}\vartheta _{n-1}^{hk},\xi ^h)-\kappa k((\theta _{n-1}^{hk})_x,\xi ^h_x)\\{} & {} \quad -a k((v_n^{hk})_{xx},\xi ^h) \quad \forall \xi ^h\in E^h. \end{aligned}$$

Using the well-known commercial code MATLAB, this algorithm was implemented on a 3.2 Ghz PC, and we note that a typical run $(h=k=0.001)$ took about 1.92 seconds of CPU time.

4.2 Numerical convergence in a simple problem

As a first example, we solve problem (1), (3) and (4) to show the accuracy of the approximations. Then, we have used the following data in the simulations:

$$\begin{aligned} \begin{array}{l} \rho =1, \quad \mu =2, \quad a=2, \quad \alpha =0.2, \quad \overline{h}=1, \quad d=10, \quad \kappa =1. \end{array} \end{aligned}$$

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

$$\begin{aligned} u^0(x)=v^0(x)=\theta ^0(x)=\vartheta ^0(x)=x^3(x-1)^3, \end{aligned}$$

considering the homogeneous Dirichlet boundary conditions (3), if we add the (artificial) supply terms for each equation of system (1), for all $(x,t)\in (0,1)\times (0,1)$,

$$\begin{aligned} \begin{array}{l} F_1(x,t)=e^t(x^6 - 3x^5 - 69x^4 + 143x^3 + (3168x^2)/5 - (3528x)/5 + 144),\\ F_2(x,t)=e^tx(11x^5 - 33x^4 + 63x^3 - 71x^2 + 36x - 6), \end{array} \end{aligned}$$

the exact solution to this slightly modified version of problem (1), (3) and (4) can be easily calculated and it has the form, for $(x,t)\in [0,1]\times [0,1]$:

$$\begin{aligned} u(x,t)=\theta (x,t)=e^tx^3(x-1)^3. \end{aligned}$$

We note that the analysis of this problem is similar to the one shown in the previous section.

Therefore, in Table 1 we show the approximation errors estimated as

$$\begin{aligned} \begin{array}{l} \displaystyle \max _{0\le n\le N}\Big \{\Vert v_n-v_n^{hk}\Vert +\Vert u_n-u_n^{hk}\Vert _V +\Vert \theta _n-\theta _n^{hk}\Vert _E+\Vert \vartheta _n-\vartheta _n^{hk}\Vert \Big \} \end{array} \end{aligned}$$

for several values of the discretization parameters h and k and, in Fig. 1, we plot the evolution of the error depending on the parameter $h+k$. As can be clearly seen, the convergence of the algorithm is observed, and the linear convergence, stated in Corollary 1, is achieved.

Table 1 Example 1: Numerical errors for some values of h and k

Fig. 1

Example 1: asymptotic constant error with respect to parameter $h+k$

If we assume that there are not supply terms, and we use the data:

$$\begin{aligned} \begin{array}{l} T=60,\quad \rho =5, \quad \mu =2, \quad a=2, \quad \alpha =0.5, \quad \overline{h}=1, \quad d=10, \quad \kappa =1. \end{array} \end{aligned}$$

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

$$\begin{aligned} u^0(x)=v^0(x)=0,\quad \theta ^0(x)=\vartheta ^0(x)=100x^3(x-1)^3, \end{aligned}$$

taking the discretization parameter $h=0.001$ and $k=0.001$, the evolution in time of the discrete energy (see Quintanilla et al. 2023) defined as:

$$\begin{aligned} E(t)=5 \Vert v_n^{hk}\Vert ^2+2 \Vert (u_n^{hk})_{xx}\Vert ^2+ \Vert \vartheta _n^{hk}\Vert ^2+\Vert (\theta _n^{hk})_{x}\Vert ^2 \end{aligned}$$

is plotted in Fig. 2 (in both natural and semi-log scales). We note that we have represented their evolution (in natural scale), on the left-hand side, until time $t=10$ to see better the shape of the curves. Even if these decays seem to be exponential, on the right-hand side (semi-log scale) we can observe that, after time $t=30$, the decay is faster for the second-order coupling corresponding to Problem 1. We recall that, for this problem, it was proved that the energy decay was exponential meanwhile it was polynomial for Problem 2.

Fig. 2

Example 1: evolution in time of the discrete energy (natural and semi-log scales)

4.3 Dependence of the solution on the coupling parameter a

In this example, we will investigate the dependence on parameter a for the solution to the thermomechanical problem defined by system (1) with initial conditions (4) and boundary conditions (3), which we name as Problem 1, and for the solution to the thermomechanical problem defined by system (2) with initial conditions (4) and boundary conditions (3), which we name as Problem 2.

In these simulations, we have used the following data:

$$\begin{aligned} \begin{array}{l} T=1,\quad \rho =10, \quad \mu =2, \quad \alpha =0.5, \quad \overline{h}=1, \quad d=10, \quad \kappa =3, \end{array} \end{aligned}$$

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

$$\begin{aligned} u^0(x)=v^0(x)=x^3(x-1)^3,\quad \theta ^0(x)=\vartheta ^0(x)=100x^2(x-1)^2. \end{aligned}$$

Taking the discretization parameters $h=k=0.001$, the solution to Problem 1 is plotted in Figs. 3 and 4 for some values of parameter a. Regarding the displacements and the velocities (see Fig. 3), we can see that both displacements and velocities have a quadratic form, being larger when parameter a increases. If we consider the temperature and the temperature speed (see Fig. 4), we observe that the shape for this largest value of the coupling parameter a is drastically different.

Fig. 3

Example 2: displacement and velocity for different values of parameter a (Problem 1)

Fig. 4

Example 2: temperature and thermal velocity for different values of parameter a (Problem 1)

By using the same data that in the previous simulations, the solution to Problem 2 is shown in Figs. 5 and 6 for those values of parameter a. We can appreciate that, for values of parameter a less than 1, the displacements and velocities (see Fig. 5) the solutions are very similar. However, for the largest value of this parameter, it seems that there is a delay on the solution. The temperature and the temperature speed are plotted in Fig. 6 and we can see that both are quadratic, being rather similar.

Fig. 5

Example 2-Problem 2: displacement and velocity for different values of parameter a (Problem 2)

Fig. 6

Example 2: temperature and thermal velocity for different values of parameter a (Problem 2)

Appendix A: Polynomial energy decay

In a recent paper (Bazarra et al. 2022), the authors considered the problem determined by system (2) with boundary conditions (3) and initial conditions (4). It was proved the existence and uniqueness of solutions and their slow decay. That is, we cannot expect the uniform exponential energy decay of the solutions. However, we are going to see that we can obtain a polynomial energy decay.

We recall that this problem can be written in the following abstract form:

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t} U(t)=\mathcal {A} U(t),\quad U(0)=U^0, \end{aligned}$$

(24)

where the unknown variable is $U=(u,v,\theta ,\vartheta )$ and so, the initial data $U^0=(u^0,v^0,\theta ^0,\vartheta ^0)$. This problem is considered in the Hilbert space

$$\begin{aligned} \mathcal {H}=(H^2(0,\ell )\cap H^1_0(0,\ell ))\times L^2(0,\ell )\times H^1_*(0,\ell )\times L^2_*(0,\ell ), \end{aligned}$$

where

$$\begin{aligned} L^2_*(0,\ell )=\{g\in L^2(0,\ell );\quad \int _0^\ell g(x)\, \textrm{d}x=0\},\quad H^1_*(0,\ell )=H^1(0,\ell )\cap L^2_*(0,\ell ). \end{aligned}$$

Moreover, the matrix operator $\mathcal {A}$ is defined as

$$\begin{aligned} \mathcal {A}=\left( \begin{array}{c c c c} 0 &{} 1 &{} 0 &{} 0\\ \displaystyle -\frac{\mu }{\rho }D^4+\frac{b}{\rho }D^2 &{} 0 &{} \displaystyle \frac{a}{\rho } D &{} \displaystyle \frac{a}{\rho }\alpha D\\ 0 &{} 0 &{} 0 &{} 1\\ 0 &{} \displaystyle \frac{a}{\overline{h}} D &{} \displaystyle \frac{h}{\overline{h}}D^2 &{} \displaystyle -\frac{d}{h} \end{array}\right) . \end{aligned}$$

It was proved that $\mathcal {A}$ generates a contractive semigroup and that zero belongs to the resolvent of the operator. We can also show that its domain is the following:

$$\begin{aligned} \begin{array}{l} \mathcal {D}(\mathcal {A})=\{(u,v,\theta ,\vartheta )\quad ; \quad v\in H^2(0,\ell )\cap H^1_0(0,\ell ),\, \vartheta \in H^1_*(0,\ell ),\, \theta \in H^2(0,\ell ),\\ \qquad \qquad \theta _x\in H^1_0(0,\ell ),\, u\in H^4(0,\ell ),\, u_{xx}\in H^1_0(0,\ell )\}, \end{array} \end{aligned}$$

which is a dense subspace.

Theorem 3

The solutions to problem (24) satisfy

$$\begin{aligned} \Vert U(t)\Vert _{\mathcal {H}}\le M t^{-1/2}\Vert U(0)\Vert _{\mathcal {D}(\mathcal {A})}, \end{aligned}$$

where M is a positive constant.

Proof

To prove this result, we recall the characterization given by Borichev and Tomilov (2009), which guarantees the proposed decay whenever the imaginary axis is contained in the resolvent of the operator and the condition

$$\begin{aligned} \overline{\lim _{|\omega |\rightarrow \infty }} \omega ^{-2}\Vert (i\omega I-\mathcal {A})^{-1}\Vert <\infty \end{aligned}$$

(25)

holds. Since it was proved in Quintanilla et al. (2023) that the imaginary axis was contained in the resolvent of the operator, we only need to show condition (25). Let us assume that it does not hold. Then, we can find a sequence of real numbers $\omega _n\rightarrow \infty$ and a sequence of unit norm vectors at the domain of the operator $(u_n,v_n,\theta _n,\vartheta _n)$ such that

$$\begin{aligned} \begin{array}{l} \omega _n^2(i\omega _n u_n-v_n)\rightarrow 0 \quad \hbox {in}\quad H^2(0,\ell ),\\ \omega _n^2(i\rho \omega _nv_n-(\mu u_{nxxxx}-bu_{nxx}-a(\theta _{nx}+d\vartheta _{nx})))\rightarrow 0 \quad \hbox {in}\quad L^2(0,\ell ),\\ \omega _n^2(i\omega _n\theta _n-\vartheta _n)\rightarrow 0\quad \hbox {in}\quad H^1(0,\ell ),\\ \omega _n^2(i\overline{h}\omega _n\vartheta _n-(d\vartheta _n-\kappa \theta _{nxx}-av_{nx}))\rightarrow 0\quad \hbox {in} \quad L^2(0,\ell ). \end{array} \end{aligned}$$

The dissipation inequality implies that

$$\begin{aligned} \omega _n\theta _{nx},\, \omega _n\vartheta _n\rightarrow 0\quad \hbox {in}\quad L^2(0,\ell ). \end{aligned}$$

Now, we see that $\omega _n^{-1}u_{nxxxx}$ is bounded by the second convergence. Therefore, we can multiply the last convergence by $\omega _n^{-1}u_{nxxxx}$ and we find that

$$\begin{aligned} \begin{array}{l} \langle \vartheta _n,u_{nxxx}\rangle =-\langle \vartheta _{nx},u_{nxx}\rangle \rightarrow \langle i\omega _n\theta _{nx},u_{nxx}\rangle \rightarrow 0,\\ \langle \theta _{nxx},\omega _n^{-1}u_{nxxx}\rangle =-\langle \theta _{nx},\omega _{n}^{-1}u_{nxxxx}\rangle \rightarrow 0. \end{array} \end{aligned}$$

Then, we conclude that

$$\begin{aligned} \Vert u_{nxx}\Vert ^2=\langle u_{nxx},u_{nxx}\rangle =-\langle u_{nx},u_{nxxx}\rangle = \langle iv_{nx},\omega _n^{-1} u_{nxxx}\rangle . \end{aligned}$$

But

$$\begin{aligned} \begin{array}{l} -\langle v_{nx},\omega _n^{-1}u_{nxxx}\rangle \rightarrow a^{-1}\left[ \langle i \overline{h}\omega _n\vartheta _n,\omega _n^{-1}u_{nxxx}\rangle -d\langle \vartheta _n,\omega _n^{-1}u_{nxxx}\rangle +\kappa \langle \theta _{nxx},\omega _n^{-1}u_{nxxx}\rangle \right] \\ \qquad \rightarrow -a^{-1}\kappa \langle \theta _{nx},\omega _{n}^{-1}u_{nxxxx}\rangle \rightarrow 0, \end{array} \end{aligned}$$

then we also find that $v_n\rightarrow 0$ in $L^2(0,\ell )$, and we arrive to a contradiction because we assumed that condition (25) did not hold. We note that, following the same argument, we can prove that the imaginary axis is contained at the resolvent. Therefore, it follows that the thesis of the theorem holds.