Derivation of Imperfect Interface Laws for Multi-Physic Composites by a Multiscale Approach: Theoretical and Numerical Studies

In the present study, we focus our attention to a speciﬁc type of composite, constituted by two media, called the adherents, bonded together with a thin interphase layer, called the adhesive. We assume that the composite constituents are made of diﬀerent multi-physic materials with highly contrasted constitutive properties. The study considers a generic multi-physic coupling in a very general framework and can be adapted to well-known multi-physic behaviors, such as piezoelectricity, thermo-elasticity, as well as to multiﬁeld microstructural theories, such as micropolar elasticity.


Introduction
Structural bonding assembly has become an important technological solution over the past few years and is increasingly replacing bolting assembly (Ascione et al, 2017) as shown in Fig. 18.1). The resulting structure has many advantages, such as weight savings or the elimination of stress concentration. Similarly, in nature there are many living or natural structures that are composed of substructures, cells or soils for example, glued together. There are also many other examples of glued structures in the field of bioengineering (Breschi et al, 2008) as illustrated in Fig. 18.2. Understanding and modelling the bonding process then becomes a necessity.  An obvious common point between all these bonded composite structures is the thinness of the adhesive compared to those of the substrates or adherents. This is true for both industrial structures and living or natural structures. This thinness will obviously lead to numerical modelling difficulties. Indeed, the mesh size of the glue will mechanically lead to computations with a very large number of degrees of freedom and therefore very expensive computations. These costs will increase further if the adhesive surface is irregular and has a high roughness. Similarly, in the presence of kinematic or behaviour non-linearities, in the presence of cracks, etc., the costs become prohibitive. "Direct" calculations are then limited to academic cases . There are two very different possibilities, either to develop suitable numerical methods (Alart and Lebon, 1995;Alart et al, 1997; Barbie et al, 2015) or to set up macroscopic models of the adhesive's behaviour. In this chapter, we will focus on the second family of solutions.
There are at least two methodologies in the literature for obtaining constitutive laws of adhesive (interface) behaviour. The most classic is the introduction of phenomenological models, usually based on experimental results, such as Coulomb models, compliances, etc. In this chapter, we will prefer to focus on deductive models. The idea of this methodology is to start from a micromechanical study of the adhesive (interphase) and to deduce, using mathematical methods, an equivalent macroscopic behaviour (interface).
In this chapter, we focus our attention to a bonded composite, constituted by two adherents, bonded together with a thin adhesive. We assume that the composite constituents are made of different multi-physic materials with highly contrasted constitutive properties. The study considers a generic multi-physic coupling in a very general framework and can be adapted to well-known multi-physic behaviors, such as piezoelectricity, thermo-elasticity, as well as to multifield microstructural theories, such as micropolar elasticity (see, e.g. Chatzigeorgiou et al, 2015). Several works have suggested a generalization of the classical interface models, including the effects of other physical (thermal, piezoelectric, etc.) interactions (dell'Isola and Romano, 1987;Chen, 2008;Wang et al, 2017;Firooz and Javili, 2019;Saeb et al, 2019), and within the framework of linear multifield theories, such as higher order continua theories (Placidi et al, 2014;Eremeyev, 2019).
The analysis has been carried out by means of the asymptotic expansions method, using the thickness as a small parameter. This technique is based on the fact that the thickness of the adhesive can be considered as a small parameter (intended to tend towards zero) and denoted by ε in the following. The asymptotic analysis has been applied to the rigorous derivation of simplified models for complex assemblies, presenting thin interphases, in the field of linear elasticity (Lebon and Rizzoni (2010); Dumont et al (2018); Rizzoni et al (2014); Serpilli and Lenci (2016)) as well as in piezoelectricity, taking into account other physical interactions, micropolar elasticity and poroelasticity (Serpilli et al (2013); Serpilli (2015Serpilli ( , 2017Serpilli ( , 2018Serpilli ( , 2019). As mentioned above, the asymptotic methods allow to replace the adhesive layer with a two-dimensional surface, the so-called imperfect interface, with non-classical transmission conditions between the two adherents. By defining the small parameter and constitutive properties of the middle layer, we perform an asymptotic analysis. We assume that the multi-physic stiffness ratios between the adherents and the adhesive depend on ε p . As proposed by Caillerie (1970), we identify three critical exponents p, corresponding to different imperfect interface models: p = 1, the soft (also called lowly-conducting) multi-physic interface model; p = 0, the hard (also called moderately-conducting) multi-physic interface model; p = −1, the rigid (also called highly-conducting) multi-physic interface model. Following the approach introduced by Rizzoni et al (2014), we characterize the order zero and the order one transmission problems. Finally, a general multi-physic interface model is developed, and numerically tested through the finite element method. In particular, in the framework of piezoelectricity, we compare the results obtained by modeling the adhesive as an interphase, having a thin finite thickness, with the results obtained with the general multi-physic interface model.

Statement of the Problem
We consider the composite assembly constituted of two solids Ω ε ± ⊂ R 3 , called the adherents, bonded together by an intermediate thin layer B ε := S × (− ε 2 , ε 2 ) of thickness ε, called the adhesive, with cross-section S ⊂ R 2 . In the following B ε and S will be called interphase and interface, respectively. Let S ε ± be the plane contact surfaces between the adhesive and the adherents and let Ω ε := Ω ε + ∪ B ε ∪ Ω ε − denote the composite system comprising the interphase and the adherents (cf. Fig.  1.3a). We suppose that the composite is constituted by a multi-physic material, in which different physical behaviors interact together, such as in piezoelectricity. Its equilibrium state is determined by a collection of order parameters s ε := (u ε 1 , . . . , u ε N , ϕ ε 1 , . . . , ϕ ε M ): N vector state variables, namely u ε i , and M scalar state variables, namely ϕ ε k . With the multi-physic state s ε , we associate its conjugated physical quantity t ε = t ε (∇ ε s ε ), where ∇ ε s ε denotes the gradient of s ε . The vector field t ε := (σ ε 1 , . . . , σ ε N , D ε 1 , . . . , D ε M ) represents a generalized stress field. We also consider the following homogeneous and linear constitutive law: where K ε is a generalized linear constitutive matrix. The constitutive tensor K ε satisfies suitable symmetry and positivity properties.
We assume that the adherents are subject to a generalized system of body Body and surface forces are neglected in adhesive We assume that on Γ lat := ∂S × (− ε 2 , ε 2 ) homogeneous Neumann boundary conditions are applied. The differential formulation of the governing equations has the following structure: represents the generalized traction vector on the boundary Γ ε g and n ε the outer normal unit vector to Γ ε g . Let us introduce the functional space The variational formulation of problem (18.1) defined on the variable domain Ω ε can be written as follows:

Method of Asymptotic Expansion
In order to perform an asymptotic analysis of problem (18.2) when ε tends to zero, we rewrite the problem on a fixed domain Ω independent of ε. By using the approach of Ciarlet (1997), we consider the change of variables π ε : where, after the change of variables, the adherents occupy 2 } denote the interfaces between B and Ω ± and Ω = Ω + ∪ Ω − ∪ B is the rescaled configuration of the composite. Lastly, Γ g and Γ u indicate the images through π ε of Γ ε g and Γ ε u (cf. Fig. 1.3b). Consequently, We assume that the constitutive coefficients of Ω ε ± are independent of ε,K ε = K, while the constitutive coefficients of B ε depend on ε,K ε = ε pK , with p ∈ {−1, 0, 1}. Finally, we assume that the forces are such that L ε (r ε ) = L(r). By virtue of the previous hypothesis, the rescaled problem can be written in the following form: andK ij denote the sub-matrices ofK, defined bŷ We can now apply the asymptotic expansions method to the rescaled problem (18.3), whose fundamental assumption relies in considering the solution s ε of the problem as a series of powers of ε: wheres ε = s ε •π ε andŝ ε = s ε •π ε . By injecting (18.4) into the rescaled problem (18.3), and by identifying the terms with identical power of ε, we obtain, as customary, a set of variational problems to be solved in order to characterize the limit multi-physic state s 0 , the first order corrector term s 1 and their associated limit problem, for p ∈ {−1, 0, 1}. Following the approach described in Rizzoni et al (2014); Dumont et al (2018), we introduce the matching conditions based on the continuity of the generalized traction t ε e 3 and multiphyisic state s ε at the interfaces S ε ± in the initial configuration and on the continuity of the traction and statet ε e 3 ,s ε ,t ε e 3 ,ŝ ε at the interfaces S ± in the rescaled configuration. Hence, one has denote, respectively, the mean value and the jump functions at the interfaces.

Multi-Physic Interface Models
In this section we present the asymptotic models for multi-physic interfaces obtained for the soft, hard and rigid cases at order 0 and order 1. For the sake of brevity, we will skip all the mathematical computations carried out in the deduction of the limit models.

The Soft Multi-Physic Interface
The transmission problems at order 0 and order 1 can be summarized as follows: • Order 0

Governing equations
The transmission problems for a soft multi-physic interface at order 0 and order 1 represent a formal generalization of the soft interface models obtained by means of the asymptotic methods in linear elasticity (see, e.g., Rizzoni et al, 2014;Dumont et al, 2018) and in other multifield frameworks, such as poroelasticity (see Serpilli, 2019). At order 0, the interface behaves from a mechanical point of view as a series of linear springs, reacting to the discontinuity of the multi-physic state between the upper and bottom faces, while the generalized traction vector remains continuous. At order 1, the interface conditions maintain a similar structure, but both the multiphysic state and the traction vector are discontinuous through the interface. Moreover, they depend on the in-plane derivatives of the jump and mean values ofs 0 , which can be considered a known source term, identified in the order 0 problem.

The Hard Multi-Physic Interface
The hard interface transmission problems at order 0 and order 1 take the following expressions: • Order 0

Governing equations
It is interesting to notice that the hard multi-physic interface problems is equivalent to the ones derived in the case of linear elasticity in Lebon and Rizzoni (2010); Rizzoni et al (2014); Dumont et al (2018). At order 0, we recover the classical continuity conditions for both the multi-physic state and generalized traction vector. Thus, the adherents are perfectly bonded together. At order 1, we encounter a mixed interface model with a jump of the state and traction vector depending on the values of the multi-physic state and traction vector at order 0. These order 0 terms, being known from the solution of the previous problem, can be viewed as external source terms.

The Rigid Multi-Physic Interface
The differential formulations of the rigid interface problems at order 0 and order 1 take the following form: • Order 0

Transmission conditions on S
The rigid multi-physic interface problems show the same features of the rigid interface asymptotic models obtained in different frameworks in Bessoud et al (2009);Serpilli (2015Serpilli ( , 2017Serpilli ( , 2018Serpilli ( , 2019. Concerning the order 0 model, we obtain a continuity of the multi-physic state at the interface level, while the traction vector is discontinuous and depends on the divergence of a generalized membrane stress vector N 0 α :=L αβ s 0 ,β . The interface behaves as a multi-physic membrane. The order 1 presents a discontinuity on both the multiphysic state and traction vector. Analogously to the order 0 model, we obtain a generalized equilibrium membrane problem defined on the plane of the interface.

The General Multi-Physic Interface
The approach of Rizzoni et al (2014) can be extended in order to obtain a condensed form of the transmission conditions summarizing both the orders 0 and 1 of the soft, hard and rigid cases in only one couple of equations in terms of the jump of the multi-physic state and generalized tractions at the interface. Therefore, we denote bys ε :=s 0 + εs 1 + ε 2s2 andt ε :=t 0 + εt 1 , two suitable approximations fors ε andt ε . Let us consider the rigid multi-physic interface conditions, as starting point. After rescaling back the constitutive coefficientsK = εK ε in B ε , we can write [s ε ] and [t ε e 3 ]. Indeed, one has An alternative expression of the above transmission conditions can be given in terms of t ε e 3 and t ε e 3 , which will be useful to write the variational formulation of the interface multi-physic problem: It is easy to prove that this interface law is general enough to describe the interface laws at order 0 and order 1 prescribing the multi-physic state jump and traction jump in the cases of the soft and hard interfaces, by choosing the following scalings for the constitutive matrices:K ε = εK, for the soft case, andK ε =K, for the hard case.
The relations (18.6) can be transformed into interface equations defined on S, by making use of the matching relations (18.5), up to higher order terms: (18.7)

Finite Element Implementation and Numerical Test
In order to derive the variational form of the multi-physic problem, which will be numerically tested through a finite element procedure, one can write the variational form of the equilibrium problem on each sub-domain Ω + and Ω − : which can be rewritten as  A standard finite element method is employed to solve this equation. In order to take into account the jumps in the displacements across the interface, a 'flat" finite element is considered on the interface S that has all nodes on S, the first ones related to Ω − , and the other ones related to Ω + . It is then possible to write a stiffness matrix of this problem that is invertible and with standard error estimates (for more details, see for example Nairn, 2007).

Numerical Study: The Piezoelectric Composite Plate
The aim of this study is to numerically test the general interface law, expressed in (18.8), comparing it to a three-dimensional analysis of the problem. As preliminary test, we consider the piezoelectric case. The piezoelectric state at the equilibrium is determined by the pair s := (u, ϕ), where u and ϕ represent the displacement field and the electric potential. The generalized stress vector is given by t := (σ, D), where σ and D denote, respectively, the Cauchy stress tensor and the electric displacement.
Let us consider a piezoelectric three-phases composite plate, which occupies a 3D domain defined by Ω = [0, L 1 ] × [0, L 2 ] × [−h/2, h/2], with L 1 = 10h and L 2 = 5h (see Fig. 1.4). The adhesive thickness is set to be ε. The adherents are constituted by PVDF (Polyvinylidene fluoride), a monoclinic piezoelectric material with poling axis e 3 , while the adhesive is made of PZT-4, a transversally isotropic piezoelectric material with poling axis e 3 . This constitutive sub-matrices (K ij ) are defined as follows:  Table 18.1. The piezoelectric composite plate is subject to surface uniform load equal to p = 1 kN/m 2 on the top face, as shown in Fig. ??. We assume that no voltage is applied on the upper and lower faces. In this case, the composite plates behaves as a sensor (see Bonaldi et al, 2017). The finite element simulations (made with GetFEM) are carried out using Q1 elements (linear on a cube), with 7280 nodes (29203 degrees of freedom) for the three phases problem and 5824 nodes (23379 degrees of freedom) for the problem with the interface law.
First, the influence of the relative thickness of the interphase ε L is investigated in order to evaluate the accuracy of the various modeling proposed in the previous sections. In particular, the quality of the solutions is evaluated considering the L 2relative error e = s ε −s model s ε , where s ε denotes the reference solution computed using the three-phases problem with a fine finite element mesh, while s model indicates the solution of the interface model (18.8). The convergence of the general interface model towards the three-phases one with respect to the thickness ratio ε L is presented in Fig. 18.5. 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 ε . The convergence rate is ε 3 . Besides, even if the relative thickness is of 1%, the relative error is equal to 7.65 · 10 −2 % for the displacement field and 9.06 · 10 −4 % for the electric potential, meaning that the general interface model can also be used for moderately thick adhesive layers. Now, let us fix the relative thickness ε L = 0.02. The numerical results for the variables are provided using the dimensionless units. We set: where we have chosen, for numerical convenience, V = 50V , E 0 = 10 9 V m −1 and C 00 = 1GP a. The results are represented in Fig. 18.6, 18.7 and 18.8. Figure 18.6 represents the trend of the displacement field and electric potential, evaluated in x = L 1 /2, y = L 2 /2, z/h ∈ [−0.5, 0.5]. The plot shows a 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 potential is 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. Figure 18.7 and Fig. 18.8 represent the trend of the jumps of the displacement and electric potential and the jumps of the stress vector and normal electric displacement along the x-axis, namely (x ∈ [0, L 1 ], y = L 2 /2, z = 0), and y-axis, namely (x = L 1 /2, y =∈ [0, L 2 ], z = 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.

Concluding Remarks
General imperfect contact conditions have been proposed, simulating the behavior of a thin interphase undergoing linear coupled multi-physic phenomena. These con- ditions link the generalized traction vector field and its jump to the multi-physic state vector field and its jump at the interface, which is the geometric limit of the interphase as its thickness parameter ε goes to zero. Three interface models (soft, hard and rigid) have been deduced by means of the asymptotic methods, by varying the rigidity ratios between the adhesive and adherents and considering the order 0 and order 1 corrector terms. Furthermore, these three different models have been condensed in one general imperfect interface model and its variational formulation has been presented. The weak formulation represents a key step towards the FEM simulation. A numerical example has been presented considering a piezoelectric composite plate, subject to an electric potential difference at the top and bottom faces. The numerical results show the convergence of the solution of the three-phases model towards the solution of the proposed model as ε tends to zero, highlighting the accuracy and "goodness" of the general interface conditions.