Analytical and numerical analyses of a viscous strain gradient problem involving type Ⅱ thermoelasticity

1.   Introduction

Many studies dealing with the time decay of solutions have been developed in the last fifty years. It is common to couple several conservative equations with other dissipative equations and try to clarify the way in which the solutions of the whole system decay (when they do it). Perhaps, it corresponds to classical thermoelasticity, where a system of conservative equations is coupled with a heat equation. It is known that, in this case, there could exist (for certain geometries of the domain) undamped solutions, but, generally, we could expect the decay (polynomial) ^[13]. For the one-dimensional case, the decay is of exponential type. Another typical example corresponds to the thermoelastic plate with a parabolic heat equation, which leads to the coupling of a conservative equation with a dissipative one (even in dimensions greater than one), and the solutions decay in an exponential way. Moreover, they are determined by using an analytic semigroup; therefore, they exhibit a strong regularity.

However, there is a family of problems of the kind considered at the beginning for which we could not find many studies. In fact, it is common is to find a conservative system that is determined by the mechanical part of the problem and a dissipative mechanism within the thermal part of the system, but very few studies have been focused on the problem in which the mechanical mechanism is dissipative (viscosity) and the thermal one is conservative (see ^[14,16,17]), as in the case of type Ⅱ thermoelastic theory (see ^{[6,7,8,9,10,11]}). In this paper, we study a problem of this type.

We will consider the problem described by a one-dimensional viscous strain gradient elastic material with type Ⅱ heat conduction. It is worth remarking that, as far as we know, this is the second contribution in the literature where the viscosity is described by using fourth-order spatial derivatives. In a recent contribution ^[21], a similar system was considered, but it is relevant to point out that, in the present contribution, we apply a different coupling among the equations describing the system and that the type of viscosity can be different. We believe that it is important to clarify the qualitative properties to describe the decay of the solutions to this new problem. In this paper, we will focus on a problem involving a strain gradient beam with type Ⅱ heat conduction (without energy dissipation), assuming that the dissipative mechanisms act on the mechanical part.

It is worth noting that we can obtain the equations by following the methods in ^[12,20], but, in order to ensure that they are self-contained, we can recall that, in the one-dimensional case, the evolution equations are given by

where $\rho$ is the mass density, $u$ represents the displacement, $t$ is the stress tensor, $S$ is the double stress tensor, $T_0$ is the reference temperature, $\eta$ is the entropy, and $q = T_0\Phi$, where $\Phi$ is the entropy flux vector.

In our case, the constitutive equations become as follows:

In this system, $u$ describes the displacement and $\alpha$ is the thermal displacement satisfying that $\dot {\alpha} = \theta$, where $\theta$ is the temperature.

If we now introduce the constitutive equations into the evolution equations, we obtain the following linear system:

with the following boundary conditions:

and the initial conditions:

The system given by (1.1) and (1.2) describes the one-dimensional strain gradient viscoelasticity with type Ⅱ heat conduction. As usual, $\rho$ is the mass density, $c$ is the heat capacity, $\mu$ is the elastic constant, $\kappa_1$ is related to the hyperelasticity, and $\kappa_2$ is a typical constant for the type Ⅱ and Ⅲ theories of Green-Naghdi, but we point out that it is not the thermal conductivity; also, $\beta$ and $m$ are two coupling terms and $\kappa_1^*$ and $\mu^*$ describe the dissipation. In this work, we assume that $\rho, c, \mu, \kappa_1, \kappa_2$ are positive and the signs of $m$ and $\beta$ are not definite, but $\kappa_1 \kappa_2 > m^2$ and $\kappa_1^*, \mu^*$ must be non-negative, where

The meaning of the assumptions regarding $\rho$ and $c$ is clear. The assumptions regarding $\kappa_1, \kappa_2$ and $\mu$ can be understood with the help of the thermoelastic stability, as well as the assumption regarding $m$, which also guarantees that the elastic energy is positive definite. The assumptions regarding $\kappa_1^*$ and $\mu^*$ are also natural to ensure the existence of dissipation.

We know that the energy of the system is defined by

and it satisfies

In this paper, we will prove the existence and uniqueness of the solutions, the exponential decay, and the lack of regularity of the solutions, which implies, in particular, that the semigroup associated with the model is not analytic or differentiable. Moreover, we will also analyze this problem from the numerical point of view. It is worth noting that we also consider the case in which $\kappa_1^* = 0$, and we obtain similar qualitative properties for the solutions. Furthermore, the numerical simulations presented in Section 6 suggest that the most efficient dissipation mechanism requires that $\kappa_1^* = 0$ but $\mu^* > 0$. This is a remarkable fact, as one could think that the case in which the dissipation is stronger ($\kappa_1^* > 0, \, \mu^* > 0$) corresponds to the most efficient dissipative mechanism. Nevertheless, we point out that similar unexpected results have been obtained recently (see ^[1]). To the best of our knowledge, there are no previous studies that deal with the so-called type Ⅱ thermal law within strain gradient theory and with dissipative mechanisms that are imposed only on the mechanical subsystem.

