# Regularized variational formulation and 3D discrete approximation with enhanced performance in thermo-elasticity and thermo-damage

Samir Suljević, Adnan Ibrahimbegovic, Samir Dolarević

► To cite this version:

Samir Suljević, Adnan Ibrahimbegovic, Samir Dolarević. Regularized variational formulation and 3D discrete approximation with enhanced performance in thermo-elasticity and thermo-damage. Computational Mechanics, 2025, 76 (4), pp.1183-1210. 10.1007/s00466-025-02645-x. hal-05382484

HAL Id: hal-05382484

https://hal.science/hal-05382484v1

Submitted on 2 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Regularized variational formulation and 3D discrete approximation with enhanced performance in thermo-elasticity and thermo-damage

Samir Suljević<fnref n="1" /> · Adnan Ibrahimbegovic<fnref n="1" />,<fnref n="3" /> · Samir Dolarević<fnref n="2" />

Received: 16 January 2025 / Accepted: 17 April 2025 / Published online: 4 June 2025

© The Author(s) 2025

## Abstract

In this paper, we present four-field variational formulation of thermo-elasticity and thermo-damage with non-symmetric stress tensor that is particularly suitable for long-term simulations in thermodynamics needed for fatigue studies. We propose the regularization of this variational formulation in order to enforce the stability of mixed formulation with respect to the Babushka-Brezzi condition. This further allows to eliminate the rotation field, and obtain hybrid-stress format for dynamics and the hybrid heat-flux for transient heat conduction. The regularized variational formulation is combined with a 3D discrete approximation based on the Raviart-Thomas vector space for both mechanic and thermal fields to ensure continuity across tetrahedron finite element faces for the stress vector and normal component of heat-flux. The proposed formulation and the chosen discrete approximations in this study provide superior accuracy for both thermal and mechanical fields compared to the classical FE methods. Both energy-conserving and energy-decaying schemes are utilized for time discretization to ensure robustness in long-term computations. The proposed method is validated through several numerical simulations in both statics and dynamics, including stationary and non-stationary heat conduction, for isotropic thermo-elastic and thermo-damage materials. The potential use in thermal fatigue studies is also demonstrated.

**Keywords** Regularized variational formulation · Hybrid stress-heat-flux discrete approximation · Raviart-Thomas vector space · Energy-conserving and energy-decaying time integration scheme · Thermo-elasticity · Thermo-damage

## 1 Introduction

There has been renewed interest in classical fatigue problems with many competing formulations that allows to predict the fatigue failure (eg. [1–3]). The ability to tackle the long-term computations is the essential factor of any such formulation, with vast majority restricted to purely mechanics framework. This work seeks to provide the novel variational formulation,

discrete approximation and time-integration that would allow dealing with the fatigue problems within the more general framework of thermodynamics.

There are many other problems in engineering practice that involve the combined effects of mechanical and thermal phenomena. These may arise, for example, when thermal stresses develop during constrained thermal deformation, or when heat is generated due to mechanical work or dissipation. The former case involves the mechanical response to thermal loads, while the latter corresponds to the thermal response to mechanical loads. An example of the latter case is the temperature rise that occurs in certain regions during the process of metal forming in materials engineering, as a result of the imposed strain rates.

In the context of thermomechanics, the coupling between mechanical and thermal fields can be considered as weak or strong. In weakly coupled thermomechanics, the mechanical field is assumed to depend on pre-computed temperature changes. This simplified approach is often used in fire safety engineering, where the behavior of structures under fire loading is considered for pre-defined thermal loading imposed

<figure><img src="image_2.png" /></figure>

Adnan Ibrahimbegovic
adnan.ibrahimbegovic@utc.fr

Samir Suljević
samir.suljevic@gf.unsa.ba

Samir Dolarević
samir.dolarevic@gf.unsa.ba

<footnote n="1">¹ Laboratoire Roberval, Université de Technologie de Compiègne - Alliance Sorbonne Université, 60200 Compiègne, France</footnote>

<footnote n="2">² Faculty of Civil Engineering, University of Sarajevo, 71000 Sarajevo, Bosnia and Herzegovina</footnote>

<footnote n="3">³ Institut Universitaire de France, IUF, 75000 Paris, France</footnote>
by design codes. Here we are interested in strongly coupled thermomechanics, considering the interdependence between deformation and temperature field. This means that we need to account for structural heating or damage dissipation as potential heat source when computing temperature field from energy balance equation and for temperature-dependent material properties in computing mechanical field from equilibrium equation.

The previous research has extensively investigated the fundamental aspects of the coupled thermomechanics approach, with notable contributions from the classical work [4]. The finite element method is used along with a staggered solution approach to solve fully coupled thermomechanics problems, as presented in [5–7]. In [8], a mixed variational statement is proposed for two-way coupled thermoelastic analysis for quasi-static applications, where stresses are assumed to be in equilibrium with field load, efficiently reducing locking phenomena. Additionally, the temperature at inter-element is independent of the temperature within the element, leading to an improved thermal response computation. Similar, staggered approach to deal with the thermomechanical coupling is investigated for inelastic behavior for plasticity and damage in [9, 10].

In this work, we deviate from such staggered solution methods and propose the simultaneous solution of thermomechanics and thermodynamics by using a second-order scheme. Our main objective is to achieve superior accuracy in stress and heat flux computations and ensure stability in long-term computations, such as those required for fatigue failure risk analysis under variable loads and in structural fire engineering under long-term thermal loading [11]. The second novelty pertains to the discrete approximation. Since the standard displacement-based finite element method’s best approximation property [12, 35] only applies to static problems, a different approach is necessary to ensure the best performance for dynamic problems. Thus, we use mixed variational formulation and hybrid stress finite elements to achieve higher accuracy in stress and heat flux computations for dynamic problems. These proposed elements provide independent interpolations for stress and displacement fields on one side and temperature and heat flux on the other side.

This is indeed the novelty with respect to most of the classical research on mixed finite elements [12–14], that consider improving the accuracy eliminating the locking of the traditional finite element method in problems with constraints, such as incompressibility, or various other approaches proposed to enhance the performance of the finite element method, by including the incompatible modes as discussed in [15, 16], as well as the interpolation of stress fields that satisfies nodal point equilibrium, as presented in [17]. In order to enhance the accuracy of results for dynamics, the finite element method has been replaced with alternative methods, such as discontinuous Galerkin [18], or more recently, iso-

geometric analysis [19]. In this paper we demonstrate that the finite element method can achieve superior accuracy in both statics and dynamics through the use of a (regularized) variational formulation and a hybrid-stress discrete approximation that enables the enforcement of stress vector continuity across element boundaries. Interpolations of this kind have proved to be effective for scalar field problems, as shown in recent work on electrostatics [20], where they are referred to as Whitney elements [21], and in studies of diffusion problems [22, 23], where they are referred to as the Raviart-Thomas discrete approximations.

This paper presents a rigorous variational formulation and corresponding regularization technique that builds upon our previous works [24–26] to provide increased stability properties for discrete approximations for thermodynamics. The regularized variational formulation is employed to eliminate the independent rotation field, which then allows for the recovery of the hybrid stress and heat flux approach for the thermodynamics. The introduction of the regularized term leads to an additional Euler-Lagrange equation, which makes it possible to condense the five-field variational formulation into a four-field variational formulation. This results with the non-symmetric stress tensor, enabling the interpolation of the stress field using either Whitney’s approximation [20] or the lowest order of the Raviart-Thomas ($RT_0$) vector space [22, 23].

The proposed approximation builds upon the traditional constant strain tetrahedron element, and extends it by adding four more nodes placed at the centroid of each element face to store the stress and heat flux components. The displacement and temperature fields are approximated using conventional shape functions for triangles [12, 35]. The stress degrees of freedom are components of traction along with the normal flux component at the corresponding element edge. The Raviart-Thomas interpolation ensures continuity of the stress vector or normal heat flux across the entire domain at a given edge of the finite element. This is particularly beneficial for dynamic problems and non-stationary heat transfer, where maintaining continuity is essential to ensure smooth propagation across element boundaries.

The final development involves the corresponding time integration schemes for long-term computations. Several such schemes have been successfully developed, and proven to remain stable over an extended period of time [29, 30] by enforcing the energy conservation in the elastic regime. Another scheme that has been shown to produce higher accuracy of stress within each time step interval [33, 34] is energy-decaying scheme. This scheme is well-suited for visco-damage analysis, where higher accuracy of stress computations is desirable compared to the backward Euler method that are typically used for operator split procedures analysis [35].
The paper outline is as follows. The variational formulation of thermodynamics using mixed functional is presented in Section 2, restricted only to small strain theory which is typical of brittle material behavior described by damage model described herein. Section 3 is devoted to numerical implementation of linear elastodynamics problem including the chosen discrete approximation of displacement, temperatures, stress and heat-flux fields. Additionally, the energy conserving time integration scheme is developed to maintain the total energy of system over computational time as well as the energy-decaying modification of such scheme. In Section 4, several illustrative simulations of the thermo-dynamics problem under static and dynamic loading scenarios are presented. In Section 5, we give the main results discussion and the conclusions.

## 2 Thermodynamics coupling variational formulation

### 2.1 General continuum thermodynamic model

The starting point is the local form of the first principle of thermodynamics on energy conserving for the case of thermo-mechanical inelastic response, for the elastic case where we assume frozen values for damage and hardening variables denoted with $\bar{D}$ and $\bar{\xi}$; (e.g. see [35]):

$$ \dot{\epsilon}(\epsilon, s, \bar{D}, \bar{\xi}) = \sigma : \dot{\epsilon} + r - \nabla \cdot \mathbf{q} \quad (1) $$

where $\epsilon$ is the internal energy density, $\epsilon$ is the strain field, $s$ is the corresponding entropy, $\sigma$ is the stress field, $r$ is the internal heat supply, and $\mathbf{q}$ is the heat flux. We note in passing that the internal energy density can also be defined with respect to infinitesimal volume with mass density stated explicitly, which is more suitable for large strains where the density can change between initial and deformed configurations in order to enforce the mass conservation; e.g. see [35].

For a more general inelastic case where the values of internal variables evolve, we assume the additive decomposition of entropy into elastic $s^e$ and damage component $s^d$ (such that $s = s^e + s^d$), and further appeal to the Legendre transformation(e.g. see [35]) to exchange the role of entropy and temperature, thus replacing the internal energy by the free energy $\psi$:

$$ \psi = e - \theta s^e \quad (2) $$

where $\theta$ denotes the absolute temperature.

For such inelastic case, the thermo-damage free energy density $\psi(\epsilon, \theta, \mathbf{D}, \xi)$ is the function of the state variables: the total strain $\epsilon$, the temperature $\theta$, the damage compliance tensor $\mathbf{D}$ and the hardening variable $\xi$. The Clausius-Duhem

inequality can be used to define dissipation rate for the damage model, which is written as

$$ \begin{aligned} D_{int} := & \sigma : \dot{\epsilon} + \theta \dot{s} - \dot{\epsilon} = \sigma : \dot{\epsilon} + \theta \dot{s} - \dot{\psi} \\ & - \dot{s}^{e} \theta - s^{e} \dot{\theta} \geq 0 \\ D_{int} := & \left( \sigma - \frac{\partial \psi}{\partial \epsilon} \right) : \dot{\epsilon} - \left( s^{e} + \frac{\partial \psi}{\partial \theta} \right) \dot{\theta} \\ & + (\theta \dot{s}^{e} - \theta \dot{s}^{e}) + \theta \dot{s}^{d} - \frac{\partial \psi}{\partial \mathbf{D}} : \dot{\mathbf{D}} - \frac{\partial \psi}{\partial \xi} \dot{\xi} \geq 0 \end{aligned} \quad (3) $$

where we neglected the dissipation by conduction (e.g. see [35]). In the elastic case, where internal variables remain frozen with $\dot{\mathbf{D}} = 0$, $\dot{\xi} = 0$ and $\dot{s}^d = 0$, the total damage dissipation rate remains equal to zero and we obtain the constitutive equations for the stress and the elastic entropy

$$ \begin{aligned} \sigma &:= \frac{\partial \psi}{\partial \epsilon} \quad (4) \\ s^e &:= -\frac{\partial \psi}{\partial \theta} \quad (5) \end{aligned} $$

By assuming that the same constitutive equations also remain valid in the inelastic case due to (much) faster rate of evolution for internal variables vs. global variables, we can obtain from (3) the (positive) inelastic dissipation rate:

$$ 0 < D_{int} = - \frac{\partial \psi}{\partial \mathbf{D}} : \dot{\mathbf{D}} - \frac{\partial \psi}{\partial \xi} \dot{\xi} + \theta \dot{s}^d \quad (6) $$

For any heat conduction with non-homogeneous temperature field, we also obtain an additional term of dissipation by conduction (e.g. see [35]). Such term confirms the Fourier law for heat transfer ($\mathbf{q} = -k\nabla\theta$, with $k > 0$ as the conductivity coefficient):

$$ 0 \le D_{cond} := - \frac{1}{\theta} \underbrace{\mathbf{q} \nabla \theta}_{\frac{1}{\theta} k(\nabla \theta)^2} \quad (7) $$

In further development, this term is neglected. Thus, the thermodynamics problem to solve is defined by (mechanics) equations that pertain to linear and angular momentum balance and definition of infinitesimal strain

$$ \rho \ddot{\mathbf{u}} = \nabla \cdot \sigma + \mathbf{b} ; \quad \sigma = \sigma^T ; \quad \epsilon = \nabla^s \mathbf{u} \quad (8) $$

where $\rho$ is mass density, $\mathbf{u}(\mathbf{x}, t)$ is displacement field, $\sigma(\mathbf{x}, t)$ is stress field, $\mathbf{b}(\mathbf{x}, t)$ is external load, and $\nabla = \frac{\partial}{\partial \mathbf{x}}$. The symmetry of the stress tensor is not enforced in the strong form, but the angular momentum balance can be recovered from the Hu-Washizu variational formulation that is further proposed in Eq. (33).

For the elastic case, with the frozen values of internal variables ($\dot{\mathbf{D}} = 0, \dot{\xi} = 0$), we also need to solve at the (thermal)
part with the reduced form of energy conservation, which can be simplified by using the constitutive equations

$$ \theta \dot{s} + \underbrace{\left(\frac{\partial \psi}{\partial \epsilon} - \boldsymbol{\sigma}\right)}_{=0} : \dot{\boldsymbol{\epsilon}} + \underbrace{\left(\frac{\partial \psi}{\partial \theta} + s\right)}_{=0} \dot{\theta} = r - \nabla \cdot \mathbf{q} \quad (9) $$

A specific approach is further proposed here for damage model by using the free energy potential written as

$$ \psi(\boldsymbol{\epsilon}, \theta, \mathbf{D}, \xi) = \frac{1}{2} [\boldsymbol{\epsilon} - \boldsymbol{\alpha}(\theta - \theta_0)] : \mathbf{D}^{-1} : [\boldsymbol{\epsilon} - \boldsymbol{\alpha}(\theta - \theta_0)] + \Xi(\xi) + \rho c \left[ (\theta - \theta_0) - \theta \ln \frac{\theta}{\theta_0} \right] \quad (10) $$

where $\boldsymbol{\alpha} = \mathbf{1}\boldsymbol{\alpha}$ with $\boldsymbol{\alpha}$ as the thermal expansion coefficient, $\theta_0$ is the reference temperature, $\Xi(\xi)$ hardening potential and $c$ is heat capacity coefficient. We note that $\mathbf{D} = \mathbf{D}^e \equiv \mathbf{C}^{-1}$ (with $\mathbf{C}^{-1}$ as the elasticity tensor) and $\xi = 0$ only at the initial stage, assuming no initial defects.

The evolution of the internal variables for inelastic case, is defined in agreement with the chosen damage criterion, written in thermodynamics framework as

$$ \phi(\boldsymbol{\sigma}, q, \theta) = \|\boldsymbol{\sigma}\|_{\mathbf{D}} - (\sigma_f(\theta) - q(\xi, \theta)) \le 0 \quad (11) \\ \|\boldsymbol{\sigma}\|_{\mathbf{D}} = [\boldsymbol{\sigma} : \mathbf{P}\boldsymbol{\sigma}]^{\frac{1}{2}} $$

where $\sigma_f$ is the limit of elasticity and $q$ is the stress-like variable which defines the threshold of fracture, both of them temperature dependent, with the simplest choice made for isotropic damage assuming $\mathbf{P} = \mathbf{D}^e$. A more general choice of the failure modes to be represented by the damage compliance tensor $\mathbf{D}$) are possible, and they can be explored with respect to different failure modes that typic of the quasi-brittle material; in that respect we refer to our recent work [54] on plasticity model for concrete that can represent full set of 3D failure modes in tension, compression and shear, which would be starting point in generalising the present approach to damage.

In agreement with a choice of damage criterion made here, we also prefer using the stress $\boldsymbol{\sigma}$ as the primary state variable, rather than strain $\boldsymbol{\epsilon}$, further postulating the quadratic form of complementary energy (rather than $\boldsymbol{\epsilon} = \mathbf{D}\boldsymbol{\sigma}$) and the quadratic form of hardening potential, such that

$$ \chi(\boldsymbol{\sigma}, \mathbf{D}, \xi) = \frac{1}{2} \boldsymbol{\sigma} : \mathbf{D}\boldsymbol{\sigma} + \Xi(\xi) \quad (12) \\ \Xi(\xi) = \frac{1}{2} \xi K \xi $$

where $K$ is hardening modulus. With the specific choice of these functions, we can further appeal to the second principle

of thermodynamics (e.g. [35]) to define the damage dissipation ($0 \le \mathcal{D}^d$)

$$ \begin{aligned} \mathcal{D}^d &:= \boldsymbol{\sigma} : \dot{\boldsymbol{\epsilon}} + \frac{\partial}{\partial t} [\chi(\boldsymbol{\sigma}, \mathbf{D}, \xi) - \boldsymbol{\sigma} \cdot \boldsymbol{\epsilon}] + \theta \dot{s}^d \\ &= \dot{\boldsymbol{\epsilon}} : (-\boldsymbol{\epsilon} + \mathbf{D}\boldsymbol{\sigma}) + \frac{1}{2} \boldsymbol{\sigma} : \dot{\mathbf{D}}\boldsymbol{\sigma} + q \dot{\xi} + \theta \dot{s}^d \quad (13) \\ q &:= -\frac{d\Xi}{d\xi} = -K\xi \end{aligned} $$

