Numerical approximation of a Bresse-Maxwell-Cattaneo system

1.   Introduction

The study of the so-called Bresse system (introduced by J. A. C. Bresse in 1856 (see ^[1])), and also known as the circular arch problem, includes the mechanical deformation of beams with natural length $\ell$. This problem can be considered as a generalization of the well-known Timoshenko problem, which is obtained in the particular case where the longitudinal displacement is not modeled and supposing zero initial curvature.

If we denote by $\varphi$, $\psi$ and $w$ the vertical displacement, the rotation angle of the cross-section and the horizontal displacement, respectively, the corresponding system of equations is written as (see, for details, the derivation provided in^[2])

where $\rho_1$ and $\rho_2$ are the mass density and the moment of mass inertia of the beam, $\kappa$ is the shear modulus of elasticity, $l$ is the initial curvature, $\kappa_0 = Eh$, where $E$ is the Young's modulus and $h$ is the cross-sectional area, and $b$ stands for the rigidity coefficient of the cross section.

There is a large number of papers dealing with different issues involving this Bresse system. In this way, just as a few examples, we can recall the work of Ma and Nunes ^[3], where the limit case corresponding to a zero arch curvature (i.e., the Timoshenko system) and the long-time dynamics are shown, the work of Baibeche et al. ^[4], who introduced the microtemperatures effect with a delay term on a boundary condition, the work of Li et al. ^[5], who studied the asymptotic behavior of this system with mass diffusion, the work of Muñoz-Rivera and Nunes ^[6], where the pointwise stabilization on an interior point is assumed to prove the exponential stability, the work of Copetti et al. ^[7], where a contact problem with the viscoelastic Bresse system is mathematically studied, the work of Bazarra et al. ^[8], where the so-called dual-phase-lag model is incorporated into the heat conduction, the works of Afilal et al. ^[9,10], where the second sound effect is also incorporated, and the work of Suzuki and Ueda ^[11], where the dissipation of a thermoelastic Bresse system is studied with different damping effects. We also recall the works of Alabau Boussouira et al. ^[12], where the lack of exponential stability is proved if the velocities of wave propagation are not the same, and the work of Almeida Júnior ^[13], where the locking phenomenon on the shear force is considered for the Timoshenko equation. Moreover, we could cite the recent numerical works of Wang et al. ^[14] and Yang et al. ^[15] for the study of some related fourth-order problems.

In this paper, we want to continue the contributions ^[16,17], studying a numerical approximation of the thermoelastic Bresse system by using the heat conduction model proposed by Maxwell and Cattaneo ^[18], where a relaxation parameter is introduced into the constitutive equations to be compatible with the causality principle and to overcome the instantaneous propagation of thermal waves ^[19].

Let us denote by $(0, \ell)$ the interval occupied by the curved thermoelastic beam of length $\ell > 0$, and let $x\in [0, \ell]$ and $t\in [0, T]$ be the spatial and time variables, respectively, where $T > 0$ denotes the final time.

The linearized version of the Bresse-Maxwell system by using the hyperbolic heat conduction proposed by Cattaneo is written as follows in $(0, \ell)\times (0, T)$ (see ^[16,17,20]):

where now $\vartheta$ is the temperature deviation from a fixed reference temperature in the vertical direction, $\xi$ is the temperature deviation from a fixed reference temperature in a horizontal direction, $p$ denotes the vertical heat flux, and $q$ stands for the horizontal heat flux. The new coefficients $\gamma$, $\varsigma$, $\varpi$, and $\rho_3$ account for the physical properties of the beam, and $\tau$ represents the delay parameter.

As it is pointed out in ^[16], these coefficients satisfy the constraints:

This system (1.1) is complemented with the following boundary conditions:

for $X = 0, \ell$, and the initial conditions:

for a.e. $x\in (0, \ell)$.

The following existence and uniqueness result was proved in ^[17].

Theorem 1. The operator that generates the solutions to problem (1.1)–(1.3) is the infinitesimal generator of a $C^0$-semigroup of contractions in a suitable Hilbert space. Therefore, problem (1.1)–(1.3) has a unique solution with the regularity:

Moreover, even if the stability of the solutions to problem (1.1)–(1.3) has been discussed by some authors, we recall the recent result obtained in ^[16].

Theorem 2. Let us define the stability numbers:

Therefore, the semigroup generated by the solutions to problem (1.1)–(1.3) is exponentially stable if and only if $\chi_{ \mbox{ $\varsigma$ }}\chi_\tau = 0$. On the contrary, it is (semiuniformly) polynomially stable with optimal decay rate $\sqrt{t}$.

From ^[17], we observe that the energy norm given by

is equivalent to the usual norm used in the analysis of these thermoelastic problems defined as

whenever we assume the general condition $l\ell\neq n\pi$ for all $n\in\mathbb{N}$.

Now, we will derive the variational formulation of problem (1.1)–(1.3). Let us define the variational spaces $V = H^1_0(0, \ell)$ and $E = H^1(0, \ell)$, and denote by $(\cdot, \cdot)$ the classical inner product in the space $Y = L^2(0, \ell)$ with corresponding norm $\|\cdot\|$. Therefore, multiplying the equations of system (1.1) by adequate test functions and taking into account the boundary conditions (1.2), we arrive at the following weak formulation written in terms of the variables $\phi = \varphi_t$, $\eta = \psi_t$, $e = w_t$, $\vartheta$, $p$, $\xi$, and $q$.

Find the vertical velocity $\phi:[0, T]\rightarrow V$, the rate of the rotation angle $\eta:[0, T]\rightarrow V$, the horizontal velocity $e:[0, T]\rightarrow V$, the vertical temperature $\vartheta:[0, T]\rightarrow V$, the horizontal temperature $\xi:[0, T]\rightarrow V$, the vertical heat flux $p:[0, T]\rightarrow E$, and the horizontal heat flux $q:[0, T]\rightarrow E$ such that, for a.e. $t\in (0, T)$ and for all $r, v, z, m, \chi\in V, \, s, \zeta\in E$,

where the vertical displacement $\varphi:[0, T]\rightarrow V$, the rotation angle $\psi:[0, T]\rightarrow V$, and the horizontal displacement $w:[0, T]\rightarrow V$ are obtained from the relations

The rest of the paper is outlined as follows. In Section 2 we study numerically an approximation of this problem (1.4)–(1.11) by using the classical finite element method to approximate the spatial variable and the implicit Euler scheme to discretize the time derivatives. A discrete stability property and some a priori error estimates are proved, from which the linear convergence is obtained under suitable regularity conditions on the continuous solution. Finally, some numerical simulations are presented in Section 3 to demonstrate the convergence of the approximation, the discrete energy decay, and the behavior of the solution.

2.   Numerical analysis of a fully discrete approximation

In this section, we will numerically study the variational problem (1.4)–(1.11). Therefore, we will provide first a numerical scheme based on the classical finite element method and the implicit Euler scheme.

2.1.   Fully discrete approximation

In this subsection, we introduce a fully discrete approximation of the weak problem introduced previously. We will do it in two steps. In order to provide the spatial approximation, let us define a uniform partition of the spatial domain denoted by $a_0 = 0 < \ldots < a_M = \ell$, and let us construct the finite element spaces:

In this definition, we have denoted by $P_1([a_i, a_{i+1}])$ the space of polynomials with degree less than or equal to one in the subinterval $[a_i, a_{i+1}]$ and, as usual, let $h = a_{i+1}-a_i = \ell/M$ be the spatial parameter.

Now, we can define the discrete initial conditions as

where $P_1^h$ and $P_2^h$ are the finite element interpolation operators over $V^h$ and $E^h$, respectively ^[21].

In order to obtain the discretization of the time derivatives, we denote by $t_0 = 0 < \ldots < t_N = T$ a uniform partition of the time interval with time step $k = t_1-t_0 = T/N$ and nodes $t_n = nk$ for $n = 0, \ldots, N$. As a usual notation, let $w_n = w(t_n)$ be the value of a continuous function $w(t)$ at time $t = t_n$, and for a sequence $\{w_n\}_{n = 0}^N$, let $\delta w_n = (w_n-w_{n-1})/k$ be its divided differences.

By using the well-known implicit Euler scheme, we have the following fully discrete problem.