The paper is structured as follows. The existence and uniqueness of solutions to the problem given by (1.1)–(1.4) are derived as shown in Section 2 by using semigroup arguments. In Section 3, the exponential energy decay is proved by applying the theory of linear semigroups. We also show that the semigroup associated with the solution to the thermoelastic problem in Section 4 is not differentiable, which implies the non-analyticity property (in particular, there is a lack of regularity). Then, a fully discrete approximation is introduced in Section 5 for a variational formulation by using the finite element method and the implicit Euler scheme. The discrete stability property and a priori error estimates are shown. Finally, some numerical simulations are presented in order to demonstrate the accuracy of the approximations, and a comparison is made among the solutions to the problems defined by the different dissipation mechanisms.

2.   Existence and uniqueness of solutions

In this section, we study the existence of solutions to the problem given by (1.1)–(1.4). It is important to recall that the total energy was defined in (1.6) and it satisfies (1.7).

Denoting $u_t = v, \alpha_t = \theta$, and $U = (u, v, \alpha, \theta)$, we can define the phase space $\mathcal{H}$ as follows:

where

This space encompasses an inner product which defines the following norm:

The previous assumptions imply that this norm is equivalent to the typical norm for the space.

Therefore, we can construct the operator $\mathcal{A}$ for a vector $U = (u, v, \alpha, \theta)$ as follows:

Let us denote by $\mathcal{A}_0$ the operator $\mathcal{A}$ with $\kappa_1^* = 0$. It is easy to verify that both are operators over the phase space $\mathcal{H}$. Therefore, the system given by (1.1)–(1.4) can be written as

or

when $\kappa_1^* > 0$ or $\kappa_1^* = 0$, respectively. Here, we have denoted $U_0 = (u_0, u_1, \alpha_0, \alpha_1).$

The corresponding domain is given by

Note that, if $\kappa_1^* = 0$, then it follows that $\kappa_1u_{xxx}+m\alpha_{xx} = g_1\in H^1 (0, \pi)$ and $mu_{xxx}+\kappa_2\alpha_{xx} = g_2\in L^2_*(0, \pi)$. By condition (1.5), we obtain that $\kappa_1\kappa_2 > m^2$; so, we can solve the two above equations for $u_{xxx}$ and $\alpha_{xx}$, concluding that $u_{xxx}, \alpha_{xx}\in L^2$. Hence, $D(\mathcal{A}_0)$ is given by

Furthermore, we see that $\mathcal{A}$ and $\mathcal{A}_0$ are dissipative operators verifying that

Using the resolvent equation $\lambda U-\mathcal{A}U = F$, we have that

Theorem 2.1. Let $\mathbf{A}$ be dissipative with $0\in\varrho(\mathbf{A}).$ If $\mathcal{H}$ is reflexive, then $\mathbf{A}$ is the infinitesimal generator of a semigroup of contractions.

Proof. Since $\varrho(\mathbf{A})$ is an open set, we have that there exists $\epsilon > 0$ such that $\epsilon \in \varrho(\mathbf {A})$. This implies that any $\lambda > 0$ belongs to $\varrho(\mathbf{A})$. In particular, we find that $R(I-\mathbf{A}) = \mathcal{H}$. Using Theorem 4.6 of ^[18], we conclude that $\overline{D(\mathbf{A})} = \mathcal{H}$; also, by applying the Lummer-Phillips theorem [18, Theorem 1.4.3], our conclusion follows. □

The method we use to show the exponential stability, as well as the the lack of differentiability, is based on the study of the resolvent system $i\lambda U-AU = F$, which, in terms of the components, is written in the following form:

Theorem 2.2. The operators $\mathcal{A}$ and $\mathcal{A}_0$ are the infinitesimal generators of a $C_0$ semigroup of contractions given by $(S(t))_{t\geq0}$ on $\mathcal{H}$. In particular, this implies the existence, uniqueness, and continuous dependence of solutions to the Cauchy problem (2.2), i.e., for any initial data $U_0\in \mathcal{H}$, there exists only one solution $U$ verifying that

Denoting by $\widetilde{\mathcal{A}}$ either the operator $\mathcal{A}$ or $\mathcal{A}_0$, we have that, if $U_0\in D(\widetilde{\mathcal{A}})$, then the solution has the following regularity:

Proof. We will prove it only for the operator $\mathcal{A}$ because the proof for the operator $\mathcal{A}_0$ is identical, taking $\kappa_1^* = 0$. Since $\mathcal{A}$ is dissipative, it is enough to show that $0\in \rho(\mathcal{A})$ since, by using Theorem 2.1, our result follows (see also ^[15]). Indeed, we show that, for any $F\in\mathcal{H}$, there exists only one $U\in D(\mathcal{A})$ verifying that $AU = F$. Let us take $F = (f_1, f_2, f_3, f_4)\in \mathcal{H}$. We must find $u$, $v$, $\alpha$, and $\theta$ that satisfy (2.10)–(2.13) for $\lambda = 0$. Thus, we have that $v = f_1$ and $\theta = f_3$, and, therefore, our problem reduces finding $u$ and $\alpha$ such that

We define the bilinear operator $a(\cdot, \cdot):V\times V \rightarrow \mathbb{R}$ over the space $V = H^2(0, \pi)\cap H_0^1(0, \pi)\times H_*^1(0, \pi)$ as follows:

and the linear form given by

It is easy to see that $T\in V^*$ and $a(\cdot, \cdot)$ is bilinear, continuous, and coercive (that is, because of conditions (1.5)). Hence, the Lax-Milgran lemma guarantees the existence and uniqueness of a solution:

Taking $\phi\in C_0^\infty(0, \pi)$ and $\psi = 0$, we get

This variational identity implies that $u$ satisfies (2.14) in the distributional sense. Moreover, since $\theta = f_3\in H^1_*$, we have that $\kappa_1 u_{xxxx} +m\alpha_{xxx}+\kappa_1^*f_{1, xxxx} = \kappa_1 u_{xxxx} +m\alpha_{xxx}+\kappa_1^*v_{xxxx}\in L^2$ or $\kappa_1 u_{xx} +m\alpha_{x}+\kappa_1^*v_{xx}\in H^2$. Similarly, taking $\psi\in C_0^\infty(0, \pi)$ and $\phi = 0$, we get

This variational identity ensures that $\alpha$ satisfies (2.15) in the distributional sense. Moreover, we also find that $-mu_{xxx}-\kappa_2\alpha_{xx}\in L^2$; so, the solution $(u, v, \alpha, \theta)\in D(\mathcal{A})$. Then, we have that $0\in \rho(\mathcal{A})$, and, by using Theorem 2.1, our result follows. Taking $\kappa_1^* = 0$ in the above procedure, we conclude that $0\in \rho(\mathcal{A}_0)$; so, $\mathcal{A}_0$ is also an infinitesimal generator of a $C_0$ semigroup of contractions. □

Despite the similarities between the operators $\mathcal{A}$ and $\mathcal{A}_0$, there is a big difference from the topological point of view, as we can see in the following lemma.

Lemma 2.1. The family of resolvent operators of $\mathcal{A}$ is not compact; however, the family of resolvent operators of $\mathcal{A}_0$ ($\mu^* = 0$) is a family of compact operators.

Proof. First, we note that the resolvent operator $(\lambda I-\mathcal{A})^{-1}:\mathcal{H}\rightarrow \mathcal{H}$ is compact if and only if $D(\mathcal{A})$ has a compact embedding in the phase space $\mathcal{H}$. In fact, let us suppose that the sequence $F_n \in \mathcal{H}$ is bounded. So, the sequence $U_n = (\lambda I-\mathcal{A})^{-1}F_n\in D(\mathcal{A})$ is bounded in $D(\mathcal{A})$. Because of the compact embedding, we conclude that there exists a subsequence of $U_n$ that converges strongly in $\mathcal{H}$.

Let us denote by $\mathcal{A}_*$ the operator $\mathcal{A}$ with $\mu^* = 0$. Since $D(\mathcal{A}) = D(\mathcal{A}_*)$, we conclude that the resolvent operator $(\lambda I-\mathcal{A})^{-1}$ is compact if and only if $(\lambda I-\mathcal{A}_*)^{-1}$ is also compact. Hence, in order to show the lack of compactness of the resolvent operators $(\lambda I-\mathcal{A})^{-1}$, we can suppose that $\mu^* = 0$.

We will show that there exists $\nu\in\mathbb{R}$, which is not an eigenvalue or an element of the resolvent set of $\mathcal{A}$. In fact, let us consider the following resolvent system for $\nu\in \mathbb{C}$:

Note that

Taking $\nu$ such that $\frac{\kappa_1\kappa_2}{m}+ \frac{\kappa_1^*\kappa_2}{m}\nu = m$, we get

We will show that $\nu$ is neither an element of the resolvent set of $\mathcal{A}$ nor an eigenvalue of $\mathcal{A}$. In fact, on the contrary, let us suppose that $\nu\in \varrho(A)$; then, the solution to the resolvent system must satisfy