It is easy to confirm that in the elastic case, with all internal variables frozen and damage dissipation equal to zero, from Eq. (13) above we can recover the strain tensor

$$ \left\{ \begin{array}{l} \dot{\mathbf{D}} = \mathbf{0} \\ \dot{\xi} = 0 \\ \dot{s}^d = 0 \\ \mathcal{D}^d = 0 \end{array} \right. \implies \boldsymbol{\epsilon} = \mathbf{D}\boldsymbol{\sigma} \quad (14) $$

Thus, for the inelastic case, we recover the positive value of damage dissipation

$$ 0 < \mathcal{D}^d := \frac{1}{2} \boldsymbol{\sigma} : \dot{\mathbf{D}}\boldsymbol{\sigma} + q \dot{\xi} + \underbrace{\theta \dot{s}^d}_{\mathcal{D}_{ther}} \quad (15) $$

The principle of maximum damage dissipation (e.g. [35]) can then be used to obtain the corresponding evolution equations of internal variables. This can be defined as the corresponding constrained evolution problem

$$ \begin{aligned} \min_{\phi(\boldsymbol{\sigma}, q, \theta) \le 0} & \quad [-\mathcal{D}^d] \Leftrightarrow \max_{\dot{\gamma} \ge 0} \min_{\forall(\boldsymbol{\sigma}, q, \theta)} L^d(\boldsymbol{\sigma}, q, \theta, \dot{\gamma}) \\ L^d(\boldsymbol{\sigma}, q, \theta, \dot{\gamma}) &:= -\mathcal{D}^d(\boldsymbol{\sigma}, q, \theta) + \dot{\gamma}\phi(\boldsymbol{\sigma}, q, \theta) \end{aligned} \quad (16) $$

The corresponding Kuhn-Tucker optimality conditions are enforced in order to obtain the corresponding evolution equations

$$ \begin{aligned} \dot{\mathbf{D}} &= \dot{\gamma} \frac{1}{\|\boldsymbol{\sigma}\|_{\mathbf{D}}} \frac{\partial \boldsymbol{\phi}}{\partial \boldsymbol{\sigma}} \otimes \frac{\partial \boldsymbol{\phi}}{\partial \boldsymbol{\sigma}} \\ \dot{\xi} &= \dot{\gamma} \frac{\partial \boldsymbol{\phi}}{\partial q} \\ \dot{s}^d &= \dot{\gamma} \frac{\partial \boldsymbol{\phi}}{\partial \theta} \end{aligned} \quad (17) $$

The expressions in Eq. (16) and (17) define rate-independent damage evolution equations. One can also obtain the rate-sensitive damage model by appealing to the penalized form of the principle of maximum damage dissipation with
$$
\begin{gathered}
\min_{\forall (\sigma^*, q^*, \theta^*)} \left[ -\mathcal{D}^d(\boldsymbol{\sigma}, q, \theta) + P(\phi(\boldsymbol{\sigma}^*, q^*, \theta^*)) \right] \\
P(\phi(\boldsymbol{\sigma}, q, \theta)) = \begin{cases} \frac{1}{2\eta} [\phi(\boldsymbol{\sigma}, q, \theta)]^2 & ; \phi(\cdot) > 0 \\ 0 & ; \phi(\cdot) \le 0 \end{cases} \tag{18}
\end{gathered}
$$

where η is the viscosity coefficient. The corresponding evo-
lution equations for internal variables can now be written
as

$$
\begin{align}
\dot{\mathbf{D}} &= \frac{\langle \phi \rangle}{\eta} \frac{1}{\|\boldsymbol{\sigma}\|_{\mathbf{D}}} \frac{\partial \phi}{\partial \boldsymbol{\sigma}} \otimes \frac{\partial \phi}{\partial \boldsymbol{\sigma}} \nonumber \\
\dot{\xi} &= \frac{\langle \phi \rangle}{\eta} \frac{\partial \phi}{\partial q} \tag{19} \\
\dot{s}^d &= \frac{\langle \phi \rangle}{\eta} \frac{\partial \phi}{\partial \theta} \nonumber
\end{align}
$$

where $\langle \cdot \rangle$ is the Macauley bracket, which is defined as

$$
\langle \phi \rangle := \frac{\phi(\cdot) + |\phi(\cdot)|}{2} = \begin{cases} \phi(\cdot); & \phi(\cdot) \ge 0 \\ 0; & \phi(\cdot) < 0 \end{cases} \quad (20)
$$

The system of local balance equations finally consists of
the force and moment balance equations and the energy con-
servation equation replacing those in Eqs. (8)$_1$ and (9)

$$
\left\{
\begin{array}{ll}
\rho \ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{b} ; & \boldsymbol{\sigma} = \boldsymbol{\sigma}^T \\
\theta \dot{s} = -\nabla \cdot \mathbf{q} + \mathcal{D}^d + r &
\end{array}
\right.
\qquad (21)
$$

Given the additive split of entropy with $s = s^e + s^d$, we can
compute the elastic entropy evolution:

$$
\begin{align}
s^e &= -\frac{\partial \psi}{\partial \theta} = \nonumber \\
    &\quad (\mathbf{D}^{-1} \cdot \boldsymbol{\alpha}) \cdot [\boldsymbol{\epsilon} - \boldsymbol{\alpha} (\theta - \theta_0)] + \rho c \ln \left(\frac{\theta}{\theta_0}\right) \tag{22} \\
    &\to \theta \dot{s}^e = \theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : (\dot{\boldsymbol{\epsilon}} - \boldsymbol{\alpha} \dot{\theta}) \nonumber \\
    &\quad - \theta \left( (\mathbf{D}^{-1} : \dot{\mathbf{D}} : \mathbf{D}^{-1}) : \boldsymbol{\alpha} \right) [\boldsymbol{\epsilon} - \boldsymbol{\alpha} (\theta - \theta_0)] \nonumber \\
    &\quad + \rho c \dot{\theta} \nonumber
\end{align}
$$

which replaced in Eq. (21)$_2$. This equation, combined with
Eq. (21)$_1$, results with the following linear and angular
momenta and energy balance equations:

$$
\left\{
\begin{array}{ll}
\rho \ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{b} ; & \boldsymbol{\sigma} = \boldsymbol{\sigma}^T \\
\underbrace{(\rho c - \theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \boldsymbol{\alpha})}_{\tilde{\rho}\tilde{c}} \dot{\theta} = -\nabla \cdot \mathbf{q} + & \\
\mathcal{D}_{\text{mech}} - \theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \dot{\boldsymbol{\epsilon}} + & \\
\theta ((\mathbf{D}^{-1} : \dot{\mathbf{D}} : \mathbf{D}^{-1}) : \boldsymbol{\alpha}) [\boldsymbol{\epsilon} - \boldsymbol{\alpha} (\theta - \theta_0)] + r. &
\end{array}
\right.
\qquad (23)
$$

<figure><img src="image_3.png" /><figcaption>Fig. 1 Thermo-mechanical framework with two-way coupling</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** The domain $\Omega \subset \mathbb{R}^3$ with boundary $\partial\Omega$ with imposed values for: a) mechanical fields and b) thermal fields</figcaption></figure>

where

$$
\tilde{\rho} \tilde{c} := \rho c - \theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \boldsymbol{\alpha} \quad (24)
$$

denotes the 'modified' heat conduction of the material, and

$$
\begin{equation}
\begin{split}
F(\epsilon, \mathbf{D}, \theta) = -\theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \dot{\boldsymbol{\epsilon}} & \\
& + \theta ((\mathbf{D}^{-1} : \dot{\mathbf{D}} : \mathbf{D}^{-1}) : \boldsymbol{\alpha}) [\epsilon - \boldsymbol{\alpha} (\theta - \theta_0)]
\end{split}
\tag{25}
\end{equation}
$$

represents the structural heating.

2.2 Variational formulation of fully coupled
thermodynamics

This section presents an extension of the regularized variational formulation [24, 26] initially proposed for mechanics to include transient heat conduction. The resulting model features both structural heating and temperature-affected mechanical fields, resulting in a two-way coupling for thermodynamics (see Fig. 1).

We assume that the solid body of interest occupies a
domain Ω ⊂ R³ with a boundary ∂Ω that can be divided
into four parts: the imposed displacement boundary ∂Ωᵢ,
the imposed stress boundary ∂Ωᵢ, the imposed temperature
boundary ∂Ωᵢ, and the imposed heat flux boundary ∂Ωⱼ. For
a well posed problem, we require that:

$$
\begin{align}
\partial \Omega &= \partial \Omega_u \cup \partial \Omega_t \quad \text{and} \quad \partial \Omega_u \cap \partial \Omega_t = \emptyset \tag{26} \\
\partial \Omega &= \partial \Omega_{\theta} \cup \partial \Omega_q \quad \text{and} \quad \partial \Omega_{\theta} \cap \partial \Omega_q = \emptyset \nonumber
\end{align}
$$

This is graphically illustrated in Fig. 2.

The state of the body is described by four fields: displace-
ment **u**(**x**, *t*), stress **σ**(**x**, *t*), temperature *θ*(**x**, *t*), and heat
flux **q**(**x**, *t*), each of which is considered independent. It is
important to note that the stress tensor is not assumed to be symmetric, but rather its symmetry is enforced through the corresponding moment equilibrium equation of the strong form of the problem. More precisely, the strong form of the central problem can be expressed $\forall \mathbf{x} \in \Omega$ as follows

(i) $\nabla \cdot \boldsymbol{\sigma} + \mathbf{b} = \rho \ddot{\mathbf{u}}$ in $\Omega$

(ii) $\tilde{\rho} \tilde{c} \dot{\theta} = -\nabla \cdot \mathbf{q} + \mathcal{D}_{\text{mech}} + F(\boldsymbol{\epsilon}, \mathbf{D}, \theta) + r$ in $\Omega$

(iii) $\text{symm } \boldsymbol{\sigma} = \mathbf{D}^{-1} \text{ symm } \nabla \mathbf{u} - (\theta - \theta_0) \boldsymbol{\beta}$ in $\Omega$

(iv) $\mathbf{q} = -k \nabla \theta = -k \mathbf{I} \nabla \theta$ in $\Omega$ \qquad (27)

(v) skew $\boldsymbol{\sigma} = \mathbf{0}$ in $\Omega$

(vi) $\mathbf{u} = \bar{\mathbf{u}}|_{\partial\Omega_u}$ and $\boldsymbol{\sigma}\mathbf{n} = \bar{\mathbf{t}}|_{\partial\Omega_t}$

(vii) $\theta = \bar{\theta}|_{\partial\Omega_{\theta}}$ and $\mathbf{q}\mathbf{n} = \bar{\mathbf{g}}|_{\partial\Omega_q}$

(viii) $\mathbf{u}(\mathbf{x}, 0) = \mathbf{u}_0(\mathbf{x}) = \bar{\mathbf{u}}$; $\theta(\mathbf{x}, 0) = \theta_0$ in $\Omega$

where $\boldsymbol{\beta} = \mathbf{D}^{-1}\boldsymbol{\alpha}$ is thermal isotropic stress tensor, $k = \text{diag}[k_x \ k_y \ k_z] = k \mathbf{I}_3$ with $k$ as the thermal conductivity coefficient derived for isotropic case and $\mathbf{I}_3$ is the identity tensor. The second-rank stress tensor can be decomposed into symmetric and skew-symmetric parts:

$$ \boldsymbol{\sigma} = \text{symm } \boldsymbol{\sigma} + \text{skew } \boldsymbol{\sigma} \qquad (28) $$

where

$$ \text{symm } \boldsymbol{\sigma} = \frac{1}{2} (\boldsymbol{\sigma} + \boldsymbol{\sigma}^T) \qquad (29) $$
$$ \text{skew } \boldsymbol{\sigma} = \frac{1}{2} (\boldsymbol{\sigma} - \boldsymbol{\sigma}^T) $$

Similarly, the displacement gradient can be split into a symmetric part, defining the infinitesimal strain tensor:

$$ \boldsymbol{\epsilon}(\mathbf{u}) := \text{symm } \nabla \mathbf{u} = \frac{1}{2} (\nabla \mathbf{u} + \nabla \mathbf{u}^T) \qquad (30) $$

and a skew-symmetric part, representing the infinitesimal rotation:

$$ \boldsymbol{\omega}(\mathbf{u}) := \text{skew } \nabla \mathbf{u} = \frac{1}{2} (\nabla \mathbf{u} - \nabla \mathbf{u}^T) \qquad (31) $$

We note that the scalar product between a symmetric second-order tensor $\boldsymbol{s}$ and a skew-symmetric second-order tensor $\boldsymbol{a}$ is equal to zero

$$ \begin{aligned} \text{skew } \boldsymbol{\sigma}^T &= - \text{skew } \boldsymbol{\sigma} \\ \text{symm } \boldsymbol{\sigma}^T &= \text{symm } \boldsymbol{\sigma} \\ \text{skew } \boldsymbol{\sigma} : \text{symm } \boldsymbol{\sigma} &= 0 \end{aligned} \qquad (32) $$

which can be exploited to simplify the variational formulation. Namely, the variational form of this boundary value problem can be defined by using the mixed thermomechanical functional given as:

$$ \begin{aligned} \Pi^{\text{Th-Mech}}(\mathbf{u}, \theta, \boldsymbol{\omega}, \boldsymbol{\sigma}, \mathbf{q}) = & -\frac{1}{2} \int_{\Omega} \text{symm } \boldsymbol{\sigma} : \mathbf{D} \text{ symm } \boldsymbol{\sigma} \, dV \\ & + \int_{\Omega} \boldsymbol{\sigma} : (\nabla \mathbf{u} - \boldsymbol{\omega}) \, dV - \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\alpha}(\theta - \theta_0) \, dV \\ & -\frac{1}{2} \int_{\Omega} \mathbf{q} \cdot \mathbf{k}^{-1} \mathbf{q} + \int_{\Omega} \mathbf{q} \cdot \nabla \theta \, dV \\ & - \Pi^{\text{ext}}(\mathbf{u}, \theta) \end{aligned} \qquad (33) $$

where $\mathbf{u} \in \mathcal{U}$, $\theta \in \mathcal{T}$, $\boldsymbol{\omega} \in \mathcal{W}$, $\boldsymbol{\sigma} \in \mathcal{S}$ and $\mathbf{q} \in \mathcal{Q}$ are respectively, trial field of displacement, rotation, stress, temperature and heat-flux. The variational statement in Eq. (33) requires that the stress tensor $\boldsymbol{\sigma}_t$, the infinitesimal rotation tensor $\boldsymbol{\omega}_t$, and the heat flux $\mathbf{q}_t$ should belong to the space of square-integrable functions over the domain $\Omega$, denoted as $L_2(\Omega)$ (which applies at any time $t$, but the time-dependence is further suppressed in order to alleviate notation). The same applies to their variation $\boldsymbol{\tau}$, $\boldsymbol{\psi}$ and $\tilde{\mathbf{q}}$ that should belong to the space of square-integrable functions $L_2(\Omega)$.

$$ \begin{aligned} \mathcal{W} &= \{\boldsymbol{\omega} \mid \boldsymbol{\omega} \in L_2(\Omega); \text{symm } \boldsymbol{\omega} = 0\} \\ \mathcal{S} &= \{\boldsymbol{\sigma} \mid \boldsymbol{\sigma} \in L_2(\Omega)\} \\ \mathcal{Q} &= \{\mathbf{q} \mid \mathbf{q} \in L_2(\Omega)\} \end{aligned} \qquad (34) $$

The set of trial displacements $\mathbf{u}$ and temperatures $\theta$ must satisfy a stronger regularity condition and belong to the Sobolev space $H^1(\Omega)$, with their generalized derivative that belong to $L_2(\Omega)$. Their variations of displacement $\mathbf{v}$ and of temperature $\vartheta$ belong to the subset denoted as $H_0^1(\Omega)$, which satisfies the homogeneous boundary condition on $\partial\Omega_u$ for displacements and $\partial\Omega_\theta$ for temperatures, which is written as follows:

$$ \begin{aligned} \mathcal{U} &:= \{\mathbf{u} \mid \mathbf{u} \in H^1(\Omega); \mathbf{u}|_{\partial\Omega_u} = \bar{\mathbf{u}}\}; \\ \mathcal{T} &:= \{\theta \mid \theta \in H^1(\Omega); \theta|_{\partial\Omega_\theta} = \bar{\theta}\}; \\ \mathcal{U}_0 &:= \{\mathbf{v} \mid \mathbf{v} \in H_0^1(\Omega); \mathbf{v}|_{\partial\Omega_u} = \mathbf{0}\}; \\ \mathcal{T}_0 &:= \{\vartheta \mid \vartheta \in H_0^1(\Omega); \vartheta|_{\partial\Omega_\theta} = 0\}; \end{aligned} \qquad (35) $$

The regularized functional is derived by adding an extra term to improve the ellipticity of the starting functional:

$$ \Pi^{\text{Th-Mech}}(\mathbf{u}, \theta, \boldsymbol{\omega}, \boldsymbol{\sigma}, \mathbf{q}) = -\frac{1}{2} \int_{\Omega} \text{symm } \boldsymbol{\sigma} : \mathbf{D} \text{ symm } \boldsymbol{\sigma} \, dV $$
$$
\begin{aligned}
&+ \int_{\Omega} \boldsymbol{\sigma} : (\nabla \mathbf{u} - \boldsymbol{\omega}) \mathrm{d}V - \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\alpha}(\theta - \theta_0) \mathrm{d}V \\
&- \frac{1}{2} \int_{\Omega} \mathbf{q} \cdot \mathbf{k}^{-1} \mathbf{q} + \int_{\Omega} \mathbf{q} \cdot \nabla \theta \mathrm{d}V \\
&- \frac{1}{2} \bar{\gamma}^{-1} \int_{\Omega} \| \text{skew } \boldsymbol{\sigma} \|^2 \mathrm{d}V - \Pi^{\text{ext}}(\mathbf{u}, \theta)
\end{aligned}
\qquad (36)
$$

