# Numerical Validation of Multiphysic Imperfect Interfaces Models

Serge Dumont, Michèle Serpilli, Raffaella Rizzoni, Frédéric Lebon

► To cite this version:

Serge Dumont, Michèle Serpilli, Raffaella Rizzoni, Frédéric Lebon. Numerical Validation of Multi-physic Imperfect Interfaces Models. Frontiers in Materials. Computational Materials Science section, 2020, 7, pp.158. 10.3389/fmats.2020.00158. [hal-02909277](https://hal.science/hal-02909277)

HAL Id: hal-02909277

https://hal.science/hal-02909277

Submitted on 2 Jun 2021

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
<figure><img src="image_2.png" /></figure>

Check for
updates

# Numerical Validation of Multiphysic Imperfect Interfaces Models

Serge Dumont¹*, Michele Serpilli², Raffaella Rizzoni³ and Frédéric C. Lebon⁴

¹IMAG CNRS UMR 5149, University of Nimes, Nimes, France, ²Department of Civil and Building Engineering and Architecture, Università Politecnica delle Marche, Ancona, Italy, ³Department of Engineering, University of Ferrara, Ferrara, Italy, ⁴Aix-Marseille Univ, CNRS, Centrale Marseille, LMA, Marseille, France

We investigate some mathematical and numerical methods based on asymptotic expansions for the modeling of bonding interfaces in the presence of linear coupled multiphysic phenomena. After reviewing new recently proposed imperfect contact conditions (Serpilli et al., 2019), we present some numerical examples designed to show the efficiency of the proposed methodology. The examples are framed within two different multiphysic theories, piezoelectricity and thermo-mechanical coupling. The numerical investigations are based on a finite element approach generalizing to multiphysic problems the procedure developed in Dumont et al. (2018).

**Keywords:** asymptotic analysis, interfaces, multiphysic materials, numerical analysis, finite element analysis

OPEN ACCESS

**Edited by:**
Chiu Tai Law,
University of Wisconsin-Milwaukee,
United States

**Reviewed by:**
Xiaowei Zeng,
University of Texas at San Antonio,
United States
Alberto Salvadori,
University of Brescia, Italy

*Correspondence:
Serge Dumont
serge.dumont@unimes.fr

**Specialty section:**
This article was submitted to
Mechanics of Materials,
a section of the journal
Frontiers in Materials

**Received:** 06 December 2019
**Accepted:** 04 May 2020
**Published:** 05 June 2020

**Citation:**
Dumont S, Serpilli M, Rizzoni R and
Lebon FC (2020) Numerical Validation
of Multiphysic Imperfect Interfaces
Models. Front. Mater. 7:158.
doi: 10.3389/fmats.2020.00158

## 1. INTRODUCTION

The design of composite materials and structures is nowadays characterized by an increasing level of complexity regarding their heterogeneous morphology and constitution in order to attain more reliable and multi-functional performances in extreme environments. Even though advances have been made in the understanding of the mechanical behavior of composite materials, their synergistic response in a multi-physical environment is not yet fully understood. Composites as well as bonded structures are obtained by joining different parts with highly contrasted mechanical properties to compose a unique assembly. The bonded joint may generally be imperfect and discontinuities in the involved physical fields can arise, drastically changing the overall mechanical response (see, e.g., Gu and He, 2011; Gu et al., 2014; Javili et al., 2014). Therefore, a rigorous modeling of the imperfect bonding plays an important role in engineering design.

In the present paper, we analyze a particular composite constituted by a three bodies: namely, two adherents separated by a thin adhesive interphase. The composite assembly domain depends on a small parameter $\epsilon$, corresponding to the constant thickness of the intermediate layer, called the adhesive. In order to take into account general multiphysic interactions among various physical behaviors, such as elasticity, magneto-electricity, thermal conduction, the constituents are made of general linear multiphysic materials.

Classically, the thin adhesive can be replaced by an interface law, as its thickness $\epsilon$ tends to zero. In the limit, the two adherents are joined by a limit surface, across which non-classical contact conditions are given to reproduce the behavior of the thin adhesive. Because the numerical treatment of a thin adhesive layers needs a very fine discretization, a 3D-2D model based on the introduction of the interface law is advantageous, requiring a much lower computation cost.

In the last years, various interface models have been proposed based on phenomenological or micromechanical approaches. On one hand, phenomenological models can be deduced directly from experiments. For example, effective non-linear interface phenomenological models have been proposed, among the others, in Fouchal et al. (2009), Lotfi and Shing (1994), and Lourenço and Rots (1997) for the analysis of masonry walls. On the other hand, interface models obtained
through the asymptotic methods (see Ciarlet, 1997 for an extensive discussion) can be regarded as micromechanical models linking the interface response to the behavior of the material layer constituting the bonding, cf. Klarbring, 1991; Benveniste and Miloh, 2001; Hashin, 2002. The asymptotic analysis is a powerful mathematical technique, employed to justify classical thin structures and layered composites model (Ciarlet, 1997; Geymonat et al., 1999, 2014; Serpilli and Lenci, 2012, 2016; Serpilli, 2017). According to these techniques, the thickness of the glue is considered as a small parameter ε and, possibly, its stiffness depends on ε^p, where p is a critical exponent. Letting ε tend to zero, the adhesive layer geometrically vanishes, reducing itself to a 2D surface, but it is accounted for by a relation linking the jumps of the traction and the displacement vectors at the interface. Another benefit of a simplified 3D-2D model based on an interface law is the possibility of incorporating multiphysic materials behavior, spanning from uncoupled phenomena, such as thermal conduction and elasticity (Lebon and Rizzoni, 2010, 2011; Bessoud et al., 2011; Rizzoni et al., 2014; Raffa et al., 2018), to multifield and multiphysic theories, such as piezoelectricity and continua with microstructure, cf. Serpilli, 2015, 2017, 2018, 2019; Serpilli et al., 2019 and the references herein.

The aim of this paper is to numerically validate the general imperfect contact conditions derived in Serpilli et al. (2019), simulating the behavior of a thin interphase undergoing linear coupled multiphysic phenomena. As already stressed in Serpilli et al. (2019), the proposed multiphysic interface law, obtained through the asymptotic methods, provides a generalized transmission condition which contains in itself the soft (also called Kapitza's or lowly-conducting), hard (perfect), and rigid (Gurtin-Murdoch's or highly-conducting) interface models. Our analysis is limited to the case of linear elastic response of the adhesive and adherents and the modeling is framed in the contest of small strains and small displacements theory. Then, some examples are given in order to show the efficiency of the proposed methodology in the case of piezoelectric and thermo-mechanical couplings. The numerical investigations are performed in the framework of the finite element method, generalizing the approach developed in Dumont et al. (2018) to multiphysic problems. In particular, finite elements are introduced to solve the fully 3D initial problem, considering a three-phase composite, and the limit 3D-2D problem with two layers with imperfect interface transmission conditions.

Finally, the current model is limited to the case of a multiphysics linear elastic response of the adhesive and the adherents. However, using the techniques proposed in Bonetti et al. (2017), Raffa et al. (2018), and Rizzoni et al. (2017), it could be extended to more general frameworks, like the non-linear regime or damaging adhesives with behavior governed by the evolution of the crack length at the micro-scale. These extensions will be the topics of further study.

## 2. THE GENERAL THEORY

In this section, some results previously obtained are recalled. We consider the assembly constituted of two solids $\Omega_{\pm}^{\epsilon} \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$, $0 < \epsilon < 1$, called the adhesive, with cross-section $S \subset \mathbb{R}^2$. In the following $B^{\epsilon}$ and $S$ will be 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.

We suppose that the composite is constituted of a material which presents a linear coupled multiphysic behavior. The state at the equilibrium of the multiphysic material is characterized by a collection of order parameters, using the multifield theory jargon (see, e.g., Mariano, 2002): $N$ vector state variables, namely $\mathbf{u}_{1}^{\epsilon}, \dots, \mathbf{u}_{N}^{\epsilon}$, and $M$ scalar state variables, namely $\phi_{1}^{\epsilon}, \dots, \phi_{M}^{\epsilon}$. Let us collect all the unknowns into a generalized vector field $\mathbf{s}^{\epsilon} := (\mathbf{u}_{1}^{\epsilon}, \dots, \mathbf{u}_{N}^{\epsilon}, \phi_{1}^{\epsilon}, \dots, \phi_{M}^{\epsilon})$, the so-called *multiphysic state*. With the multiphysic state $\mathbf{s}^{\epsilon}$, we associate its conjugated physical quantity $\mathbf{t}^{\epsilon} = \mathbf{t}^{\epsilon}(\nabla^{\epsilon}\mathbf{s}^{\epsilon})$, depending on the gradient of $\mathbf{s}^{\epsilon}$, defined by $(\nabla^{\epsilon}\mathbf{s}^{\epsilon})_{i} := \mathbf{s}_{i}^{\epsilon} = (\mathbf{u}_{i,1}^{\epsilon}, \dots, \mathbf{u}_{i,N}^{\epsilon}, \phi_{1}^{\epsilon}, \dots, \phi_{M,i}^{\epsilon})$. The vector field $\mathbf{t}^{\epsilon} := (\sigma_{1}^{\epsilon}, \dots, \sigma_{N}^{\epsilon}, \mathbf{D}_{1}^{\epsilon}, \dots, \mathbf{D}_{M}^{\epsilon})$ represents a generalized stress field.

We consider the following homogeneous and linear constitutive law:

$$ \mathbf{t}^{\epsilon} = \mathbb{K}^{\epsilon} \nabla^{\epsilon} \mathbf{s}^{\epsilon}, $$

where $\mathbb{K}^{\epsilon}$ is a generalized linear constitutive matrix, that we can write, component-wise,

$$ \left\{ \begin{array}{l} \sigma_{J}^{\epsilon} = \mathbf{C}_{JK}^{\epsilon} \nabla^{\epsilon} \mathbf{u}_{K}^{\epsilon} + \mathbf{P}_{JL}^{\epsilon} \nabla^{\epsilon} \varphi_{I}^{\epsilon}, \\ \mathbf{D}_{L}^{\epsilon} = \mathbf{R}_{LK}^{\epsilon} \nabla^{\epsilon} \mathbf{u}_{K}^{\epsilon} + \mathbf{H}_{LI}^{\epsilon} \nabla^{\epsilon} \varphi_{I}^{\epsilon}, \end{array} \right. $$

with $J, K = 1, \dots, N$ and $I, L = 1, \dots, M$ and using Einstein's notations. The constitutive tensor $\mathbb{K}^{\epsilon}$ satisfies the classical symmetry and positivity properties.

A straight-forward example of such behavior is represented by piezoelectricity, in which an electric field is generated by a mechanical strain (direct effect) and, viceversa, a mechanical strain is produced as a result of an electric field (converse effect). The piezoeletric state consists of a pair $\mathbf{s} = (\mathbf{u}, \varphi)$, namely the displacement field $\mathbf{u} = (u_i)$ and the electric potential $\varphi$. The constitutive law takes the form takes the form:

$$ \left\{ \begin{array}{l} \sigma = \mathbf{Ce}(\mathbf{u}) - \mathbf{PE}, \\ \mathbf{D} = \mathbf{P}^{T}\mathbf{e}(\mathbf{u}) + \mathbf{HE}, \end{array} \right. \qquad (1) $$

where $\mathbf{e}(\mathbf{u}) := \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ is the linearized strain tensor, $\mathbf{D} = (D_i)$ represent the electric displacement field, $\mathbf{E} := -\nabla\varphi$ represent the electric field, while $\mathbf{C} = (C_{ijkl})$, $\mathbf{P} = (P_{ijk})$ and $\mathbf{H} = (H_{ij})$ denote, respectively, the elasticity tensor, the piezoelectric coupling tensor and the dielectric tensor. In this case, tensor $\mathbb{K}$ reduces to

$$ \mathbb{K} = \begin{bmatrix} \mathbf{C} & \mathbf{P} \\ \mathbf{P}^T & -\mathbf{H} \end{bmatrix}. $$

Any other multiphysic behavior can be generalized by adding other possible couplings within the constitutive law, coming from different physical interactions, as in magneto-electro-thermo-elastic materials (see, e.g., Serpilli, 2017, 2018).
It is assumed that the adherents are subject to a generalized system of body forces $F^\epsilon: \Omega_\pm^\epsilon \to \mathbb{R}^{3N \times M}$ and surface forces $G^\epsilon: \Gamma_g^\epsilon \to \mathbb{R}^{3N \times M}$, where $\Gamma_g^\epsilon \subset (\partial\Omega_+^\epsilon \setminus S_+^\epsilon) \cup (\partial\Omega_-^\epsilon \setminus S_-^\epsilon)$. Body and surface forces are neglected in adhesive layer. On $\Gamma_u^\epsilon := (\partial\Omega_+^\epsilon \setminus S_+^\epsilon) \cup (\partial\Omega_-^\epsilon \setminus S_-^\epsilon) \setminus \Gamma_g^\epsilon$ homogeneous boundary conditions are prescribed, so that $\mathbf{s}^\epsilon = \mathbf{0}$ on $\Gamma_u^\epsilon$. We assume that everywhere, near the interphase boundary $\Gamma_{lat} := \partial S \times (-\frac{\epsilon}{2}, \frac{\epsilon}{2})$, homogeneous Neumann boundary conditions are applied. It is assumed that the adhesive and the adherents are perfectly bonded in order to ensure the continuity of the multiphysic state and generalized stress vector fields across $S_\pm^\epsilon$. The differential (strong) formulation of the governing equations has the following structure:

$$ \left\{ \begin{array}{ll} -\mathrm{div} \, \mathbf{t}^{\epsilon} = \mathbf{F}^{\epsilon} & \text{in } \Omega^{\epsilon}, \\ \mathbf{t}^{\epsilon} \mathbf{n}^{\epsilon} = \mathbf{G}^{\epsilon} & \text{on } \Gamma_{g}^{\epsilon}, \\ \mathbf{s}^{\epsilon} = \mathbf{0} & \text{on } \Gamma_{u}^{\epsilon}, \end{array} \right. \qquad (2) $$

where $\mathbf{t}^\epsilon \mathbf{n}^\epsilon := (\sigma_1^\epsilon \mathbf{n}^\epsilon, \dots, \sigma_N^\epsilon \mathbf{n}^\epsilon, \mathbf{D}_1^\epsilon \cdot \mathbf{n}^\epsilon, \dots, \mathbf{D}_M^\epsilon \cdot \mathbf{n}^\epsilon)$ represents the generalized traction vector on the boundary $\Gamma_g^\epsilon$ and $\mathbf{n}^\epsilon$ the outer normal unit vector to $\Gamma_g^\epsilon$. We take $V(\Omega^\epsilon)$ to denote the space of state variables admissible fields. The weak formulation of problem (2) defined on the variable domain $\Omega^\epsilon$ is:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{s}^\epsilon \in V(\Omega^\epsilon) \text{ such that} \\ \bar{A}_-^\epsilon(\mathbf{s}^\epsilon, \mathbf{r}^\epsilon) + \bar{A}_+^\epsilon(\mathbf{s}^\epsilon, \mathbf{r}^\epsilon) + \hat{A}^\epsilon(\mathbf{s}^\epsilon, \mathbf{r}^\epsilon) = L^\epsilon(\mathbf{r}^\epsilon), \end{array} \right. \qquad (3) $$

for all $\mathbf{r}^\epsilon := (\mathbf{v}_1^\epsilon, \dots, \mathbf{v}_N^\epsilon, \psi_1^\epsilon, \dots, \psi_M^\epsilon) \in V(\Omega^\epsilon)$, where the bilinear forms $\bar{A}_\pm^\epsilon(\cdot, \cdot)$, and $\hat{A}^\epsilon(\cdot, \cdot)$ and the linear form $L^\epsilon(\cdot)$ are defined by

$$ \begin{aligned} \bar{A}_\pm^\epsilon(\mathbf{s}^\epsilon, \mathbf{r}^\epsilon) &:= \int_{\Omega_\pm^\epsilon} \bar{\mathbf{K}}^\epsilon \nabla^\epsilon \mathbf{s}^\epsilon \cdot \nabla^\epsilon \mathbf{r}^\epsilon dx^\epsilon, \\ \hat{A}^\epsilon(\mathbf{s}^\epsilon, \mathbf{r}^\epsilon) &:= \int_{B^\epsilon} \hat{\mathbf{K}}^\epsilon \nabla^\epsilon \mathbf{s}^\epsilon \cdot \nabla^\epsilon \mathbf{r}^\epsilon dx^\epsilon, \\ L^\epsilon(\mathbf{r}^\epsilon) &:= \int_{\Omega_\pm^\epsilon} \mathbf{F}^\epsilon \cdot \mathbf{r}^\epsilon dx^\epsilon + \int_{\Gamma_g^\epsilon} \mathbf{G}^\epsilon \cdot \mathbf{r}^\epsilon d\Gamma^\epsilon. \end{aligned} $$

In the sequel, the following notations are used:

$$ \begin{aligned} \langle f \rangle(\tilde{x}) &:= \frac{1}{2}(f(\tilde{x}, (1/2)^+) + f(\tilde{x}, -(1/2)^-)), \quad \tilde{x} := (x_\alpha) \in S, \\ [f](\tilde{x}) &:= f(\tilde{x}, (1/2)^+) - f(\tilde{x}, -(1/2)^-), \\ \langle \langle f \rangle \rangle(\tilde{x}) &:= \frac{1}{2}(f(\tilde{x}, 0^+) + f(\tilde{x}, 0^-)), \\ [f \rangle](\tilde{x}) &:= f(\tilde{x}, 0^+) - f(\tilde{x}, 0^-). \end{aligned} $$

The asymptotic behavior of problem (3), i.e., when $\epsilon \to 0$ is studied. The first step of the approach given by Ciarlet (1997) is to introduce a rescaling 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} $$

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_g$ and $\Gamma_u$ indicate the images through $\pi^\epsilon$ of $\Gamma_g^\epsilon$ and $\Gamma_u^\epsilon$. 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{1}{\epsilon} \frac{\partial}{\partial x_3}$ in $B$. It is assumed that the constitutive coefficients of $\Omega_\pm^\epsilon$ are independent of $\epsilon$, so that $\bar{\mathbf{K}}^\epsilon = \bar{\mathbf{K}}$, while the constitutive coefficients of $B^\epsilon$ present the following dependences on $\epsilon$, so that $\hat{\mathbf{K}}^\epsilon = \epsilon^p \bar{\mathbf{K}}$, with $p \in \{-1, 0, 1\}$. Three different limit behaviors will be characterized according to the choice of the exponent $p$: in the case of $p = -1$, we derive a model for a rigid interface; when $p = 0$, we derive a model for a hard interface; by choosing $p = 1$, we deduce a model for a soft interface. By virtue of classical hypothesis, the rescaled problem can be written in the following form:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{s}^\epsilon \in V(\Omega), \text{ such that} \\ \bar{A}_-(\mathbf{s}^\epsilon, \mathbf{r}) + \bar{A}_+(\mathbf{s}^\epsilon, \mathbf{r}) + \epsilon^{p-1}\hat{a}(\mathbf{s}^\epsilon, \mathbf{r}) + \epsilon^p \hat{b}(\mathbf{s}^\epsilon, \mathbf{r}) + \epsilon^{p+1}\hat{c}(\mathbf{s}^\epsilon, \mathbf{r}) = L(\mathbf{r}), \end{array} \right. \qquad (4) $$

for all $\mathbf{r} \in V(\Omega)$, the set of rescaled admissible fields, where the bilinear forms $\bar{A}_\pm(\cdot, \cdot)$, $\hat{a}(\cdot, \cdot)$, $\hat{b}(\cdot, \cdot)$ and $\hat{c}(\cdot, \cdot)$ are defined by

$$ \begin{aligned} \bar{A}_\pm(\mathbf{s}^\epsilon, \mathbf{r}) &:= \int_{\Omega_\pm} \bar{\mathbf{K}} \nabla^\epsilon \mathbf{s}^\epsilon \cdot \nabla^\epsilon \mathbf{r} dx^\epsilon, \\ \hat{a}(\mathbf{s}^\epsilon, \mathbf{r}) &:= \int_B \hat{\mathbf{K}}_{33} \mathbf{s}_3^\epsilon \cdot \mathbf{r}_{,3} dx^\epsilon, \\ \hat{b}(\mathbf{s}^\epsilon, \mathbf{r}) &:= \int_B \{ \hat{\mathbf{K}}_{3\alpha} \mathbf{s}_3^\epsilon \cdot \mathbf{r}_{,\alpha} + \hat{\mathbf{K}}_{\alpha3} \mathbf{s}_{,\alpha}^\epsilon \cdot \mathbf{r}_{,3} \} dx^\epsilon, \\ \hat{c}(\mathbf{s}^\epsilon, \mathbf{r}) &:= \int_B \hat{\mathbf{K}}_{\alpha\beta} \mathbf{s}_{,\beta}^\epsilon \cdot \mathbf{r}_{,\alpha} dx^\epsilon, \end{aligned} $$

and $\hat{\mathbf{K}}_{ij}$ denote the sub-matrices of $\hat{\mathbf{K}}$, defined by

$$ \hat{\mathbf{K}} = \begin{bmatrix} \hat{\mathbf{K}}_{\alpha\beta} & \hat{\mathbf{K}}_{\alpha3} \\ \hat{\mathbf{K}}_{3\alpha} & \hat{\mathbf{K}}_{33} \end{bmatrix}, \quad (\hat{\mathbf{K}}_{ij})^T = \hat{\mathbf{K}}_{ji}. $$

At this level, one can perform an asymptotic analysis of the rescaled problem (4), introducing the asymptotic expansions of the solution as a series of powers of $\epsilon$:

$$ \begin{aligned} \mathbf{s}^\epsilon &= \mathbf{s}^0 + \epsilon \mathbf{s}^1 + \epsilon^2 \mathbf{s}^2 + \dots, \\ \bar{\mathbf{s}}^\epsilon &= \bar{\mathbf{s}}^0 + \epsilon \bar{\mathbf{s}}^1 + \epsilon^2 \bar{\mathbf{s}}^2 + \dots, \\ \hat{\mathbf{s}}^\epsilon &= \hat{\mathbf{s}}^0 + \epsilon \hat{\mathbf{s}}^1 + \epsilon^2 \hat{\mathbf{s}}^2 + \dots \end{aligned} \qquad (5) $$

where $\bar{\mathbf{s}}^\epsilon = \mathbf{s}^\epsilon \circ \bar{\pi}^\epsilon$ and $\hat{\mathbf{s}}^\epsilon = \mathbf{s}^\epsilon \circ \hat{\pi}^\epsilon$, by substituting (5) into the rescaled problem (4), and by identifying the terms with identical power of $\epsilon$, we can recover the limit problems at order 0 and order 1. Finally, matching conditions are introduced based on the continuity of the generalized traction $\mathbf{t}^\epsilon e_3$ and multiphysic state $\mathbf{s}^\epsilon$ at the interfaces $S_\pm^\epsilon$ in the initial configuration and on the continuity of the traction and state $\bar{\mathbf{t}}^\epsilon e_3$, $\bar{\mathbf{s}}^\epsilon$, $\hat{\mathbf{t}}^\epsilon e_3$, $\hat{\mathbf{s}}^\epsilon$ at the interfaces $S_\pm$ in the rescaled configuration (see Rizzoni et al., 2014; Dumont et al., 2018).

## 2.1. The Soft Multiphysic Interface Model

In the limit case $p = 1$, the adhesive is weaker than the adherents and the transmission problems at order 0 and order 1 can be summarized as follows (Serpilli et al., 2019):
* Order 0

**Governing equations**

$$ egin{cases} -\operatorname{div} \bar{\mathbf{t}}^0 = \mathbf{F} & \text{in } \Omega_{\pm}, \\ \bar{\mathbf{t}}^0 \mathbf{n} = \mathbf{G} & \text{on } \Gamma_g, \\ \bar{\mathbf{s}}^0 = \mathbf{0} & \text{on } \Gamma_u, \end{cases} $$

**Transmission conditions on $S_\pm$**

$$ egin{cases} [\bar{s}^0] = (\hat{\mathbf{K}}_{33})^{-1} \langle \bar{\mathbf{t}}^0 \mathbf{e}_3 \rangle, \\ [\bar{\mathbf{t}}^0 \mathbf{e}_3] = \mathbf{0}. \end{cases} $$

* Order 1

**Governing equations**

$$ egin{cases} -\operatorname{div} \bar{\mathbf{t}}^1 = \mathbf{0} & \text{in } \Omega_{\pm}, \\ \bar{\mathbf{t}}^1 \mathbf{n} = \mathbf{0} & \text{on } \Gamma_g, \\ \bar{\mathbf{s}}^1 = \mathbf{0} & \text{on } \Gamma_u, \end{cases} $$

**Transmission conditions on $S_\pm$**

$$ egin{cases} [\bar{s}^1] = (\hat{\mathbf{K}}_{33})^{-1} (\langle \bar{\mathbf{t}}^1 \mathbf{e}_3 \rangle - \hat{\mathbf{K}}_{\alpha 3} \langle \bar{s}^0 \rangle_{,\alpha}), \\ [\bar{\mathbf{t}}^1 \mathbf{e}_3] = - \hat{\mathbf{K}}_{3\alpha} [\bar{s}^0]_{,\alpha}. \end{cases} $$

## 2.2. The Hard Multiphysic Interface Model

When $p = 0$, the limit model corresponds to a hard multiphysic interface. The order 0 and order 1 interface transmission problems are presented in the sequel (Serpilli et al., 2019):

* Order 0

**Governing equations**

$$ egin{cases} -\operatorname{div} \bar{\mathbf{t}}^0 = \mathbf{F} & \text{in } \Omega_{\pm}, \\ \bar{\mathbf{t}}^0 \mathbf{n} = \mathbf{G} & \text{on } \Gamma_g, \\ \bar{\mathbf{s}}^0 = \mathbf{0} & \text{on } \Gamma_u, \end{cases} $$

**Transmission conditions on $S_\pm$**

$$ egin{cases} [\bar{s}^0] = \mathbf{0}, \\ [\bar{\mathbf{t}}^0 \mathbf{e}_3] = \mathbf{0}. \end{cases} $$

* Order 1

**Governing equations**

$$ egin{cases} -\operatorname{div} \bar{\mathbf{t}}^1 = \mathbf{0} & \text{in } \Omega_{\pm}, \\ \bar{\mathbf{t}}^1 \mathbf{n} = \mathbf{0} & \text{on } \Gamma_g, \\ \bar{\mathbf{s}}^1 = \mathbf{0} & \text{on } \Gamma_u, \end{cases} $$

**Transmission conditions on $S_\pm$**

$$ egin{cases} [\bar{s}^1] = (\hat{\mathbf{K}}_{33})^{-1} (\langle \bar{\mathbf{t}}^0 \mathbf{e}_3 \rangle - \hat{\mathbf{K}}_{\alpha 3} \langle \bar{s}^0 \rangle_{,\alpha}), \\ [\bar{\mathbf{t}}^1 \mathbf{e}_3] = - (\hat{\mathbf{K}}_{3\alpha} [\bar{s}^1]_{,\alpha} + \hat{\mathbf{K}}_{\alpha\beta} \langle \bar{s}^0 \rangle_{,\alpha\beta}). \end{cases} $$

## 2.3. The Rigid Multiphysic Interface Model

A rigid multiphysic interface can be obtained for $p = -1$. The order 0 and order 1 limit interface models are the following ones (Serpilli et al., 2019):

* Order 0

**Governing equations**

$$ egin{cases} -\operatorname{div} \bar{\mathbf{t}}^0 = \mathbf{F} & \text{in } \Omega_{\pm}, \\ \bar{\mathbf{t}}^0 \mathbf{n} = \mathbf{G} & \text{on } \Gamma_g, \\ \bar{\mathbf{s}}^0 = \mathbf{0} & \text{on } \Gamma_u, \end{cases} $$

**Transmission conditions on $S_\pm$**

$$ egin{cases} [\bar{s}^0] = \mathbf{0}, \\ [\bar{\mathbf{t}}^0 \mathbf{e}_3] = - \hat{\mathcal{L}}_{\alpha\beta} \langle \bar{s}^0 \rangle_{,\alpha\beta}. \end{cases} $$

* Order 1

**Governing equations**

$$ egin{cases} -\operatorname{div} \bar{\mathbf{t}}^1 = \mathbf{0} & \text{in } \Omega_{\pm}, \\ \bar{\mathbf{t}}^1 \mathbf{n} = \mathbf{0} & \text{on } \Gamma_g, \\ \bar{\mathbf{s}}^1 = \mathbf{0} & \text{on } \Gamma_u, \end{cases} $$

**Transmission conditions on $S_\pm$**

$$ egin{cases} [\bar{s}^1] = -(\hat{\mathbf{K}}_{33})^{-1} \hat{\mathbf{K}}_{\alpha 3} \langle \bar{s}^0 \rangle_{,\alpha}, \\ [\bar{\mathbf{t}}^1 \mathbf{e}_3] = - \hat{\mathbf{K}}_{3\alpha} (\hat{\mathbf{K}}_{33})^{-1} \langle \bar{\mathbf{t}}^0 \mathbf{e}_3 \rangle_{,\alpha} - \hat{\mathcal{L}}_{\alpha\beta} \langle \bar{s}^1 \rangle_{,\alpha\beta}. \end{cases} $$

Note that it is possible to obtain a condensed form of transmission 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 (Rizzoni et al., 2014). To this end, we denote by $\bar{s}^\epsilon := \bar{s}^0 + \epsilon \bar{s}^1 + \epsilon^2 \bar{s}^2$ and $\bar{\mathbf{t}}^\epsilon := \bar{\mathbf{t}}^0 + \epsilon \bar{\mathbf{t}}^1$, two suitable approximations for $\bar{s}^\epsilon$ and $\bar{\mathbf{t}}^\epsilon$. The rigid multiphysic interface conditions are considered. Let $\hat{\mathbf{K}} = \epsilon \hat{\mathbf{K}}^\epsilon$ in $B^\epsilon$, we can write $[\bar{s}^\epsilon]$ and $[\bar{\mathbf{t}}^\epsilon \mathbf{e}_3]$, as

$$ egin{cases} [\bar{s}^\epsilon] = -\epsilon (\hat{\mathbf{K}}_{33})^{-1} (\hat{\mathbf{K}}_{\alpha 3}^\epsilon \langle \bar{s}^\epsilon \rangle_{,\alpha} - (\bar{\mathbf{t}}^\epsilon \mathbf{e}_3)) + o(\epsilon^2), \\ [\bar{\mathbf{t}}^\epsilon \mathbf{e}_3] = -\epsilon \hat{\mathbf{K}}_{3\alpha}^\epsilon (\hat{\mathbf{K}}_{33})^{-1} (\bar{\mathbf{t}}^\epsilon \mathbf{e}_3)_{,\alpha} - \epsilon \hat{\mathcal{L}}_{\alpha\beta}^\epsilon \langle \bar{s}^\epsilon \rangle_{,\alpha\beta} + o(\epsilon^2). \end{cases} $$

An alternative expression of the above transmission conditions can be given in terms of $\langle \bar{\mathbf{t}}^\epsilon \mathbf{e}_3 \rangle$ and $[\bar{\mathbf{t}}^\epsilon \mathbf{e}_3]$, which will be useful to write the variational formulation of the interface multiphysic problem:

$$ egin{cases} \langle \bar{\mathbf{t}}^\epsilon \mathbf{e}_3 \rangle = \frac{1}{\epsilon} \hat{\mathbf{K}}_{33}^\epsilon [\bar{s}^\epsilon] + \hat{\mathbf{K}}_{\alpha 3}^\epsilon \langle \bar{s}^\epsilon \rangle_{,\alpha} + o(\epsilon^2), \\ [\bar{\mathbf{t}}^\epsilon \mathbf{e}_3] = - \hat{\mathbf{K}}_{3\alpha}^\epsilon [\bar{s}^\epsilon]_{,\alpha} - \epsilon \hat{\mathbf{K}}_{\alpha\beta}^\epsilon \langle \bar{s}^\epsilon \rangle_{,\alpha\beta} + o(\epsilon^2). \end{cases} $$

The variational formulation of the general multiphysic interface problem is given by Serpilli et al. (2019):

$$ egin{cases} \text{Find } \mathbf{s} \in W(\tilde{\Omega}), \text{ such that} \\ \tilde{A}_-(\mathbf{s}, \mathbf{r}) + \tilde{A}_+(\mathbf{s}, \mathbf{r}) + \mathcal{A}(\mathbf{s}, \mathbf{r}) = \mathcal{L}(\mathbf{r}), \end{cases} \quad (6) $$

for all $\mathbf{r} \in W(\tilde{\Omega})$, where $W(\tilde{\Omega})$ is the set of admissible fields, with $\tilde{\Omega} := \Omega_+ \cup S \cup \Omega_-$ and

$$ egin{aligned} \mathcal{A}(\mathbf{s}, \mathbf{r}) := & \int_S \left( \frac{1}{\epsilon} \hat{\mathbf{K}}_{33} [[\mathbf{s}]] \cdot [[\mathbf{r}]] + \hat{\mathbf{K}}_{\alpha 3} \langle \langle \mathbf{s} \rangle \rangle_{,\alpha} \cdot [[\mathbf{r}]] \right. \\ & \left. + \hat{\mathbf{K}}_{3\alpha} [[\mathbf{s}]] \cdot \langle \langle \mathbf{r} \rangle \rangle_{,\alpha} + \epsilon \hat{\mathbf{K}}_{\alpha\beta} \langle \langle \mathbf{s} \rangle \rangle_{,\alpha} \cdot \langle \langle \mathbf{r} \rangle \rangle_{,\beta} \right) d\tilde{x}, \\ \mathcal{L}(\mathbf{r}) := & \int_{\partial \Omega_{\pm}} \mathbf{F} \cdot \mathbf{r} dx + \int_{\Gamma_g} \mathbf{G} \cdot \mathbf{r} d\Gamma + \int_{\partial S} \mathbf{f} \cdot \langle \langle \mathbf{r} \rangle \rangle d\gamma, \end{aligned} $$

where $\mathbf{f} := (\hat{\mathbf{K}}_{3\alpha} [[\mathbf{s}]] + \epsilon \hat{\mathbf{K}}_{\alpha\beta} \langle \langle \mathbf{s} \rangle \rangle_{,\beta}) v_\alpha$ denotes the load on the lateral boundary of the interface, which can be evaluated. By virtue of the Lax-Milgram lemma, we can infer that the interface variational problem (6) admits one and only one solution (see Serpilli et al., 2019, for more details).

## 3. NUMERICAL VALIDATION OF THE GENERALIZED INTERFACE CONDITIONS

This Section is devoted to the numerical implementation of the generalized interface conditions, identified in Serpilli et al. (2019) and expressed by means of the transmission problem (6). The numerical simulations will help to prove the numerical performance of the proposed asymptotic approach. The above problem is numerically solved through the finite element method by considering a multiphysic composite constituted by two parallelepiped adherents separated by a thin adhesive, subject to general loading conditions. The analysis is carried out by choosing two particular multiphysic constitutive laws,
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | The 3D geometry of the piezoelectric laminated plate represented in the plane $(x_1, x_3)$.</figcaption></figure>

corresponding to the case of piezoelectricity and thermo-elasticity. Finally, the numerical solution of the three-phase model (two adherents and adhesive) is compared with the numerical solution provided by the generalized interface problem (6), obtained by means of the asymptotic methods. These simple examples can be employed to model and understand the mechanical behavior of real laminated structures presenting multiphysic couplings such as piezoelectric stack actuators, piezo-patch sensor glued on a controlled substructure (see, for instance, Geis et al., 2004).

The numerical simulations hereinafter are performed using a finite element method to solve the weak problem (6). For that purpose, standard piecewise linear finite element are considered. As in the description of the interface model, the multiphysic variables are considered using a generalized variable $s$, leading to a fully coupled problem.

Finally, the problem is solved employing the software GetFem++ (see Renard and Pommier, 2002; Geuzaine and Remacle, 2009 for more details), with a standard linear solver (conjugate gradient).

## 3.1. The Piezoelectric Case

A piezoelectric material represents one of the most peculiar multiphysic material, combining the linear elastic behavior with the electric counterpart. As mentioned in section 2, the piezoeletric state reduces to a pair $s = (u, \varphi)$, constituted by the displacement field $u = (u_i)$ and the electric potential $\varphi$. The constitutive law takes the form (1).

In what follows, we define the geometrical and physical settings of the numerical simulations. Let us consider a piezoelectric laminated plate occupying a 3D domain defined by $\Omega = [0, L_1] \times [0, L_2] \times [-h/2, h/2]$, with $L_1 = 10h$ and $L_2 = 5h$ (see **Figure 1**).

The adherents (Piezoelectric 1) are constituted by PVDF (Polyvinylidene fluoride), a monoclinic piezoelectric material with poling axis $e_3$, while the adhesive (Piezoelectric 2) is made of PZT-4, a transversally isotropic piezoelectric material with poling axis $e_3$, whose mechanical properties are shown in **Table 1** (Fernades and Pouget, 2002). The constitutive sub-matrices ($K_{ij}$) are defined as follows:

<table><caption>**TABLE 1** | Piezoelectric material properties for PZT-4 and PVDF, Fernades and Pouget (2002).</caption><thead><tr><th>Moduli</th><th>PZT-4</th><th>PVDF</th></tr></thead><tbody><tr><td>c<sub>11</sub>, GPa</td><td>139</td><td>238.24</td></tr><tr><td>c<sub>22</sub>, GPa</td><td>139</td><td>23.6</td></tr><tr><td>c<sub>33</sub>, GPa</td><td>115</td><td>10.64</td></tr><tr><td>c<sub>12</sub>, GPa</td><td>77.8</td><td>3.98</td></tr><tr><td>c<sub>13</sub>, GPa</td><td>74.3</td><td>2.19</td></tr><tr><td>c<sub>23</sub>, GPa</td><td>74.3</td><td>1.92</td></tr><tr><td>2c<sub>44</sub>, GPa</td><td>25.6</td><td>2.15</td></tr><tr><td>2c<sub>55</sub>, GPa</td><td>25.6</td><td>4.4</td></tr><tr><td>2c<sub>66</sub>, GPa</td><td>30.6</td><td>6.43</td></tr><tr><td>e<sub>31</sub>, C/m<sup>2</sup></td><td>-5.2</td><td>-0.13</td></tr><tr><td>e<sub>32</sub>, C/m<sup>2</sup></td><td>-5.2</td><td>-0.145</td></tr><tr><td>e<sub>33</sub>, C/m<sup>2</sup></td><td>15.1</td><td>-0.276</td></tr><tr><td>e<sub>24</sub>, C/m<sup>2</sup></td><td>12.7</td><td>-0.009</td></tr><tr><td>e<sub>15</sub>, C/m<sup>2</sup></td><td>12.7</td><td>-0.135</td></tr><tr><td>H<sub>11</sub>, nF/m</td><td>13.06</td><td>0.111</td></tr><tr><td>H<sub>22</sub>, nF/m</td><td>13.06</td><td>0.106</td></tr><tr><td>H<sub>33</sub>, nF/m</td><td>11.51</td><td>0.106</td></tr></tbody></table>

$$K_{33} = \begin{pmatrix} 2c_{55} & 0 & 0 & 0 \\ 0 & 2c_{44} & 0 & 0 \\ 0 & 0 & c_{33} & e_{33} \\ 0 & 0 & -e_{33} & -H_{33} \end{pmatrix}, K_{12} = \begin{pmatrix} 0 & 2c_{66} & 0 & 0 \\ c_{12} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix},$$

$$K_{13} = \begin{pmatrix} 0 & 0 & 2c_{55} & e_{15} \\ 0 & 0 & 0 & 0 \\ c_{13} & 0 & 0 & 0 \\ -e_{31} & 0 & 0 & 0 \end{pmatrix}, K_{23} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 2c_{44} & e_{24} \\ 0 & c_{23} & 0 & 0 \\ 0 & -e_{32} & 0 & 0 \end{pmatrix},$$

$$K_{11} = \begin{pmatrix} c_{11} & 0 & 0 & 0 \\ 0 & 2c_{66} & 0 & 0 \\ 0 & 0 & 2c_{55} & e_{15} \\ 0 & 0 & -e_{15} & -H_{11} \end{pmatrix}, K_{22} = \begin{pmatrix} 2c_{66} & 0 & 0 & 0 \\ 0 & c_{22} & 0 & 0 \\ 0 & 0 & 2c_{44} & e_{24} \\ 0 & 0 & -e_{24} & -H_{22} \end{pmatrix}.$$

For a transversally isotropic material with poling axis $e_3$, one has $c_{11} = c_{22}$, $c_{13} = c_{23}$, $c_{55} = c_{44}$, $c_{66} = (c_{11} - c_{12})/2$, $e_{15} = e_{24}$, $e_{31} = e_{32}$ and $H_{11} = H_{22}$.

### 3.1.1. Numerical Results: The Applied Electric Potentials

The FEM discretization is carried out using P1 finite elements (linear on thetrahedrons), with 6552 nodes (27187 degrees of freedom) for the three-phases problem and 5824 nodes (24275 degrees of freedom) for the problem with the generalized interface law.

In this section, we study the case of two electric potentials $V_0 = \pm 25V$, applied on both the upper and the lower faces of the composite plate (see **Figure 1**). No mechanical loading is considered in this example. The piezoelectric laminated plates behaves as a piezoelectric actuator.

Following the ideas proposed in Fernades and Pouget (2002), the numerical results for the variables are provided using the
dimensionless units. For an applied electric potential $V_0$, we set:

$$ (U_i, \Phi_i) = \frac{E_0}{V_0}(u_i, \frac{\varphi_i}{E_0}) \qquad (T_{ij}, D_k) = \frac{hE_0}{C_{00}V_0}(\sigma_{ij}, E_0 D_k) $$

where, for numerical convenience, $E_0 = 10^9 \text{Vm}^{-1}$ and $C_{00} = 1 \text{GPa}$.

First, the influence of the relative thickness of the interphase $\epsilon_h$ is investigated in order to evaluate the accuracy of the asymptotic modeling. In particular, the quality of the solutions is evaluated considering the $L^2$-relative error $\frac{\|s^\epsilon - s_{model}\|}{\|s^\epsilon\|}$, where $s^\epsilon$ denotes the reference solution computed using the three-phases problem with a sufficiently fine finite element mesh to ensure a good

approximation of the solution, while $s_{model}$ indicates the solution of the interface model (6). Let us notice that the mesh does not depend on the thickness $\epsilon$ for the interface model, on the contrary of the computation of the reference solution, where the meshing of the interphase is necessary. The convergence of the general interface model toward the three-phases one with respect to the thickness ratio $\epsilon_h$ is presented in Figure 2.

From the plot, it can be observed that, by reducing the thickness of the adhesive, the relative error has a drastic reduction and so, the proposed general interface model provides an acceptable solution and it is able to correctly approximate the solution $s^\epsilon$ (see Bonnet et al., 2016). The convergence rate is $(\epsilon/h)^3$ and is consistent with expectancy. Besides, even if the

<figure><img src="image_4.png" /><figcaption>FIGURE 2 | Convergence results with respect to the relative thickness $\epsilon/h$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3 | Displacement field and electric potential on a section along the $x_3$ axis, with $\epsilon/h = 0.1$.</figcaption></figure>
relative thickness is of 1%, the relative error is close to $7.05 \cdot 10^{-4}$ for the displacement field and close to $2.42 \cdot 10^{-5}$ for the electric potential, meaning that the general interface model can also be used for moderately thick adhesive layers.

In **Figure 2**, we also compared the solutions of the generalized interface conditions with some other classical interface models, namely, the hard and soft interfaces. Firstly, we consider the model, referenced as hard model, where the displacements and the electrical potential fields, and the traction vector and normal electric displacement at the interface are considered as continuous. Hence, the hard interface can be considered perfect. This case is obtained with asymptotic expansions method at order 0, when the materials properties in the interphase are independent of the parameter $\epsilon$, and thus they have similar rigidities with respect to the adherents (see the equations at order 0 in section 2.2). In addition, we present the results obtained by the so-called soft model, obtained with asymptotic expansions at order 0 when the interphase material properties depend linearly on the parameter $\epsilon$. In the soft case, it can be shown that the interface behaves as a linear multiphysic spring, where only the first term of the bilinear form $A(s, r) = \int_S \frac{1}{\epsilon} \hat{\mathbf{K}}_{33} [[s]] \cdot [[r]] dx$

is considered in the interface model (6) (see, also, equations at order 0 in section 2.1). **Figure 2** shows that the present modeling significantly improves the quality of the approximation, even for large $\epsilon$: the trend of the $L^2$-relative error convergence rate, both for the displacements and the electric potential field, is passing from $\epsilon/h$ to $(\epsilon/h)^3$. This means that the solution of the three-phases model converges more quickly toward the solution of the present reduced model, as $\epsilon$ tends to zero, then the solution of the classical hard or soft interface models. Moreover, it is important to stress that the current generalized multiphysic interface model has been built in order to adapt itself to the any kind of adhesive, since it comprises and also improve both the hard and soft cases, taking into account higher order terms, see Serpilli et al., 2019 for a detailed discussion.

We now present a comparison between the solutions obtained on the three-phases problem and on the asymptotic

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Jumps of the displacement field and the electric potential across the interface on a line along the $x_1$ axis, with $\epsilon/h = 0.1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | Jumps of the displacement field and the electric potential across the interface on a line along the $x_2$ axis, with $\epsilon/h = 0.1$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | Jumps of the traction vector and normal electric displacement across the interface on a line along the $x_1$ axis, with $\epsilon/h = 0.1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 7** | Jumps of the traction vector and normal electric displacement across the interface on a line along the $x_2$ axis, with $\epsilon/h = 0.1$.</figcaption></figure>
approximations. In this case, the relative thickness of the interface is equal to 0.1, and the relative error is equal to 0.04% in terms of displacements and is equal to 10⁻³% in terms of electric potential. In the sequel, the plots of the solutions ($U_i, \Phi_i$) are represented on the orthogonal line to the interface mid-plane defined by ($x_1 = L_1/2, x_2 = L_2/2, x_3/h \in [-0.5, 0.5]$), while the jumps of the various physical quantities are plotted on the lines parallel to

the interface mid-plane, namely ($x_1 \in [0, L_1]$, $x_2 = L_2/2$, $x_3 = 0$) and ($x_1 = L_1/2$, $x_2 \in [0, L_2]$, $x_3 = 0$).

**Figure 3** represents the trend of the displacement field and electric potential, evaluated on the central orthogonal fiber to the mid-plane of the interface. The plot shows a very good agreement between the solution of the general interface problem (dotted line) and the solution of the three-phases problem (solid line). The composite plate behaves mostly as a Kirchhoff-Love single-layer plate, taking also into account the transversal deformation of the adhesive. From the electric point of view, the electric

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

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

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 11 | Jumps of the displacement field and the electric potential across the interface on a line along the x₂ axis, with ε/h = 0.1.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12 | Jumps of the traction vector and normal electric displacement across the interface on a line along the x₁ axis, with ε/h = 0.1.</figcaption></figure>

potential is linear through the adherents and constant through the adhesive: this is due to the fact that the intermediate layer (PZT-4) has a higher electrical conductivity with respect to upper and lower bodies (PVDF), see Table 1, and, hence, it behaves as a highly conducting interface.

Figures 4, 6 represent the trend of the jumps of the displacement and electric potential and the jumps of the traction vector and normal electric displacement along the x₁-axis, namely (x₁ ∈ [0, L₁], x₂ = L₂/2, x₃ = 0), and Figures 5, 7 along the x₂-axis, namely (x₁ = L₁/2, x₂ = ∈ [0, L₂], x₃ = 0). The numerical simulations highlight that the proposed model is able to describe the mechanical behavior of the composite. Few solution oscillations can be found close to the lateral boundaries, due to the presence of edges, which produce expected stress concentrations and singularities.

### 3.1.2. Numerical Results: The Applied Surface Loads

In this section, we study the case of a constant surface load equal to p = 1Pa, applied on the top face of the plate, with no electrical

<figure><img src="image_15.png" /><figcaption>FIGURE 13 | Jumps of the traction vector and normal electric displacement across the interface on a line along the x₂ axis, with ε/h = 0.1.</figcaption></figure>

loadings (see Figure 1). The composite laminated plate behaves as a sensor. As previously shown, the numerical results for the unknowns are provided using the dimensionless units. For an applied pressure p, we set:

$$ (U_i, \Phi_i) = \frac{C_{00}}{hp} (u_i, \frac{\varphi_i}{E_0}) \quad (T_{ij}, D_k) = \frac{1}{p} (\sigma_{ij}, E_0 D_k) $$

where E₀ = 10⁹Vm⁻¹ and C₀₀ = 1GPa.

In Figure 8 one can observe, as before, that the rate of convergence of the solution of the approximated interface problem to the solution of the initial three phases problem is of order 3. As a consequence, even if the relative thickness of the interphase is large, the asymptotic approximation provides accurate results. More precisely, in the case of a relative thickness of the interphase equal to 0.1, the relative difference is equal to 0.07% for the displacement field and is equal to 0.1% for the electric potential.

Figure 9 presents the trend of the displacement field and electric potential on the orthogonal line to the interface mid-plane defined by (x = L₁/2, y = L₂/2, z/h ∈ [-0.5, 0.5]), showing a good numerical fit between the three-phases problem solution (solid line) and our generalized interface model solution (dotted line). As expected, the composite laminated plate, being more elongated on the x₁ axis, shows a Kirchhoff-Love kinematics with respect to U₁, while the displacement U₂ almost vanishes along the fiber. The electric potential presents a quadratic trend which coherent with the results proposed in Fernades and Pouget (2002). Moreover, since the PZT-4 adhesive layer is much more stiffer and has higher electric conductivity with respect to the PVDF adherents, it behaves as a rigid highly conducting interface, as shown in Figures 10–13, according to which the jump of the piezoelectric state almost vanishes at the interface, while the jump of the traction vector can be different from zero (see [T₃₃]).
## 3.2. The Thermo-Elastic Case

The interface problem (6) can be easily adapted in the case of linear elasticity with thermal effect. The thermo-elastic state is

defined by the pair $\mathbf{s} = (\mathbf{u}, \theta)$, constituted by the displacement field $\mathbf{u} = (u_i)$ and the variation of temperature $\theta := T - T_0$, where $T$ and $T_0$ denote the temperature field and a reference temperature, respectively. The corresponding constitutive law takes the following form:

<table><caption>TABLE 2 | Thermo-elastic material properties.</caption><thead><tr><th>Moduli</th><th>Upper adherent</th><th>Lower adherent</th><th>Adhesive</th></tr></thead><tbody><tr><td>Young's modulus <i>E</i>, GPa</td><td>14.53</td><td>200</td><td>2</td></tr><tr><td>Poisson's ratio <i>ν</i></td><td>0.33</td><td>0.33</td><td>0.2</td></tr><tr><td>Thermal expansion <i>β</i>, K<sup>-1</sup></td><td>0.8·10<sup>-6</sup></td><td>12·10<sup>-6</sup></td><td>76·10<sup>-6</sup></td></tr><tr><td>Thermal conductivity <i>k</i>, Wm<sup>-1</sup>K<sup>-1</sup></td><td>37</td><td>40</td><td>20</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>FIGURE 14 | Convergence with respect to the relative thickness $\varepsilon/h$.</figcaption></figure>

$$ \begin{cases} \sigma = C\mathrm{e}(\mathbf{u}) - \mathbf{X}\theta, \\ \mathbf{q} = -K\nabla\theta, \end{cases} $$

where $\mathbf{q} = (q_i)$ denotes the thermal flow, $\mathbf{X} = (X_{ij})$ and $\mathbf{K} = (K_{ij})$ represent, respectively, the thermal expansion tensor and the thermal conductivity tensor.

Adapting the same method used for the case of piezoelectricity, the variational problem for thermo-elasticity can be easily derived, where the bilinear form defined on the interface $S$ is given by

$$ A(\mathbf{s}, \mathbf{r}) = \int_S \left( \frac{1}{\varepsilon} \hat{\mathbf{K}}_{33} [[\mathbf{s}]] \cdot [[\mathbf{r}]] + \hat{\mathbf{K}}_{\alpha 3} \langle\langle \mathbf{u} \rangle\rangle_{,\alpha} \cdot [[\mathbf{v}]] \right. \\ \left. + \hat{\mathbf{K}}_{3\alpha} [[\mathbf{u}]] \cdot \langle\langle \mathbf{v} \rangle\rangle_{,\alpha} + \varepsilon \hat{\mathbf{K}}_{\alpha\beta} \langle\langle \mathbf{s} \rangle\rangle_{,\alpha} \cdot \langle\langle \mathbf{r} \rangle\rangle_{,\beta} - \right. \\ \left. - \langle\langle \theta \rangle\rangle (\hat{\mathbf{X}}^3 \cdot [[\mathbf{v}]] + \varepsilon \hat{\mathbf{X}}^\alpha \cdot \langle\langle \mathbf{v} \rangle\rangle_{,\alpha}) \right) d\tilde{\mathbf{x}}, $$

where $\mathbf{v}$ denotes the displacement test function and $\hat{\mathbf{X}}^k = (\hat{\mathbf{X}}_i^k) := (\hat{\mathbf{X}}_{ik})$.

In what follows, we consider the same geometrical setting of Figure 1, namely a thermo-elastic laminated plate occupying a 3D domain defined by $\Omega = [0, L_1] \times [0, L_2] \times [-h/2, h/2]$, with $L_1 = 10h$ and $L_2 = 5h$. The boundary and loading conditions are slightly different with respect to the piezoelectric cases. We suppose that:

* The union of the lateral boundaries $\Gamma_1$ is a free boundary;

<figure><img src="image_17.png" /><figcaption>FIGURE 15 | Displacement field and temperature on a section along the $x_3$ axis, with $\varepsilon/h = 0.04$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**FIGURE 16** | Jumps of the displacement field, temperature (**left**), traction vector, and thermal flux (**right**) across the interface on a line along the $x_1$ axis, with $\varepsilon/h = 0.04$.</figcaption></figure>

* A portion of (or the whole) top face $\Gamma_2$ is subject to a variation of temperature $\theta_1 = 1$ K around a given temperature $T_0$ and is mechanically free;

* The bottom face $\Gamma_3$ is vertically clamped and subject to a variation of temperature $\theta = 0$ K around a given temperature $T_0$;

* No mechanical volume or surface loads are applied.

The reference temperature $T_0$ is set at 310 K. We assume that both the adherents and the adhesive are made of thermo-elastic isotropic materials, such that $C_{ijkl} = \frac{E}{1+\nu}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{E\nu}{(1+\nu)(1-2\nu)}\delta_{ij}\delta_{kl}$, $K_{ij} = k\delta_{ij}$, $X_{ij} = \beta\delta_{ij}$, where $\delta_{ij}$ represents the Kronecker's delta. The constitutive coefficients are defined in **Table 2**.

**Remark.** These data comes from Rakotomanana and Ramaniraka (2004), devoted to the bonding of metallic implant on bones using cement. The upper body models the bone, the lower body represents the metallic implant, while the adhesive is made of classical cement used in this type of surgery. Let us notice that the rigidity of the interphase is very smaller than the rigidity of the implant (the ratio between the Young's modulus is equal to 1%). In that case, the displacements in the interphase is large compared to the displacements in the adherents. This case is generally called a soft interface problem, and it is often considered that the rigidity of the interphase depends linearly of the thickness in the asymptotic analysis (see, e.g., Dumont et al., 2018). It is worth noticing that in the present work the interphase is considered as a hard interface.

### 3.2.1. Numerical Results: The Applied Variation of Temperature

In all the numerical simulations, the problem is treated as a full coupled problem, with a variable **s** = (**u**, **θ**), with values in $\mathbb{R}^4$. The Finite Elements are piecewise linear on the tetrahedrons (commonly called P1 Finite Elements). The number of nodes of the mesh is equal to 5824 (resp. 6552) for the interface problem (resp. the three phases problem), leading to 24,275

<figure><img src="image_19.png" /><figcaption>**FIGURE 17** | Jumps of the displacement field, temperature (**left**), traction vector, and thermal flux (**right**) across the interface on a line along the $x_1$ axis, with $\varepsilon/h = 0.04$.</figcaption></figure>

(resp. 27,187) degrees of freedom. In this case, we did not use dimensionless unit.

In **Figure 14** is shown the trend of convergence of the solution of the interface problem to the solution of the three phases problem. Even in the less favorable case of a relative thickness $\varepsilon/h$ of $10^{-1}$, corresponding to comparable thickness for the adhesive and the adherents, the relative $L^2$ error is below 1% for the displacement field and $10^{-2}$% for the temperature field.

In **Figure 15** is reported the displacement and temperature fields on a section along the $x_3$ axis for a relative thickness of the interphase equal to 0.04. Continuous curves correspond to the numerical solution of the original three-dimensional problem, while the dotted curves to the data calculated by numerically implementing the transmission conditions. We obtain a relative difference of 0.06% in terms of displacement, and of $9 \cdot 10^{-4}$% in terms of temperatures. The applied temperature difference between the upper and lower faces provides an elongation along $x_3$ of the composite. The adhesive layer is characterized by weaker elastic material properties with respect to the adherents, while the thermal conductivities are slightly similar (see **Table 2**). Thus the interphase behaves as a soft interface from the mechanical point of view and almost hard interface in terms of the heat conduction, providing a jump on the displacement field and a quasi-constant temperature variation across the interface.

In **Figures 16, 17**, we present the jumps of the displacement, temperature, traction vector and normal heat flow across the interface on a line along the $x_1$ axis. Even if the jumps are large relatively to the displacements and the temperatures, one can conclude that the interface approximation is able to correctly reproduce their distributions. The soft elastic interface behavior is also corroborated by the vanishing of the traction vector jump, which is typical of weak interface models.

## 4. CONCLUSIONS

The validity of multiphysic interface transmission conditions as an approximated model of a thin multiphysic thin layer has
been discussed through numerical examples. The investigation
has be framed within the context of two different multiphysical
context: piezoelectricity and thermo-mechanical coupling. The
numerical simulation have been developed by using a finite
element approach which generalizes an analogous methodology
already proposed in Dumont et al. (2018) in the framework
of linearized elasticity. In the numerical examples, a composite
made of two phases joined by an adhesive has been studied in two
different configurations. In the first configuration, the adhesive
is a three-dimensional thin endowed with material parameters
taken to be softer or harder than the material parameters of
the two adherents. In the second configuration, the adhesive
is a material interface, whose transmission conditions are the
multiphysic interface laws obtained in Serpilli et al. (2019). The
relevant fields of the two configurations (displacement, electric
potential, temperature) are then compared to test the validity
of the proposed interface laws. The latter have been obtained
for several cases of material behavior of the adhesive: soft, hard
and rigid. These case are numerically compared with the three-
dimensional solution (the first configuration). A good agreement
has been found even if the elastic moduli of some layers of

the adhesive are of the same order as that of the adherents.
Moreover, we show that the rate of convergence of the L²-
relative error of the current model increases when compared
to the one obtained with two classical interface models, namely
the hard and soft interfaces. These findings clearly indicate
that the approach of substituting the multiphysic interphase
with the proposed interface law provides a robust modeling
for the composite and improves the numerical results with
respect to classical interface laws. This fact, intuitive but not
a priori obvious, represents the main original result of the
present paper.

DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to
the corresponding author.

AUTHOR CONTRIBUTIONS

SD: numerical analysis. MS and RR: asymptotic theory. FL: mechanical modeling.

REFERENCES

Benveniste, Y., and Miloh, T. (2001). Imperfect soft and stiff interfaces in two-dimensional elasticity. *Mech. Mater.* 33, 309–323. doi: 10.1016/S0167-6636(01)00055-2

Bessoud, A.-L., Krasucki, F., and Serpilli, M. (2011). Asymptotic analysis of shell-like inclusions with high rigidity. *J. Elastic*. 103, 153–172. doi: 10.1007/s10659-010-9278-1

Bonetti, E., Bonfanti, G., Lebon, F., and Rizzoni, R. (2017). A model of imperfect interface with damage. *Meccanica* 52, 1911–1922. doi: 10.1007/s11012-016-0520-1

Bonnet, M., Burel, A., Duruflé, M., and Joly, P. (2016). Effective transmission conditions for thin-layer transmission problems in elastodynamics. The case of a planar layer model. *ESAIM: M2AN* 50, 43–75. doi: 10.1051/m2an/2015030

Ciarlet, P. (1997). *Mathematical Elasticity, Vol. II: Theory of Plates*. Amsterdam: Elsevier.

Dumont, S., Rizzoni, R., Lebon, F., and Sacco, E. (2018). Soft and hard interface models for bonded elements. *Composit. B Eng.* 153, 480–490. doi: 10.1016/j.compositesb.2018.08.076

Fernades, A., and Pouget, J. (2002). An accurate modelling of piezoelectric multi-layer plates. *Eur. J. Mech. A Solids* 21, 629–651. doi: 10.1016/S0997-7538(02)01224-X

Fouchal, F., Lebon, F., and Titeux, I. (2009). Contribution to the modelling of interfaces in masonry construction. *Construct. Build. Mater.* 23, 2428–2441. doi: 10.1016/j.conbuildmat.2008.10.011

Geis, W., Mishuris, G., and Sandig, A. (2004). *Asymptotic Models for Piezoelectric Stack Actuators With Thin Metal Inclusions*. Berichte aus dem Institut fur Angewandte Analysis und Numerische Simulation 2004/001, Preprint IANS.

Geuzaine, C., and Remacle, J.-F. (2009). GMSH: a 3-D finite element mesh generator with built-in pre- and post-processing facilities. *Int. J. Num. Methods Eng.* 79, 1309–1331. doi: 10.1002/nme.2579

Geymonat, G., Hendili, S., Krasucki, F., Serpilli, M., and Vidrascu, M. (2014). “Asymptotic expansions and domain decomposition,” in *Domain Decomposition Methods in Science and Engineering XXI*. Lecture Notes in Computational Science and Engineering, eds J. Erhel, M. Gander, L. Halpern, G. Pichot, T. Sassi, and O. Widlund (Cham: Springer), 749–758. doi: 10.1007/978-3-319-05789-7_72

Geymonat, G., Krasucki, F., and Lenci, S. (1999). Mathematical analysis of a bonded joint with a soft thin adhesive. *Math. Mech. Solids* 16, 201–225. doi: 10.1177/108128659900400204

Gu, S. T., and He, Q. C. (2011). Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces. *J. Mech. Phys. Solids* 59, 1413–1426. doi: 10.1016/j.jmps.2011.04.004

Gu, S. T., Liu, J.-T., and He, Q. C. (2014). The strong and weak forms of a general imperfect interface model for linear coupled multifield phenomena. *Int. J. Eng. Sci.* 85, 31–46. doi: 10.1016/j.ijengsci.2014.07.007

Hashin, Z. (2002). Thin interphase/imperfect interface in elasticity with application to coated fiber composites. *J. Mech. Phys. Solids* 50, 2509–2537. doi: 10.1016/S0022-5096(02)00050-9

Javili, A., Kaessmair, S., and Steinmann, P. (2014). General imperfect interfaces. *Comput. Methods Appl. Mech. Eng.* 275, 76–97. doi: 10.1016/j.cma.2014.02.022

Klarbring, A. (1991). Derivation of a model of adhesively bonded joints by the asymptotic expansion method. *Int. J. Eng. Sci.* 29, 493–512. doi: 10.1016/0020-7225(91)90090-P

Lebon, F., and Rizzoni, R. (2010). Asymptotic analysis of a thin interface: the case involving similar rigidity. *Int. J. Eng. Sci.* 48, 473–486. doi: 10.1016/j.ijengsci.2009.12.001

Lebon, F., and Rizzoni, R. (2011). Asymptotic behavior of a hard thin linear interphase: an energy approach. *Int. J. Solids Struct.* 48, 441–449. doi: 10.1016/j.ijsolstr.2010.10.006

Lotfi, H., and Shing, P. (1994). Interface model applied to fracture of masonry structures. *J. Struct. Eng.* 120, 63–80. doi: 10.1061/(ASCE)0733-9445(1994)120:1(63)

Lourenço, P., and Rots, J. (1997). Multisurface interface model for analysis of masonry structures. *J. Eng. Mech.* 123, 660–668. doi: 10.1061/(ASCE)0733-9399(1997)123:7(660)

Mariano, P. (2002). Multifield theories in mechanics of solids. *Adv. Appl. Mech.* 38, 1–93. doi: 10.1016/S0065-2156(02)80102-8

Raffa, M., Lebon, F., and Rizzoni, R. (2018). Derivation of a model of imperfect interface with finite strains and damage by asymptotic techniques: an application to masonry structures. *Meccanica* 53, 1645–1660. doi: 10.1007/s11012-017-0765-3

Rakotomanana, L. R., and Ramaniraka, N. A. (2004). “Thermomechanics of bone-
implant interfaces after cemented joint arthroplasty: theoretical models and
computational aspects,” in *Computational Bioengineering*, eds M. Cerrolaza, G. Martinez, M. Doblare, and B. Calvo (London: Imperial College Press), 185–214. doi: 10.1142/9781860945403_0009

Renard, Y., and Pommier, J. (2002). *GetFem++*. *An Open Source Generic c++ Library for Finite Element Methods*. Technical report, INSA, Lyon.

Rizzoni, R., Dumont, S., and Lebon, F. (2017). On Saint Venant-Kirchhoff imperfect interfaces. *Int. J. Nonlinear Mech.* 89, 101–115. doi: 10.1016/j.ijnonlinmec.2016.12.002

Rizzoni, R., Dumont, S., Lebon, F., and Sacco, E. (2014). Higher order model for soft and hard elastic interfaces. *Int. J. Solids Struct.* 51, 4137–4148. doi: 10.1016/j.ijsolstr.2014.08.005

Serpilli, M. (2015). Mathematical modeling of weak and strong piezoelectric interfaces. *J. Elastic*. 121, 235–254. doi: 10.1007/s10659-015-9526-5

Serpilli, M. (2017). Asymptotic interface models in magneto-electro-thermo-elastic composites. *Meccanica* 52, 1407–1424. doi: 10.1007/s11012-016-0481-4

Serpilli, M. (2018). On modeling interfaces in linear micropolar composites. *Math. Mech. Solids* 23, 667–685. doi: 10.1177/1081286517692391

Serpilli, M. (2019). Classical and higher order interface conditions in poroelasticity. *Ann. Solid Struct. Mech.* 11, 1–10. doi: 10.1007/s12356-019-00052-5

Serpilli, M., and Lenci, S. (2012). Asymptotic modelling of the linear dynamics of laminated beams. *Int. J.*

Solids Struct. 49, 1147–1157. doi: 10.1016/j.ijsolstr.2012.01.012

Serpilli, M., and Lenci, S. (2016). An overview of different asymptotic models for anisotropic three-layer plates with soft adhesive. *Int. J. Solids Struct.* 81, 130–140. doi: 10.1016/j.ijsolstr.2015.11.020

Serpilli, M., Rizzoni, R., Lebon, F., and Dumont, S. (2019). An asymptotic derivation of a general imperfect interface law for linear multiphysics composites. *Int. J. Solids Struct.* 180–181:97–107. doi: 10.1016/j.ijsolstr.2019.07.014

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Dumont, Serpilli, Rizzoni and Lebon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.