Taking $f_1$ such that $f_1\in H^2(0, \pi)$ but $f_1\notin H^3(0, \pi)$, we get that the right hand side of (2.20) belongs to $H^{-2}(0, \pi)$, which is a contradiction because $\kappa_1 u_{xxxx}+\kappa_1^*v_{xxxx}+m\alpha_{xxx}$ must be in $L^2(0, \pi)$.

Finally, we will show that $\nu$ is not an eigenvalue. In fact, multiplying (2.17) by $u$ and (2.19) by $\alpha$ and recalling that $f_1 = f_2 = f_3 = f_4 = 0$, we get

Recalling the definition of $\nu$, we have

which implies that

Therefore, we conclude that $U = (u, v, \alpha, \theta)^\top = 0$; so, $\nu$ is not an eigenvalue. This implies that the family of the resolvent operators of $\mathcal{A}$ is not compact.

Finally, it is not difficult to see that $D(\mathcal{A}_0)$ has a compact embedding over the phase space $\mathcal{H}$. In fact, let us suppose that a sequence $F_n$ is bounded in $\mathcal{H}$. Then, $(i\lambda I-\mathcal{A}_0)^{-1}F_n = U_n = (u_n, v_n, \alpha_n, \theta_n)$ is bounded in $D(\mathcal{A}_0)$. In particular, we have that $(u_n, v_n, \alpha_n, \theta_n)$ is bounded in $H^3\times H^1\times H^2 \times H^1\subset \mathcal{H}$ with a compact embedding. Thus, there exists a subsequence of $U_n$ that can be denoted in the same way, and it converges strongly in $\mathcal{H}$. Hence, the resolvent operator $(\lambda I-\mathcal{A}_0)^{-1}$ is compact. The proof is now complete. □

3.   Exponential decay

In this section, we will prove that the solutions to the problems studied in the previous section decay in an exponential way. To this end, we assume that $m\neq 0.$ The main tool we use is the following result according to Prüss ^[19].

Theorem 3.1. Let $S(t)$ be a contraction $C_{0}$ semigroup, generated by ${\mathbb A}$ over a Hilbert space ${\cal H}$. Then, there exists $C, \gamma > 0$ verifying that

Lemma 3.1. The operator $\mathcal{A}$ ($\mathcal{A}_0$) defined by (2.1) verifies that $i \mathbb{R}\subset \varrho(\mathcal{A})$ ($i \mathbb{R}\subset \varrho(\mathcal{A}_0)$).

Proof. In the case that $\mu^* = 0$, from Lemma 2.1, the operator $\mathcal{A}_0$ has a compact family of resolvent operators. So, in order to show that $i \mathbb{R}\subset \varrho(\mathcal{A}_0)$, it is enough to prove that there are no imaginary eigenvalues. We will proceed by contradiction. Let us suppose that there exists $U\ne0$ such that $\mathcal{A}_0U = i\lambda U$. Since $F = 0$, the relation (2.9) yields that $v_x = 0$, and, by using (2.10), we get that $u = v = 0$. From (2.11), we get that

Taking the real and imaginary parts in the implication above yields that $\alpha = 0$, and, by using (2.12), we get that $\theta = 0$. Therefore, we conclude that $U = 0$, which is a contradiction.

Finally, we consider the case in which $\mu^* > 0$. Since $0 \in \varrho (\mathcal{A})$, the set

is not empty. Setting $\sigma = \sup \mathcal{N}$, if $\sigma = +\infty$, there is nothing to prove. By contradiction, if we suppose that $\sigma < \infty$, then there exist a sequence $\lambda_n \in \mathbb{R}$ such that $\lambda_n\to \sigma < \infty$ and a sequence of unit norm vectors $\{ \widetilde{F}_n \}_n$ of elements of $\mathcal{H}$ such that

Letting $U_n = \frac{W_n}{ \Vert W_n \Vert_{\mathcal{H}}}$ and $F_n = \frac{\widetilde{F}_n}{ \Vert W_n \Vert_{\mathcal{H}}}$, we have

Note that the sequence $\mathcal{A}U_n$ is bounded in $\mathcal{H}$. From (2.8) and (2.10), we get

On the other hand, applying (2.11) with the sequence $(u_n, v_n, \alpha_n, \theta_n)^\top$ and multiplying the result by $\int_0^x \overline{\alpha_n} d\tau$, we find that

Since

in view of the strong convergence of $u_n$ and $v_n$, we conclude that

From (3.3) and (3.4), it follows that $U_n$ converges strongly to $0$. This contradicts the fact that $\left\| U_n\right\|_{\mathcal{H}} = 1$ for all $n\in\mathbb{N}$.

The proof is now complete. □

Theorem 3.2. The semigroup generated by the operators $\mathcal{A}$ and $\mathcal{A}_0$ defined by (2.1) are exponentially stable; that is, there exist two positive constants $M$ and $\omega$ such that