where $\bar{\gamma}$ is a regularization parameter that can be taken equal to the shear modulus [32]. This kind of regularization in the variational statement allows for flexibility in choosing the finite element discrete approximation which is the most appropriate for the problem on hand. To fully describe the proposed boundary value problem, it is necessary to consider the external virtual work, $G^{\text{ext}}(\mathbf{v})$ and $G(\theta)$ which account for thermo-plastic coupling and the thermal energy rate. The external virtual work, $G^{\text{ext}}(\mathbf{v})$, is derived in accordance with d'Alembert's principle, accounting for the inertial forces term that arise in dynamic problems

$$ G^{\text{ext}}(\mathbf{v}) = - \int_{\Omega} \rho \ddot{\mathbf{u}} \cdot \mathbf{v} \mathrm{d}V + \int_{\Omega} \mathbf{b} \cdot \mathbf{v} \mathrm{d}V \\ + \int_{\partial \Omega_t} \mathbf{t} \cdot \mathbf{v} \mathrm{d}A \qquad (37) $$

The term $G(\vartheta)$ can be written as

$$ G(\vartheta) = \int_{\Omega} \rho c \dot{\vartheta} \vartheta \mathrm{d}V - \int_{\Omega} r \vartheta \mathrm{d}V - \int_{\Omega} \mathcal{D}_{\text{mech}} \vartheta \mathrm{d}V \\ - \int_{\Omega} F(\epsilon, \mathbf{D}, \theta) \vartheta \mathrm{d}V \qquad (38) $$

The weak form equations with respect to the displacement, stress, rotation, temperature, and heat-flux fields to obtain are derived by taking variations of the regularized variational functional.

$$
\begin{aligned}
G_{\mathbf{u}}(\mathbf{u}, \theta, \omega, \sigma, \mathbf{q}; \mathbf{v}) &:= \int_{\partial \Omega_t} \mathbf{v} \cdot (\boldsymbol{\sigma} \cdot \mathbf{n} - \bar{\mathbf{t}}) \mathrm{d}A \\
&\quad - \int_{\Omega} \mathbf{v} \cdot (\nabla \cdot \boldsymbol{\sigma} + \mathbf{b} - \rho \ddot{\mathbf{u}}) \mathrm{d}V = 0 \\
G_{\theta}(\mathbf{u}, \theta, \omega, \sigma, \mathbf{q}; \vartheta) &:= - \int_{\Omega} \vartheta \boldsymbol{\sigma} : \mathbf{1} \boldsymbol{\alpha} \mathrm{d}V \\
&\quad + \int_{\Omega} \vartheta \rho c \dot{\theta} \mathrm{d}V - \int_{\Omega} \vartheta \theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \boldsymbol{\alpha} \dot{\theta} \mathrm{d}V \\
&\quad - \int_{\Omega} \vartheta \mathcal{D}_{\text{mech}} \mathrm{d}V + \int_{\Omega} \vartheta \theta (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \dot{\boldsymbol{\epsilon}} \mathrm{d}V \\
&\quad + \int_{\Omega} \vartheta \theta \left( (\mathbf{D}^{-1} : \dot{\mathbf{D}} : \mathbf{D}^{-1}) : \boldsymbol{\alpha} \right) : [\boldsymbol{\epsilon} - \boldsymbol{\alpha}(\theta - \theta_0)] \mathrm{d}V \\
&\quad + \int_{\Omega} \nabla \vartheta \mathbf{q} \mathrm{d}V - \int_{\Omega} r \vartheta \mathrm{d}V \\
&\quad + \int_{\partial \Omega_q} \vartheta (\mathbf{q} \cdot \mathbf{n} - \bar{\mathbf{g}}) \mathrm{d}A = 0
\end{aligned}
$$

$$
\begin{aligned}
G_{\omega}(\mathbf{u}, \theta, \omega, \sigma, \mathbf{q}; \psi) &:= \int_{\Omega} \psi : \text{skew } \boldsymbol{\sigma} \mathrm{d}V = 0 \\
G_{\sigma}(\mathbf{u}, \theta, \omega, \sigma, \mathbf{q}; \tau) &:= \\
&\quad - \int_{\Omega} \text{symm } \boldsymbol{\tau} : \mathbf{D} \text{symm } \boldsymbol{\sigma} \mathrm{d}V \\
&\quad - \int_{\Omega} \text{symm } \boldsymbol{\tau} : (-\text{symm } \nabla \mathbf{u} + \mathbf{alpha}(\theta - \theta_0)) \mathrm{d}V \\
&\quad + \int_{\Omega} \boldsymbol{\tau} : (\text{skew } \nabla \mathbf{u} - \omega - \bar{\gamma}^{-1} \text{skew } \boldsymbol{\sigma}) \mathrm{d}V = 0 \\
G_{q}(\mathbf{u}, \theta, \omega, \sigma, \mathbf{q}; \tilde{\mathbf{q}}) &:= - \int_{\Omega} \tilde{\mathbf{q}} \cdot (\mathbf{k}^{-1} \mathbf{q}) \mathrm{d}V \\
&\quad + \int_{\Omega} \tilde{\mathbf{q}} \cdot \nabla \theta \mathrm{d}V = 0
\end{aligned}
\qquad (39)
$$

where the result in Eq. (39) is obtained by applying the Gauss divergence theorem. By considering the arbitrary variations of displacement $\mathbf{v}$, temperature $\vartheta$, stress $\boldsymbol{\tau}$, rotation $\psi$ and heat flux $\tilde{\mathbf{q}}$, we can obtain all the Euler-Lagrange equations corresponding to the strong form stated previously. In addition, the regularization also allows us to recover an extra relationship from Eq. (39)(4), showing that:

$$ \omega = \text{skew } \nabla \mathbf{u} - \bar{\gamma}^{-1} \text{skew } \boldsymbol{\sigma} \qquad (40) $$

We can further use the best result in Eq. (40) to eliminate the rotation field in the mixed variational formulation in Eq. (36), and thus to recover a four-field hybrid-stress-heat-flux variational formulation with displacement, temperature, stress, and heat-flux as independent fields:

$$
\begin{aligned}
\Pi^{\text{Th-Mech}}(\mathbf{u}, \theta, \omega, \sigma, \mathbf{q}) = & -\frac{1}{2} \int_{\Omega} \text{symm } \boldsymbol{\sigma} : \mathbf{D} \text{symm } \boldsymbol{\sigma} \mathrm{d}V \\
& + \int_{\Omega} \text{symm } \boldsymbol{\sigma} : \text{symm } \nabla \mathbf{u} \mathrm{d}V - \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\alpha}(\theta - \theta_0) \mathrm{d}V \qquad (41) \\
& +\frac{1}{2} \bar{\gamma}^{-1} \int_{\Omega} \| \text{skew } \boldsymbol{\sigma} \|^2 \mathrm{d}V - \frac{1}{2} \int_{\Omega} \mathbf{q} \cdot \mathbf{k}^{-1} \mathbf{q} \\
& + \int_{\Omega} \mathbf{q} \cdot \nabla \theta \mathrm{d}V - \Pi^{\text{ext}}(\mathbf{u}, \theta)
\end{aligned}
$$

The corresponding variational equation $\Pi_{\gamma}^{\text{Th-Mech}}(\mathbf{u}, \theta, \sigma, \mathbf{q})$ with respect to the trial displacement field $\mathbf{u}$ can then be expressed

$$
\begin{aligned}
G_{\mathbf{u}}(\mathbf{u}, \theta, \sigma, \mathbf{q}; \mathbf{v}) &:= \int_{\Omega} \boldsymbol{\sigma} : \text{symm } \nabla \mathbf{v} \mathrm{d}V \\
&\quad - \int_{\Omega} (\mathbf{b} - \rho \ddot{\mathbf{u}}) \cdot \mathbf{v} \mathrm{d}V - \int_{\partial \Omega_t} \bar{\mathbf{t}} \cdot \mathbf{v} \mathrm{d}A = 0
\end{aligned}
\qquad (42)
$$

Similarly, by using Eq. (32) we can obtain the variational equation of the functional $\Pi_{\gamma}^{\text{Th-Mech}}(\mathbf{u}, \theta, \sigma, \mathbf{q})$ with respect to the trial stress field $\boldsymbol{\tau}$
$$
\begin{align}
G_{\sigma}(\mathbf{u}, \theta, \sigma, \mathbf{q}; \boldsymbol{\tau}) &= \int_{\Omega} [\text{symm } \boldsymbol{\tau} : (\text{symm } \nabla \mathbf{u} - \boldsymbol{\alpha}(\theta - \theta_0))] dV \nonumber \\
&\quad - \int_{\Omega} [\text{symm } \boldsymbol{\tau} : (\mathbf{D} \text{ symm } \boldsymbol{\sigma})] dV \nonumber \\
&\quad + \text{skew } \boldsymbol{\tau} : \frac{1}{\bar{\gamma}} \text{ skew } \boldsymbol{\sigma} dV = 0 \tag{43} \\
&= \int_{\Omega} \boldsymbol{\tau} : [\text{symm } \nabla \mathbf{u} - \boldsymbol{\alpha}(\theta - \theta_0)] dV \nonumber \\
&\quad - \int_{\Omega} \boldsymbol{\tau} : \left[ \underbrace{\left( \mathbf{D} \text{ symm } \boldsymbol{\sigma} - \frac{1}{\bar{\gamma}} \text{ skew } \boldsymbol{\sigma} \right)}_{\mathbf{D}\boldsymbol{\sigma}} \right] dV = 0 \nonumber
\end{align}
$$

Here an optimal choice $\bar{\gamma} = 2\mu$, where $\mu$ is the Lamé parameter defined as a function of the elastic modulus and Poisson's ratio: $\frac{E}{2(1+\nu)}$, is used.

### 3 Discrete approximations enforcing stress and heat flux continuity

Here, we propose a special discrete approximation for the displacement, temperature, stress, and heat-flux fields using a conventional four-node tetrahedron element geared towards non-stationary evolution. Namely, the interpolation for the heat flux and stress fields are chosen in the lowest order Raviart-Thomas vector space to ensure continuity of fields across element boundaries. This requires added nodes to the faces of the element, resulting in an enriched element.

The development of the enriched element was primarily driven by the need for improved performance in dynamic simulations, which require fully discretizing the problem in a dynamic framework with a time-integration scheme. To achieve superior performance, the proposed scheme preserves the total energy of the system in free vibrations for the elastic case. This involves evaluating all fields at the time step $t_{n+1/2}$ for linear thermo-elastodynamics that satisfy the force and thermal balance equations. The same scheme can be extended to energy-decaying with the corresponding evolution equations for damage behavior, making it well-suited for a wide range of applications.

#### 3.1 Interpolation of displacement and temperature fields

For 3D problems, we adopt a discrete approximation method using a conventional four-node tetrahedron element as the starting point to approximate the displacement and temperature fields, see Fig. 3. The lowest-order linear shape functions

<figure><img src="image_5.png" /><figcaption>Fig. 3 4-node isoparametric tetrahedron element</figcaption></figure>

<table><thead><tr><th>Nodes</th></tr></thead><tbody><tr><td>1 (0,0,0)</td></tr><tr><td>2 (1,0,0)</td></tr><tr><td>3 (0,1,0)</td></tr><tr><td>4 (0,0,1)</td></tr></tbody></table>

N are utilized and expressed in natural coordinates, as illustrated in A, in Fig. 24.

$$
\begin{align}
N_1 &= 1 - \xi - \eta - \zeta \nonumber \\
N_2 &= \xi \nonumber \\
N_3 &= \eta \nonumber \\
N_4 &= \zeta \nonumber \\
0 &\le \xi, \eta, \zeta \le 1 \tag{44}
\end{align}
$$

The discrete approximation for the displacement field is written in terms of nodal displacement

$$
\begin{align}
\mathbf{d} \approx \mathbf{d}^h (\xi, \eta, \zeta) &= \begin{bmatrix} \mathbf{d}_x^h \\ \mathbf{d}_y^h \\ \mathbf{d}_z^h \end{bmatrix} = \sum_{i=1}^{4} N_i (\xi, \eta, \zeta) \begin{bmatrix} d_x^i \\ d_y^i \\ d_z^i \end{bmatrix} \tag{45} \\
\mathbf{d} &= \mathbf{N}_d d_d
\end{align}
$$

The strain components can be expressed as the symmetric part of the displacement gradients as:

$$
\begin{align}
\epsilon \approx \epsilon^h (\xi, \eta, \zeta) &= \sum_{i=1}^{4} \mathbf{B}_{d,i}(\xi, \eta, \zeta) \begin{bmatrix} d_x^i \\ d_y^i \\ d_z^i \end{bmatrix} \tag{46} \\
\epsilon &= \mathbf{B}_d d_d
\end{align}
$$

where

$$
\mathbf{B}_{d,i}(\xi, \eta, \zeta) = \begin{bmatrix} \frac{\partial N_i}{\partial x} & 0 & 0 \\ 0 & \frac{\partial N_i}{\partial y} & 0 \\ 0 & 0 & \frac{\partial N_i}{\partial z} \\ \frac{\partial N_i}{\partial y} & \frac{\partial N_i}{\partial x} & 0 \\ 0 & \frac{\partial N_i}{\partial z} & \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z} & 0 & \frac{\partial N_i}{\partial x} \end{bmatrix} \quad (47)
$$
In similar way, we can derive

$$ 
\begin{aligned}
\dot{\epsilon} &\cong \dot{\epsilon}^h (\xi, \eta, \zeta) = \sum_{i=1}^{4} \mathbf{B}_{d,i} (\xi, \eta, \zeta) \begin{bmatrix} v_x^i \\ v_y^i \\ v_z^i \end{bmatrix} \\
\epsilon &= \mathbf{B}_d v_d
\end{aligned} 
\quad (48) $$

The corresponding acceleration field required for dynamics can be obtained from the nodal values of acceleration, by using the same shape functions as those used for the displacement field

$$ 
\begin{aligned}
\ddot{\mathbf{d}} &\cong \ddot{\mathbf{d}}^h (\xi, \eta, \zeta) = \begin{bmatrix} \ddot{\mathbf{d}}_x^h \\ \ddot{\mathbf{d}}_y^h \\ \ddot{\mathbf{d}}_z^h \end{bmatrix} = \sum_{i=1}^{4} N_i (\xi, \eta, \zeta) \begin{bmatrix} a_x^i \\ a_y^i \\ a_z^i \end{bmatrix} \\
\ddot{\mathbf{d}} &= \mathbf{N}_d a_d
\end{aligned} 
\quad (49) $$

The discrete approximation of temperature at any point within the element is expressed as a linear combination of nodal temperature values

$$ \theta = \sum_{i=1}^{4} N_i (\xi, \eta, \zeta) \theta_i \mapsto \theta = \mathbf{N}_{\theta} \theta_d \quad (50) $$

Using the nodal values of temperature, the approximation of temperature gradients can be expressed as the following equation:

$$ 
\begin{aligned}
\nabla \theta &\cong \nabla \theta^h (\xi, \eta, \zeta) = \begin{bmatrix} \nabla \theta_x^h \\ \nabla \theta_y^h \\ \nabla \theta_z^h \end{bmatrix} = \sum_{i=1}^{4} \mathbf{B}_{\theta,i} (\xi, \eta, \zeta) \theta_i \\
\nabla \theta &= \mathbf{B}_{\theta} \theta_d
\end{aligned} 
\quad (51) $$

where

$$ \mathbf{B}_{\theta,i} (\xi, \eta, \zeta) = \begin{bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z} \end{bmatrix} \quad (52) $$

## 3.2 Semi-discrete approximations enforcing stress and heat flux continuity

The nodal displacement $\mathbf{d} = (d_x^i, d_y^i, d_z^i)^T$ and temperature $\theta_i$ are assigned to the degrees of freedom of the element's corner nodes 1, 2, 3 and 4. In order to represent the stress and heat-flux fields, we use the lowest-order Raviart-Thomas vector space ($RT_0$) and introduce on 4 element's faces 4 extra nodes labeled as 11, 12, 13 and 14 in Fig. 4(a). These nodes are used to store the degrees of freedom of the stress vector components $\tau = (\tau_x^i, \tau_y^i, \tau_z^i)$ and the normal component of

<figure><img src="image_6.png" /><figcaption>(a) Element nodes</figcaption></figure>

<table><thead><tr><th>Corner nodes</th><th>Edge nodes</th></tr></thead><tbody><tr><td>1 (0,0,0)</td><td>5 (1/2,0,0)</td></tr><tr><td>2 (1,0,0)</td><td>6 (1/2,1/2,0)</td></tr><tr><td>3 (0,1,0)</td><td>7 (0,1/2,0)</td></tr><tr><td>4 (0,0,1)</td><td>8 (0,0,1/2)</td></tr><tr><td></td><td>9 (1/2,0,1/2)</td></tr><tr><td></td><td>10 (0,1/2,1/2)</td></tr></tbody><thead><tr><th>Face nodes</th><th>Interior nodes</th></tr></thead><tbody><tr><td>11 (1/3,0,1/3)</td><td>15 (1/3,1/3,1/3)</td></tr><tr><td>12 (1/3,1/3,1/3)</td><td></td></tr><tr><td>13 (0,1/3,1/3)</td><td></td></tr><tr><td>14 (1/3,1/3,0)</td><td></td></tr></tbody></table>

(b) Description of nodes and its coordinates

<caption>**Fig. 4** Enhanced tetrahedron element and the corresponding degrees of freedom for displacement and temperature at corner nodes, and stress and heat-flux at element's faces</caption>

the heat flux $g_i$. This setup is designed to utilize the geometry of a 15-node tetrahedral element, implemented in the research version of finite element software FEAP, developed by RL Taylor [39]. The coordinates of all the nodes in the parent element are shown in Fig. 4b.

Since the edge nodes and the interior node do not store any degrees of freedom for the proposed boundary value problem, these nodes are deactivated. Consequently, the analysis focuses on a hybrid element with standard corner nodes and face nodes, with the degrees of freedom graphically illustrated in Fig. 5.

The basis functions of the Raviart-Thomas vector space, which are illustrated in B, in Figure 25, can be defined as follows:

$$ \hat{\mathbf{\Phi}}^1 (\xi, \eta, \zeta) = \begin{bmatrix} \phi_1^1 \\ \phi_2^1 \\ \phi_3^1 \end{bmatrix} = \sqrt{3} \begin{bmatrix} \xi \\ \eta \\ \zeta \end{bmatrix}; $$
<figure><img src="image_7.png" /><figcaption>Fig. 5 Degrees of freedom for the proposed hybrid tetrahedron model</figcaption></figure>