Find the discrete vertical velocity $\{\phi_n^{hk}\}_{n = 0}^N\subset V^h$, the discrete rate of the rotation angle $\{\eta_n^{hk}\}_{n = 0}^N\subset V^h$, the discrete horizontal velocity $\{e_n^{hk}\}_{n = 0}^N\subset V^h$, the discrete vertical temperature $\{\vartheta_n^{hk}\}_{n = 0}^N\subset V^h$, the discrete horizontal temperature $\{\xi_n^{hk}\}_{n = 0}^N\subset V^h$, the discrete vertical heat flux $\{p_n^{hk}\}_{n = 0}^N\subset E^h$, and the discrete horizontal heat flux $\{q_n^{hk}\}_{n = 0}^N\subset E^h$ such that, for $n = 1, \ldots, n$ and for all $r^h, v^h, z^h, m^h, \chi^h\in V^h, \, s^h, \zeta^h\in E^h$,

where the discrete vertical displacement $\{\varphi_n^{hk}\}_{n = 0}^N\subset V^h$, the discrete rotation angle $\{\psi_n^{hk}\}_{n = 0}^N\subset V^h$, and the discrete horizontal displacement $\{w_n^{hk}\}_{n = 0}^N\subset V^h$ are obtained from the relations

Remark 3. We note that the above fully discrete problem can be written in terms of the unknowns $\phi_n^{hk}$, $\eta_n^{hk}$, $e_n^{hk}$, $\vartheta_n^{hk}$, $p_n^{hk}$, $\xi_n^{hk}$, and $q_n^{hk}$ by using the definition of the divided differences (we recall the notation $\delta w_n = (w_n-w_{n-1})/k$); see the description of the numerical resolution given in Section 3. Therefore, it leads to a linear system written in terms of a product variable $U_n^{hk} = (\phi_n^{hk}, \eta_n^{hk}, e_n^{hk}, \vartheta_n^{hk}, p_n^{hk}, \xi_n^{hk}, q_n^{hk})$, where the associated matrix is positive definite thanks to the required assumptions on the constitutive coefficients. Therefore, applying Lax-Milgram lemma, it is straightforward to show that the fully discrete problem (2.1)–(2.8) has a unique solution.

2.2.   Discrete stability and a priori error analysis

In the rest of the section, we will show the numerical analysis of the above problem, providing its discrete stability and a main a priori error estimate result.

First, the discrete stability is summarized in the following.

Lemma 4. There exists a positive constant $C$, which is independent of the discretization parameters $h$ and $k$, such that, for $n = 1, \ldots, N$,

where $\|\cdot\|_X$ stands for the usual norm in the Hilbert space $X$.

Proof.

First, we obtain the estimates for the discrete vertical velocity. Taking as a test function $r^h = \phi_n^{hk}$ in the discrete variational equation (2.1), we find that

Taking into account that $\phi_n^{hk} = \delta \varphi_n^{hk}$ and the estimates:

using Cauchy-Schwarz inequality and Cauchy's inequality

we have

where, here and in what follows, $C > 0$ is a positive constant that does not depend on the discretization parameters $h$ and $k$.

In the same way, we can obtain the estimates for the discrete rate of the rotation angle and the discrete horizontal velocity (we omit here the details for the sake of clarity in the presentation). So, we find that

Now, we derive the estimates for the discrete vertical temperature and the discrete vertical heat flux. By using the test functions $m^h = \vartheta_n^{hk}$ and $s^h = p_n^{hk}$ in discrete variational equations (2.4) and (2.5), respectively, we find that

Keeping in mind the estimates

using several times Cauchy-Schwarz inequality and Cauchy inequality (2.9), we have

Proceeding in an analogous way with the discrete horizontal temperature and the discrete horizontal heat flux, it follows that

Combining all these estimates, we find that

Multiplying these estimates by $k$ and summing up to $n$, we have

Finally, applying a discrete version of Gronwall's lemma (see ^[22] for details), we obtain the desired stability estimates.

□

Now, we will obtain some a priori error estimates. First, subtracting the variational equation (1.4) at time $t = t_n$ and for a test function $r = r^h\in V^h\subset V$ and discrete variational equation (2.1), it follows that

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