Proof. Thanks to Theorem 3.1 and Lemma 3.1, we only need to show that

Using inequality (2.9), we get

Multiplying (2.11) by $\overline{u}$, we have

Using (2.10) and (2.12), we obtain

which gives

where

Using relation (3.6), we get

On the other hand, multiplying (2.19) by $\overline{u}_x$, we find that

Then, we have

where

Using the above procedure, we also conclude that

Solving (3.7) and (3.8) for $\int_0^\pi |u_{xx}|^2$ and $\int_0^\pi\alpha_{x}\overline{u}_{xx}$, we find that

Multiplying (2.11) by $\int_0^x \overline{\alpha}$, we get

Using (2.12), we obtain

Therefore, we have

where

Keeping in mind (3.6) and (2.12), we conclude that

Hence, taking the real and imaginary parts in (3.11), we get

Multiplying (2.13) by $\overline{\alpha}$ and performing an integration by parts, we have

Using (2.12), we find that

and, taking into account (3.9) and (3.12), we get

Therefore, we obtain

which implies that $\|U\|_{\mathcal{H}}\leq c\|F\|_{\mathcal{H}}$. Therefore, the exponential stability follows.

Finally, in the case in which $\mu^* > 0$, relation (2.8) implies that

Multiplying (2.11) by $\int_0^x \overline{\alpha}$, and by using the same procedure as that for inequality (3.10), we get

Then, we have

where

Using (3.15) and (2.12), we conclude that

Hence, taking the real and imaginary parts in (3.12), we get

Finally, using (3.14), we obtain

Therefore, we conclude that

It implies that $\|U\|_{\mathcal{H}}\leq c\|F\|_{\mathcal{H}}$; so, the proof is now complete.

Remark 3.1. We have proved the existence and the exponential stability of the solutions to the problem given by (1.1)–(1.4). In this case, we have assumed homogeneous Neumann boundary conditions on the thermal displacement. It is clear that, using the standard modifications, we can obtain the same kind of results when we consider homogeneous Dirichlet boundary conditions on the thermal displacement. In fact, in Sections 5 and 6, we obtain similar results when these boundary conditions are imposed, but from a numerical point of view.

4.   Lack of differentiability of the semigroup

In this section, we will show that the semigroup associated with the system given by (1.1)–(1.4) is not differentiable ^[18] (and not immediately differentiable ^[5]). To see this, we recall the following results.

Theorem 4.1. Let $S = (S(t))_{t \geq 0}$ be an immediately differentiable semigroup on the Banach space $X$; then, $S(t)$ is an immediately norm-continuous semigroup (see ^[5], Definition 4.17, page 112).

Proof. If $S(t)$ is immediately differentiable, then $S(t)$ is immediately differentiable with the uniform norm of $\mathcal{L}(X)$ for any $t > 0$. This implies that the semigroup is immediately norm-continuous. □

Theorem 4.2. If $\mathcal{A}$ is the generator of an immediately norm-continuous exponentially stable semigroup, then

Proof. See ^[5], Corollary 4.19, page 114. □

Theorem 4.3. The semigroup $S = (S(t))_{t \geq 0}$ defined by the system given by (1.1)–(1.4) is not differentiable.

Proof. From Theorems 4.1 and 4.2, it is enough to show that there exist a sequence $(\lambda_n)_n$ of real numbers and a bounded sequence $(F_n)_n$ in $\mathcal{H}$ with $\lim\limits_{n\to\infty} \lambda_n = \infty$ such that

In fact, for each $n\in \mathbb{N}$, we consider $F_n = (0, 0, 0, \sqrt{\dfrac{2}{\pi}}\cos(nx))$, which is bounded in $\mathcal{H}$. Let $(u_n, v_n, \alpha_n, \theta_n)\in D(\mathcal{A})$ be the unique solution to the following resolvent equation:

Now, let us fix any constant $\mu^* > 0$. We consider two cases:

Case 1. $\kappa_1^*\neq 0.$ In this case, system (4.1) can be broken down into its components as follows:

Given by boundary conditions described by (1.3), the solutions to the above system have the following forms:

Substituting this into the previous system, we get that

Then, the following two equations must be satisfied:

Then, from (4.7) and (4.8), we have

We let $\kappa_2n^2-\lambda_n^2 c = 0$. Note that our choice is correct because

given the conditions on $\kappa_2$ and $c$, with $\gamma = \sqrt{\frac{\kappa_2}{c}}$. Under the above conditions, we have

Hence, we find that the $L^2(0, \pi)$-norm of $\theta_n$ is given by

Hence, we conclude that

Since $U_n = (i\lambda_n I-\mathcal{A})^{-1}F_n$, our conclusion follows.