$$
\begin{align}
\hat{\mathbf{\Phi}}^2(\xi, \eta, \zeta) &= \begin{bmatrix} \phi_1^2 \\ \phi_2^2 \\ \phi_3^2 \end{bmatrix} = \begin{bmatrix} \xi - 1 \\ \eta \\ \zeta \end{bmatrix}; \\
\hat{\mathbf{\Phi}}^3(\xi, \eta, \zeta) &= \begin{bmatrix} \phi_1^3 \\ \phi_2^3 \\ \phi_3^3 \end{bmatrix} = \begin{bmatrix} \xi \\ \eta - 1 \\ \zeta \end{bmatrix}; \\
\hat{\mathbf{\Phi}}^4(\xi, \eta, \zeta) &= \begin{bmatrix} \phi_1^4 \\ \phi_2^4 \\ \phi_3^4 \end{bmatrix} = \begin{bmatrix} \xi \\ \eta \\ \zeta - 1 \end{bmatrix}; \tag{53}
\end{align}
$$

One can show that the RT₀ normal basis functions possess the independence property for any midpoint $\hat{\mathbf{p}}_i$ on face $j \in [1, 2, 3, 4]$. Hence, the normal basis function $\hat{\mathbf{\Phi}}^i$ related

to face $i$ remains unaffected by the normal basis functions related to the other faces, which can be written as

$$ \hat{\mathbf{\Phi}}^i(\hat{\mathbf{p}}_j) \cdot \hat{\mathbf{n}}_i = \begin{cases} 1, & \text{if } j = i \\ 0, & \text{if } j \neq i \end{cases}; \quad i, j \in [1, 2, 3, 4] \quad (54) $$

where $\hat{\mathbf{n}}_i$ represents a unit vector of exterior normal at the face $i \in [1, 2, 3, 4]$, placed opposite node $i$ as shown in Fig. 6.

The Piola transformation is employed to transform the Raviart-Thomas normal basis functions $\hat{\mathbf{\Phi}}^i$ from the local coordinate system to the global coordinate system, which is represented by an affine function $\mathcal{F}$ and defined as follows:

$$
\begin{align}
\hat{\mathbf{\Phi}}^i &\mapsto \hat{\mathbf{\Phi}}^i(x) = \mathcal{F}(\hat{\mathbf{\Phi}}^i)(x) := \frac{|e_k|}{|\hat{e}_k|} \frac{\mathbf{J}}{|\mathbf{J}|} \hat{\mathbf{\Phi}}^i(\hat{x}) \tag{55} \\
\mathbf{\Phi}^i(\mathbf{p}_j) \cdot \mathbf{n}_i &= \begin{cases} 1, & \text{if } j = i \\ 0, & \text{if } j \neq i \end{cases} \nonumber
\end{align}
$$

Here, the terms $|e_k|$ and $|\hat{e}_k|$ are respectively the length of considered edges in the physical and parent space, and $\mathbf{J}$ is the Jacobian matrix defined by

$$ \mathbf{J} = \begin{bmatrix} \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial x} & \frac{\partial N_4}{\partial x} \\ \frac{\partial N_1}{\partial y} & \frac{\partial N_2}{\partial y} & \frac{\partial N_3}{\partial y} & \frac{\partial N_4}{\partial y} \\ \frac{\partial N_1}{\partial z} & \frac{\partial N_2}{\partial z} & \frac{\partial N_3}{\partial z} & \frac{\partial N_4}{\partial z} \end{bmatrix} \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ x_4 & y_4 & z_4 \end{bmatrix} \quad (56) $$

while $|\mathbf{J}|$ is its determinant. The independence property of the Raviart-Thomas basis functions in parent space is maintained in the physical space by the Piola transformation.

The Raviart-Thomas basis functions can be used to express the stress field as a linear combination of stress vector components entry of each element face.

$$
\begin{align}
\begin{bmatrix} \sigma_{11}^h \\ \sigma_{12}^h \\ \sigma_{13}^h \end{bmatrix} &= \sigma_1^h|_{\Omega^e} = \tau_1^1 \mathbf{\Phi}^1 + \tau_1^2 \mathbf{\Phi}^2 + \tau_1^3 \mathbf{\Phi}^3 + \tau_1^4 \mathbf{\Phi}^4 \nonumber \\
\begin{bmatrix} \sigma_{21}^h \\ \sigma_{22}^h \\ \sigma_{23}^h \end{bmatrix} &= \sigma_2^h|_{\Omega^e} = \tau_2^1 \mathbf{\Phi}^1 + \tau_2^2 \mathbf{\Phi}^2 + \tau_2^3 \mathbf{\Phi}^3 + \tau_2^4 \mathbf{\Phi}^4 \tag{57} \\
\begin{bmatrix} \sigma_{31}^h \\ \sigma_{32}^h \\ \sigma_{33}^h \end{bmatrix} &= \sigma_3^h|_{\Omega^e} = \tau_3^1 \mathbf{\Phi}^1 + \tau_3^2 \mathbf{\Phi}^2 + \tau_3^3 \mathbf{\Phi}^3 + \tau_3^4 \mathbf{\Phi}^4 \nonumber
\end{align}
$$

By rearranging the terms on the right-hand side, the last result can be expressed in a more compact form as follows:

$$ \sigma_{n+1/2}^h|_{\Omega^e} = S_{\tau} \tau_{n+1/2}^h|_{\Omega^e} \quad (58) $$

where $S_{\tau}$ denotes Raviart-Thomas matrix for stress interpolation. The fully developed form of this matrix is written in
<figure><img src="image_8.png" /><figcaption>**Fig. 6** The mapping $\mathcal{F}$ between a local isoparametric element and a global tetrahedron element</figcaption></figure>

C. By using Eq. (55) and (57), we can write the projection of the stress tensor onto a unit normal vector $\mathbf{n}_i$ at face $i$ as follows:

$$ \mathbf{t}_i^h = \boldsymbol{\sigma}^h \cdot \mathbf{n}_i^h = \begin{bmatrix} \boldsymbol{\sigma}_1^h \cdot \mathbf{n}_i^h \\ \boldsymbol{\sigma}_2^h \cdot \mathbf{n}_i^h \\ \boldsymbol{\sigma}_3^h \cdot \mathbf{n}_i^h \end{bmatrix} = \begin{bmatrix} \tau_1^i \\ \tau_2^i \\ \tau_3^i \end{bmatrix} \quad (59) $$

This projection defines the stress vector $\boldsymbol{\tau} = (\tau_x^i, \tau_y^i, \tau_z^i)$ at the face $i$ and the nodal degrees of freedom of stress interpolations. Since the degrees of freedom are shared by adjacent tetrahedron elements, the chosen discrete approximation of the stress field ensures continuity of the stress vector across the proposed enhanced tetrahedron element boundaries. This continuity is crucial in dynamics to prevent undesirable stress and heat flux jumps, which can compromise smooth wave propagation.

The same principle of ensuring boundary continuity can also be applied to the heat flux normal component $g^i$, by using

$$ \mathbf{q}^h|_{\Omega^e} = \begin{bmatrix} q_x^h \\ q_y^h \\ q_z^h \end{bmatrix} = g^1 \boldsymbol{\Phi}^1 + g^2 \boldsymbol{\Phi}^2 + g^3 \boldsymbol{\Phi}^3 + g^4 \boldsymbol{\Phi}^4 \quad (60) $$

which can be recast

$$ \mathbf{q}_{n+1/2}^h|_{\Omega^e} = \mathbf{S}_g \mathbf{g}_{n+1/2}^h|_{\Omega^e} \quad (61) $$

where $\mathbf{S}_g$ is Raviart-Thomas matrix for heat-flux interpolation and is written explicitly in terms of basis functions of the Raviart-Thomas vector space in C. The projection of the

heat flux field onto a normal vector $\mathbf{n}_i$ at face $i$ is heat flux degree of freedom, yielding the relation

$$ g_i^h = \mathbf{q}^h \cdot \mathbf{n}_i^h \quad (62) $$

The nodal degrees of freedom $g_i$ are identical on the corresponding face $i$, which is shared by a neighbouring enhanced tetrahedron elements. Thus, it follows that the lowest order Raviart-Thomas interpolation ($RT_0$) also guarantees the continuity of heat flux field projection $g^h$ across the faces of each tetrahedron element.

### 3.3 Discrete approximation of the weak form

The discrete weak form for a typical element $\Omega^e$ can be expressed in matrix notation by using the discrete approximations chosen for displacement, stress, temperature, and heat flux fields, as shown below:

$$ \begin{align*} G_{\mathbf{u}}^{e}(\mathbf{u}, \boldsymbol{\sigma}; \boldsymbol{\nu}) &:= \delta\boldsymbol{\nu}^{T} \int_{\Omega^{e}} \rho \mathbf{N}_{d}^{T} \mathbf{N}_{d} \mathbf{a} \, dV \\ &\quad + \delta\boldsymbol{\nu}^{T} \int_{\Omega^{e}} \mathbf{B}_{d}^{T} S_{\boldsymbol{\tau}} \boldsymbol{\tau} \, dV - \delta\boldsymbol{\nu}^{T} \int_{\Omega^{e}} \mathbf{N}_{d}^{T} \mathbf{b} \, dV \\ &\quad - \delta\boldsymbol{\nu}^{T} \int_{\partial\Omega_{i}^{e}} \mathbf{N}_{d}^{T} \boldsymbol{t} \, dA = 0 \\ G_{\theta}^{e}(\mathbf{u}, \theta, \boldsymbol{\sigma}, \mathbf{q}; \vartheta) &:= \delta\vartheta^{T} \int_{\Omega^{e}} \rho c \mathbf{N}_{\theta}^{T} \mathbf{N}_{\theta} \hat{\boldsymbol{\theta}} \, dV \\ &\quad - \delta\vartheta^{T} \int_{\Omega^{e}} \mathbf{N}_{\theta}^{T} \left[ (\mathbf{N}_{\theta} \hat{\boldsymbol{\theta}}) (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \boldsymbol{\alpha} \right] \mathbf{N}_{\theta} \hat{\boldsymbol{\theta}} \, dV \\ &\quad + \delta\vartheta^{T} \int_{\Omega^{e}} \mathbf{N}_{\theta}^{T} (\mathbf{D}^{-1} : \boldsymbol{\alpha}) (\mathbf{B}_{d} \mathbf{v}) \mathbf{N}_{\theta} \hat{\boldsymbol{\theta}} \, dV \end{align*} $$
$$
\begin{align*}
& -\delta\vartheta^T \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\hat{\mathbf{D}} : \mathbf{alpha}) (\mathbf{B}_d \mathbf{d}) \mathbf{N}_{\theta} \hat{\theta} \, \mathrm{d}V \right] \\
& +\delta\vartheta^T \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\hat{\mathbf{D}} : \boldsymbol{\alpha}) (\boldsymbol{\alpha} (\mathbf{N}_{\theta} \hat{\theta} - \theta_0)) \right] \mathbf{N}_{\theta} \hat{\theta} \, \mathrm{d}V \\
& -\delta\vartheta^T \int_{\Omega^e} \mathbf{N}_{\theta}^T \boldsymbol{\alpha} \mathbf{S}_{\tau} \boldsymbol{\tau} \, \mathrm{d}V + \delta\vartheta^T \int_{\Omega^e} \mathbf{B}_{\theta}^T \mathbf{S}_g \mathbf{g} \, \mathrm{d}V \\
& -\delta\vartheta^T \int_{\Omega^e} \mathbf{N}_{\theta}^T r \, \mathrm{d}V - \delta\vartheta^T \int_{\Omega^e} \mathbf{N}_{\theta}^T \mathcal{D}_{\mathrm{mech}} \, \mathrm{d}V \\
& -\delta\vartheta^T \int_{\partial\Omega_q^e} \mathbf{N}_{\theta}^T \bar{\mathbf{g}} \, \mathrm{d}A = 0
\end{align*}
$$

$$
\begin{align*}
G_{\sigma}^{e}(\mathbf{u}, \theta, \boldsymbol{\sigma}; \boldsymbol{\iota}) &:= \delta \boldsymbol{\iota}^{T} \int_{\Omega^{e}} S_{\boldsymbol{\tau}}^{T} \mathbf{B}_{d} \mathbf{d} \, \mathrm{d}V \\
&\quad - \delta \boldsymbol{\iota}^{T} \int_{\Omega^{e}} S_{\boldsymbol{\tau}}^{T} \boldsymbol{\alpha} (\mathbf{N}_{\theta} \hat{\theta} - \theta_{0}) \, \mathrm{d}V \\
&\quad - \delta \boldsymbol{\iota}^{T} \int_{\Omega^{e}} S_{\boldsymbol{\tau}}^{T} \mathbf{D} \mathbf{S}_{\boldsymbol{\tau}} \boldsymbol{\tau} \, \mathrm{d}V = 0
\end{align*}
$$

$$
\begin{align}
G_{\mathbf{q}}^{e}(\theta, \mathbf{q}; \tilde{\mathbf{q}}) &:= \delta \tilde{\mathbf{q}}^{T} \int_{\Omega^{e}} S_{\mathbf{g}}^{T} \mathbf{B}_{\theta} \hat{\theta} \, \mathrm{d}V \nonumber \\
&\quad - \delta \tilde{\mathbf{q}}^{T} \int_{\Omega^{e}} S_{\mathbf{g}}^{T} \mathbf{k}^{-1} S_{\mathbf{g}} \mathbf{g} \, \mathrm{d}V = 0 \tag{63}
\end{align}
$$

Here we introduce: $\hat{\mathbf{D}} = \mathbf{D}^{-1} : \hat{\mathbf{D}} : \mathbf{D}^{-1}$. By considering that the nodal values of virtual fields can be picked up arbitrarily, it is possible to obtain from Eq. (63) the final set of non-linear residual equations that need to be solved:

$$
\begin{align}
\mathcal{R}_{\mathbf{u}} &= \int_{\Omega^e} \rho \mathbf{N}_d^T \mathbf{N}_d \mathbf{a} \, \mathrm{d}V + \int_{\Omega^e} \mathbf{B}_d^T S_{\boldsymbol{\tau}} \boldsymbol{\tau} \, \mathrm{d}V \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_d^T \mathbf{f} \, \mathrm{d}V - \int_{\partial \Omega_t^e} \mathbf{N}_d^T \mathbf{t} \, \mathrm{d}A \nonumber \\
\mathcal{R}_{\theta} &= \int_{\Omega^e} \rho c \mathbf{N}_{\theta}^T \mathbf{N}_{\theta} \dot{\hat{\theta}} \, \mathrm{d}V \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\mathbf{N}_{\theta} \hat{\theta}) (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \boldsymbol{\alpha} \right] \mathbf{N}_{\theta} \dot{\hat{\theta}} \, \mathrm{d}V \nonumber \\
&\quad + \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\mathbf{D}^{-1} : \boldsymbol{\alpha}) (\mathbf{B}_d \mathbf{v}) \right] \mathbf{N}_{\theta} \hat{\theta} \, \mathrm{d}V \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\hat{\mathbf{D}} : \boldsymbol{\alpha}) (\mathbf{B}_d \mathbf{d}) \right] \mathbf{N}_{\theta} \hat{\theta} \, \mathrm{d}V \nonumber \\
&\quad + \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\hat{\mathbf{D}} : \boldsymbol{\alpha}) (\boldsymbol{\alpha} (\mathbf{N}_{\theta} \hat{\theta} - \theta_0)) \right] \mathbf{N}_{\theta} \hat{\theta} \, \mathrm{d}V \nonumber \\
&\quad + \left( - \int_{\Omega^e} \mathbf{N}_{\theta}^T \boldsymbol{\alpha} \mathbf{S}_{\boldsymbol{\tau}} \boldsymbol{\tau} \, \mathrm{d}V + \int_{\Omega^e} \mathbf{B}_{\theta}^T \mathbf{S}_g \mathbf{g} \, \mathrm{d}V \right) \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T r \, \mathrm{d}V - \int_{\Omega^e} \mathbf{N}_{\theta}^T \mathcal{D}_{\mathrm{mech}} \, \mathrm{d}V - \int_{\partial \Omega_q^e} \mathbf{N}_{\theta}^T \bar{\mathbf{g}} \, \mathrm{d}A \tag{64}
\end{align}
$$

$$
\begin{align*}
\mathcal{R}_{\sigma} &= \int_{\Omega^e} S_{\tau}^T \mathbf{B}_d \mathbf{d} \, \mathrm{d}V - \int_{\Omega^e} S_{\tau}^T \boldsymbol{\alpha} (\mathbf{N}_{\theta} \hat{\theta} - \theta_0) \, \mathrm{d}V \\
&\quad - \int_{\Omega^e} S_{\tau}^T \mathbf{D} \mathbf{S}_{\tau} \boldsymbol{\tau} \, \mathrm{d}V
\end{align*}
$$

$$
\mathcal{R}_{\mathbf{q}} = \int_{\Omega^e} S_{\mathbf{g}}^{T} \mathbf{B}_{\theta} \hat{\theta} \, \mathrm{d}V - \int_{\Omega^e} S_{\mathbf{g}}^{T} \mathbf{k}^{-1} S_{\mathbf{g}} \mathbf{g} \, \mathrm{d}V
$$

To solve such a nonlinear problem, we use Newton’s iterative method where at each iteration *i* (with (*i*) as iteration counter) we perform the consistent linearization of the residual leading to:

$$
\begin{align}
\mathcal{R}_{\mathbf{u},n+1/2}^{e,(i+1)}(\mathbf{d}, \boldsymbol{\tau}) &:= \mathcal{R}_{\mathbf{u},n+1/2}^{e,(i)}(\mathbf{d}, \boldsymbol{\tau}) \nonumber \\
&\quad + \Delta \mathcal{R}_{\mathbf{u},n+1/2}^{e,(i)}(\Delta \mathbf{d}, \Delta \boldsymbol{\tau}) := \mathbf{0}; \nonumber \\
\mathcal{R}_{\mathbf{u},n+1/2}^{e,(i)}(\cdot) &:= \frac{\int_{\Omega^e} \rho \mathbf{N}_d^T \mathbf{N}_d \mathrm{d}V \mathbf{a}_{n+1/2}^{e,(i)} + \int_{\Omega^e} \mathbf{B}_d^T \mathbf{S}_{\tau} \mathrm{d}V \boldsymbol{\tau}_{n+1/2}^{e,(i)}}{\mathbf{f}_{u,int}^{e,(i)}} \tag{65} \\
&\quad - \underbrace{\left( \int_{\Omega^e} \mathbf{N}_d^T \mathbf{b}_{n+1/2}^e \mathrm{d}V + \int_{\partial \Omega_2^e} \mathbf{N}_d^T \bar{\boldsymbol{\tau}}_{n+1/2}^e \mathrm{d}A \right)}_{\mathbf{f}_{u ext}^{e,(i)}}, \nonumber \\
\Delta \mathcal{R}_{\mathbf{u},n+1/2}^{e,(i)}(\cdot) &:= \chi_{dd} \mathbf{M}_{dd}^{e,(i)} \Delta \mathbf{d}_{n+1/2}^{e,(i)} + \mathbf{F}_{d\tau}^{e,(i)} \Delta \boldsymbol{\tau}_{n+1/2}^{e,(i)} \nonumber
\end{align}
$$

where we denote the matrices as:

$$
\begin{align}
\mathbf{M}_{dd}^{e} &:= \int_{\Omega^e} \rho \mathbf{N}_d^T \mathbf{N}_d \, \mathrm{d}V \nonumber \\
\mathbf{F}_{d\tau}^{e} &:= \int_{\Omega^e} \mathbf{B}_d^T \mathbf{S}_{\tau} \, \mathrm{d}V \tag{66} \\
\chi_{dd} &:= \partial \mathbf{a}_{n+1/2} / \partial \mathbf{d}_{n+1/2} \nonumber
\end{align}
$$

Linearizing the second equation in weak form, given by Eq. (63), leads to the following expression:

$$
\begin{align}
\mathcal{R}_{\theta,n+1/2}^{e,(i+1)}(\mathbf{d}, \hat{\theta}, \boldsymbol{\tau}, \mathbf{g}) &:= \mathcal{R}_{\theta,n+1/2}^{e,(i)}(\mathbf{d}, \hat{\theta}, \boldsymbol{\tau}, \mathbf{g}) \nonumber \\
&\quad + \Delta \mathcal{R}_{\theta,n+1/2}^{e,(i)}(\Delta \mathbf{d}, \Delta \hat{\theta}, \Delta \boldsymbol{\tau}, \Delta \mathbf{g}) = \mathbf{0}; \nonumber \\
\mathcal{R}_{\theta,n+1/2}^{e,(i)}(\cdot) &:= \int_{\Omega^e} \rho c \mathbf{N}_{\theta}^T \mathbf{N}_{\theta} \, \mathrm{d}V \hat{\theta}_{n+1/2}^{e,(i)} \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\mathbf{N}_{\theta} \hat{\theta}) (\mathbf{D}^{-1} : \boldsymbol{\alpha}) : \boldsymbol{\alpha} \right] \mathbf{N}_{\theta} \, \mathrm{d}V \hat{\theta}_{n+1/2}^{e,(i)} \nonumber \\
&\quad + \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\mathbf{D}^{-1} : \boldsymbol{\alpha}) (\mathbf{B}_d \mathbf{v}_{n+1/2}^{e,(i)}) \right] \mathbf{N}_{\theta} \hat{\theta}_{n+1/2}^{e,(i)} \, \mathrm{d}V \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\hat{\mathbf{D}} : \boldsymbol{\alpha}) (\mathbf{B}_d \mathbf{d}_{n+1/2}^{e,(i)}) \right] \mathbf{N}_{\theta} \hat{\theta}_{n+1/2}^{e,(i)} \, \mathrm{d}V \hat{\theta}_{n+1/2}^{e,(i)} \tag{67} \\
&\quad + \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\hat{\mathbf{D}} : \boldsymbol{\alpha}) (\boldsymbol{\alpha} (\mathbf{N}_{\theta} \hat{\theta}_{n+1/2}^{e,(i)} - \theta_0)) \right] \mathbf{N}_{\theta} \, \mathrm{d}V \hat{\theta}_{n+1/2}^{e,(i)} \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T \boldsymbol{\alpha} \mathbf{S}_{\tau} \, \mathrm{d}V \boldsymbol{\tau}_{n+1/2}^{e,(i)} + \int_{\Omega^e} \mathbf{B}_{\theta}^T \mathbf{S}_g \, \mathrm{d}V \mathbf{g}_{n+1/2}^{e,(i)} \nonumber \\
&\quad - \int_{\Omega^e} \mathbf{N}_{\theta}^T r \, \mathrm{d}V - \int_{\Omega^e} \mathbf{N}_{\theta}^T \mathcal{D}_{\mathrm{mech}} \, \mathrm{d}V - \int_{\partial \Omega_q^e} \mathbf{N}_{\theta}^T \bar{\mathbf{g}} \, \mathrm{d}A \nonumber \\
&\qquad \underbrace{\mathbf{f}_{\theta,ext}}_{\mathbf{f}_{\theta,ext}^{e,(i)}} \nonumber \\
&\quad \Delta \mathcal{R}_{\theta,n+1/2}^{e,(i)}(\cdot) := (\chi_{\theta d} \mathbf{D}_{\theta d}^e - \mathbf{F}_{\theta d}^e) \Delta \mathbf{d}_{n+1/2}^{e,(i)} \nonumber \\
&\quad + (\chi_{\theta \theta} \mathbf{D}_{\theta \theta}^e + \mathbf{F}_{\theta \theta}^e) \Delta \hat{\theta}_{n+1/2}^{e,(i)} - \mathbf{F}_{\theta \tau}^e \Delta \boldsymbol{\tau}_{n+1/2}^{e,(i)} + \mathbf{F}_{\theta g}^e \Delta \mathbf{g}_{n+1/2}^{e,(i)} \nonumber
\end{align}
$$

where we introduced
$$ \mathbf{D}_{\theta d}^e := \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\mathbf{N}_{\theta} \theta) \left( \mathbf{D}^{-1} : \boldsymbol{\alpha} \right) \right] \mathbf{B}_d \, dV; $$

$$ \mathbf{F}_{\theta d}^e := \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ (\mathbf{N}_{\theta} \theta) \left( \left( \mathbf{D}^{-1} : \dot{\mathbf{D}} : \mathbf{D}^{-1} \right) : \boldsymbol{\alpha} \right) \right] \mathbf{B}_d \, dV $$

$$ \mathbf{D}_{\theta\theta}^e := \mathbf{D}_{\theta\theta}^{e,1} - \mathbf{D}_{\theta\theta}^{e,2}; $$

$$ \mathbf{D}_{\theta\theta}^{e,1} := \int_{\Omega^e} \rho c \mathbf{N}_{\theta}^T \mathbf{N}_{\theta} \, dV; $$

$$ \mathbf{D}_{\theta\theta}^{e,2} := \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ \left( \mathbf{N}_{\theta} \hat{\theta} \right) \left( \mathbf{D}^{-1} : \boldsymbol{\alpha} \right) : \boldsymbol{\alpha} \right] \mathbf{N}_{\theta} \, dV; $$

$$ \mathbf{F}_{\theta\theta}^e := \mathbf{F}_{\theta\theta}^{e,1} - \mathbf{F}_{\theta\theta}^{e,2} + \mathbf{F}_{\theta\theta}^{e,3}; $$

$$ \mathbf{F}_{\theta\theta}^{e,1} := \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ \left( \mathbf{D}^{-1} : \boldsymbol{\alpha} \right) (\mathbf{B}_d \mathbf{v}) \right] \mathbf{N}_{\theta} \, dV; \quad (68) $$

$$ \mathbf{F}_{\theta\theta}^{e,2} := \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ \left( \hat{\mathbf{D}} : \boldsymbol{\alpha} \right) (\mathbf{B}_d \mathbf{d}) \right] \mathbf{N}_{\theta} \, dV; $$

$$ \mathbf{F}_{\theta\theta}^{e,3} := \int_{\Omega^e} \mathbf{N}_{\theta}^T \left[ \left( \hat{\mathbf{D}} : \boldsymbol{\alpha} \right) \left( \boldsymbol{\alpha} \left( \mathbf{N}_{\theta} \hat{\theta} - \theta_0 \right) \right) \right] \mathbf{N}_{\theta} \, dV; $$

$$ \mathbf{F}_{\theta\tau}^e := \int_{\Omega^e} \mathbf{N}_{\theta}^T \boldsymbol{\alpha} \mathbf{S}_{\tau} \, dV; $$

$$ \mathbf{F}_{\theta g}^e := \int_{\Omega^e} \mathbf{B}_{\theta}^T \mathbf{S}_g \, dV; $$

$$ \chi_{\theta d} = \partial \mathbf{v}_{n+1/2} / \partial \mathbf{d}_{n+1/2}; $$

$$ \chi_{\theta\theta} = \partial \dot{\theta}_{n+1/2} / \partial \dot{\theta}_{n+1/2} $$

Upon linearizing the third equation in the Eq. (63), we arrive at the following expression:

$$ \mathbf{R}_{\sigma,n+1/2}^{e,(i+1)}(\mathbf{d}, \hat{\theta}, \boldsymbol{\tau}) := \mathbf{R}_{\sigma,n+1/2}^{e,(i)}(\mathbf{d}, \hat{\theta}, \boldsymbol{\tau}) + \Delta \mathbf{R}_{\sigma,n+1/2}^{e,(i)}(\Delta \mathbf{d}, \Delta \hat{\theta}, \Delta \boldsymbol{\tau}) = \mathbf{0}; $$

$$ \begin{aligned} \mathbf{R}_{\sigma,n+1/2}^{e,(i)}(\cdot) := & \int_{\Omega^e} \mathbf{S}_{\tau}^T \mathbf{B}_d \, dV \, \mathbf{d}_{n+1/2}^{e,(i)} \\ & - \int_{\Omega^e} \mathbf{S}_{\tau}^T \boldsymbol{\alpha} \left( \mathbf{N}_{\theta} \hat{\theta}_{n+1/2}^{e,(i)} - \theta_0 \right) \, dV \\ & - \int_{\Omega^e} \mathbf{S}_{\tau}^T \mathbf{D} \mathbf{S}_{\tau} \, dV \, \boldsymbol{\tau}_{n+1/2}^{e,(i)}; \end{aligned} $$

$$ \begin{aligned} \mathbf{R}_{\sigma,n+1/2}^{e,(i)}(\cdot) := & \int_{\Omega^e} \mathbf{S}_{\tau}^T \mathbf{B}_d \, dV \, \mathbf{d}_{n+1/2}^{e,(i)} \\ & - \int_{\Omega^e} \mathbf{S}_{\tau}^T \boldsymbol{\alpha} \mathbf{N}_{\theta} \, dV \, \hat{\theta}_{n+1/2}^{e,(i)} \\ & - \int_{\Omega^e} \mathbf{S}_{\tau}^T \mathbf{D} \mathbf{S}_{\tau} \, dV \, \boldsymbol{\tau}_{n+1/2}^{e,(i)} \\ & - \underbrace{\left( - \int_{\Omega^e} \mathbf{S}_{\tau}^T \boldsymbol{\alpha} \theta_0 \, dV \right)}_{\mathbf{f}_{\sigma,ext}}; \end{aligned} \quad (69) $$

$$ \begin{aligned} \Delta \mathbf{R}_{\sigma,n+1/2}^{e,(i)}(\cdot) := & \mathbf{F}_{\tau d}^e \Delta \mathbf{d}_{n+1/2}^{e,(i)} \\ & - \mathbf{F}_{\tau \theta}^e \Delta \hat{\theta}_{n+1/2}^{e,(i)} - \mathbf{F}_{\tau \tau}^e \Delta \boldsymbol{\tau}_{n+1/2}^{e,(i)} \end{aligned} $$

where the corresponding matrices are

$$ \begin{aligned} \mathbf{F}_{\tau d}^e &:= \int_{\Omega^e} \mathbf{S}_{\tau}^T \mathbf{B}_d \, dV = \mathbf{F}_{d\tau}^{e,T}; \\ \mathbf{F}_{\tau \theta}^e &:= \int_{\Omega^e} \mathbf{S}_{\tau}^T \boldsymbol{\alpha} \mathbf{N}_{\theta} \, dV = \mathbf{F}_{\theta \tau}^{e,T}; \\ \mathbf{F}_{\tau \tau}^e &:= \int_{\Omega^e} \mathbf{S}_{\tau}^T \mathbf{D} \mathbf{S}_{\tau} \, dV \end{aligned} \quad (70) $$

Finally, we can express the linearized form of the last expression in Eq. (63) as follows:

$$ \begin{aligned} \mathbf{R}_{\mathbf{q},n+1/2}^{e,(i+1)}(\hat{\theta}, \mathbf{g}) &= \mathbf{R}_{\mathbf{q},n+1/2}^{e,(i)}(\hat{\theta}, \mathbf{g}) \\ &\quad + \Delta \mathbf{R}_{\mathbf{q},n+1/2}^{e,(i)}(\Delta \hat{\theta}, \Delta \mathbf{g}) = \mathbf{0}; \end{aligned} \quad (71) $$
$$ \begin{aligned} \mathbf{R}_{\mathbf{q},n+1/2}^{e,(i)}(\cdot) := & \int_{\Omega^e} \mathbf{S}_g^T \mathbf{B}_\theta \, dV \, \hat{\theta}_{n+1/2}^{e,(i)} \\ & - \int_{\Omega^e} \mathbf{S}_g^T \mathbf{k}^{-1} \mathbf{S}_g \, dV \, \mathbf{g}_{n+1/2}^{e,(i)}; \end{aligned} \quad (72) $$
$$ \begin{aligned} \Delta \mathbf{R}_{\mathbf{q},n+1/2}^{e,(i)}(\cdot) := & \mathbf{F}_{g\theta}^e \Delta \hat{\theta}_{n+1/2}^{e,(i)} \\ & - \mathbf{F}_{gg}^e \Delta \mathbf{g}_{n+1/2}^{e,(i)} \end{aligned} $$

where we define the matrices as:

$$ \begin{aligned} \mathbf{F}_{g\theta}^e &:= \int_{\Omega^e} \mathbf{S}_g^T \mathbf{B}_\theta \, dV = \mathbf{F}_{\theta g}^{e,T}; \\ \mathbf{F}_{gg}^e &:= \int_{\Omega^e} \mathbf{S}_g^T \mathbf{k}^{-1} \mathbf{S}_g \, dV \end{aligned} \quad (72) $$

By combining the linearization results for the weak form equations given by Eq. (63), we can formulate a set of discrete equations in matrix notation. The unknowns in these equations are the nodal displacement, stress, temperature, and heat flux fields, which need to be computed at each iteration.

$$ \mathbf{A}_{e=1}^{nel} [K^e] \begin{bmatrix} \Delta \mathbf{d}_{n+1/2}^{e,(i)} \\ \Delta \hat{\theta}_{n+1/2}^{e,(i)} \\ \Delta \boldsymbol{\tau}_{n+1/2}^{e,(i)} \\ \Delta \mathbf{g}_{n+1/2}^{e,(i)} \end{bmatrix} = - \mathbf{A}_{e=1}^{nel} \begin{bmatrix} \mathbf{r}_{u,n+1/2}^{e,(i)}(\cdot) \\ \mathbf{r}_{\theta,n+1/2}^{e,(i)}(\cdot) \\ \mathbf{r}_{\sigma,n+1/2}^{e,(i)}(\cdot) \\ \mathbf{r}_{q,n+1/2}^{e,(i)}(\cdot) \end{bmatrix} \quad (73) $$

$$ [K^e] = \begin{bmatrix} \chi_{dd} \mathbf{M}_{dd}^e & 0 & \mathbf{F}_{d\tau}^e & 0 \\ \chi_{\theta d} \mathbf{D}_{\theta d}^e - \mathbf{F}_{\theta d}^e & \chi_{\theta \theta} \mathbf{D}_{\theta \theta}^e + \mathbf{F}_{\theta \theta}^e & -\mathbf{F}_{\theta \tau}^e & \mathbf{F}_{\theta g}^e \\ \mathbf{F}_{\tau d}^e & -\mathbf{F}_{\tau \theta}^e & -\mathbf{F}_{\tau \tau}^e & 0 \\ 0 & \mathbf{F}_{g\theta}^e & 0 & -\mathbf{F}_{gg}^e \end{bmatrix} \quad (74) $$

After obtaining the nodal values through the solution of the linearized discrete equations, we can then update the values of stress, displacement, temperature, and heat-flux accordingly.
$$
\begin{align}
\mathbf{d}_{n+1/2}^{(i+1)} &= \mathbf{d}_{n+1/2}^{(i)} + \Delta\mathbf{d}_{n+1/2}^{(i)} \nonumber \\
\hat{\theta}_{n+1/2}^{(i+1)} &= \hat{\theta}_{n+1/2}^{(i)} + \Delta\hat{\theta}_{n+1/2}^{(i)} \nonumber \\
\boldsymbol{\tau}_{n+1/2}^{(i+1)} &= \boldsymbol{\tau}_{n+1/2}^{(i)} + \Delta\boldsymbol{\tau}_{n+1/2}^{(i)} \nonumber \\
\mathbf{g}_{n+1/2}^{(i+1)} &= \mathbf{g}_{n+1/2}^{(i)} + \Delta\mathbf{g}_{n+1/2}^{(i)} \tag{75}
\end{align}
$$

3.4 Time-stepping integration schemes enforcing energy-conservation and energy-decaying of high-frequency modes

Here, we propose an energy-conserving method for address-
ing linear thermo-elastodynamics problems. To achieve
second-order accuracy, we use the mid-point integration
scheme [36]. When applied to damage, the scheme involves
the simultaneous integration of balance equations along with
the damage compliance and hardening variable evolution
equations. We note in passing that the latter implies turn-
ing towards the viscous regularisation of damage, in order to
avoid the operator split procedure typical of inelastic consti-
tutive models (e.g. see [35]).

At the time $t_{n+1/2}$, we can obtain the nodal values of veloc-
ity and acceleration by computing the difference between
their values at two consecutive time steps.