Now, keeping in mind that

we obtain the estimates:

where $C > 0$ is again a positive constant that does not depend on the discretization parameters $h$ and $k$, but it depends now on the continuous solution.

Similarly, we derive the following estimates for the rate of the rotation angle and the horizontal velocity. As we pointed out in the proof of the discrete stability, we omit the details for the sake of clarity in the presentation. So, we obtain

Secondly, we obtain the estimates for the vertical temperature and the vertical heat flux. Thus, we subtract the variational equations (1.7) and (1.8) at time $t = t_n$ and for test functions $m = m^h\in V^h\subset V$ and $s = s^h\in E^h\subset E$, respectively, and the discrete variational equations (2.4) and (2.5) to obtain

and so, we conclude that, for all $m^h\in V^h$ and $s^h\in E^h,$

By using the following estimates:

applying several times Cauchy-Schwarz inequality and Cauchy's inequality (2.9), it follows that, for all $m^h\in V^h$ and $s^h\in E^h$,

Proceeding in an analogous way with the discrete horizontal temperature and the discrete horizontal heat flux, we obtain, for all $\chi^h\in V^h$ and $\zeta^h\in E^h$,

Combining all the previous estimates, we have, for all $r^h, v^h, z^h, m^h, \chi^h\in V^h, \, s^h, \zeta^h\in E^h$,

Multiplying the above estimates by $k$ and summing up to $n$, we find that, for all $\{r_j^h\}_{j = 1}^n, \{v_j^h\}_{j = 1}^n, \{z_j^h\}_{j = 1}^n, \{m_j^h\}_{j = 1}^n, \{\chi_j^h\}_{j = 1}^n\subset V^h, \, \{s_j^h\}_{j = 1}^n, \{\zeta_j^h\}_{j = 1}^n\subset E^h$,

Now, we observe that

where similar relations can be found for the remaining variables.

Finally, applying a discrete version of Gronwall's lemma (see again ^[22]), we conclude the following a priori error estimates result.

Theorem 5. Under the assumptions imposed on the constitutive coefficients, let $(\phi, \varphi, \eta, \psi, e, w, \vartheta, p, \xi, q)$ and $(\phi^{hk}, \varphi^{hk}, \eta^{hk}, \psi^{hk}, e^{hk}, w^{hk}, \vartheta^{hk}, p^{hk}, \xi^{hk}, q^{hk})$ be the solutions to problems (1.4)–(1.11) and (2.1)–(2.8), respectively. Therefore, it follows that, for all $\{r_n^h\}_{n = 0}^N, \{v_n^h\}_{n = 0}^N, \{z_n^h\}_{n = 0}^N, \{m_n^h\}_{n = 0}^N, \{\chi_n^h\}_{n = 0}^N\subset V^h, \, \{s_n^h\}_{n = 0}^N, \{\zeta_n^h\}_{n = 0}^N\subset E^h$,

where $C$ is a positive constant that depends on the continuous solution but it is independent of the discretization parameters.

The a priori error estimates shown above can be used to derive the convergence order. As an example, if we assume that the continuous solution has the regularity

then we approximations provided by the fully discrete problem (2.1)–(2.8) converge linearly; that is, there exists a positive constant $C > 0$, independent of the discretization parameters $h$ and $k$, such that

3.   Numerical results

In this final section, we describe the numerical scheme implemented in MATLAB for solving the fully discrete problem (2.1)–(2.8), and we show some numerical examples to demonstrate the accuracy of the approximations and the decay of the discrete energy.

3.1.   Numerical scheme

As a first step, given the solution $\varphi_{n-1}^{hk}, \psi_{n-1}^{hk}, w_{n-1}^{hk}, \vartheta_{n-1}^{hk}, p_{n-1}^{hk}, \xi_{n-1}^{hk}$ and $q_{n-1}^{hk}$ at time $t_{n-1}$, variables $\phi_{n}^{hk}, \eta_{n}^{hk}, e_{n}^{hk}, \vartheta_{n}^{hk}, p_{n}^{hk}, \xi_{n}^{hk}$ and $q_{n}^{hk}$ are obtained by solving the discrete linear system, for all $r^h, v^h, z^h, m^h, \chi^h\in V^h, \, s^h, \zeta^h\in E^h$,