Case 2. $\kappa_1^* = 0.$ In this case, we are solving the system given by (4.2)–(4.5), and by using the solutions detailed in (4.6), we get

Then, the following two equations must be satisfied:

Then, from (4.7) and (4.8), we get

We let $\kappa_2n^2-\lambda_n^2 c = \frac{m^2}{\kappa_1}n^2$, which is correct because

given the conditions on $\kappa_1, \kappa_2$, and $m$, with $\gamma = \sqrt{\frac{\kappa_1\kappa_2-m^2}{c\kappa_1}}$. Under the above conditions, we have

Hence, we find that the $L^2(0, \pi)$-norm of $\theta_n$ is given by

Therefore, we find that

So, our conclusion follows.

Remark 4.1. The above result implies that the semigroup defined by system given by (1.1)–(1.4) is not differentiable and, hence, not analytic. In particular, it implies that the solution does not have the smoothing effect property.

5.   A fully discrete approximation: a priori error estimates

In this section, we will provide an a priori error analysis of the thermoelastic problems described in the previous sections. However, for the sake of generality, we will assume that the length of the beam is $\ell > 0$ instead of $\pi$, and we will study its deformation over a finite time interval $[0, T]$, with $T > 0$. We recall the assumptions required for the constitutive coefficients of (1.5). Now, we shall derive the weak form of the thermomechanical problem defined by system given by (1.1) and (1.2) with the initial conditions presented in (1.4), as well as the following boundary conditions:

Thus, let us denote $Y = L^2(0, \ell)$, $E = H^1_0(0, \ell)$, and $V = H^2_0(0, \ell)$, and denote by $(\cdot, \cdot)$ the inner product of $L^2(0, \ell)$.

Therefore, by applying integration by parts to the constitutive equations (1.1) and (1.2), and by taking into account the new boundary conditions of (5.1), we obtain the following variational formulation written in terms of the velocity $v = u_t$ and the temperature $\theta = \alpha_t$.

Find the velocity field $v:[0, T]\rightarrow V$ and the temperature $\theta:[0, T]\rightarrow E$ such that $v(0) = u_1$ and $\theta(0) = \alpha_1$, and, for all $w\in V$ and $r\in E$,

where the displacement field $u$ and the thermal displacement $\alpha$ are recovered via the following relations:

In what follows, we provide an a priori error analysis of a fully discrete scheme that can be obtained by using the classical finite element method and the implicit Euler scheme.

Now, we introduce the approximation of problem given by (5.2)–(5.4). This is done first from the spatial point of view. So, we assume that the interval $[0, \ell]$ is split into $M$ subintervals $a_0 = 0 < a_1 < \ldots < a_M = \ell$ of uniform length $h = a_{i+1}-a_i = \ell/M$. Therefore, we approximate the variational spaces $E$ and $V$ by the finite element spaces $E^h\subset E$ and $V^h\subset V$, defined as follows:

In these definitions, we have represented by $P_q([a_i, a_{i+1}])$ the space of polynomials of degree less or equal to $q$ in the subinterval $[a_i, a_{i+1}]$. Moreover, let us define the discrete initial conditions $u^{0h}$, $v^{0h}$, $\alpha^{0h}$ and $\theta^{0h}$ as follows:

Here, $\mathcal{P}_1^h$ and $\mathcal{P}_2^h$ denote the respective projection operators over the finite element spaces $E^h$ and $V^h$ (see ^[3]).

Now, we discretize the time derivatives. Therefore, we use a uniform partition of the time interval $[0, T]$, denoted by $0 = t_0 < t_1 < \ldots < t_N = T$, with time step size $k = T/N$ and nodes $t_n = n\, k$ for $n = 0, 1, \ldots, N$. As usual, for a continuous function $z(t)$, we denote $z_n = z(t_n)$ and, given a sequence $\{z_n\}_{n = 0}^N$, let $\delta z_n = (z_n-z_{n-1})/k$ be its divided differences.

By using the well-known implicit Euler scheme, a fully discrete approximation of the problem given by (5.2)–(5.4) can be obtained as follows.

Find the discrete velocity $\{v^{hk}_n\}_{n = 0}^N \subset V^h$ and the discrete temperature $\{\theta_n^{hk}\}_{n = 0}^N \subset E^h$ such that $v_0^{hk} = v^{0h}$, $\theta_0^{hk} = \theta^{0h}$, and, for all $w^h\in V^h, \, r^h\in E^h$, and $n = 1, \ldots, N$,

where the discrete displacements denoted by $u^{hk}_n$ and the discrete thermal displacement denoted by $\alpha^{hk}_n$ are now recovered from the following relations:

We note that it is easy to prove that the above fully discrete problem admits a unique solution. It is obtained by using the classical Lax-Milgram lemma and the conditions defined in (1.5).