$$
\begin{align}
\mathbf{v}_{n+1/2} &= (\mathbf{d}_{n+1} - \mathbf{d}_n) / \Delta t; \tag{76} \\
\mathbf{a}_{n+1/2} &= (\mathbf{v}_{n+1} - \mathbf{v}_n) / \Delta t \nonumber
\end{align}
$$

We can express the incremental displacement by using the
last result:

$$
\begin{align}
\mathbf{d}_{n+1} - \mathbf{d}_n &= \frac{1}{2} \Delta t (\mathbf{v}_{n+1} + \mathbf{v}_n) \tag{77} \\
\Rightarrow \quad \mathbf{d}_{n+1} &= \mathbf{d}_n + \mathbf{u}; \notag
\end{align}
$$

To ensure energy conservation, the stress value at the time
$t_{n+1/2}$ is chosen as follows:

$$
\tau_{n+1/2}^{\text{avg}} := \frac{1}{2} (\tau_{n+1} + \tau_n) \qquad (78)
$$

By selecting a test displacement vector $\delta \mathbf{d} = \mathbf{d}_{n+1} - \mathbf{d}_n$, we can demonstrate that the work done by internal $f_{\text{int}}$ and external forces $f_{\text{ext}}$ from the weak form would be balanced, as shown below:

$$
(\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{f}_{\text{int}} = (\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{f}_{\text{ext}} \quad (79)
$$

We can reformulate the relevant balance equation as fol-
lows:

$$
\begin{equation}
\begin{split}
(\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{f}_{\text{int}} ={}& (\mathbf{d}_{n+1} - \mathbf{d}_n)^T \int_{\Omega} \rho \mathbf{N}^T \mathbf{N} \mathrm{d}V \mathbf{a}_{n+1/2} \\
& + (\mathbf{d}_{n+1} - \mathbf{d}_n)^T \int_{\Omega} \mathbf{B}^T \mathbf{S}_{\tau} \mathrm{d}V \boldsymbol{\tau}_{n+1/2}
\end{split}
\tag{80}
\end{equation}
$$

In this context, the work done by the inertia force, which
contributes to the external load in the spirit of d'Alembert,
can be simplified by using the symmetry of the mass matrix
$\mathbf{M}_{dd}$, which leads to the expressions

$$
\begin{align}
(\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{M}_{dd} \mathbf{a}_{n+1/2} &= \frac{1}{2} (\mathbf{v}_{n+1} + \mathbf{v}_n)^T \mathbf{M}_{dd} (\mathbf{v}_{n+1} - \mathbf{v}_n) \nonumber \\
&= \frac{1}{2} \mathbf{v}_{n+1}^T \mathbf{M}_{dd} \mathbf{v}_{n+1} - \frac{1}{2} \mathbf{v}_n^T \mathbf{M}_{dd} \mathbf{v}_n \nonumber \\
&= K_{n+1} - K_n \tag{81}
\end{align}
$$

where $K_{n+1}$ and $K_n$ represents the kinetic energy of the sys-
tem at time $t_{n+1}$ and $t_n$, respectively. The work done by the
internal force can be expressed as a linear combination of
two terms: the first term involves values at time step $t_{n+1}$,
while the second term is composed of field variables from
two consecutive time steps.

$$
\begin{align}
& (\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{F}_{d\tau} \boldsymbol{\tau}_{n+1/2} \notag \\
&= \frac{1}{2} (\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{F}_{d\tau} (\boldsymbol{\tau}_{n+1} + \boldsymbol{\tau}_n) \notag \\
&= \left( \frac{1}{2} \mathbf{d}_{n+1}^T \mathbf{F}_{d\tau} \boldsymbol{\tau}_{n+1} - \frac{1}{2} \mathbf{d}_n^T \mathbf{F}_{d\tau} \boldsymbol{\tau}_n \right) \notag \\
&\qquad + \underbrace{\left( \frac{1}{2} \mathbf{d}_{n+1}^T \mathbf{F}_{d\tau} \boldsymbol{\tau}_n - \frac{1}{2} \mathbf{d}_n^T \mathbf{F}_{d\tau} \boldsymbol{\tau}_{n+1} \right)}_{L_1} \tag{82}
\end{align}
$$

The second term $L_1$ can be further simplified by selecting
a test stress function $\delta \boldsymbol{\tau} = \boldsymbol{\tau}_{n+1} - \boldsymbol{\tau}_n$ and utilizing the result
from Eq. (69). By exploiting the symmetry of the matrix $\mathbf{F}_{\tau\tau}$
and the relation $\mathbf{F}_{d\tau} = \mathbf{F}_{\tau d}^T$, we obtain the expression:

$$
\begin{align}
& (\boldsymbol{\tau}_{n+1} - \boldsymbol{\tau}_n)^T \nonumber \\
& \quad [\mathbf{F}_{\tau d} \mathbf{d}_{n+1/2} - \mathbf{F}_{\tau \theta} \theta_{n+1/2} - \mathbf{F}_{\tau \tau} \boldsymbol{\tau}_{n+1/2}] = 0 \nonumber \\
& + (\boldsymbol{\tau}_{n+1} - \boldsymbol{\tau}_n)^T \mathbf{F}_{\tau d} \left( \frac{\mathbf{d}_{n+1} + \mathbf{d}_n}{2} \right) \nonumber \\
& - (\boldsymbol{\tau}_{n+1} - \boldsymbol{\tau}_n)^T \mathbf{F}_{\tau \theta} \left( \frac{\theta_{n+1} + \theta_n}{2} \right) \nonumber \\
& - (\boldsymbol{\tau}_{n+1} - \boldsymbol{\tau}_n)^T \mathbf{F}_{\tau \tau} \left( \frac{\boldsymbol{\tau}_{n+1} + \boldsymbol{\tau}_n}{2} \right) = 0 \nonumber \\
& \Rightarrow \underbrace{\left( \frac{1}{2} \mathbf{d}_{n+1}^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\tau \theta}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2} \boldsymbol{\tau}_{n+1}^T \mathbf{F}_{\tau \tau}^T \boldsymbol{\tau}_{n+1} \right)}_{L_{1,n+1}} \tag{83} \\
& - \underbrace{\left( \frac{1}{2} \mathbf{d}_n^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_n - \frac{1}{2} \theta_n^T \mathbf{F}_{\tau \theta}^T \boldsymbol{\tau}_n - \frac{1}{2} \boldsymbol{\tau}_n^T \mathbf{F}_{\tau \tau}^T \boldsymbol{\tau}_n \right)}_{L_{1,n}} \nonumber \\
& + \underbrace{\left( \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\tau \theta}^T \boldsymbol{\tau}_n - \frac{1}{2} \theta_n^T \mathbf{F}_{\tau \theta}^T \boldsymbol{\tau}_{n+1} \right)}_{L_{1,comb}} \nonumber \\
& = \underbrace{\frac{1}{2} \mathbf{d}_{n+1}^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_n - \frac{1}{2} \mathbf{d}_n^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_{n+1}}_{L_1} \nonumber
\end{align}
$$
By utilizing the previous result and the one presented in Eq. (82), we can rewrite the work done by internal forces as follows:

$$ 
\begin{aligned}
(\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{F}_{d\tau} \boldsymbol{\tau}_{n+1/2} &= \\
& \quad \left(\frac{1}{2}\mathbf{d}_{n+1}^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2}\mathbf{d}_n^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_n\right) + \\
& \quad \underbrace{\left(\frac{1}{2}\mathbf{d}_{n+1}^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2}\theta_{n+1}^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2}\boldsymbol{\tau}_{n+1}^T \mathbf{F}_{\tau\tau} \boldsymbol{\tau}_{n+1}\right)}_{L_{1,n+1}} \\
& \quad - \underbrace{\left(\frac{1}{2}\mathbf{d}_n^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_n - \frac{1}{2}\theta_n^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_n - \frac{1}{2}\boldsymbol{\tau}_n^T \mathbf{F}_{\tau\tau} \boldsymbol{\tau}_n\right)}_{L_{1,n}} \\
& \quad + \underbrace{\left(\frac{1}{2}\theta_{n+1}^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_n - \frac{1}{2}\theta_n^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_{n+1}\right)}_{L_{1,comb}} = \\
&= \underbrace{\left(\mathbf{d}_{n+1}^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2}\boldsymbol{\tau}_{n+1}^T \mathbf{F}_{\tau\tau} \boldsymbol{\tau}_{n+1}\right)}_{P_{n+1}} \\
& \quad - \underbrace{\left(\mathbf{d}_n^T \mathbf{F}_{\tau d}^T \boldsymbol{\tau}_n - \frac{1}{2}\boldsymbol{\tau}_n^T \mathbf{F}_{\tau\tau} \boldsymbol{\tau}_n\right)}_{P_n} \\
& \quad - \underbrace{\left(\frac{1}{2}\theta_{n+1}^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_{n+1} - \frac{1}{2}\theta_n^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_n\right)}_{\mathcal{D}_{\theta\tau}} + \\
& \quad \underbrace{\left(\frac{1}{2}\theta_{n+1}^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_n - \frac{1}{2}\theta_n^T \mathbf{F}_{\tau\theta}^T \boldsymbol{\tau}_{n+1}\right)}_{L_2} \\
&= P_{n+1} - P_n + \mathcal{D}_{\theta\tau} + L_2
\end{aligned}
\quad (84) $$

In the above equation, $P_n$ and $P_{n+1}$ denote the potential energies of the system at times $t_n$ and $t_{n+1}$, respectively. The term $\mathcal{D}_{\theta\tau}$ denotes the dissipation that arises from the combined temperature and stress fields, while the scalar $L_2$ is a combination of values from two consecutive time steps.

Furthermore, we can derive the following identity by selecting a test temperature function $\delta\theta = \theta_{n+1} - \theta_n$ and utilizing the relation $\mathbf{F}_{\theta\tau} = \mathbf{F}_{\tau\theta}^T$:

$$ 
\begin{aligned}
& + (\theta_{n+1} - \theta_n)^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_{n+1/2} \\
& + (\theta_{n+1} - \theta_n)^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_{n+1/2} \\
& - (\theta_{n+1} - \theta_n)^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_{n+1/2} \\
& + (\theta_{n+1} - \theta_n)^T \mathbf{F}_{\theta g} \mathbf{g}_{n+1/2} = 0 \\
& (\theta_{n+1} - \theta_n)^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \left(\frac{\mathbf{d}_{n+1} + \mathbf{d}_n}{2}\right) \\
& + (\theta_{n+1} - \theta_n)^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \left(\frac{\theta_{n+1} + \theta_n}{2}\right)
\end{aligned}
$$

$$ 
\begin{aligned}
& - (\theta_{n+1} - \theta_n)^T \mathbf{F}_{\theta\tau} \left(\frac{\boldsymbol{\tau}_{n+1} + \boldsymbol{\tau}_n}{2}\right) \\
& + (\theta_{n+1} - \theta_n)^T \mathbf{F}_{\theta g} \left(\frac{\mathbf{g}_{n+1} + \mathbf{g}_n}{2}\right) = 0 \\
\Rightarrow \quad & + \frac{1}{2} \theta_{n+1}^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_{n+1} \\
& + \frac{1}{2} \theta_{n+1}^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_n \\
& + \frac{1}{2} \theta_{n+1}^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_{n+1} \\
& + \frac{1}{2} \theta_{n+1}^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_n \\
& - \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_{n+1} - \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_n \\
& + \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta g} \mathbf{g}_{n+1} + \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta g} \mathbf{g}_n \\
& - \frac{1}{2} \theta_n^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_{n+1} \\
& - \frac{1}{2} \theta_n^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_n \\
& - \frac{1}{2} \theta_n^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_{n+1} - \frac{1}{2} \theta_n^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_n \\
& + \frac{1}{2} \theta_n^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_{n+1} + \frac{1}{2} \theta_n^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_n \\
& - \frac{1}{2} \theta_n^T \mathbf{F}_{\theta g} \mathbf{g}_{n+1} - \frac{1}{2} \theta_n^T \mathbf{F}_{\theta g} \mathbf{g}_n = 0
\end{aligned}
\quad (85) $$

This yields the following:

$$ 
\begin{aligned}
& \frac{1}{2} \theta_{n+1}^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_{n+1} \\
& + \frac{1}{2} \theta_{n+1}^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_{n+1} \\
& - \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_{n+1} + \frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta g} \mathbf{g}_{n+1} \\
& - \left(\frac{1}{2} \theta_n^T (\chi_{\theta d} \mathbf{D}_{\theta d} - \mathbf{F}_{\theta d}) \mathbf{d}_n\right) \\
& - \left(\frac{1}{2} \theta_n^T (\chi_{\theta\theta} \mathbf{D}_{\theta\theta} + \mathbf{F}_{\theta\theta}) \theta_n\right) \\
& - \left(-\frac{1}{2} \theta_n^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_n + \frac{1}{2} \theta_n^T \mathbf{F}_{\theta g} \mathbf{g}_n\right) \\
& - \underbrace{\left(\frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta g} \mathbf{g}_n - \frac{1}{2} \theta_n^T \mathbf{F}_{\theta g} \mathbf{g}_{n+1}\right)}_{L_3} \\
& = \underbrace{\frac{1}{2} \theta_{n+1}^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_n - \frac{1}{2} \theta_n^T \mathbf{F}_{\theta\tau} \boldsymbol{\tau}_{n+1}}_{L_2}
\end{aligned}
\quad (86) $$

The term $\mathcal{D}_{\theta}$ represents the dissipation that arises from the temperature field. To further elaborate on the term $L_2$, we can use the result in Eq. (71) and select a test flux function
<table><caption>**Table 1** Adopted material parameters for the proposed 3D bending test</caption><thead><tr><th>E<br>[N/m<sup>2</sup>]</th><th>ν<br>[-]</th><th>ρ<br>[Ns<sup>2</sup>m<sup>-4</sup>]</th><th>c<br>[m<sup>2</sup>s<sup>-2</sup>K<sup>-1</sup>]</th><th>k<br>[Ns<sup>-1</sup>K<sup>-1</sup>]</th><th>α<br>[K<sup>-1</sup>]</th><th>θ<sub>0</sub><br>[K]</th></tr></thead><tbody><tr><td>3 · 10<sup>10</sup></td><td>0.2</td><td>2400</td><td>880</td><td>2.0</td><td>1 · 10<sup>-5</sup></td><td>273.15</td></tr></tbody></table>

$δg = g_{n+1} - g_n$. This yields the following expression:

$$ 
\begin{aligned}
& (g_{n+1} - g_n)^T [\mathbf{F}_{g\theta}\theta_{n+1/2} - \mathbf{F}_{gg}g_{n+1/2}] = 0 \\
& (g_{n+1} - g_n)^T \left[ \mathbf{F}_{g\theta} \left( \frac{\theta_{n+1} + \theta_n}{2} \right) - \mathbf{F}_{gg} \left( \frac{g_{n+1} + g_n}{2} \right) \right] = 0 \\
& \Rightarrow \underbrace{\left( \frac{1}{2}\theta_{n+1}^T \mathbf{F}_{g\theta}^T \mathbf{g}_{n+1} - \frac{1}{2}\mathbf{g}_{n+1}^T \mathbf{F}_{gg}^T \mathbf{g}_{n+1} \right)}_{\mathcal{D}_{\theta_g}^{n+1}} \\
& \quad - \underbrace{\left( \frac{1}{2}\theta_n^T \mathbf{F}_{g\theta}^T \mathbf{g}_n - \frac{1}{2}\mathbf{g}_n^T \mathbf{F}_{gg}^T \mathbf{g}_n \right)}_{\mathcal{D}_{\theta_g}^n} \\
& = \underbrace{\frac{1}{2}\theta_{n+1}^T \mathbf{F}_{g\theta}^T \mathbf{g}_n - \frac{1}{2}\theta_n^T \mathbf{F}_{g\theta}^T \mathbf{g}_{n+1}}_{L_3}
\end{aligned}
\quad (87) $$

The above equation is derived by utilizing the symmetry of $\mathbf{F}_{gg}$ and $\mathbf{F}_{\theta_g} = \mathbf{F}_{g\theta}^T$. By combining the results from the previous equations, we can add the contributions of each term in Eq. (79) to confirm the balance of internal energy and external work.

