# Interface Models in Coupled Thermoelasticity

Michele Serpilli, Serge Dumont, Raffaella Rizzoni, Frédéric Lebon

► To cite this version:

Michele Serpilli, Serge Dumont, Raffaella Rizzoni, Frédéric Lebon. Interface Models in Coupled Thermoelasticity. Technologies, 2021, 9 (1), pp.17. 10.3390/technologies9010017. hal-03354120

HAL Id: hal-03354120

https://hal.science/hal-03354120

Submitted on 22 Jun 2023

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>
Article

# Interface Models in Coupled Thermoelasticity

Michele Serpilli <fnref n="1" />,*, Serge Dumont <fnref n="2" />, Raffaella Rizzoni <fnref n="3" /> and Frédéric Lebon <fnref n="4" />

**Abstract:** This work proposes new interface conditions between the layers of a three-dimensional composite structure in the framework of coupled thermoelasticity. More precisely, the mechanical behavior of two linear isotropic thermoelastic solids, bonded together by a thin layer, constituted of a linear isotropic thermoelastic material, is studied by means of an asymptotic analysis. After defining a small parameter $\epsilon$, which tends to zero, associated with the thickness and constitutive coefficients of the intermediate layer, two different limit models and their associated limit problems, the so-called *soft* and *hard* thermoelastic interface models, are characterized. The asymptotic expansion method is reviewed by taking into account the effect of higher-order terms and defining a generalized thermoelastic interface law which comprises the above aforementioned models, as presented previously. A numerical example is presented to show the efficiency of the proposed methodology, based on a finite element approach developed previously.

**Keywords:** interfaces; asymptotic analysis; coupled thermoelasticity

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

check for
updates

**Citation:** Serpilli, M.; Dumont, S.; Rizzoni, R.; Lebon, F. Interface Models in Coupled Thermoelasticity. *Technologies* 2021, *9*, 17.
https://doi.org/10.3390/technologies9010017

Received: 12 February 2021

Accepted: 2 March 2021

Published: 4 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. License MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The use of composite structures, obtained by bonding together simpler structural members, has spread in all fields of engineering in the last decades. On the one hand, the structural assembly presents a significant improvement of the mechanical properties and an enhancement of its performances. On the other hand, the bonded joints among the composite constituents may cause a jump of the physical fields at the interface level and radically modify the global mechanical response. Thus, the correct modeling of composite interfaces is crucial in the understanding and design of complex structures.

From a theoretical point of view, the bonded region is considered as a thin interphase between two adjacent parts. By letting the thickness of this layer tend to zero, the interphase is reduced into a two-dimensional surface, called imperfect interface, where ad-hoc transmission conditions in terms of the representative physical fields are prescribed. The contact laws can be derived by means of classical variational tools and more refined mathematical techniques, in different physical frameworks, involving uncoupled (thermal conduction and elasticity) and coupled (piezoelectricity and multiphysics) phenomena.

Concerning the thermal (or electrical) conduction case, two main interface laws have been formulated: the lowly-conducting (LC) or Kapitza's model and highly-conducting (HC) model. The LC model provides a discontinuity of the temperature field (electric potential) and a continuity of the normal heat flow (electric displacement) across the interface (see, e.g., [1–3]). The HC model gives rise to two-dimensional Young–Laplace equation, defined on the interface, depending on the jump of the normal heat flow (electric displacement) and maintaining the temperature (electric potential) continuous (see, e.g., [4,5]). A unifying approach of a general imperfect interface model, involving the concurrent jump of both the temperature field and the normal heat flow, recovering both the LC and HC

<footnote n="1">¹ Department of Civil and Building Engineering, and Architecture, Università Politecnica delle Marche, 60121 Ancona, Italy</footnote>

<footnote n="2">² IMAG CNRS UMR 5149, University of Nîmes, 30000 Nîmes, France; serge.dumont@unimes.fr</footnote>

<footnote n="3">³ Department of Engineering, University of Ferrara, 44122 Ferrara, Italy; raffaella.rizzoni@unife.it</footnote>

<footnote n="4">⁴ CNRS, Centrale Marseille, Laboratoire de Mécanique et d'Acoustique, Aix-Marseille University, 13453 Marseille, France; lebon@lma.cnrs-mrs.fr</footnote>

<footnote>* Correspondence: m.serpilli@univpm.it</footnote>
models, was proposed by [6,7]. Concerning the linear elastic case, three types of imperfect interfaces have been proposed: the spring-layer interface model (SL) (soft interface), the coherent interface (CI) (rigid interface), and the general imperfect interface. The SL models considers that the traction vector is continuous across the interface, while the displacement presents a jump linearly proportional to the traction vector (see, e.g., [8,9]). The CI model has been developed for continuum theories with surface effects and nano-sized materials (see, e.g., [10–12]): the traction vector suffers a jump, while the displacement field is continuous across the interface. Finally, in the general imperfect model, both the displacement and normal traction fields are discontinuous across the interface [13,14].

The asymptotic expansions method and convergence approaches represent mathematical tools, usually employed in the derivation and justification of classical thin structures and layered plates [15–19]. These methodologies are based on the behavior of the problem solution, when a small parameter $\epsilon$, related to the thickness of the interphase, tends to zero. Considering that the material properties of the intermediate layer depend on $\epsilon^p$, different limit behaviors can be derived by means of the asymptotic analysis: for $p = 1$, an SL interface model can be recovered (see, e.g., [20,21]); and, for $p = -1$, the CI interface model is mathematically justified by means of strong convergence arguments in [22,23]. Within the framework of a higher-order theory, assuming the interphase elastic constants are independent of the small thickness ($p = 0$), the asymptotic analysis yields to a general stiff imperfect interface condition, prescribing both the jumps of the displacement and traction vector fields and recovering as a particular case the perfect contact conditions at the zeroth-order [24–28]. The above transmission conditions have been generalized by considering some multiphysics and multifield couplings, such as in piezoelectricity and magneto-electro-thermo-elasticity [29,30], poroelasticity [31], and micropolar elasticity [32].

The goal of the present work is to identify the interface limit models of a composite constituted by a thin thermoelastic layer surrounded by two thermoelastic bodies in the framework of dynamic coupled thermoelasticity. Different situations are analyzed by varying the stiffnesses ratios between the middle layer and the adherents: namely, the *soft* thermoelastic lowly conducting interface, where the intermediate material coefficients have the order of magnitude $\epsilon$ with respect to those of the surrounding bodies, and the *hard* thermoelastic moderately conducting interface, where the constitutive parameters have the same order of magnitude. Following the ideas of [33], a generalized interface law is derived, comprising the aforementioned behaviors. A numerical investigation was performed in the framework of the finite element method (FEM), employing the approach developed in [34] for multiphysics problems, in order to assess the validity of the asymptotic models. Convergence results and a comparison between the full 3D model and the generalized interface problem are given.

## 2. Position of the Problem

In the sequel, Greek indices range in the set {1, 2}, Latin indices range in the set {1, 2, 3}, and the Einstein's summation convention with respect to the repeated indices is adopted. Let us consider a three-dimensional Euclidian space identified by $\mathbb{R}^3$ and such that the three vectors $e_i$ form an orthonormal basis. Let $M^n$ be the space of $n \times n$ square matrices. We introduce the following notations for the inner products: $\mathbf{a} \cdot \mathbf{b} := a_i b_i$, for all vectors $\mathbf{a} = (a_i)$ and $\mathbf{b} = (b_i)$ in $\mathbb{R}^3$ and $\mathbf{A} : \mathbf{B} := A_{ij} B_{ij}$, for all $\mathbf{A} = (A_{ij})$ and $\mathbf{B} = (B_{ij})$ in $M^3$.

Let us define a small parameter $0 < \epsilon < 1$. We consider the assembly constituted of two solids $\Omega_{\pm} \subset \mathbb{R}^3$, called the adherents, bonded together by an intermediate thin layer $B^{\epsilon} := S \times (-\frac{\epsilon}{2}, \frac{\epsilon}{2})$ of thickness $\epsilon$, called the adhesive, with cross-section $S \subset \mathbb{R}^2$. In the following, $B^{\epsilon}$ and $S$ are called interphase and interface, respectively. Let $S_{\pm}^{\epsilon}$ be the plane interfaces between the interphase and the adherents and let $\Omega^{\epsilon} := \Omega_{+}^{\epsilon} \cup B^{\epsilon} \cup \Omega_{-}^{\epsilon}$ denote the composite system comprising the interphase and the adherents (cf. Figure 1a).
<figure><img src="image_4.png" /><figcaption>Figure 1. Initial (a); rescaled (b); and limit (c) configurations of the composite.</figcaption></figure>

Let us assume that $\Omega_{\pm}^{\epsilon}$ and $B^{\epsilon}$ are constituted by three homogeneous linear isotropic thermoelastic materials, whose constitutive laws are defined as follows:

$$ egin{cases} 
\sigma_{ij}^{\epsilon} = \lambda^{\epsilon} e_{pp}^{\epsilon} \delta_{ij} + 2\mu^{\epsilon} e_{ij}^{\epsilon} - \beta^{\epsilon} \theta^{\epsilon} \delta_{ij}, \\ 
S^{\epsilon} = c_{v}^{\epsilon} \theta^{\epsilon} + \beta^{\epsilon} e_{pp}^{\epsilon}, \\ 
q_{i}^{\epsilon} = -k^{\epsilon} \theta_{,i}^{\epsilon} 
\end{cases} 
\quad (1) $$