Now, we will confirm the discrete stability property, which is summarized as follows.

Lemma 5.1. Under conditions defined in (1.5), it follows that the discrete solution to the problem given by (5.8)–(5.10) $\{u_n^{hk}, v_n^{hk}, \alpha_n^{hk}, \theta_n^{hk}\}_{n = 0}^N$ satisfies the following estimates:

where the positive constant $C$ is independent of the discretization parameters $h$ and $k$, and, from now on, let $\|\cdot\|_{\mathcal{H}}$ be the typical norm in the Hilbert space $\mathcal{H}$.

Proof. Taking $w^h = v_n^{hk}$ as a test function in discrete variational equation (5.8), we find that

Now, keeping in mind that (given the conditions of (1.5))

we obtain

Proceeding in a similar manner as with the discrete variational equation (5.9), for a test function $r^h = \theta_n^{hk}$, we obtain

Combining the previous estimates, it follows that

Observing that

and taking into account that, under the condition that $\kappa_1\kappa_2 > m^2$, we have

Taking the above estimates, multiplying them by $k$, and summing up to $n$, we find that

where $C$ is a positive constant which does not depend on $h$ or $k$.

By again applying the condition that $\kappa_1\kappa_2 > m^2$, we have verified the desired stability property. □

In the rest of the section, we will derive some estimates on the numerical errors $v_n-v_n^{hk}$ and $\theta_n-\theta_n^{hk}$. First, we will obtain the estimates for the velocity field. So, taking (5.2) with a test function $w = w^h\in V^h$ at time $t = t_n$, and subtracting it by discrete variational equation (5.8), it follows that, for all $w^h \in V^h$,

Keeping in mind that

and by applying Cauchy's inequality $ab\leq \varepsilon a^2+\frac{1}{4\varepsilon}b^2$ for $a, b, \varepsilon \in \mathbb{R}$ with $\varepsilon > 0$, it follows that, for all $w^h\in V^h$,

Proceeding in a similar manner, we obtain the following estimates for the temperature field for all $r^h\in E^h$:

After easy algebraic manipulations, we have

Observing that, thanks again to the condition that $\kappa_1\kappa_2 > m^2$, we find that

and, summing up to $n$ the previous estimates on the velocity and temperatures, it follows that

Now, taking into account that

where we have again applied the conditions of (1.5) as well as the estimates

with a similar expression for the temperature, after an application of a discrete version of Gronwall's lemma (see ^[4]), we obtain the following a priori error estimate result.

Theorem 5.1. Under the conditions of (1.5), let us denote by $(u, v, \alpha, \theta)$ the solution to the problem given by (5.2)–(5.4), and denote by $(u^{hk}, v^{hk}, \alpha^{hk}, \theta^{hk})$ the solution to the problem given by (5.8)–(5.10); then, we have the following a priori error estimates for all $w^h = \{w_j^{h}\}_{j = 0}^N \subset V^h$: $r^h = \{r_j^{h}\}_{j = 0}^N\subset E^h$:

where $C$ is a positive constant which does not depend on the discretization parameters $h$ and $k$, and we recall that $\|\cdot\|_{\mathcal{H}}$ denotes the norm of the Hilbert space $\mathcal{H}$.

From the above result, we can derive the convergence order of the approximations for the discrete problem given by (5.8)–(5.10). As an example, if we assume the following additional regularity:

we conclude that the convergence of the algorithm is linear. It can be proved by applying some well-known results on the approximation through the use of finite elements (see, e.g., ^[2]) and some estimates already used in ^[4]. That is, it follows that there exists a positive constant $C > 0$ such that

6.   Numerical examples

In this final section, the numerical scheme for solving the problem given by (5.8)–(5.10) is described and some numerical simulations are presented to show the convergence of the approximations and the behavior of the discrete energy.

First, we note that, given the solutions $u_{n-1}^{hk}$ and $\alpha_{n-1}^{hk}$ at time $t_{n-1}$, the discrete velocity $v_{n}^{hk}$ and the discrete temperature $\theta_{n}^{hk}$ are obtained by solving the following linear system for all $w^h\in V^h, \, r^h\in E^h$:

This numerical scheme was implemented on a 3.2 GHz PC by using MATLAB, and a typical run (using parameters $h = k = 0.001$) took about 0.52s of CPU time.

As an example, in order to show the accuracy of the approximations and the behavior of the solution, we solved problem given by (1.1), (1.2), (1.4), and (5.1) with the following data:

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

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

the exact solution to this problem can be obtained, and it has the following form for $(x, t)\in [0, 1]\times[0, 1]$:

We note that the analysis performed in the previous section and the numerical scheme presented above should be modified to include these supply terms, but the corresponding changes are really minor.

Thus, taking some discretization parameters $h$ and $k$,