$$ 
\begin{aligned}
& \underbrace{(K_{n+1} + P_{n+1})}_{E_{n+1}} - \underbrace{(K_n + P_n)}_{E_n} + \\
& \underbrace{\left[ (\mathcal{D}_{\theta}^{n+1} + \mathcal{D}_{\theta_g}^{n+1}) - (\mathcal{D}_{\theta}^n + \mathcal{D}_{\theta_g}^n) + \mathcal{D}_{\theta\tau} \right]}_{\mathcal{D}_{\theta,\text{total}}} = \\
& \underbrace{(\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{f}_{\text{ext}}}_{\Delta W} \\
& \Leftrightarrow E_{n+1} - E_n + \mathcal{D}_{\theta,\text{total}} = \Delta W
\end{aligned}
\quad (88) $$

### 3.5 Numerical simulations for thermo-elastostatics

In this section we present the performance of the proposed numerical framework for pure statics and thermo-mechanics. We analyzed the rectangular beam with dimensions $b/h/l = 0.15/0.15/1.20$ m (where the general form is $a/a/8a$), with material properties given in Table 1.

Boundary conditions are chosen to cause dominantly bending deformation mode, for both mechanical and thermo-mechanical case defined by the corresponding boundary conditions illustrated in Fig. 7 and 8. In the proposed mechan-

<figure><img src="image_9.png" /><figcaption>(a) Mechanical loading</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Boundary conditions</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(a) Mechanical loading</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Boundary conditions</figcaption></figure>

<caption>**Fig. 8** Numerical setup for the thermo-mechanical bending test in 3D setting</caption>

ical test, the beam is loaded with a distributed load at the top and bottom edges to induce bending deformation, as shown in Fig. 7. The resultant load is equal to $\bar{p} = 60$ kN/m.

The mechanical test was extended to include thermal loading. The bottom of the beam was thermally insulated, with the boundary condition $\partial\Omega_{\theta}: \theta|_{z=0} = 0$ K, while the top of the beam was subjected to a temperature of $\partial\Omega_{\theta}: \theta|_{z=a} = \bar{\theta} = 373.15$ K (i.e 100°C), as shown in Fig. 8.

Several numerical simulations were conducted at various levels of mesh refinement to evaluate the convergence performance of the proposed hybrid element. We started with
<figure><img src="image_13.png" /><figcaption>(a) Mesh size 8x1x1</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Mesh size 16x2x2</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Mesh size 32x4x4</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Mesh size 64x8x8</figcaption></figure>

<caption>**Fig. 9** Different mesh types used in analysis</caption>

a coarse mesh of 8x1x1 and advanced to a finest mesh of 64x8x8, as illustrated in Fig. 9.

The evaluation of the convergence and quality of the finite element was based on the estimation of the Z-displacement value of the upper beam edge at the middle-length intersection. The Z-displacement field values, which were interpolated using standard shape functions for the tetrahedron element, are displayed in a colormap in Fig. 10 for the mechanical test and thermo-mechanical test from a finite element mesh type 2: 16x2x2.

The displacement values show gradual convergence towards the exact solution, highlighting the effectiveness of the proposed approach. To achieve the exact solution, the domain was discretized with a very fine mesh of 8-node brick elements, using *h*-refinement technique to improve the accuracy of the finite element method results. Figure 11 reveals that the proposed approach with an enhanced formulation, resulting in a hybrid tetrahedral element, demonstrates significantly better performance compared to its counterpart, the 4-node tetrahedron, and 8-node brick under both compressible and incompressible conditions.

<figure><img src="image_17.png" /><figcaption>(c) Displacements $d_z$ [m]</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Displacements $d_z$ [m]</figcaption></figure>

<caption>**Fig. 10** Displacement results for bending tests: mechanical and thermo-mechanical cases</caption>

The traction (or stress vector) components at the corresponding nodes on the faces of the tetrahedral elements are shown. Since the focus is on the bending test, we present the results for the traction components in the X-direction, as they are the most significant, along with the total stress vectors for both the mechanical and thermo-mechanical tests. For mechanical test, the primary traction component in the X-direction dominates in the upper zone, suggesting that the load is mainly due to compression rather than bending of the beam itself as illustrated in Fig. 12.

For thermo-mechanical test, the contribution of other two traction components in *Y* and *Z* directions are slightly more emphasized, appearing mainly in the support zones due to disturbances in those regions. In this way, the total traction component is more different from the Traction *X*-component as illustrated in Fig. 13.
<figure><img src="image_19.png" /><figcaption>(a) Case 1: Mechanical loading</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Case 2: Thermo-mechanical
loading</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(c) Case 1: Mechanical loading</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) Case 2: Thermo-mechanical
loading</figcaption></figure>

<caption>Fig. 11 Convergence properties of numerical model for: (a), (b) compressible; (c), (d) incompressible case</caption>

**3.6 Numerical simulations for non-stationary heat transfer**

This section briefly discusses the thermal energy conservation properties within the proposed numerical framework. Using the thermal properties provided in Table 2, we analyzed heat propagation and energy storage within the cube specimen, as introduced in the Section 4. The boundary conditions and the selected tetrahedral mesh division are depicted in Fig. 14.

**Remark 1** For the Case 1 with *k* = 1 in Fig. 15, energy conservation is not achieved since the rate of heat flux still exist due to thermal inertia and slow heat-wave propagation. Other cases are capable of representing thermal energy conservation.

**3.7 Numerical simulations for thermo-elastodynamics**

The results of the numerical simulations presented in this section are obtained by a research version of computer code FEAP, developed by RL Taylor [39]. In addition to demonstrating the performance of the hybrid model for statics, which is achieved through the use of special interpolations for the stress and heat flux fields, the performance of the proposed formulation is also presented for dynamic case. The satisfying numerical model’s performance under dynamic loading conditions is attributed to the time-stepping integration scheme, which conserves the total energy of the system during the free-vibration phase for the elasticity.

All numerical simulations were performed on a cube model with dimensions of a=0.15, made of a thermally

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>(a) Traction $\tau_x$ [N/m²]</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b) Total traction $\tau$ [N/m²]</figcaption></figure>

<caption>**Fig. 12** Illustration of traction components for bending mechanical test in 3D setting</caption>

isotropic material with elastic properties presented in Table 3.

We first analyzed the behavior of the model under mechanical and thermomechanical loading, which consisted of imposed tractions on the top face and boundary heat flux on the top face, as shown in Fig. 16(a,b). The displacement, temperature, traction, and heat flux values were extracted for characteristic nodes at the top boundary of the cube, as illustrated in Fig. 16c.

The dynamic loading was in the form of a symmetric triangular pulse with varying duration and time step used in the analysis, as illustrated in Table 3 for the elastodynamics case. First, the thermo-elastodynamic case was analyzed by considering the values of the relevant fields and the energy of the system. Since the displacements and temperatures were interpolated using standard shape functions, the validation of the results was performed by comparing them with those

<figure><img src="image_26.png" /><figcaption>(a) Traction $\tau_x$ [N/m²]</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) Total traction $\tau$ [N/m²]</figcaption></figure>

<caption>**Fig. 13** Illustration of traction components for bending thermomechanical test in 3D setting</caption>

<table><caption>**Table 2** Adopted material parameters for the thermal test</caption><thead><tr><th rowspan="2">Material parameter</th><th colspan="2">Thermal parameters</th><th rowspan="2">k</th></tr><tr><th>ρ</th><th>c</th></tr></thead><tbody><tr><td>Case 1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Case 2</td><td>1</td><td>1</td><td>10</td></tr><tr><td>Case 3</td><td>1</td><td>1</td><td>100</td></tr><tr><td>Case 4</td><td>1</td><td>1</td><td>1000</td></tr></tbody></table>

obtained from a standard 4-node tetrahedral element. The values were compared for node 63, which is located at the center of the exposed face of the cube specimen (0.075, 0.075, 0.15).

The numerical results for the displacement field are in perfect agreement with each other for Case 1 as shown in Fig. 17a, while there is a small difference in behavior for
<figure><img src="image_28.png" /><figcaption>Fig. 14 Numerical model setup for the proposed thermal test</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 15 Energy conservation for the thermal test</figcaption></figure>

Case 2 described by diagram in Fig. 17b. Overall, the peak value of displacement is practically the same for both case.

Considering the temperature field, the adopted values of thermal parameters are such that Case 1 describes a stationary heat flow, given the identical shape of the imposed boundary heat flux load diagram and the time-history temperature dia-

<figure><img src="image_30.png" /><figcaption>Fig. 16 Cube specimen with: (a) boundary conditions for the thermo-mechanical test, (b) boundary conditions for the mechanical test, (c) nodes for diagrams</figcaption></figure>

gram for Case 1 (see Fig. 18a). For Case 2, the peak values are smaller than those for Case 1, implying that a non-stationary

<table><caption>**Table 3** Adopted material parameters for the proposed 3D dynamic loading test</caption><thead><tr><th rowspan="2">Material parameter</th><th colspan="6">Thermo-elastic properties</th><th colspan="3">Triangular pulse</th></tr><tr><th><i>E</i></th><th><i>ν</i></th><th><i>ρ</i></th><th><i>c</i></th><th><i>k</i></th><th><i>α</i></th><th><i>p̄</i></th><th><i>q̄</i></th><th>Δ<i>t</i>, (dt)</th></tr></thead><tbody><tr><td>Case 1</td><td>200</td><td>0.2</td><td>1</td><td>1</td><td>1</td><td>0.0001</td><td>1</td><td>1</td><td>(8 · 10<sup>-1</sup>), (1 · 10<sup>-2</sup>)</td></tr><tr><td>Case 2</td><td>200</td><td>0.2</td><td>1</td><td>1</td><td>1</td><td>0.0001</td><td>1</td><td>1</td><td>(8 · 10<sup>-3</sup>), (1 · 10<sup>-4</sup>)</td></tr></tbody></table>
<figure><img src="image_31.png" /><figcaption>Case 1: Displacement $d_z$</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Case 2: Displacement $d_z$</figcaption></figure>

<caption>**Fig. 17** Displacement comparison for loading cases 1 and 2 for thermo-elastodynamics in 3D setting</caption>

heat flow. An excellent agreement between the results for both cases is illustrated in Fig. 18b.

In order to illustrate and verify the quality of the proposed element, comparative diagrams were presented for two additional cases of dynamic loading, with parameters: $[\Delta t, (dt) = 80, (1)]$, and then $[\Delta t, (dt) = (8 \cdot 10^{-2}), (1 \cdot 10^{-3})]$.