where $F_i$, for $i = 1, \ldots, 7,$ represents artificial supply terms that are introduced in order to make the system more general.

The numerical scheme was implemented on a 2.9 GHz PC using MATLAB, and a typical run (using parameters $h = k = 0.01$) took about 0.12 sec of CPU time.

3.2.   First example: numerical convergence

As an academic example, in order to show the accuracy of the approximations, we have solved an example with the following data:

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

considering homogeneous Dirichlet boundary conditions for the variables $\varphi, \psi, w, \vartheta, \xi$ and the (artificial) supply terms, for all $(x, t)\in [0, 1]\times [0, 1]$,

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

Here, if we approximate the numerical errors by

they are depicted in Table 1 using different values of $h$ and $k$. Furthermore, in Figure 1 we show the evolution of the above numerical errors depending on the parameter $h+k$. It is clear that the convergence is found, and the linear convergence (that is, first order), stated in the above section, is achieved.

3.3.   Second example: Energy decay

In this last example, we are going to study the energy decay of the discrete problem. Therefore, we assume that there are no supply terms, and we use the following data:

and the initial conditions, for all $x\in [0, 1]$,

With this data, we note that the stability numbers defined in Theorem 2 are both equal to zero. Thus, the semigroup generated by the solutions of the problem is exponentially stable because the condition $\chi_{ \mbox{ $\varsigma$ }}\chi_\tau = 0$ is satisfied. Taking the discretization parameters $h = 0.001$ and $k = 0.0001$, the evolution in time of the discrete energy of the problem, defined as

is plotted in Figure 2 (in both natural and semi-log scales).

As can be seen, it converges to zero, and the exponential decay seems to be achieved.

Taking the data

and the same initial conditions as before, for all $x\in [0, 1]$. With the above data, the stability numbers defined in Theorem 2 are both non equal to zero ($\chi_{ \mbox{ $\varsigma$ }} = \chi_\tau = 10$), and so the semigroup generated by the solutions to this problem is not exponentially stable. In fact, the energy decay is (semiuniformly) polynomially stable with optimal decay rate $\sqrt{t}$. Taking the discretization parameters $h = 10^{-3}$ and $k = 10^{-5}$, the evolution in time of the discrete energy of the problem is plotted in Figure 3 (in both natural and semi-log scales).

Finally, in order to verify if the theoretical polynomial energy decay is achieved, we compare the above discrete energy with the function $-C(\log{t})/2$, where $C$ is a positive constant to be fixed. After an adjustment procedure, we have found the value $C = 4$, and so in Figure 4 the evolution of the discrete energy (in semi-log scale) is plotted against the evolution of the function $-2\log{t}$. As can be seen, a similar behavior of both graphs is obtained.

4.   Conclusions

In this work, we have studied, from the numerical point of view, a new thermoelastic Bresse system where the heat conduction has been modeled by using the Cattaneo-Maxwell law. This leads to a hyperbolic law, and it was considered, from the analytical point of view, by Dell'Oro ^[16] and Lima and Fernández-Sare ^[17]. Here, we have focused on a fully discrete approximation based on the finite element method to approximate the spatial variable and the implicit Euler scheme to discretize the time derivatives. The stability of the discrete solution has been proved, and an a priori error analysis has been provided. As a particular case, the linear convergence of the approximations has been deduced assuming some regularity conditions on the continuous solution. Finally, we have shown some numerical results, including an example to demonstrate the numerical convergence (Example 1) and another one to analyze the discrete energy decay (Example 2). In the latter one, we have seen how it depends on the condition proved in ^[17] for the so-called stability numbers $\chi_{ \mbox{ $\varsigma$ }}$ and $\chi_\tau$.

Use of AI tools declaration

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

Acknowledgements

This paper is part of the project "Qualitative and numerical analyses of some thermomechanical problems (ACUANUTER)" (Ref. PID2024-156827NB-I00), which is currently under evaluation by the Spanish Ministry of Science, Innovation and Universities.

The authors would like to thank the anonymous reviewers for their helpful and interesting comments, which have improved the quality of this work.

Conflict of interest

The authors declare there are no conflicts of interest.