was used to obtain the estimated approximation errors as are shown in Table 1. Moreover, the evolution of the error according to the parameter $h+k$ was determined and plotted in Figure 1. We notice that the convergence of the algorithm can be clearly observed, and that the linear convergence, as stated in the previous section, is achieved.

Now, we are going to compare the energy decay of the two problems corresponding to the dissipation mechanisms studied in this work. So, let us denote as "Problem 1" the problem resulting from applying $\mu^* = 1$ and $\kappa_1^* = 0$ (that is, with a dissipation mechanism of second order), denote as "Problem 2" the problem associated with applying $\kappa_1^* = 1$ and $\mu^* = 0$ (i.e., a dissipation mechanism of fourth order), and, finally, denote as "Problem 3" the problem associated with applying $\kappa_1^* = 1$ and $\mu^* = 1$ (that is, with the previous two dissipation mechanisms).

In these simulations, we have assumed that there are no supply terms, and we have applied the following data:

By applying the initial conditions for all $x\in (0, 1)$:

and, by taking the discretization parameters $h = k = 0.001$, the evolution in time of the discrete energy of the two problems given by

was obtained as in Figure 2 (in both natural and semi-log scales) for the three cases described above. As can be seen, the three discrete energies converge to zero and an exponential decay seems to be achieved. Moreover, as expected, the decay is faster for the second-order mechanism (Problem 1); meanwhile, the case with the two dissipation mechanisms together has a worse rate decay, which is similar to that of the fourth-order dissipation mechanism. This is an amazing result.

Finally, our aim will be to numerically demonstrate the lack of differentiability of the solutions to the problem given by (1.1)–(1.4). Hence, following the convergences described by (4.2)–(4.5), we will solve the following coupled discrete linear system:

where the function $F_n(x) = \sqrt{\frac{2}{\pi}}\cos{nx}$ for all $x\in [0, \pi]$. We note that, in this case, the problem is static; so, the approximation is done only in space. Moreover, the subscript $n$ represents the discrete solution obtained for each value $\lambda_n$.

In this example, we have used the following data:

Our aim here will be to show that the discrete solution to the above system $U_n^{h} = (u_n^{h}, v_n^{h}, \alpha_n^{h}, \theta_n^{h})$ has a norm defined by

which is greater than a positive constant when $\lambda_n\rightarrow \infty$.

Therefore, taking as a discretization parameter $h = \pi/1000$, in Figure 3, we show the evolution of the norm $\|U_n^{h}\|^2$ when $\lambda_n$ increases to $10^5$. As we can see on the left-hand side, the norm of these discrete solutions seems to stabilize near $0.663$; this can be clearly checked if we zoom in near value $10^6$ (see the right-hand side). Thus, proceeding as in the proofs in Section 4, we could conclude that the semigroup is not differentiable.

7.   Conclusions

In this work, we studied a one-dimensional viscous strain gradient problem with several dissipation mechanisms from the analytical and numerical points of view. It is remarkable that we have assumed mechanical dissipation of strain gradient type under the condition of a conservative heat conduction.

From the analytical point of view, we have proved that the problem has a unique solution, and that the energy decays in an exponential way. Moreover, we have also shown the lack of differentiability of the semigroup.

From the numerical point of view, we have introduced a fully discrete approximation of a weak form of the above thermal problem, and we have proved the existence of a discrete stability property and a main a priori error estimate result. Through the numerical example, we have demonstrated the linear convergence and analyzed the convergence rate according to the number of dissipation mechanisms considered. We have found that the energy decays faster when the second-order dissipation is used. Finally, we have shown numerically that the discrete solutions have a norm which is greater than a positive value for any value $\lambda_n$ introduced into the problem; so, we have concluded that the semigroup is not differentiable.

Author contributions

J. E. Muñoz-Rivera, E. Ochoa and R. Quintanilla: Conceptualization, Methodology; J. R. Fernández, J. E. Muñoz-Rivera, E. Ochoa and R. Quintanilla: Formal analysis; J. R. Fernández, E. Ochoa and R. Quintanilla: Writing-original draft; B. Noelia and J. R. Fernández: Writing-review & editing, Visualization. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The work of N. Bazarra, J. R. Fernández, and R. Quintanilla was supported by the project "Qualitative and numerical analyses of problems in Thermomechanics" (Ref. PID2023-146359NB-I00), which is currently under evaluation in the Spanish Ministry of Science, Innovation and Universities.

The work of J. E. Muñoz-Rivera was supported by CNPq project 307947/2022-0 and Fondecyt project 1230914.

The work of E. Ochoa was supported by CONICYT-PFCHA/doctorado nacional/2020-21200268.

Conflict of interest

Ramón Quintanilla is an editorial board member and was not involved in the editorial review or the decision to publish this article. The authors declare that they have no conflict of interest.