where $\sigma^{\epsilon} = (\sigma_{ij}^{\epsilon})$ is the Cauchy stress tensor, associated with the linearized strain tensor $e^{\epsilon} = (e_{ij}^{\epsilon}) := \frac{1}{2}(u_{ij}^{\epsilon} + u_{ji}^{\epsilon})$, $S^{\epsilon}$ represents the thermodynamic entropy and $q^{\epsilon} = (q_{i}^{\epsilon})$ is the heat flow field. Constants $\lambda^{\epsilon}, \mu^{\epsilon}, \beta^{\epsilon}, c_{v}^{\epsilon}$, and $k^{\epsilon}$ represent the Lamé's constants, the thermal stress coefficient, the calorific capacity, and the thermal conductivity, respectively.

The thermoelastic state is defined by the couple $s^{\epsilon} := (u^{\epsilon}, \theta^{\epsilon})$, where $u^{\epsilon}$ and $\theta^{\epsilon}$ are the displacement field and variation of temperature, respectively. The thermoelastic composite is subject to body forces $f^{\epsilon} = (f_{i}^{\epsilon})$: $\Omega_{\pm}^{\epsilon} \times (0, T) \to \mathbb{R}^3$ and heat source $h^{\epsilon}: \Omega_{\pm}^{\epsilon} \times (0, T) \to \mathbb{R}$, applied on the top and bottom bodies, while all thermo-mechanical loadings and inertia forces are neglected in the intermediate layer $B^{\epsilon}$. The thermoelastic state $s^{\epsilon}$ verifies the following coupled thermoelasticity system:

$$ egin{cases} 
\rho^{\epsilon} \ddot{u}^{\epsilon} - \operatorname{div}^{\epsilon} \sigma^{\epsilon} = \mathbf{f}^{\epsilon} & \text{in } \Omega_{\pm}^{\epsilon} \times (0, T), \\ 
\dot{S}^{\epsilon} + \frac{1}{T_0} \operatorname{div}^{\epsilon} \mathbf{q}^{\epsilon} = h^{\epsilon} & \text{in } \Omega_{\pm}^{\epsilon} \times (0, T), 
\end{cases} 
\quad 
\begin{cases} 
\operatorname{div}^{\epsilon} \sigma^{\epsilon} = \mathbf{0} & \text{in } B^{\epsilon} \times (0, T), \\ 
\dot{S}^{\epsilon} + \frac{1}{T_0} \operatorname{div}^{\epsilon} \mathbf{q}^{\epsilon} = 0 & \text{in } B^{\epsilon} \times (0, T), 
\end{cases} 
\quad (2) $$

where $\dot{f} = \partial_t f$ denotes the time derivative of $f$ and $T_0$ represents a reference temperature. The transmission conditions across the interfaces $S^{+},^{\epsilon}$ and $S^{-},^{\epsilon}$ implies the continuity of the state $s^{\epsilon}$ and of its normal dual counterpart with respect to $S^{\pm},^{\epsilon}$, meaning that $[u^{\epsilon}] = \mathbf{0}$, $[\theta^{\epsilon}] = 0$, $[\sigma^{\epsilon}e_3] = \mathbf{0}$, $[q^{\epsilon} \cdot e_3] = 0$ on $S^{\pm},^{\epsilon} \times (0, T)$, where $[f]$ stands for the jump function evaluated at the interface $S^{\pm},^{\epsilon}$. The boundary conditions are posed on $\Gamma^{\epsilon} \times (0, T)$, with $\Gamma^{\epsilon} := \Gamma^{+},^{\epsilon} \cup \Gamma^{-},^{\epsilon}$; we recall that $\Gamma^{\epsilon} = \Gamma_{g}^{\epsilon} \cup \Gamma_{u}^{\epsilon}$. For simplicity, we assume homogeneous boundary conditions on $\Gamma_{u}^{\epsilon} \times (0, T)$, concerning displacements and temperature, and non-homogeneous boundary conditions on $\Gamma_{g}^{\epsilon} \times (0, T)$, concerning surface forces $g^{\epsilon} = (g_{i}^{\epsilon})$ and surface heat flow $q^{\epsilon}$. Hence, one has: $\sigma^{\epsilon}n^{\epsilon} = g^{\epsilon}$ and $-q^{\epsilon} \cdot n^{\epsilon} = q^{\epsilon}$ on $\Gamma_{g}^{\epsilon} \times (0, T)$, and $u^{\epsilon} = \mathbf{0}$ and $\theta^{\epsilon} = 0$ on $\Gamma_{u}^{\epsilon} \times (0, T)$, where $n^{\epsilon} = (n_{i}^{\epsilon})$ is the outer unit normal vector to $\partial\Omega^{\epsilon}$. The initial conditions are posed in $\Omega^{\epsilon}$. Let $\theta_{in}^{\epsilon}, u_{in}^{\epsilon}$, and $\dot{u}_{in}^{\epsilon}$ be, respectively, the variation of temperature, the displacement, and velocity fields at time $t = 0$; one has $\theta^{\epsilon}(x^{\epsilon}, 0) = \theta^{\epsilon}(0) = \theta_{in}^{\epsilon} u^{\epsilon}(x^{\epsilon}, 0) = u^{\epsilon}(0) = u_{in}^{\epsilon}$ and $\dot{u}^{\epsilon}(x^{\epsilon}, 0) = \dot{u}^{\epsilon}(0) = \dot{u}_{in}^{\epsilon}$ in $\Omega^{\epsilon}$.

Let us introduce the functional spaces $V(\Omega^{\epsilon}) := \{v^{\epsilon} \in H^1(\Omega^{\epsilon}); v^{\epsilon} = 0 \text{ on } \Gamma_{u}^{\epsilon}\}$ and $V(\Omega^{\epsilon}) := [V(\Omega^{\epsilon})]^3$. Given a certain state $s^{\epsilon} := (u^{\epsilon}, \theta^{\epsilon}) \in V(\Omega^{\epsilon}) := V(\Omega^{\epsilon}) \times V(\Omega^{\epsilon})$, for all
test functions $r^\epsilon = (\mathbf{v}^\epsilon, \xi^\epsilon) \in \mathbb{V}(\Omega^\epsilon)$ and for any fixed $t \in (0, T)$, we introduce the following bilinear and linear forms:

$$A^\epsilon(s^\epsilon, r^\epsilon) := \int_{\Omega^\epsilon} \left\{ \rho^\epsilon \ddot{\mathbf{u}}^\epsilon \cdot \mathbf{v}^\epsilon + \sigma^\epsilon : \mathbf{e}^\epsilon(\mathbf{v}^\epsilon) + \dot{S}^\epsilon \xi^\epsilon - \frac{1}{T_0} \mathbf{q}^\epsilon \cdot \nabla^\epsilon \xi^\epsilon \right\} dx^\epsilon, \quad (3)$$

$$L^\epsilon(r^\epsilon) := \int_{\Omega_\pm^\epsilon} \{\mathbf{f}^\epsilon \cdot \mathbf{v}^\epsilon + h^\epsilon \xi^\epsilon\} dx^\epsilon + \int_{\Gamma_\pm^\epsilon} \{\mathbf{g}^\epsilon \cdot \mathbf{v}^\epsilon + q^\epsilon \xi^\epsilon\} d\Gamma^\epsilon. \quad (4)$$

The variational form of the coupled thermoelastic system defined on the variable domain $\Omega^\epsilon$ reads as follows:

$$\left\{ \begin{array}{l} \text{Find } s^\epsilon(t) \in \mathbb{V}(\Omega^\epsilon), t \in (0, T), \text{ such that} \\ \bar{A}_-^\epsilon(s^\epsilon, r^\epsilon) + \bar{A}_+^\epsilon(s^\epsilon, r^\epsilon) + \hat{A}^\epsilon(s^\epsilon, r^\epsilon) = L^\epsilon(r^\epsilon), \end{array} \right. \quad (5)$$

for all $r^\epsilon \in \mathbb{V}(\Omega^\epsilon)$, with initial condition $\theta_{in}^\epsilon$, $\mathbf{u}_{in}^\epsilon$, and $\dot{\mathbf{u}}_{in}^\epsilon$. The coupled hyperbolic–parabolic equations associated with variation problem (5) imply a degenerate system. Hence, the standard existence theorems are not applicable. For instance, in [35,36], by applying the pseudo-monotone theory, a weak solution is provided for a sufficiently small thermal stress coupling coefficient. In [37,38], a solution to the implicit evolution equation is derived after time-differentiation of the equilibrium equation provided by sufficiently smooth data of the problem. Under suitable regularity properties of the initial data, source and boundary values, and constitutive parameters, the well-posedness of thermo-electro-elastic evolution problem is extensively discussed in [39]: the proof of existence, uniqueness, and regularity of the solution has been obtained through the Faedo–Galerkin method. The existence and uniqueness theorems have also been extended to the thermo-electro-magneto-elastic case [40] and can be easily adapted to the present coupled thermoelastic problem.

## Rescaling

To study the asymptotic behavior of the solution of problem (5) when $\epsilon$ tends to zero, we rewrite the problem on a fixed domain $\Omega$ independent of $\epsilon$. By using the approach of [15], we consider the bijection $\pi^\epsilon: x \in \bar{\Omega} \mapsto x^\epsilon \in \bar{\Omega}^\epsilon$ given by

$$\pi^\epsilon: \begin{cases} \bar{\pi}^\epsilon(x_1, x_2, x_3) = (x_1, x_2, x_3 \mp \frac{1}{2}(1-\epsilon)), & \text{for all } x \in \bar{\Omega}_\pm, \\ \hat{\pi}^\epsilon(x_1, x_2, x_3) = (x_1, x_2, \epsilon x_3), & \text{for all } x \in \bar{B}, \end{cases} \quad (6)$$

where, after the change of variables, the adherents occupy $\Omega_\pm := \Omega_\pm^\epsilon \pm \frac{1}{2}(1-\epsilon)e_3$ and the interphase $B = \{x \in \mathbb{R}^3 : (x_1, x_2) \in S, |x_3| < \frac{1}{2}\}$. The sets $S_\pm = \{x \in \mathbb{R}^3 : (x_1, x_2) \in S, x_3 = \pm \frac{1}{2}\}$ denote the interfaces between $B$ and $\Omega_\pm$ and $\Omega = \Omega_+ \cup \Omega_- \cup B$ is the rescaled configuration of the composite. Lastly, $\Gamma_u$ and $\Gamma_g$ indicate the images through $\pi^\epsilon$ of $\Gamma_u^\epsilon$ and $\Gamma_g^\epsilon$ (cf. Figure 1b). Consequently, $\frac{\partial}{\partial x_\alpha^\epsilon} = \frac{\partial}{\partial x_\alpha}$ and $\frac{\partial}{\partial x_3^\epsilon} = \frac{\partial}{\partial x_3}$ in $\Omega_\pm$, and $\frac{\partial}{\partial x_\alpha^\epsilon} = \frac{\partial}{\partial x_\alpha}$ and $\frac{\partial}{\partial x_3^\epsilon} = \frac{1}{\epsilon} \frac{\partial}{\partial x_3}$ in $B$. In the sequel, only if necessary, $\mathfrak{s}^\epsilon = (\bar{\mathbf{u}}^\epsilon, \bar{\theta}^\epsilon)$ and $\mathfrak{s}^\epsilon = (\hat{\mathbf{u}}^\epsilon, \hat{\theta}^\epsilon)$ denote the restrictions of functions $s^\epsilon = (\mathbf{u}^\epsilon, \theta^\epsilon)$ to $\Omega_\pm$ and $B$.

The constitutive coefficients of $\Omega_\pm^\epsilon$ are assumed to be independent of $\epsilon$, while the constitutive coefficients of $B^\epsilon$ present the following dependences on $\epsilon$: $\hat{\lambda}^\epsilon = \epsilon^p \hat{\lambda}$, $\hat{\mu}^\epsilon = \epsilon^p \hat{\mu}$, $\hat{\beta}^\epsilon = \epsilon^p \hat{\beta}$, $\hat{c}_v^\epsilon = \epsilon^p \hat{c}_v$, and $\hat{k}^{m,\epsilon} = \epsilon^p \hat{k}$, with $p \in \{0, 1\}$. Two different limit behaviors are characterized according to the choice of the exponent $p$: by choosing $p = 1$, a model for a soft thermoelastic interface with low conductivity is deduced; and, when $p = 0$, a model for a hard thermoelastic interface with moderate conductivity is obtained. Finally, the data, unknowns, and test functions verify the following scaling assumptions: $s^\epsilon(x^\epsilon) = s^\epsilon(x)$, $r^\epsilon(x^\epsilon) = r(x)$, $x \in \Omega$, $\mathbf{f}^\epsilon(x^\epsilon) = \mathbf{f}(x)$, $h^\epsilon(x^\epsilon) = h(x)$, $x \in \Omega_\pm$, $\mathbf{g}^\epsilon(x^\epsilon) = \mathbf{g}(x)$, $q^\epsilon(x^\epsilon) = q(x)$, $x \in \Gamma_g$. Thus, $L^\epsilon(r^\epsilon) = L(r)$.
According to the previous hypothesis, problem (5) can be reformulated on a fixed domain Ω independent of ε. Thus, the following rescaled problem (in the sequel, we omit the explicit dependences on time t of the unknowns and data) is obtained:

$$ \begin{cases} \text{Find } s^{\epsilon} \in \mathbb{V}(\Omega), t \in (0, T), \text{ such that} \\ \bar{A}_{-}(s^{\epsilon}, r) + \bar{A}_{+}(s^{\epsilon}, r) + \epsilon^{p+1}\hat{A}(s^{\epsilon}, r) = L(r), \end{cases} \quad (7) $$

for all $r \in \mathbb{V}(\Omega)$, $p \in \{0, 1\}$, with initial condition $\theta_{in}$, $\mathbf{u}_{in}$, and $\dot{\mathbf{u}}_{in}$, where

$$ \bar{A}_{\pm}(s^{\epsilon}, r) := \int_{\Omega_{\pm}} \left\{ \rho^{\epsilon} \ddot{\mathbf{u}}^{\epsilon} \cdot \mathbf{v}^{\epsilon} + \sigma^{\epsilon} : \mathbf{e}(\mathbf{v}) + \dot{S}^{\epsilon} \xi - \frac{1}{T_0} \mathbf{q}^{\epsilon} \cdot \nabla \xi \right\} dx, \quad (8) $$

$$ \hat{A}(s^{\epsilon}, r) := \frac{1}{\epsilon^2} a_0(s^{\epsilon}, r) + \frac{1}{\epsilon} a_1(s^{\epsilon}, r) + a_2(s^{\epsilon}, r), \quad (9) $$

where

$$ a_0(s^{\epsilon}, r) := \int_B \left\{ \hat{\mathbf{K}} \mathbf{u}_{,3}^{\epsilon} \cdot \mathbf{v}_{,3} + \frac{\hat{k}}{T_0} \theta_{,3}^{\epsilon} \xi_{,3} \right\} dx, \quad (10) $$

$$ a_1(s^{\epsilon}, r) := \int_B \left\{ (\hat{\mathbf{K}}^{\alpha})^T \mathbf{u}_{,\alpha}^{\epsilon} \cdot \mathbf{v}_{,3} + \hat{\mathbf{K}}^{\alpha} \mathbf{u}_{,3}^{\epsilon} \cdot \mathbf{v}_{,\alpha} - \hat{\beta} \theta^{\epsilon} v_{3,3} + \hat{\beta} \dot{u}_{3,3}^{\epsilon} \xi \right\} dx, \quad (11) $$

$$ a_2(s^{\epsilon}, r) := \int_B \left\{ \hat{\mathbf{K}}^{\alpha\beta} \mathbf{u}_{,\beta}^{\epsilon} \cdot \mathbf{v}_{,\alpha} - \hat{\beta}^{\epsilon} \theta^{\epsilon} v_{\tau,\tau} + \frac{\hat{k}}{T_0} \theta_{,\alpha}^{\epsilon} \xi_{,\alpha} + (\hat{c}_v \dot{\theta}^{\epsilon} + \hat{\beta} \dot{u}_{\alpha,\alpha}^{\epsilon}) \xi \right\} dx \quad (12) $$

and

$$ \hat{\mathbf{K}} := \begin{bmatrix} \hat{\mu} & 0 & 0 \\ 0 & \hat{\mu} & 0 \\ 0 & 0 & 2\hat{\mu} + \hat{\lambda} \end{bmatrix}, \quad \hat{\mathbf{K}}^1 := \begin{bmatrix} 0 & 0 & \hat{\lambda} \\ 0 & 0 & 0 \\ \hat{\mu} & 0 & 0 \end{bmatrix}, \quad \hat{\mathbf{K}}^2 := \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & \hat{\lambda} \\ 0 & \hat{\mu} & 0 \end{bmatrix}, \quad (13) $$

$$ \hat{\mathbf{K}}^{11} := \begin{bmatrix} 2\hat{\mu} + \hat{\lambda} & 0 & 0 \\ 0 & \hat{\mu} & 0 \\ 0 & 0 & \hat{\mu} \end{bmatrix}, \quad \hat{\mathbf{K}}^{22} := \begin{bmatrix} \hat{\mu} & 0 & 0 \\ 0 & 2\hat{\mu} + \hat{\lambda} & 0 \\ 0 & 0 & \hat{\mu} \end{bmatrix}, \quad (14) $$

$$ \hat{\mathbf{K}}^{12} := \begin{bmatrix} 0 & \hat{\lambda} & 0 \\ \hat{\mu} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \hat{\mathbf{K}}^{21} = (\hat{\mathbf{K}}^{12})^T. \quad (15) $$

Now, an asymptotic analysis of the rescaled problem (7) can be performed. Since the rescaled problem (7) has a polynomial structure with respect to the small parameter $\epsilon$, we can look for the solution $s^{\epsilon}$ of the problem as a series of powers of $\epsilon$:

$$ s^{\epsilon} = s^0 + \epsilon s^1 + \epsilon^2 s^2 + \dots, \quad \bar{s}^{\epsilon} = \bar{s}^0 + \epsilon \bar{s}^1 + \epsilon^2 \bar{s}^2 + \dots, \quad \hat{s}^{\epsilon} = \hat{s}^0 + \epsilon \hat{s}^1 + \epsilon^2 \hat{s}^2 + \dots \quad (16) $$

where $\bar{s}^{\epsilon} = s^{\epsilon} \circ \bar{\pi}^{\epsilon}$ and $\hat{s}^{\epsilon} = s^{\epsilon} \circ \hat{\pi}^{\epsilon}$. By substituting (16) into the rescaled problem (7), and by identifying the terms with identical power of $\epsilon$, as customary, a set of variational problems is obtained to be solved in order to characterize the limit thermoelastic state $s^0$, the first-order corrector term $s^1$ and their associated limit problem, for $p \in \{0, 1\}$. The order 1 can be considered as a corrector term of the order 0, giving a better approximation of the initial model.

### 3. The Soft Thermoelastic Interface Model

In this section, the limit model for a soft thermoelastic interface model, corresponding to an adhesive which is weaker with respect to the adherents, is derived. By choosing $p = 1$ and injecting (16) into (7), the following set of variational problems $\mathcal{P}_q$ is obtained:

$$ \left\{ \begin{array}{ll} \mathcal{P}_0: & \bar{A}_-(s^0, r) + \bar{A}_+(s^0, r) + a_0(s^0, r) = L(r), \\ \mathcal{P}_1: & \bar{A}_-(s^1, r) + \bar{A}_+(s^1, r) + a_0(s^1, r) + a_1(s^0, r) = 0, \\ \mathcal{P}_q: & \bar{A}_-(s^q, r) + \bar{A}_+(s^q, r) + a_0(s^q, r) + a_1(s^{q-1}, r) + a_2(s^{q-2}, r) = 0, \quad q \ge 2 \end{array} \right. \quad (17) $$
In the sequel, the limit problems at order 0 and order 1 are presented, by skipping all the mathematical technicalities involved in the solution of problems $\mathcal{P}_q$ (see [33] for a detailed description of the asymptotic analysis).

* Order 0 model

$$
\begin{equation}
\begin{array}{ll}
\text{Governing equations} & \text{Transmission conditions on } S_{\pm} \\
\left\{
\begin{array}{ll}
\rho \ddot{\bar{\mathbf{u}}}^0 - \operatorname{div} \bar{\sigma}^0 = \mathbf{f} & \text{in } \Omega_{\pm}, \\
\dot{\bar{S}}^0 + \frac{1}{T_0} \operatorname{div} \bar{\mathbf{q}}^0 = h & \text{in } \Omega_{\pm}, \\
\bar{\sigma}^0 \mathbf{n} = \mathbf{g} & \text{on } \Gamma_g \\
-\bar{\mathbf{q}}^0 \cdot \mathbf{n} = q & \text{on } \Gamma_g, \\
s^0 = 0 & \text{on } \Gamma_u,
\end{array}
\right.
&
\left\{
\begin{array}{ll}
[\bar{\mathbf{u}}^0] = \hat{\mathbf{K}}^{-1} \langle \bar{\sigma}^0 \mathbf{e}_3 \rangle, \\
[\bar{\theta}^0] = -\frac{T_0}{k} \langle \bar{\mathbf{q}}^0 \cdot \mathbf{e}_3 \rangle, \\
[\bar{\sigma}^0 \mathbf{e}_3] = \mathbf{0}, \\
[\bar{\mathbf{q}}^0 \cdot \mathbf{e}_3] = 0.
\end{array}
\right.
\tag{18}
\end{array}
\end{equation}
$$

* Order 1 model

$$
\begin{equation}
\begin{array}{ll}
\text{Governing equations} & \text{Transmission conditions on } S_{\pm} \\
\left\{
\begin{array}{ll}
\rho \ddot{\bar{\mathbf{u}}}^1 - \operatorname{div} \bar{\sigma}^1 = \mathbf{0} & \text{in } \Omega_{\pm}, \\
\dot{\bar{S}}^1 + \frac{1}{T_0} \operatorname{div} \bar{\mathbf{q}}^1 = 0 & \text{in } \Omega_{\pm}, \\
\bar{\sigma}^1 \mathbf{n} = \mathbf{0} & \text{on } \Gamma_g, \\
\bar{\mathbf{q}}^1 \cdot \mathbf{n} = 0 & \text{on } \Gamma_g, \\
s^1 = 0 & \text{on } \Gamma_u,
\end{array}
\right.
&
\left\{
\begin{array}{ll}
[\bar{\mathbf{u}}^1] = \hat{\mathbf{K}}^{-1} \{ \langle \bar{\sigma}^1 \mathbf{e}_3 \rangle - (\hat{\mathbf{K}}^{\alpha})^T \langle \bar{\mathbf{u}}^0 \rangle_{,\alpha} + \hat{\beta} \langle \bar{\theta}^0 \rangle \mathbf{e}_3 \}, & \\
[\bar{\theta}^1] = -\frac{T_0}{k} \langle \bar{\mathbf{q}}^1 \cdot \mathbf{e}_3 \rangle, & \\
[\bar{\sigma}^1 \mathbf{e}_3] = -\mathbf{K}^{\alpha} [\bar{\mathbf{u}}^0]_{,\alpha}, & \\
[\bar{\mathbf{q}}^1 \cdot \mathbf{e}_3] = -\hat{\beta} [\dot{u}_3^0], &
\end{array}
\right.
\tag{19}
\end{array}
\end{equation}
$$

where $\langle f \rangle := \frac{1}{2}(f(\tilde{x}, + (1/2)^+) + f(\tilde{x}, - (1/2)^-))$ and $[f] := f(\tilde{x}, + (1/2)^+) - f(\tilde{x}, - (1/2)^-)$, $\tilde{x} := (x_\alpha) \in S$ denote, respectively, the mean value and the jump functions at the interfaces. The soft thermoelastic interface models at order 0 and order 1 present various similarities, compared with the linear elastic case [27]. At order 0, from a mechanical point of view, the interface behaves as linear springs reacting to the jump between the top and bottom displacements and temperature, while the traction vector and normal heat flow are remains continuous. The order 1 model provides a mixed contact law, expressed by a concurrent discontinuity in terms of thermoelastic state, traction vector, and normal heat flow. The order 1 transmission conditions can be also rewritten in terms of $\langle \bar{\mathbf{q}}^1 \cdot \mathbf{e}_3 \rangle$ and $\langle \bar{\sigma}^1 \mathbf{e}_3 \rangle$, as follows:

$$
\left\{
\begin{array}{l}
\langle \bar{\sigma}^1 \mathbf{e}_3 \rangle = \hat{\mathbf{K}}[\bar{\mathbf{u}}^1] + (\hat{\mathbf{K}}^\alpha)^T \langle \bar{\mathbf{u}}^0 \rangle_{,\alpha} - \hat{\beta} \langle \bar{\theta}^0 \rangle \mathbf{e}_3, \\
\langle \bar{\mathbf{q}}^1 \cdot \mathbf{e}_3 \rangle = -\frac{\hat{k}}{T_0} [\bar{\theta}^1].
\end{array}
\right.
\tag{20}
$$

The jump and mean values of the traction vector and normal heat flow at the interface depend on $s^0$ and are analogous to those obtained for the soft elastic case in [26]. It is interesting to notice that, at order 1, the jump of the heat flow at the interface inside the intermediate layer depend on the variation in time of the normal displacement $u_3$.

4. The Hard Thermoelastic Interface Model

In this section, the limit model for a hard thermoelastic interface, corresponding to an intermediate layer having the same rigidities of the top and bottom bodies, is derived. Let $p = 0$, the asymptotic expansion (16) is inserted in (7), and the following set of variational problems $\mathcal{P}_q$ is obtained:

$$
\left\{
\begin{array}{ll}
\mathcal{P}_{-1}: & a_0(s^0, r) = 0, \\
\mathcal{P}_0: & \bar{A}_-(s^0, r) + \bar{A}_+(s^0, r) + a_0(s^1, r) + a_1(s^0, r) = L(r), \\
\mathcal{P}_1: & \bar{A}_-(s^1, r) + \bar{A}_+(s^1, r) + a_0(s^2, r) + a_1(s^1, r) + a_2(s^0, r) = 0, \\
\mathcal{P}_q: & \bar{A}_-(s^q, r) + \bar{A}_+(s^q, r) + a_0(s^{q+1}, r) + a_1(s^q, r) + a_2(s^{q-1}, r) = 0, \quad q \ge 2
\end{array}
\right.
\tag{21}
$$

A detailed equivalent analysis on the solution of the variational problems $\mathcal{P}_q$ can be found in [33]. In the sequel, the limit problems at order 0 and order 1 are presented.
• Order 0 model

$$
\begin{equation}
\begin{aligned}
& \text{Governing equations} \qquad \text{Transmission conditions on } S_{\pm} \\
& \left\{
  \begin{array}{ll}
    \rho \ddot{\bar{\mathbf{u}}^0} - \operatorname{div} \bar{\sigma}^0 = \mathbf{f} & \text{in } \Omega_{\pm}, \\
    \dot{\bar{S}}^0 + \frac{1}{T_0} \operatorname{div} \bar{\mathbf{q}}^0 = h & \text{in } \Omega_{\pm}, \\
    \bar{\sigma}^0 \mathbf{n} = \mathbf{g} & \text{on } \Gamma_g, \\
    -\bar{\mathbf{q}}^0 \cdot \mathbf{n} = q & \text{on } \Gamma_g, \\
    s^0 = 0 & \text{on } \Gamma_u,
  \end{array}
\right.
\qquad
\left\{
  \begin{array}{l}
    [\bar{\mathbf{u}}^0] = \mathbf{0}, \\
    [\bar{\theta}^0] = 0, \\
    [\bar{\sigma}^0 \mathbf{e}_3] = \mathbf{0}, \\
    [\bar{\mathbf{q}}^0 \cdot \mathbf{e}_3] = 0.
  \end{array}
\right.
\tag{22}
\end{aligned}
\end{equation}
$$

• Order 1 model

$$
\begin{equation}
\begin{aligned}
& \text{Governing equations} \qquad \text{Transmission conditions on } S_{\pm} \\
& \left\{
  \begin{array}{ll}
    \rho \ddot{\bar{\mathbf{u}}}^1 - \operatorname{div} \bar{\sigma}^1 = \mathbf{0} & \text{in } \Omega_{\pm}, \\
    \dot{\bar{S}}^1 + \frac{1}{T_0} \operatorname{div} \bar{\mathbf{q}}^1 = \mathbf{0} & \text{in } \Omega_{\pm}, \\
    \bar{\sigma}^1 \mathbf{n} = \mathbf{0} & \text{on } \Gamma_g, \\
    \bar{\mathbf{q}}^1 \cdot \mathbf{n} = \mathbf{0} & \text{on } \Gamma_g, \\
    s^1 = \mathbf{0} & \text{on } \Gamma_u,
  \end{array}
\right.
\qquad
\left\{
  \begin{aligned}
    [\bar{\mathbf{u}}^1] &= \hat{\mathbf{K}}^{-1}\{\langle\bar{\sigma}^0 \mathbf{e}_3\rangle - (\hat{\mathbf{K}}^\alpha)^T \bar{\mathbf{u}}^0_{,\alpha} + \hat{\beta} \bar{\theta}^0 \mathbf{e}_3\}, \\
    [\bar{\theta}^1] &= -\frac{T_0}{\hat{k}} \langle\bar{\mathbf{q}}^0 \cdot \mathbf{e}_3\rangle, \\
    [\bar{\sigma}^1 \mathbf{e}_3] &= -(\hat{\mathbf{K}}^\alpha \hat{\mathbf{K}}^{-1} \langle\bar{\sigma}^0 \mathbf{e}_3\rangle_{,\alpha} + \hat{\mathbf{L}}^{\alpha\beta} \bar{\mathbf{u}}^0_{,\alpha\beta} \\
    & \qquad + \hat{\beta}(\hat{\mathbf{K}}^\alpha \hat{\mathbf{K}}^{-1} \bar{\theta}^0_{,\alpha} \mathbf{e}_3 - \bar{\theta}^0_{,\tau} \mathbf{e}_\tau)), \\
    [\bar{\mathbf{q}}^1 \cdot \mathbf{e}_3] &= -\left(\frac{\hat{\beta}}{\lambda+2\hat{\mu}} \langle\bar{\theta}^0_{33}\rangle + \hat{\Sigma}^0 - \frac{\hat{k}}{T_0} \Delta_s \bar{\theta}^0\right),
  \end{aligned}
\right.
\tag{23}
\end{aligned}
\end{equation}
$$

where $\hat{\mathbf{L}}^{\alpha\beta} := \hat{\mathbf{K}}^{\alpha\beta} - \hat{\mathbf{K}}^{\beta}\hat{\mathbf{K}}^{-1}(\hat{\mathbf{K}}^{\alpha})^{T}$, $\hat{\Sigma}^{0} := \tilde{\beta}\bar{\mathbf{u}}^{0}_{,\alpha} + \tilde{\mathcal{C}}_{v}\bar{\theta}^{0}$, with $\tilde{\beta} := \frac{2\hat{\mu}\hat{\beta}}{\lambda+2\hat{\mu}}$ and $\tilde{\mathcal{C}}_{v} := \hat{\mathcal{C}}_{v} + \frac{\hat{\beta}^{2}}{\lambda+2\hat{\mu}}$,
$\Delta_s$ denotes the two-dimensional Laplacian operator. Note that, in this case, $\langle\bar{\theta}^{0}\rangle = \bar{\theta}^{0}$ and $\langle\bar{\mathbf{u}}^{0}\rangle = \bar{\mathbf{u}}^{0}$.

The hard thermoelastic interface problems above present the same structures of the
analogous linear elastic hard interface models [24–26]. Concerning the order 0, the trans-
mission conditions provide a continuity of the thermoelastic state and of its conjugated
counterpart, which is typical for adhesives having the same rigidity properties of the
adherents. In this case, the upper and lower bodies are perfectly bonded together. At
order 1, a mixed interface model is obtained, characterized by a jump of the state and
traction vector depending on the values of the thermoelastic state and traction vector at
order 0. These order 0 terms are known since they have been determined in the previous
problem and they appear in the formulation as source terms. The interface conditions at
order 1 can be interpreted as the two-dimensional coupled thermoelastic problem defined
on the plane of the interface.

**5. Generalized Interface Transmission Conditions**

In [26,33], it has been shown that it is possible to obtain a condensed form of trans-
mission conditions summarizing both the orders 0 and 1 of the soft and hard cases in
only one couple of equations in terms of the jump of the displacement field and tractions
at the interface. Equivalently, it is possible to define an implicit general thermoelastic
interface law starting from the hard case, comprising the order 0 and order 1 soft and hard
thermoelastic interface models.

To this end, by denoting by $\bar{s}^\epsilon := \bar{s}^0 + \epsilon\bar{s}^1$, $\bar{\sigma}^\epsilon := \bar{\sigma}^0 + \epsilon\bar{\sigma}^1$ and $\bar{\mathbf{q}}^\epsilon := \bar{\mathbf{q}}^0 + \epsilon\bar{\mathbf{q}}^1$, suitable approximations of $s^\epsilon$, $\bar{\sigma}^\epsilon$ and $\bar{\mathbf{q}}^\epsilon$, respectively, and following the approach developed in [26,33], one can obtain the implicit form of the transmission conditions:

$$
\left\{
\begin{aligned}
\langle \bar{\sigma}^{\epsilon} \mathbf{e}_3 \rangle &= \frac{1}{\epsilon} \hat{\mathbf{K}}[\bar{\mathbf{u}}^{\epsilon}] + (\hat{\mathbf{K}}^{\alpha})^{T} \langle \bar{\mathbf{u}}^{\epsilon} \rangle_{,\alpha} - \hat{\beta} \langle \bar{\theta}^{\epsilon} \rangle \mathbf{e}_3, \\
\langle \bar{\mathbf{q}}^{\epsilon} \cdot \mathbf{e}_3 \rangle &= -\frac{\hat{k}}{\epsilon T_0} [\bar{\theta}^{\epsilon}], \\
[\bar{\sigma}^{\epsilon} \mathbf{e}_3] &= -\hat{\mathbf{K}}^{\alpha}[\bar{\mathbf{u}}^{\epsilon}]_{,\alpha} - \epsilon \hat{\mathbf{K}}^{\alpha\beta} \langle \bar{\mathbf{u}}^{\epsilon} \rangle_{,\alpha\beta} + \epsilon \hat{\beta} \langle \bar{\theta}^{\epsilon} \rangle_{,\alpha} \mathbf{e}_{\alpha}, \\
[\bar{\mathbf{q}}^{\epsilon} \cdot \mathbf{e}_3] &= -\left(\hat{\beta}[\bar{\mathbf{u}}^{\epsilon}]_3 + \epsilon \langle \hat{\Sigma}^{\epsilon} \rangle - \epsilon \frac{\hat{k}}{T_0} \Delta_s \langle \bar{\theta}^{\epsilon} \rangle\right),
\end{aligned}
\right.
\tag{24}
$$
with $\langle \Sigma^\epsilon \rangle := \hat{c}_v \langle \tilde{\theta}^\epsilon \rangle + \hat{\beta} \langle \tilde{u}_\alpha^\epsilon \rangle_\alpha$.

To write the variational formulation of the general coupled thermoelastic interface problem, the expression of the general transmission conditions presented in (24) is employed. In what follows, for the sake of simplicity, the indices $\epsilon$ and symbol (·) are omitted. Let us write the variational form of the equilibrium equations on each sub-domain $\Omega_+$ and $\Omega_-$.

The sum of the two equations leads to

$$ \int_{\Omega_{\pm}} \left\{ \rho \ddot{\mathbf{u}} \cdot \mathbf{v} + \sigma : \mathbf{e}(\mathbf{v}) + \dot{\mathcal{S}} \tilde{\xi} - \frac{1}{T_0} \mathbf{q} \cdot \nabla \tilde{\xi} \right\} dx - \int_S \left\{ \sigma(\tilde{x}, 0^+) \mathbf{n}(\tilde{x}, 0^+) \cdot \mathbf{v} + \sigma(\tilde{x}, 0^-) \mathbf{n}(\tilde{x}, 0^-) \cdot \mathbf{v} \right\} d\Gamma + \int_S \left\{ \mathbf{q}(\tilde{x}, 0^+) \cdot \mathbf{n}(\tilde{x}, 0^+) \tilde{\xi} + \mathbf{q}(\tilde{x}, 0^-) \cdot \mathbf{n}(\tilde{x}, 0^-) \tilde{\xi} \right\} d\Gamma = L(r), \quad (25) $$

which can be written

$$ \int_{\Omega_{\pm}} \left\{ \rho \ddot{\mathbf{u}} \cdot \mathbf{v} + \sigma : \mathbf{e}(\mathbf{v}) + \dot{\mathcal{S}} \tilde{\xi} - \frac{1}{T_0} \mathbf{q} \cdot \nabla \tilde{\xi} \right\} dx + \int_S [\sigma \mathbf{e}_3 \cdot \mathbf{v}] - [\mathbf{q} \cdot \mathbf{e}_3 \tilde{\xi}] d\tilde{x} = L(r), \quad (26) $$

letting $\mathbf{e}_3 = \mathbf{n}(\tilde{x}, 0^-) = -\mathbf{n}(\tilde{x}, 0^+)$ and $d\Gamma = d\tilde{x}$. Then, using the property $[ab] = \langle a \rangle [b] + [a] \langle b \rangle$ and relations (24), and after an integration by parts, one has

$$ \left\{ \begin{array}{l} \text{Find } s \in W(\tilde{\Omega}), \tilde{\Omega} := \Omega_+ \cup S \cup \Omega_-, \ t \in (0, T), \text{ such that} \\ \tilde{A}_-(s,r) + \tilde{A}_+(s,r) + \mathcal{A}(s,r) = \mathcal{L}(r), \end{array} \right. \quad (27) $$

for all $r \in W(\tilde{\Omega})$, where $W(\tilde{\Omega}) := W(\tilde{\Omega}) \times W(\tilde{\Omega})$, with $W(\tilde{\Omega}) := \{r \in H^1(\tilde{\Omega}), r|_S \in H^1(S), r=0 \text{ on } \Gamma_u\}, W(\tilde{\Omega}) := [W(\tilde{\Omega})]^3$, and

$$ \mathcal{A}(s,r) := \int_S \left\{ \frac{1}{\epsilon} \hat{\mathbf{K}}[\mathbf{u}] \cdot [\mathbf{v}] + (\hat{\mathbf{K}}^\alpha)^T \langle \mathbf{u} \rangle_\alpha \cdot [\mathbf{v}] + \hat{\mathbf{K}}^\alpha[\mathbf{u}] \cdot \langle \mathbf{v} \rangle_\alpha + \epsilon \hat{\mathbf{K}}^{\alpha\beta} \langle \mathbf{u} \rangle_\beta \cdot \langle \mathbf{v} \rangle_\alpha - \right. \\ \left. - \hat{\beta} \langle \theta \rangle [v_3] - \epsilon \hat{\beta} \langle \theta \rangle \langle v_\alpha \rangle_\alpha + \frac{1}{\epsilon} \frac{\hat{k}}{T_0} [\theta][\tilde{\xi}] + \epsilon \frac{\hat{k}}{T_0} \langle \theta \rangle_\alpha \langle \tilde{\xi} \rangle_\alpha + \hat{\beta}[i\mathbf{u}_3] \langle \tilde{\xi} \rangle + \epsilon \langle \dot{\Sigma} \rangle \langle \tilde{\xi} \rangle \right\} d\tilde{x}, \quad (28) $$

$$ \mathcal{L}(r) := \int_{\Omega_{\pm}} \{\mathbf{f} \cdot \mathbf{v} + h\tilde{\xi}\}dx + \int_{\Gamma_g} \{\mathbf{g} \cdot \mathbf{v} + q\tilde{\xi}\}d\Gamma + \int_{\partial S} \{\mathbf{F} \cdot \langle \mathbf{v} \rangle + \mathcal{H}(\tilde{\xi})\}d\gamma, \quad (29) $$

where $\langle \Sigma \rangle := \hat{c}_v \langle \theta \rangle + \hat{\beta} \langle u_\alpha \rangle_\alpha$, $(\hat{\mathbf{K}}^\alpha[\mathbf{u}] + \epsilon \hat{\mathbf{K}}^{\alpha\beta} \langle \mathbf{u} \rangle_\beta - \epsilon \hat{\beta} \langle \theta \rangle \mathbf{e}_\alpha) v_\alpha := \mathbf{F}$ and $\epsilon \frac{\hat{k}}{T_0} \langle \theta \rangle_\alpha v_\alpha := \mathcal{H}$ denote the loads on the lateral boundary of the interface $\partial S$, with outer unit normal vector $(v_\alpha)$ (see [33]).

## 6. FEM Implementation

The numerical simulations were carried out by means of the finite element method, discretizing the variational problem (27). This helped validate the proposed asymptotic approach. The FEM analysis was performed considering the coupled dynamic thermoelastic problem and comparing the solution of the three-phase model (two adherents and adhesive) with the generalized interface (two adherents + interface). The problem was solved employing the software GetFem++ (see [41,42] for more details), with a standard linear solver (conjugate gradient). For that purpose, standard piecewise linear finite elements were considered.

Let us consider a thermoelastic laminated plate occupying a 3D domain defined by $\Omega^\epsilon = [0, L_1] \times [0, L_2] \times [0, 2h + \epsilon]$, with $h = 1$ cm, $L_1/h = 10$, $L_2/h = 5$. (see Figure 2). Clearly, with self-explanatory notation, $x_1 = x$, $x_2 = y$, and $x_3 = z$. The adherents are made of Material 1, while the adhesive is constituted by Material 2.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The 3D geometry of the thermoelastic laminated plate represented in the plane $(x_1, x_3)$.</figcaption></figure>

Simply supported boundary conditions are considered on the bottom edges of the composite plate. The plate is subject to a thermal shock $q \cdot n|_{\Gamma_{top}} = q(t) = ate^{-bt}$ on the top face, with $a = 30$ and $b = 0.8$, whose plot is illustrated in Figure 3. The bottom face is thermally insulated $q \cdot n|_{\Gamma_{bottom}} = 0$. No volume or surface mechanical loads were applied $f = g = 0$.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Applied heat flow shock.</figcaption></figure>

The FEM discretization was carried out using piecewise linear finite elements on hexahedrons, with 7280 nodes (29,203 degrees of freedom) for the three-phase problem and 5824 nodes (23,635 degrees of freedom) for the problem with the generalized interface law. The time discretization was realized using a Newmark-beta scheme with $\beta = 0.25$ and $\gamma = 0.5$ for the second time derivative of the displacements, and a Crank-Nicholson scheme for the first derivative of the thermodynamic entropy. It is worth noting that both methods are unconditionally stable and of order two in time. The time step of discretization is equal to $\delta t = 0.1$ s.

The numerical example considers a composite plate, in which the adherents and the adhesive have very different thermo-mechanical properties. Material 1 is aluminum (Al), while Material 2 is a polyvinyl chloride (PVC) foam. The constitutive parameters are listed in Table 1.

<table><caption>**Table 1.** Thermoelastic material properties for Al and PVC Foam.</caption><thead><tr><th colspan="3">Material 1: Al</th><th colspan="3">Material 2: PVC Foam</th></tr></thead><tbody><tr><td>ρ<sub>1</sub></td><td>2700</td><td>[kg/m<sup>3</sup>]</td><td>ρ<sub>2</sub></td><td>250</td><td>[kg/m<sup>3</sup>]</td></tr><tr><td>E<sub>1</sub></td><td>72.4</td><td>[GPa]</td><td>E<sub>2</sub></td><td>0.28</td><td>[GPa]</td></tr><tr><td>ν<sub>1</sub></td><td>0.32</td><td></td><td>ν<sub>2</sub></td><td>0.40</td><td></td></tr><tr><td>α<sub>1</sub></td><td>40.0</td><td>[µm/m K]</td><td>α<sub>2</sub></td><td>22.4</td><td>[µm/m K]</td></tr><tr><td>k<sub>1</sub></td><td>122.2</td><td>[W/mK]</td><td>k<sub>2</sub></td><td>0.05</td><td>[W/mK]</td></tr><tr><td>c<sup>1</sup><sub>v</sub></td><td>900</td><td>[J/kg K]</td><td>c<sup>2</sup><sub>v</sub></td><td>1900</td><td>[J/kg K]</td></tr></tbody></table>
To evaluate the accuracy of the asymptotic analysis, the influence of the relative thickness $\epsilon/h$, for fixed time instants, on the $L^2$-relative error was investigated. The $L^2$-relative errors $\frac{\|u^\epsilon - u\|_{L^2}}{\|u\|_{L^2}}$ and $\frac{\|\theta^\epsilon - \theta\|_{L^2}}{\|\theta\|_{L^2}}$ was computed taking into account the solution $(u^\epsilon, \theta^\epsilon)$ of the initial three-phase problem, discretized with a FE mesh, and the solution $(u, \theta)$ of the interface problem (27). Tables 2 and 3 report the relative error values for increasing time and vanishing relative thickness.

<table><caption>**Table 2.** Relative error $\frac{\|u^\epsilon - u\|_{L^2}}{\|u\|_{L^2}}$.</caption><thead><tr><th>ε/h / t</th><th>5</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td>0.1</td><td>4.20 × 10<sup>-3</sup></td><td>1.89 × 10<sup>-3</sup></td><td>1.58 × 10<sup>-3</sup></td><td>1.60 × 10<sup>-3</sup></td></tr><tr><td>0.05</td><td>7.55 × 10<sup>-4</sup></td><td>3.44 × 10<sup>-4</sup></td><td>2.76 × 10<sup>-4</sup></td><td>3.45 × 10<sup>-4</sup></td></tr><tr><td>0.01</td><td>2.25 × 10<sup>-5</sup></td><td>3.48 × 10<sup>-6</sup></td><td>4.52 × 10<sup>-6</sup></td><td>4.56 × 10<sup>-6</sup></td></tr></tbody></table>

<table><caption>**Table 3.** Relative error $\frac{\|\theta^\epsilon - \theta\|_{L^2}}{\|\theta\|_{L^2}}$.</caption><thead><tr><th>ε/h / t</th><th>5</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td>0.1</td><td>1.56 × 10<sup>-3</sup></td><td>5.16 × 10<sup>-4</sup></td><td>7.57 × 10<sup>-4</sup></td><td>4.97 × 10<sup>-4</sup></td></tr><tr><td>0.05</td><td>1.07 × 10<sup>-4</sup></td><td>2.83 × 10<sup>-4</sup></td><td>9.94 × 10<sup>-5</sup></td><td>2.27 × 10<sup>-5</sup></td></tr><tr><td>0.01</td><td>9.88 × 10<sup>-8</sup></td><td>5.58 × 10<sup>-7</sup></td><td>1.78 × 10<sup>-8</sup></td><td>7.43 × 10<sup>-10</sup></td></tr></tbody></table>

The convergence diagrams of the relative $L^2$-norms of the displacements and temperatures, obtained with the three-phase problem and the reduced interface problem, are plotted in Figure 4, as the thickness ratio $\epsilon/h$ tends to zero, at time $t = 1$ s and $t = 10$ s, respectively. Moreover, the evolution in time of the $L^2$-relative error is reported in Figure 5.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Convergence diagrams with respect to $\epsilon/h$ for: $t = 1$ s (**left**); and $t = 10$ s (**right**).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5. Evolution of the relative error with respect to the time (ε/h = 0.01).</figcaption></figure>

From the results in Tables 2 and 3 and, especially, Figure 4, it can be noticed that, by decreasing the thickness ratio ε/h, the relative errors present an immediate reduction for fixed times. The convergence rate is of the order (ε/h)² and remains constant for increasing time instants. As illustrated in Figure 5, the evolution in time of the L²-relative error, for fixed ε/h, becomes approximately steady after t = 10 s for the displacement field, while it presents a decreasing trend concerning the temperature field. Besides, even for a relative thickness ε/h = 0.1, at time t = 10 s, the relative error is close to about 1.89 × 10⁻³, for the displacement field, and about 5.16 × 10⁻⁴, for the temperature field. Hence, the proposed general thermoelastic interface model provides an acceptable solution and it is able to correctly approximate the solution (u^e, θ^e) of the three-phase problem. Moreover, the reduced model can also be employed for moderately thick adhesives.

In the sequel, the numerical results obtained by solving the general interface model are presented, considering a relative thickness of ε/h = 0.01. Following the approach by [43], hereinafter, the results are provided using dimensionless variables:

* $U(X_1, X_2, X_3, t) := \frac{1-\nu_1}{\ell(1+\nu_1)\alpha_1 T_0} u(x_1, x_2, x_3, t)$,

* $\Theta(X_1, X_2, X_3, t) := \frac{\theta(x_1, x_2, x_3, t) - T_0}{T_0}$,

* $\Sigma_{ij}(X_1, X_2, X_3, t) := \frac{1}{\rho_1 V} \sigma_{ij}(x_1, x_2, x_3, t)$

where $X_i = x_i/\ell$, $t = V/\ell$, and $\ell$ and $V$ are defined by

$$V = \sqrt{\frac{E_1(1-\nu_1)}{(1+\nu_1)(1-2\nu_1)\rho_1}}, \quad \ell = \frac{k_1}{\rho_1 c_v^1 V} \quad (30)$$

Let us notice that the domain Ω^e is chosen such that $X_1 \in [0; 10]$, $X_2 \in [0; 10]$ and $X_3 \in [0; 1]$.

Figure 6 represents the trend of the displacement $U_3$ and temperature $\Theta$, evaluated along $X_3$ on the orthogonal fiber to the mid-plane of the interface at point $(\bar{X}_1 = 6, \bar{X}_2 = 6)$, for given times. The plot shows that, after the thermal shock, the displacement $U_3$ evolves in opposite directions within the adherents: the composite laminated plate tends to expand and contract itself along the through-the-thickness axis. On the other hand, the temperature field $\Theta$ remains constant along the $X_3$-axis within the adherents, for given times, reaching a steady value after a certain time interval. As expected, the plots also report a jump of the state fields ($U_3, \Theta$) in correspondence of the intermediate layer, and, thus, the adhesive behaves as a soft thermoelastic interface. This is mainly due to the material properties of the adhesive, which are smaller with respect to those of the adherents.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Displacement $U_3$ and temperature $\Theta$ along the $X_3$-axis, on a fiber $(\bar{X}_1 = 6, \bar{X}_2 = 6, X_3)$, for given times.</figcaption></figure>

Figure 7 illustrates the evolution of the displacement field **U** and temperature **Θ** with respect to the time *t*, at a given point X̄ = (6.5, 3.5, 0), placed on the bottom face of the composite plate. As expected, the thermal shock induces an oscillatory trend concerning the displacements. Conversely, the temperature evolves to a steady state, corresponding to a constant value, after a sudden increase related to the thermal shock application.

Figure 8 represents the trend of the stresses $\Sigma_{33}$ and $\Sigma_{13}$, evaluated along $X_3$ on the orthogonal fiber to the mid-plane of the interface at point $(\bar{X}_1 = 6, \bar{X}_2 = 6)$, for given times. The plot shows that, after the thermal shock, the stress $\Sigma_{33}$ remains constant along the $X_3$-axis within the adherents. In this particular case, the thermal contribution to $\Sigma_{33}$ is predominant with respect to the elastic one, i.e., $\Sigma_{33} \approx -B(\Theta + 1)$: indeed, their diagrams present analogous trends and differ for a constant of proportionality $B := \frac{\beta T_0}{\rho_1 V}$ (see Figures 6 and 8). The stress $\Sigma_{13}$ presents an oscillating behavior along $X_3$ inside the adherents, but its contribution is negligible compared with $\Sigma_{33}$. Moreover, the normal $\Sigma_{33}$ and shear stresses $\Sigma_{13}$, evaluated at the top and bottom faces of the intermediate layer, are very similar and, thus, their jump almost vanishes. This is typical of soft interface models, in which the thermoelastic state presents a discontinuity, while its conjugated quantities (traction vector and normal heat flow) are continuous across the interface.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Displacement field $U = (U_i)$ and temperature $Θ$ versus time $t$, at a given point $\bar{X} = (6.5, 3.5, 0)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Stresses $\Sigma_{33}$ and $\Sigma_{13}$ along the $X_3$-axis, on a fiber ($\bar{X}_1 = 6, \bar{X}_2 = 6, X_3$), for given times.</figcaption></figure>
Figure 9 shows a comparison between the evolution in time of $U_3(\bar{X}, t)$, at a given point $\bar{X} = (6.5, 3.5, 0.)$, of a homogeneous three-layer plate, made of aluminum, and the Al/PVC composite plate. Although the thickness ratio is small ($\varepsilon/h = 0.01$), the effect of the adhesive becomes relevant concerning the response of the plate to the thermal shock. The homogeneous plate appears to be stiffer with respect to the composite one, which manifests a significant amplitude and period increase of the $U_3$ motion.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Evolution in time of $U_3(\bar{X}, t)$ at a given point $\bar{X} = (6.5, 3.5, 0.)$ for a homogeneous Al-plate and an Al/PVC composite plate.</figcaption></figure>

**7. Concluding Remarks**

General imperfect interface conditions are proposed in the framework of coupled thermoelasticity, simulating the thermomechanical behavior of a thin-bonded joint. The approach used to obtain the transmission conditions is based on the asymptotic expansions method. Zero- and higher-order interface models are derived for soft and hard interphases. Following [33], a general transmission law, comprising the two regimes (soft and hard) at the various order, is derived. To assess the validity of the previous asymptotic approach, numerical simulations were developed using a finite element method, which generalizes an analogous methodology to dynamical coupled thermoelasticity, already proposed in [34] in the framework of piezoelectricity. The numerical example consisted of a thermoelastic composite three-layer aluminum plate, with a PVC adhesive, subject to a thermal shock. Two different configurations were considered: the first one consisted of an initial three-phase problem, while the second one took into account the FE discretized form of interface problem (27). The most significant fields (displacement and temperature) and their $L^2$-relative errors were then computed and compared to test the validity of the proposed interface laws and the accuracy of the asymptotic model. The proposed general thermoelastic interface model provides an acceptable solution and it is able to correctly approximate the solution of the three-phase problem. These findings clearly indicate that the approach of substituting the interphase with the proposed interface law provides a robust modeling for the composite.

**Author Contributions:** The following statements should be used "Conceptualization, M.S., S.D., R.R. and F.L.; methodology, M.S., S.D., R.R. and F.L.; software, S.D.; validation, M.S. and S.D.; formal analysis, M.S., S.D. and R.R.; investigation, M.S.; data curation, M.S., S.D. and R.R.; writing—original draft preparation, M.S.; writing—review and editing, M.S., S.D., R.R. and F.L.; funding acquisition, R.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
## References

1. Benveniste Y. The effective conductivity of composites with imperfect thermal contact at constituent interfaces. *Int. J. Eng. Sci* **1986**, *24*,1537–1552. [CrossRef]

2. Benveniste Y. Effective thermal-conductivity of composites with a thermal contact resistance between the constituents-nondilute case. *J. Appl. Phys.* **1987**, *61*, 2840–2843. [CrossRef]

3. Kapitza P.L. *Collected Papers of P. L. Kapitza*; ter Haar, D., Ed.; reprinted 1965; Pergamon: Oxford, UK, 1941.

4. Pham Huy H.; Sanchez-Palencia, E. Phenomenes de transmission à travers des couches minces de conductivité élevée. *J. Math. Anal. Appl.* **1974**, *47*, 284–309.

5. Miloh, T.; Benveniste, Y. On the effective conductivity of composites with ellipsoidal inhomogeneities and highly conducting interfaces. *Proc. R. Soc. Lond. A* **1999**, *455*, 2687–2706. [CrossRef]

6. Hashin, Z. Thin interphase/imperfect interface in conduction. *J. App. Phys.* **2001**, *89*, 2261–2267. [CrossRef]

7. Javili, A.; Kaessmair, S.; Steinmann, P. General imperfect interfaces. *Comput. Methods Appl. Mech. Engrg.* **2014**, *275*, 76–97. [CrossRef]

8. Benveniste, Y. Effective mechanical behaviour of composite materials with imperfect contact between the constituents. *Mech. Mat.* **1985**, *4*, 197–208. [CrossRef]

9. Hashin, Z. Thin interphase/imperfect interface in elasticity with application to coated fiber composites. *J. Mech. Phys. Solids* **2002**, *50*, 2509–2537. [CrossRef]

10. Gurtin, M.E.; Murdoch, A.I. A continuum theory of elastic material surfaces. *Arch. Ration. Mech. Anal.* **1975**, *57*, 291–323. [CrossRef]

11. Yvonnet, J.; Le Quang, H.; He, Q.C. An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites. *Comput. Mech.* **2008**, *42*, 119–131. [CrossRef]

12. Yvonnet, J.; Mitrushchenkov, A.; Chambaud, G.; He, Q.C. Finite element model of ionic nanowires with size-dependent mechanical properties determined by ab initio calculations. *Comput. Methods Appl. Mech. Eng.* **2011**, *200*, 614–625. [CrossRef]

13. Benveniste, Y.; Miloh, T. Imperfect soft and stiff interfaces in two-dimensional elasticity. *Mech. Mat.* **2001**, *33* 309–323. [CrossRef]

14. Benveniste, Y. A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media. *J. Mech. Phys. Solids* **2006**, *54*, 708–734. [CrossRef]

15. Ciarlet, P.G. *Mathematical Elasticity*, Vol. II: *Theory of Plates*; North-Holland: Amsterdam, The Netherlands, 1997.

16. Serpilli, M.; Lenci, S. Asymptotic modelling of the linear dynamics of laminated beams. *Int. J. Solids Struct.* **2012**, *49*, 1147–1157. [CrossRef]

17. Serpilli, M.; Lenci, S. An overview of different asymptotic models for anisotropic three-layer plates with soft adhesive. *Int. J. Solids Struct.* **2016**, *81* 130–140. [CrossRef]

18. Furtsev, A.; Itou, H.; Rudoy, E. Modeling of bonded elastic structures by a variational method: Theoretical analysis and numerical simulation. *Int. J. Solids Struct.* **2020**, *182–183*, 100–111. [CrossRef]

19. Furtsev, A.; Rudoy, E. Variational approach to modeling soft and stiff interfaces in the Kirchhoff-Love theory of plates. *Int. J. Solids Struct.* **2020**, *202*, 562–574. [CrossRef]

20. Klarbring, A. Derivation of the adhesively bonded joints by the asymptotic expansion method. *Int. J. Eng. Sci.* **1991**, *29*, 493–512. [CrossRef]

21. Geymonat, G.; Krasucki, F.; Lenci, S. Mathematical Analysis of a bonded joint with a soft thin adhesive. *Math. Mech. Solids* **1999**, *16*, 201–225. [CrossRef]

22. Bessoud, A.L.; Krasucki, F.; Serpilli, M. Asymptotic analysis of shell-like inclusions with high rigidity. *J. Elast.* **2011**, *103*, 153–172. [CrossRef]

23. Ljulj, M.; Tambaca, J. 3D structure – 2D plate interaction model. *Math. Mech. Solids* **2019**, *4*, 3354–3377. [CrossRef]

24. Lebon, F.; Rizzoni, R. Asymptotic analysis of a thin interface: The case involving similar rigidity. *Int. J. Eng. Sci.* **2010**, *48*, 473–486. [CrossRef]

25. Lebon, F.; Rizzoni, R. Asymptotic behavior of a hard thin linear interphase: An energy approach. *Int. J. Solids Struct.* **2011**, *48*, 441–449. [CrossRef]

26. Rizzoni, R.; Dumont, S.; Lebon, F.; Sacco, E. Higher order model for soft and hard elastic interfaces. *Int. J. Solids Struct.* **2014**, *51*, 4137–4148. [CrossRef]

27. Dumont, S.; Rizzoni, R.; Lebon, F.; Sacco, E. Soft and hard interface models for bonded elements. *Compos. Part B Eng.* **2018**, *153*, 480–490. [CrossRef]

28. Lebon, F.; Rizzoni, R. Higher order interfacial effects for elastic waves in one dimensional phononic crystals via the Lagrange-Hamilton’s principle. *Eur. J. Mech. A Solids* **2018**, *67*, 57–70. [CrossRef]

29. Serpilli, M. Mathematical modeling of weak and strong piezoelectric interfaces. *J. Elast.* **2015**, *121*, 235–254. [CrossRef]

30. Serpilli, M. Asymptotic interface models in magneto-electro-thermo-elastic composites. *Meccanica* **2017**, *52*, 1407–1424. [CrossRef]

31. Serpilli, M. Classical and higher order interface conditions in poroelasticity. *Ann. Solid Struct. Mech.* **2019**, *11*, 1–10. [CrossRef]

32. Serpilli, M. On modeling interfaces in linear micropolar composites. *Math. Mech. Solids* **2018**, *23*, 667–685. [CrossRef]

33. Serpilli, M.; Rizzoni, R.; Lebon, F.; Dumont, S. An asymptotic derivation of a general imperfect interface law for linear multiphysics composites. *Int. J. Solids Struct* **2019**, *180–181*, 97–107. [CrossRef]
34. Dumont, S.; Serpilli, M.; Rizzoni, R.; Lebon, F. Numerical Validation of Multiphysic Imperfect Interfaces Models. *Front. Mater.* **2020**, *7*, 158. [CrossRef]

35. Shi, P.; Shillor, M. Existence of a solution to the N dimensional problem of thermoelastic contact. *Indiana Univ. Math. J.* **1992**, *17*, 1597–1618.

36. Khludnev, M.; Kovtunenko, V.A. *Analysis of Cracks in Solids*; WIT-Press: Southampton, UK; Boston, MA, USA, 2000.

37. Xu, X. The N-dimensional quasistatic problem of thermoelastic contact with Barbers heat exchange conditions. *Adv. Math. Sci. Appl.* **1996**, *6*, 559–587.

38. Showalter, R.E. Diffusion in poro-elastic media. *J. Math. Anal. Appl.* **2000**, *251*, 310–340. [CrossRef]

39. Miara, B.; Suarez, J.S. Asymptotic pyroelectricity and pyroelasticity in thermopiezoelectric plates. *Asympt. Anal.* **2013** *81*, 211–250.

40. Bonaldi, F.; Geymonat, G.; Krasucki, F.; Serpilli, M. An asymptotic plate model for magneto-electro-thermo-elastic sensors and actuators. *Math. Mech. Solids* **2017**, *22*, 798–822. [CrossRef]

41. Renard, Y.; Pommier, J. *Getfem++*. *An Open Source Generic C++ Library for Finite Element Methods*; Tech. Rep.; INSA Lyon: Villeurbanne, France, 2002.

42. Geuzaine, C.; Remacle, J.-F. Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities. *Int. J. Numer. Methods Eng.* **2009**, *79*, 1309–1331. [CrossRef]

43. Jafarinezhad, M.R.; Eslami, M.R. Coupled thermoelasticity of FGM annular plate under lateral thermal shock. *Compos. Struct.* **2017**, *168*, 758–771. [CrossRef]