This difference in behavior can be attributed to the different stiffness of the element, given the different degrees of freedom, which results in slightly different behavior for dynamic loads in Case 2 and the additional analyzed case $[\Delta t, (dt) = (8 \cdot 10^{-2}), (1 \cdot 10^{-3})]$, as illustrated in Figs. 17(b and 19b.

<figure><img src="image_33.png" /><figcaption>Case 1: Temperature $\theta$</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Case 2: Temperature $\theta$</figcaption></figure>

<caption>**Fig. 18** Temperature comparison for loading cases 1 and 2 for thermo-elastodynamics in 3D setting</caption>

The numerical results for the displacement field show perfect agreement for Case 1, while a small discrepancy is observed in the behavior for Case 2. However, the peak displacement value remains nearly identical in both cases.

By using the proposed formulation, the results of elasto-dynamic numerical simulation in terms of total mechanical energy of the system are shown in Fig. 20(a,b). Obviously, the total mechanical energy is conserved during the free vibration phase in both cases. In Case 1, the contribution of kinetic energy to the total energy is practically negligible, leading to the conclusion that this is a practically static load. In Case 2, both kinetic and potential energy are of a similar order of magnitude, leading to the conservation of total energy as the sum of these terms.
<figure><img src="image_35.png" /></figure>

$$[\Delta t, (dt) = 80, (1)]$$

<figure><img src="image_36.png" /></figure>

$$[\Delta t, (dt) = (8 \cdot 10^{-2}, (1 \cdot 10^{-3})]$$

<figure><img src="image_37.png" /><figcaption>Case 1: Energy diagram</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Case 2: Energy diagram</figcaption></figure>

<caption>Fig. 20 Mechanical energy time-history plot under dynamic pulse loading</caption>

### 3.7.1 Numerical simulations for dynamics: Thermo-damage case

The proposed approach for the thermo-elastic case has been extended to include the modeling of the material’s inelastic behavior. A simple isotropic damage model with viscous regularisation (see [35]) is used to describe the inelastic behavior. The model requires two additional material parameters to characterize the damage: the damage stress $\sigma_f$ and the linear hardening modulus $K_{h,lin}$.

The results of a series of numerical simulations on the same cube sample with a side length of $a = 0.15$ are illustrated for the thermo-mechanical test and pure mechanical test with the load values given in Table 4.

Since the dynamic load has the same duration across all analyzed cases, maintaining its previously defined symmetric triangular pulse shape, the goal of the thermo-damage analysis was only to analyze the influence of different values of damage stress on the intensity and time-history diagram of traction components, as well as the total energy evolution of the system.

Diagrams of the total mechanical energy for the proposed rate-sensitive thermo-damage numerical model, defined by the parameters in Table 4, are shown in Fig. 21a. As expected, the total mechanical energy is conserved in the elastic case. However, for other cases, due to damage dissipation, the energy is no longer conserved. This effect is more
<table><caption>**Table 4** Load data and parameters for 3D Dynamic Analysis: Thermo-damage case</caption><thead><tr><th>Load character</th><th>Case</th><th>p&#773;</th><th>q&#773;</th><th>&Delta;t, (dt)</th><th>&sigma;<sub>f</sub></th><th>K<sub>h,lin</sub></th><th>&omega;</th></tr></thead><tbody><tr><td>Triangular pulse</td><td>1</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>0.0</td></tr><tr><td>Triangular pulse</td><td>2</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.5</td><td>100</td><td>0.0</td></tr><tr><td>Triangular pulse</td><td>3</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.6</td><td>100</td><td>0.0</td></tr><tr><td>Triangular pulse</td><td>4</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.7</td><td>100</td><td>0.0</td></tr><tr><td>Triangular pulse</td><td>5</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>&infin;</td><td>&mdash;</td><td>&mdash;</td></tr></tbody></table>

<figure><img src="image_39.png" /><figcaption>Energy diagram ($\eta = 1 \cdot 10^{-8}$)</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Accumulated damage dissipation ($\eta = 1 \cdot 10^{-8}$)</figcaption></figure>

<caption>**Fig. 21** Energy and damage dissipation time-history plots</caption>

pronounced for lower yield stress values, implying higher dissipation in these cases, as shown in Fig. 21b.

When the temperature-dependent yield parameters are considered, according to the data in Table 5, thermal dissipation is activated upon reaching the yield stress and triggering internal variables. Notably, larger values of the parameter $\omega$ lead to lower yield stress, resulting in higher thermal dissipation, as illustrated in Fig. 22a.

In all these cases, the viscous parameter $\eta$ is assumed constant, with a value of $1 \cdot 10^{-8}$. To evaluate the effect of varying viscous parameters, additional numerical simulations were performed using the data in Table 6 on the same cube specimen. This approach allowed us to compare the viscous damage dissipation for different viscous parameter values, as shown in Fig. 22b.

A convergence study was conducted for all three analyzed material behavior cases. The results are presented in Tables 7, 8, and 9, respectively, for the thermo-damage case, the thermo-damage case with temperature-dependent yield stress, and the rate-sensitive thermo-damage case, for different load cases as previously done. A time step of 100 was taken as the reference, where the material already exhibits damage, and internal variables have been activated.

Additionally, the traction component values are provided for node 191, located at coordinates (0.05, 0.05, 0.15). The time-history diagrams of these traction components are shown in Fig. 23a. The z-component values of traction are of primary importance here, given the direction of the applied load, while the other two components are less important, especially for Case 1, which resembles the static loading case. Given that the diagrams for the traction components $\tau_x$ and $\tau_y$ are practically identical, we can confirm that the proposed model perfectly describes the isotropy of the material. In any of these cases, the peak value of the traction component exceeds the damage stress $\sigma_f$, causing damage

<table><caption>**Table 5** Load data and parameters for 3D Dynamic Analysis: Thermo-damage case with temperature-dependent yield stress</caption><thead><tr><th>Load character</th><th>Case</th><th>p&#773;</th><th>q&#773;</th><th>&Delta;t, (dt)</th><th>&sigma;<sub>f</sub></th><th>K<sub>h,lin</sub></th><th>&omega;</th><th>&eta;</th></tr></thead><tbody><tr><td>Triangular pulse</td><td>1a</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>10</td><td>1 &middot; 10<sup>-8</sup></td></tr><tr><td>Triangular pulse</td><td>1b</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>20</td><td>1 &middot; 10<sup>-8</sup></td></tr><tr><td>Triangular pulse</td><td>1c</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>50</td><td>1 &middot; 10<sup>-8</sup></td></tr><tr><td>Triangular pulse</td><td>1d</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>), (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>100</td><td>1 &middot; 10<sup>-8</sup></td></tr></tbody></table>
<figure><img src="image_41.png" /><figcaption>Thermal dissipation ($\eta = 1 \cdot 10^{-8}$)</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Viscous damage dissipation</figcaption></figure>

<caption>**Fig. 22** Thermal and viscous dissipation time-history plots</caption>

deformation within the element and resulting in a decrease in the total energy of the system due to damage dissipation (see Fig. 23b). The chosen value of $K_{h,lin} = 0.5E$ allows us to illustrate the contribution of damage to energy and dissipation time history.

<table><caption>**Table 6** Load data and parameters for 3D Dynamic Analysis: Rate-sensitive thermo-damage case</caption><thead><tr><th>Load character</th><th>Case</th><th>p&#773;</th><th>q&#773;</th><th>&Delta;t, (dt)</th><th>&sigma;<sub>f</sub></th><th>K<sub>h,lin</sub></th><th>&omega;</th><th>&eta;</th></tr></thead><tbody><tr><td>Triangular pulse</td><td>1c</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>) , (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>50</td><td>1 &middot; 10<sup>-8</sup></td></tr><tr><td>Triangular pulse</td><td>1c*</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>) , (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>50</td><td>1 &middot; 10<sup>-7</sup></td></tr><tr><td>Triangular pulse</td><td>1c**</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>) , (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>50</td><td>1 &middot; 10<sup>-6</sup></td></tr><tr><td>Triangular pulse</td><td>1c***</td><td>1</td><td>0.1</td><td>(8 &middot; 10<sup>-3</sup>) , (1 &middot; 10<sup>-4</sup>)</td><td>0.4</td><td>100</td><td>50</td><td>1 &middot; 10<sup>-5</sup></td></tr></tbody></table>

## 4 Conclusion

We presented regularized variational formulation that allows to construct 3D discrete approximation using an enriched tetrahedral finite element that shows exceptional performance in static and dynamic applications under mechanical and thermo-mechanical loads. The hybrid finite element was created with 4 degrees of freedom for displacements and temperature at the corner nodes, and 4 degrees of freedom for stress vector components and normal heat flux component at the nodes located at element faces. With the choice of the Raviart-Thomas vector field interpolations for stress fields and heat flux we can enforce continuity of stress vector and normal component of heat-flux across element boundary. Such a discrete approximation can be chosen due to the regularization of the variational formulation.

The hybrid-stress-heat-flux variational formulation and chosen discretization provide what appears to be optimal result delivering high accuracy in stress and heat-flux computations. The performance of the proposed element has been demonstrated with a number of numerical simulations results for static problems (both compressible and incompressible case), with proposed hybrid element typically outperforms the standard finite elements based on the Galerkin discrete approximation.

The enhanced accuracy in computing the time-evolution for dynamic simulations was achieved by utilizing the mid-point integration scheme, which is enforcing energy-conservation in elasticity during the free vibration phase or energy-decaying either due to damage or to reducing the contribution of high frequency modes. This enables stable calculations under long-term dynamic and thermal loading, such as in fire-resistance problems. Furthermore, the continuity of the stress and heat flux fields enforced by the Raviart-Thomas interpolation is also useful for wave-propagation problems in dynamics, as it prevents undesirable jumps in these fields on the shared faces of adjacent tetrahedral elements that are typical of the standard displacement-based elements.

The performance of the proposed hybrid tetrahedron element was evaluated for both static and dynamic loading conditions, and also accounting for inelastic behavior. For illustration of element's potential for long-term simulation in inelasticity and accumulation of damage, we employed here a simple isotropic damage model, which is defined
<table><caption>**Table 7** Convergence study for 3D Dynamic Analysis: Thermo-damage case</caption><thead><tr><th rowspan="2">Iter (i)</th><th colspan="5">Residual norm</th></tr><tr><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 4</th><th>Case 5</th></tr></thead><tbody><tr><td>1</td><td>4.03E-01</td><td>4.15E-01</td><td>4.35E-01</td><td>4.58E-01</td><td>4.96E-01</td></tr><tr><td>2</td><td>1.73E-03</td><td>1.68E-03</td><td>1.70E-03</td><td>1.63E-03</td><td>1.36E-03</td></tr><tr><td>3</td><td>2.69E-04</td><td>3.69E-04</td><td>4.49E-04</td><td>5.49E-04</td><td>7.97E-04</td></tr><tr><td>4</td><td>1.21E-06</td><td>1.26E-06</td><td>1.35E-06</td><td>1.48E-06</td><td>3.38E-06</td></tr><tr><td>5</td><td>7.13E-10</td><td>6.61E-10</td><td>7.01E-10</td><td>6.13E-10</td><td>3.16E-10</td></tr></tbody></table>

<table><caption>**Table 8** Convergence study for 3D Dynamic Analysis: Thermo-damage case with temperature-dependent yield stress</caption><thead><tr><th rowspan="2">Iter (i)</th><th colspan="4">Residual norm</th></tr><tr><th>Case 1a</th><th>Case 1b</th><th>Case 1c</th><th>Case 1d</th></tr></thead><tbody><tr><td>1</td><td>4.02E-01</td><td>4.02E-01</td><td>4.03E-01</td><td>4.16E-01</td></tr><tr><td>2</td><td>1.73E-03</td><td>1.72E-03</td><td>1.70E-03</td><td>1.66E-03</td></tr><tr><td>3</td><td>2.60E-04</td><td>2.50E-04</td><td>2.25E-04</td><td>1.93E-04</td></tr><tr><td>4</td><td>1.24E-06</td><td>1.21E-06</td><td>1.33E-06</td><td>1.44E-06</td></tr><tr><td>5</td><td>4.32E-10</td><td>6.10E-10</td><td>5.31E-10</td><td>4.31E-10</td></tr></tbody></table>

<table><caption>**Table 9** Convergence study for 3D Dynamic Analysis: Rate-sensitive thermo-damage case</caption><thead><tr><th rowspan="2">Iter (i)</th><th colspan="4">Residual norm</th></tr><tr><th>Case 1c</th><th>Case 1c*</th><th>Case 1c**</th><th>Case 1c***</th></tr></thead><tbody><tr><td>1</td><td>4.03E-01</td><td>4.03E-01</td><td>4.03E-01</td><td>4.03E-01</td></tr><tr><td>2</td><td>1.70E-03</td><td>1.71E-03</td><td>1.71E-03</td><td>1.71E-03</td></tr><tr><td>3</td><td>2.25E-04</td><td>2.25E-04</td><td>2.25E-04</td><td>2.25E-04</td></tr><tr><td>4</td><td>1.33E-06</td><td>1.33E-06</td><td>1.33E-06</td><td>1.33E-06</td></tr><tr><td>5</td><td>3.43E-10</td><td>8.16E-10</td><td>5.15E-10</td><td>4.91E-10</td></tr></tbody></table>

with two material parameters, the damage stress $\sigma_f$ and linear hardening modulus $K_{h,lin}$. For dynamic loads, the element is capable of describing damage dissipation induced by entering the hardening phase of material behavior. The computational accuracy is enhanced with viscous regularization of the damage model, which allows to simultaneously integrate global balance equations and internal variable evolution equations, and thus provide the second order accuracy. The simulation results have confirmed the proposed hybrid tetrahedral element’s high performance and effectiveness in a number of illustrative numerical simulations for statics, dynamics, stationary and non-stationary heat transfer.

<figure><img src="image_43.png" /><figcaption>Case 2: Traction components</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Case 2: Energy diagram</figcaption></figure>

<caption>**Fig. 23** Case 3 ($\sigma_f$ = 0.6): Traction components and energy time-history plot for the proposed test</caption>

## Appendix A: Shape functions for standard 4-node tetrahedron element

## Appendix B: Lowest-order Raviart-Thomas shape functions
<figure><img src="image_45.png" /><figcaption>**Fig. 24** Shape functions for displacement and temperature fields</figcaption></figure>

## Appendix C: Raviart-Thomas interpolation

### a) Stress field

* Stress vector components

$$ \tau = [\tau_1^1, \tau_2^1, \tau_3^1, \tau_1^2, \tau_2^2, \tau_3^2, \tau_1^3, \tau_2^3, \tau_3^3, \tau_1^4, \tau_2^4, \tau_3^4]^T \quad (89) $$

* Discrete stress vector in Voigt notation

$$ \sigma = [\sigma_{11}^h, \sigma_{22}^h, \sigma_{33}^h, \sigma_{12}^h, \sigma_{23}^h, \sigma_{13}^h, \sigma_{21}^h, \sigma_{32}^h, \sigma_{31}^h]^T \quad (90) $$

* Raviart-Thomas matrix for stress field

$$ S_{\tau} = \begin{bmatrix} \varphi_1^1 & 0 & 0 & \varphi_1^2 & 0 & 0 & \varphi_1^3 & 0 & 0 & \varphi_1^4 & 0 & 0 \\ 0 & \varphi_2^1 & 0 & 0 & \varphi_2^2 & 0 & 0 & \varphi_2^3 & 0 & 0 & \varphi_2^4 & 0 \\ 0 & 0 & \varphi_3^1 & 0 & 0 & \varphi_3^2 & 0 & 0 & \varphi_3^3 & 0 & 0 & \varphi_3^4 \\ \varphi_2^1 & 0 & 0 & \varphi_2^2 & 0 & 0 & \varphi_2^3 & 0 & 0 & \varphi_2^4 & 0 & 0 \\ 0 & \varphi_3^1 & 0 & 0 & \varphi_3^2 & 0 & 0 & \varphi_3^3 & 0 & 0 & \varphi_3^4 & 0 \\ \varphi_3^1 & 0 & 0 & \varphi_3^2 & 0 & 0 & \varphi_3^3 & 0 & 0 & \varphi_3^4 & 0 & 0 \\ 0 & \varphi_1^1 & 0 & 0 & \varphi_1^2 & 0 & 0 & \varphi_1^3 & 0 & 0 & \varphi_1^4 & 0 \\ 0 & 0 & \varphi_2^1 & 0 & 0 & \varphi_2^2 & 0 & 0 & \varphi_2^3 & 0 & 0 & \varphi_2^4 \\ 0 & 0 & \varphi_1^1 & 0 & 0 & \varphi_1^2 & 0 & 0 & \varphi_1^3 & 0 & 0 & \varphi_1^4 \end{bmatrix} \quad (91) $$

<figure><img src="image_46.png" /><figcaption>**Fig. 25** Lowest order Raviart-Thomas shape functions</figcaption></figure>
## b) Heat-flux field

* Vector form of heat-flux components

$$g = [g^1, g^2, g^3, g^4]^T \qquad (92)$$

* Discrete heat-flux vector

$$q = [q_x, q_y, q_z]^T \qquad (93)$$

* Raviart-Thomas matrix for heat-flux field

$$S_g = \begin{bmatrix} \varphi_1^1 & \varphi_1^2 & \varphi_1^3 & \varphi_1^4 \\ \varphi_2^1 & \varphi_2^2 & \varphi_2^3 & \varphi_2^4 \\ \varphi_3^1 & \varphi_3^2 & \varphi_3^3 & \varphi_3^4 \end{bmatrix} \qquad (94)$$

**Acknowledgements** This work was supported by funding from MEAE-Ministry of Europe and Foreign Affaires (project CESPA), IUF-Institut Universitaire de France (project MS1479), MONKS-Ministry of Science for Kanton Sarajevo (project STEMS) and doctoral scholarship funding by French Embassy in Sarajevo, Bosnia and Herzegovina. All this support is gratefully acknowledged.

**Funding** Open access funding provided by Université de Technologie de Compiègne.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Chudoba R, Vorechovský M, Aguilar M, Baktheer A (2022) Coupled sliding-decohesion-compression model for a consistent description of monotonic and fatigue behavior of material interfaces. Comput. Methods Appl. Mech. Engrg. 398:115259

2. Seleš K, Aldakheel F, Tonković Z, Sorić J, Wriggers P (2021) A general phase-field model for fatigue failure in brittle and ductile solids. Comput Mech 67:1431–145

3. Krüger C, Curosu V, Loehnert S (2024) An extended phase-field approach for the efficient simulation of fatigue fracture processes. Int J Numer Methods Eng. 125:e7422

4. Biot MA (1956) Thermoelasticity and irreversible thermodynamics. J Appl Phys 27:250–253

5. Niekamp R, Ibrahimbegovic A, Matthies HG (2014) Finite element analysis of coupled thermoelasticity. Coupled Systems Mechanics 3(1):1–25

6. Dunic V, Busarac N, Slavkovic V, Rosic B, Niekamp R, Matthies H, Slavkovic R, Zivkovic M (2016) A thermo-mechanically coupled

finite strain model considering inelastic heat generation. Continuum Mech Thermodyn 28:993–1007

7. Farhat C, Park KC, Pelerin YD (1991) An Unconditionally Stable Staggered Algorithm for Transient Finite Element Analysis of Coupled Thermoelastic Problems. Comput. Methods Appl. Mech. Engrg. 85:349–365

8. Cannarozzi AA (2001) Francesco Ubertini, A mixed variational method for linear coupled thermoelastic analysis. Int J Solids Struct 38:717–739

9. Ngo VM, Ibrahimbegovic A, Brancherie D (2013) Model for localized failure with thermo-plastic coupling. Theoretical formulation and ED-FEM implementation. Comput Struct 127:2–18

10. Ngo M, Ibrahimbegovic A, Brancherie D (2014) Softening behavior of quasi brittle material under full thermo-mechanical coupling condition: Theoretical formulation and finite element implementation. Comput. Methods Appl. Mech. Engrg. 281:1–28

11. Ibrahimbegovic A, Colliat JB, Davenne L (2005) Thermomechanical coupling in folded plates and non-smooth shells. Comput. Methods Appl. Mech. Engrg. 194:2686–2707

12. Hughes TJR (1987) The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall,

13. Boffi D, Brezzi F, Fortin M, et al (2013) Mixed Finite Element Methods and Applications, volume 44, Springer,

14. Brezzi F, Fortin M (2012) Mixed and Hybrid Finite Element Methods, volume 15, Springer Science & Business Media,

15. Hughes TJR, Masud A, Harari I (1995) Numerical assessment of some membrane elements with drilling degrees of freedom. Comput Struct 55(2):297–314

16. Ibrahimbegovic A, Frey F (1993) Geometrically non-linear method of incompatible modes in application to finite elasticity with independent rotations. Internat. J. Numer. Methods Engrg. 36(24):4185–4200

17. Stein E, Ahmad R (1977) An equilibrium method for stress calculation using finite element displacement models. Comput. Methods Appl. Mech. Engrg. 10(2):175–198

18. Johnson C (1987) Numerical Solution of Partial Differential Equations By the Finite Element Method, Cambridge University Press,

19. Hughes TJR, Cottrell JA, Bazilevs Y (2005) Isogeometric analysis: Cad, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Engrg. 194(39–41):4135–4195

20. Moreno-Navarro P, Ibrahimbegovic A, Ospina A (2020) Multi-field variational formulations and mixed finite element approximations for electrostatics and magnetostatics. Comput Mech 65(1):41–59

21. Bossavit A (1988) Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism. IEE Proc. A 135(8):493–500

22. Dubois F (2010) Dual Raviart-Thomas mixed finite elements, arXiv preprint arXiv:1012.1691

23. Ervin VJ (2012) Computational bases for RTk and BDMk on triangles. Comput Math Appl 64(8):2765–2774

24. Hughes TJR, Brezzi F (1989) On drilling degrees of freedom. Comput. Methods Appl. Mech. Engrg. 72(1):105–121

25. Ibrahimbegovic A (1993) Mixed finite element with drilling rotations for plane problems in finite elasticity. Comput. Methods Appl. Mech. Engrg. 107(1–2):225–238

26. Ibrahimbegovic A, Taylor RL, Wilson EL (1990) A robust quadrilateral membrane finite element with drilling degrees of freedom. Internat. J. Numer. Methods Engrg. 30(3):445–457

27. Braess D, Klaas O, Niekamp R, Stein E, Wobschal F (1995) Error indicators for mixed finite elements in 2-dimensional linear elasticity. Comput. Methods Appl. Mech. Engrg. 127(1–4):345–356

28. Bathe KJ, Baig MMI (2005) On a composite implicit time integration procedure for nonlinear dynamics. Comput Struct 83(31–32):2513–2524
29. Armero F (2006) Energy-dissipative momentum-conserving time-
stepping algorithms for finite strain multiplicative plasticity. Com-
put. Methods Appl. Mech. Engrg. 195(37–40):4862–4889

30. Leyendecker S, Betsch P, Steinmann P (2006) Objective energy?momentum conserving integration for the constrained dynamics of geometrically exact beams. Comput. Methods Appl. Mech. Engrg. 195(19–22):2313–2333

31. Nguyen CU, Ibrahimbegovic A (2020) Visco-plasticity stress-
based solid dynamics formulation and time-stepping algorithms
for stiff case. Int J Solids Struct 196–197:154–170

32. Nguyen CU, Ibrahimbegovic A (2020) Hybrid-stress triangular finite element with enhanced performance for statics and dynamics. Comput. Methods Appl. Mech. Engrg. 372:113381

33. Nguyen CU, Ibrahimbegovic A (2022) Regularized variational for-
mulation for nonlinear dynamics of viscoplastic plates. Comput.
Methods Appl. Mech. Engrg. 397:115146

34. Ibrahimbegovic A, Mamouri S (2002) Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations. Comput. Methods Appl. Mech. Engrg. 191(37–38):4241–4258

35. Ibrahimbegovic A (2009) Nonlinear Solid Mechanics: Theoretical Formulations and Finite Element Solution Methods, volume 160, Springer Science & Business Media,

36. Artioli E, Auricchio F, da Veiga LB (2007) Second-order accurate integration algorithms for von-Mises plasticity with a nonlinear kinematic hardening mechanism. Comput. Methods Appl. Mech. Engrg. 196(9–12):1827–1846

37. Hughes TJR, Taylor RL (1978) Unconditionally stable algorithms for quasi-static elasto visco-plastic finite element analysis. Comput Struct 8(2):169–173

38. Zienkiewicz OC, Taylor RL, Zhu JZ (2005) The Finite Element Method: Its Basis and Fundamentals, Elsevier,

39. Taylor RL (2014) FEAP-Finite Element Analysis Program, University of California, Berkeley, URL http://projects.ce.berkeley.edu/feap/

40. Auricchio F, Da Veiga LB, Hughes TJ, Reali A, Sangalli G (2012) Isogeometric collocation for elastostatics and explicit dynamics. Comput. Methods Appl. Mech. Engrg. 249:2–14

41. Timoshenko S, Timoshenko S, Goodier JN (1951) Theory of Elasticity. By S. Timoshenko and JN Goodier, McGraw-Hill book Company

42. Cook RD (1974) Improved two-dimensional finite element. J Struct Div 100(9):1851–1865

43. Fissolo A, Amiable S, Ancelet O, Mermaz F, Stelmaszyk JM et al (2009) Crack initiation under thermal fatigue: An overview of CEA experience. part I: Thermal fatigue appears to be more damaging than uniaxial isothermal fatigue. International Journal of Fatigue 31(3):587–600 ff10.1016/j.ijfatigue.2008.03.038ff. ffhal-00438553ff

44. Le J-L, Bažant Z (2022) Spontaneous Collapse Mechanism of World Trade Center Twin Towers and Progressive Collapse in General. J Struct Eng 148(6):04022065

45. Bugrimov AL, Bobrova IA, Chukalovskaya EM (2020) Analysis of the thermodynamic criteria for long-term strength of materials. J Phys: Conf Ser 1560:012015

46. Gong F, Ueda T, Wang Y, Zhang D, Wang Z (2017) Mesoscale simulation of fatigue behavior of concrete materials damaged by freeze-thaw cycles. Constr Build Mater 144:702–716

47. Abdoh DA, Yin BB, Liew KM (2021) A phase-field thermomechanical framework for modeling failure and crack evolution in glass panes under fire. Comput. Methods Appl. Mech. Engrg. 385:114068

48. Luo C, Sanaviab L, De Lorenzis L (2023) Phase-field modeling of drying-induced cracks: Choice of coupling and study of homogeneous and localized damage. Comput. Methods Appl. Mech. Engrg. 410:115962

49. Ba K, Gakwaya A (2018) Thermomechanical total Lagrangian SPH formulation for solid mechanics in large deformation problems. Comput. Methods Appl. Mech. Engrg. 342:458–473

50. Yuan H, Tang W, He J, Guan X (2024) An efficient phase-field model for fatigue fracture in viscoelastic solids using cyclic load decomposition and damage superposition. Comput Methods Appl Mech Eng 418:116469

51. Koniorczyka M, Gawina D, Schrefler BA (2015) Modeling evolution of frost damage in fully saturated porous materials exposed to variable hygro-thermal conditions. Comput. Methods Appl. Mech. Engrg. 297:38–61

52. Andréa D, Levrauta B, Tessier-Doyenb N, Hugera M (2017) A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expansion mismatch of two-phase materials. Comput. Methods Appl. Mech. Engrg. 318:898–916

53. Bie Y, Ren H, Rabczuk T, Bui TQ, Wei Y (2024) The fully coupled thermo mechanical dual-horizon peridynamic correspondence damage model for homogeneous and heterogeneous materials. Comput Methods Appl Mech Eng 420:116730

54. Karavelic E, Ibrahimbegovic A, Dolarevic S (2019) Multi-surface plasticity model for concrete with 3D hardening/softening failure modes for tension, compression and shear. Comput Struct 221:74–90

**Publisher's Note** Springer Nature remains neutral with regard to juris-
dictional claims in published maps and institutional affiliations.