Raul Guinovart-Diaz

Reinaldo Rodriguez-Ramos
email: reinaldo@matcom.uh.cu

Juan Carlos López-Realpozo

Julian Bravo-Castillero

José Antonio Otero
email: j.a.otero@itesm.mx

Frederico J Sabina

Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Serge Dumont
email: serge.dumont@u-picardie.fr

R Guinovart-Díaz

R Rodríguez-Ramos

J C López-Realpozo

F J R In

Guinovart-Díaz

Analysis of fibrous elastic composites with nonuniform imperfect adhesion

à 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.

matrix-fiber contact which is called nonuniform spring imperfect contact. In this case, the contrast or jump of the displacements in the boundary of each phase is proportional to the corresponding component of the tension in the interface in terms of a parameter given by a certain function that depends on the position. The constituents of the composites exhibit transversely isotropic properties. A doubly periodic parallelogram array of cylindrical inclusions under longitudinal shear is considered. The three-phase model is validated by the Finite Element Method and the AHM both approaches applied to two-phase composites with nonuniform spring imperfect contact. Comparisons with theoretical and experimental results verified that the present model is efficient for the analysis of composites with presence of nonuniform imperfect interface and parallelogram cell. The effect of the nonuniform imperfection on the shear effective property is observed. The present method can provide benchmark results for other numerical and approximate methods.

Introduction

Composite materials containing microstructural constituents are very extensively used in industrial applications. However, defects induced due to manufacturing processes have been recognized as a general problem in this type of composites. These defects can be found in the fiber (i.e., fiber misalignment, irregular fiber distribution and broken fibers), in the matrix (i.e., voids and non-homogeneity) and/or in the interface regions (i.e., debonding, delimitations and imperfect interface) [START_REF] Huang | Effects of void geometry on elastic properties of unidirectional fiber reinforced composites[END_REF][START_REF] Mansonwakeman | Composite processing and manufacturing: an overview[END_REF]. The presence of such defects affects the overall properties of composites and nonhomogeneous materials [START_REF] Kushch | Transverse conductivity of unidirectional fibrous composite with interface arc cracks[END_REF][START_REF] Prakash | Significance of defects in the fatigue failure of carbon fibre reinforced[END_REF]. An imperfect interfacial bonding is a common defect in composites, which can significantly affect the mechanical properties of the material [START_REF] Theocaris | The Mesophase Concept in Composites[END_REF].

In most composites, the fiber-matrix adhesion is imperfect; the continuity conditions for stresses and displacement are not satisfied. Thus, various approaches have been used, in which the bond between the reinforcement and the matrix is modeled by an interphase with specified thickness [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF][START_REF] Rodríguez-Ramos | Closed-form thermo-elastic moduli of a periodic three-phase fiber-reinforced composite[END_REF]. Another approach consists in considering that the thickness and the elasticity coefficients of the bond are small parameters compared to those of the other components of the composite. Then, using matched asymptotic expansions, it can be shown that in that case, the behavior of the interphase can be accurately approximated by an interface with an appropriate contact law, where the jump in the displacements is proportional to the stress, which is supposed to be continuous (see, for example, [START_REF] Abdelmoula | Comportement asymptotique d'une interface mince[END_REF][START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF][START_REF] Klarbring | Derivation of a model of adhesively bonded joints by the asymptotic expansion method[END_REF][START_REF] Rizzoni | Asymptotic analysis of an elastic thin interphase with mismatch strain[END_REF]).

This kind of imperfect contact (spring type) in the modeling of composites was investigated by Benveniste and Miloh [START_REF] Benveniste | Imperfect soft and stiff interfaces in two dimensional elasticity[END_REF] among others and has been used, for instance, by Achenbach and Zhu [START_REF] Achenbach | Effect of interphases on micro and macro mechanical behaviour of hexagonal-array fiber composites[END_REF] and Hashin [START_REF] Hashin | Thermoelastic properties of fiber composites with imperfect interface[END_REF][START_REF] Hashin | The spherical inclusion with imperfect interface[END_REF][START_REF] Hashin | Thermoelastic properties of particulate composites with imperfect interface[END_REF]. Some previous studies also reported the elastic effective coefficients for two-phase fibrous composites with rhombic array of periodic cells, both with perfect and with imperfect contact conditions, analytically and numerically [START_REF] López-Realpozo | Effective elastic shear stiffness of a periodic fibrous composite with non-uniform imperfect contact between the matrix and the fibers[END_REF][START_REF] Molkov | Effective characteristic of fibrous unidirectional composite with periodic structure[END_REF][START_REF] Abolfathi | A micromechanical characterization of angular bidirectional fibrous composites[END_REF][START_REF] Jiang | A rigorous analytical method for doubly periodic cylindrical inclusions under longitudinal shear and its application[END_REF].

In this work, a micromechanical analysis method is applied to a periodic composite with unidirectional fibers and parallelogram cells, in particular, rhombic periodic cells. The analytical expressions of the homogenized elastic properties are calculated for a two-phase composite with interphase. In this study, coated fibers or a thin mesophase is considered which often contains the contrast or jump of the displacements in the boundary of matrix and fibers [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF]. The mesophase is assumed as a function that depends on the position (nonuniform imperfect adherence). The asymptotic homogenization method (AHM), for three-phase fibrous periodic composites with nonuniform mesophase and oblique cell, is used for the calculation of shear elastic effective coefficients. This contribution is an extension of previous works [START_REF] Rodríguez-Ramos | Two analytical models for the study of periodic fibrous elastic composite with different unit cells[END_REF][START_REF] Guinovart-Díaz | Influence of parallelogram cells in the axial behaviour of fibrous composite[END_REF][START_REF] López-Realpozo | Transport properties in fibrous elastic rhombic composite with imperfect contact condition[END_REF] where the imperfection on the interface is a constant function (uniform). The results in this paper are mainly focused on the impact of the arrangement of the fibers and the mechanical nonuniform imperfection at the interface on the stiffness properties. The novelty of the present work is that the imperfection of the interface in the composite with parallelogram cell is studied by two different models (three-phase and spring models) where the interface is in general nonuniform, which makes the composite more complex due to the presence of an heterogeneous interface. Moreover, both theoretical approaches are validated by the Finite Element Method.

General considerations for three-phase composites: contact at the interface

Let us consider a composite of unidirectional periodic fibers embedded in a matrix as shown in Fig. 1. The angle of the cell μ is assumed to remain constant so that a parallelogram cell with periods w 1 and w 2 can be defined. The periodicity of the microstructure determines the geometry of the periodic cell S. Thus, a three-phase periodic composite is considered here which consists of a parallelogram array of identical parallel circular concentric cylinders embedded in a homogeneous medium (Fig. 2). As a unidirectional fibrous composite, it is assumed that the microstructure of the composite along the third direction (perpendicular to the plane of cross section) remains constant. The central fibers are all assumed straight and of circular cross section with radius R, and the mesophase has a thickness h (see Fig. 2). The material properties of each phase belong to the crystal symmetry class 6 mm, where the axes of material and geometric symmetry are parallel. The composite is not well bonded at the contact between the matrix and the fiber. Therefore, imperfect contact conditions at the interface ε are considered, where ε = l/L is a small geometrical parameter relating the distance l between the centers of two neighboring cylinders and the diameter L of the composite.

A family of boundary value problems depending on ε for the linear elastic heterogeneous media can be formulated in general form for a three-phase composite as follows:

∇ • C ε : ∇u ε = 0, x ∈ , (1) T(u ε ) = 0, x ∈ ε , (2) [[u ε ]] = 0, x ∈ ε , (3) u = u 0 on ∂ , ( 4
)
where C ε (x) = C(y) is the S-periodic stiffness elastic tensor and y = x/ε is the local variable. Moreover, ∇ is the gradient operator and u ε = (u ε 1 , u ε 2 , u ε 3 ) the mechanical displacement vector referring to the Cartesian coordinate system. u 0 is a prescribed displacement at the boundary ∂ . T is the traction vector, and [[ f ]] γ (γ = 1, 2) indicates the jump in the quantity at the common interface ε between the matrix and the mesophase, and the mesophase with the fiber, i.e, (2) .

[[ f ]] 1 = f (1) -f (I ) , [[ f ]] 2 = f (I ) -f
Since a periodic composite with nonuniform interphase is studied, three distinct phases, occupying S 1 , S I and S 2 (Fig. 2), are assumed. The region S I characterizes the non-perfect contact along the interfaces 1 = z = Re iθ , 0 ≤ θ ≤ 2π and 2 = z = (R + h)e iθ , 0 ≤ θ ≤ 2π of each cylinder of radius R and R 1 = R + h, respectively, where h denotes the thickness of the mesophase S I . The effect of mechanical imperfection in the effective properties of the composite is incorporated through the presence of the mesophase S I . In the forthcoming sections, some summary ideas of the antiplane problem and a brief comment about the asymptotic homogenization method (AHM) are given.

2.1 Asymptotic homogenization method (AHM) for the antiplane problems α3 L It is well known that the problem (1)-( 4) is decoupled into two independent systems for unidirectional fibrous composites with components of transversely isotropy symmetry. One of them involves u ε 1 , u ε 2 , i.e., it is a state of in-plane mechanical deformation field. The other state, which is of particular interest in this work, is characterized by an out-of-plane mechanical displacement u ε 3 which is a function of the variables x 1 and x 2 . The main aim of this section is the determination of effective properties in three-phase composites for the out of plane using the homogenization method, considering that the mesophase S I is not a homogeneous medium with its properties being functions of the angle θ with vertices in the center of the fiber (Fig. 3). In this case, the relevant constitutive relations are

σ 13 = C 1313 u 3,1 , σ 23 = C 2323 u 3,2 , ( 5
)
where σ 13 and σ 23 are the components of out-of-plane mechanical stress, C 1313 and C 2323 are the elastic moduli, and the comma denotes partial differentiation.

Here, we are interested to consider the antiplane case, and the BVP (1)-( 4) is transformed into the problem

C 3η3β (y)u 3,β ,η = 0, in , (6)
[[C ε 1313 (u ε 3,1 n 1 + u ε 3,2 n 2 )]] = 0, on ε , (7) [[u ε 3 ]] = 0, on ε , ( 8
) u ε 3 = u 0 , on ∂ , ( 9
)
Fig. 3 Three-phase composite with nonuniform interface where the Greek indices run from 1 to 2 and ε = 1 ∪ 2 . The usual convention of index summation is used. Equation ( 6) has S-periodic coefficients and they are rapidly oscillating (see, Fig. 1). In order to obtain the homogenized problem and the corresponding effective coefficients, the solution of ( 6)-( 9) is sought using the method of two scales by the ansatz:

u ε 3 (x) = v 0 (x) + εv 1 (x, y) + O(ε 2 ), x = (x 1 , x 2 ), y = (y 1 , y 2 ), (10)
with v 1 being an S-periodic function with respect to y. Substituting ( 10) into ( 6), applying the chain rule considering that x and y are independent and equating the terms of orders ε -1 , ε 0 to zero, one can obtain that

v 1 (x, y) = α3 N (y)v 0,α (x), where v 0 (x) is the solution of the homogenized equation C * 3α3β v 0,αβ = 0, where C * 3α3β = C 3α3β (y) + C 3α3η (y) α3 N ,η (y)
are the so-called effective coefficients and α3 N (y) is an S-periodic harmonic function of the α3 L local problem. The comma will also be used to denote a derivative with respect to the local variable. The angular brackets define the volume average per unit length over the unit cell, that is, F = 1

V S F(y)dy. The index α = 1, 2 denotes two different problems over the periodic cell which are to be solved. It is a well-known derivation whose details can be found elsewhere and is omitted here (see, for instance, [START_REF] Pobedrya | Mechanics of composite materials[END_REF][START_REF] Bakhvalov | Homogenisation: Averaging Processes in Periodic Media[END_REF]).

The local problems are referred to as 13 L and 23 L with solutions 13 N (y) and 23 N (y), respectively. The solution of such local problems on the periodic cell described in Fig. 2, taking into account the perfect contact conditions ( 7) and ( 8), is one of the main objectives of the present work to obtain the overall properties C * 55 , C * 54 , C * 45 and C * 44 with the short index notation. The alternative form is found as follows:

C * 55 = p 1 V 1 + p I V I + p 2 V 2 + p 13 N ,1 , (11)
C * 54 = p 13 N ,2 , C * 45 = p 23 N ,1 , (12)
C * 44 = p 1 V 1 + p I V I + p 2 V 2 + p 23 N ,2 , (13)
where

p γ = C (γ )
1313 with (γ = 1, I, 2 ) are the shear moduli of the matrix, mesophase and fibers, respectively.

V 2 = π R 2 /V, V I = π 2h R + h 2 /V and V 1 = 1 -V I -V 2 are
the percentages of concentrations of fiber, mesophase and matrix, respectively, and

V = |w 1 | |w 2 | sin μ is the area of periodic cell.
The mathematical statement for the two local problems α3 L(α = 1, 2) in the case of three-phase composites consist of finding doubly periodic functions α3 N = α3 N (γ ) (y) if y ∈ S γ , that satisfy the following Laplace equations with the contact condition in each phase (γ = 1, I, 2):

∇ 2 α3 N (γ ) = 0 in S γ , ( 14
)
where the continuity condition for the stress and displacement at the interfaces 1 and 2 depends on the angle θ of the point z = re iθ ∈ S I , which describes a particular type of imperfect contact between the matrix and fiber with nonhomogeneous mesophase. The interface ε = 1 ∪ 2 by the continuity of displacement, therefore, yields

[[ α3 N ]] = 0 on ε = 1 ∪ 2 , ( 15
)
and the continuity condition for the stress

[[ p γ ( α3 N (γ ) ,1 n 1 + α3 N (γ ) ,2 n 2 )]] = -[[ p γ ]]n α on ε = 1 ∪ 2 , ( 16
)
and the condition for uniqueness

α3 N = 0 in S = S 1 ∪ S I ∪ S 2 . ( 17
)
The outward unit normal vector to the interface is n = (n 1 , n 2 ). The pre-index used in ( 14)-( 17) for the local function will be omitted for simplicity. The problem ( 14)-( 17) should be converted into dimensionless problems making the appropriate change. Indeed, defining the dimensionless variable ξ = y/l (this is not a small parameter!), then

N (γ ) ,β = u (γ ) ,β and N (γ ) ,ββ = u (γ )
,ββ /l where u (γ ) ≡ N (γ ) /l and the derivative u ,β is with respect to the variable ξ β . R 0 is the true radius of the fibers in the composite. Also the dimensionless parameter R = R 0 /l is introduced.

Solution of local problems α3 L (α = 1, 2) for three-phase composites with nonuniform imperfect mesophase

The well-developed theory of analytical functions in Muskhelishvili [START_REF] Muskhelishvili | Some Basic Problems in the Mathematical Theory of Elasticity[END_REF] can be applied to solve the problem ( 14)- [START_REF] López-Realpozo | Effective elastic shear stiffness of a periodic fibrous composite with non-uniform imperfect contact between the matrix and the fibers[END_REF]. Thus, it is necessary to solve the local problems α3 L. Doubly periodic harmonic functions in the matrix (γ = 1), mesophase (γ = I ) and fiber (γ = 2) region are to be found for the α3 L local problems in terms of harmonic functions ϕ γ (z) analogous to Lopez-Realpozo et al. [START_REF] López-Realpozo | Transport properties in fibrous elastic rhombic composite with imperfect contact condition[END_REF],

u (γ ) = Re ϕ γ (z) . ( 18
)
Since the composite exhibits a nonuniform mesophase between the matrix and fiber, we consider a partition of the cell S γ = n j=1 S γ j , where

S γ j (θ ) = {z = re iθ : z ∈ S γ , θ j-1 < θ < θ j }. This partition divides the interfaces β = n j=1 β j , βi ∩ β j = ∅, i = j, β = 1, 2, into n arcs.
Over the mesophase S I , the material property p I is a piecewise constant function p I (θ ) = p I j , θ j-1 < θ < θ j with j = 1, . . . , n, θ 0 ≡ 0, θ n ≡ 2π and θ j denotes the angle of each arc (Fig. 3).

In each region S γ j , we consider

ϕ 1 j = a 0 j z R 1 + ∞ k=1 o ζ (k-1) (z/R 1 ) (k -1)! a k j , ( 19
)
φ I j (z) = ∞ p=1 o b pj R 1 z p + ∞ p=1 o b -pj z R p , ( 20
)
ϕ 2 j (z) = ∞ p=1 o z R p c pj . ( 21
)
ζ(z) is the Weierstrass quasi-periodic function, which is defined as

ζ(z) = 1 z + ∞ m,n 1 z-P nm + 1 P nm + z P 2
nm , P nm = nw 1 + mw 2 , for m, n ∈ Z , and the prime over the summation symbol means that the pair (m, n) = (0, 0) is excluded. The Laurent's expansion of function ϕ j 1 (z) is given by the following expression:

ϕ 1 j (z) = z R a 0 j + ∞ p=1 o R 1 z p a pj - ∞ p=1 o ∞ k=1 o z R 1 p k p w kp a k j , ( 22
)
where

w kp = (k+ p-1)! (k-1)!( p-1)! S k+ p R k+ p 1 √ kp , S k+ p = m,n (mw 1 + nw 2 ) -(k+ p) , m 2 + n 2 = 0, k + p ≥ 2, S 2 = 0. The constants a 0 j , a pj , b pj , b -pj ,
c pj and z = ξ 1 + iξ 2 are complex numbers; the overbar denotes the complex conjugate and the superscript "o" on the summation indicates that the sum is carried out only over odd indices. w 1 , w 2 are the periods and

a 0 j = -R 2 1 H 1 ā1 j -R 2 1 H 2 a 1 j , where H 1 = δ1 w2 -δ2 w1 w 1 w2 -w 2 w1 , H 2 = δ 1 w2 -δ 2 w1 w 1 w2 -w 2 w1 , with δ i = ζ(z + w i ) -ζ(z).
Replacing ( 18)-( 22) into ( 15)-( 16) and performing some algebraic manipulations, the following system of equations is obtained in order to find the complex unknown coefficients a k j for each arc by the system:

āpj + χ 1 j R 2 1 H 1 δ 1 p ā1 j + χ 1 j R 2 1 H 2 δ 1 p a 1 j + χ pj ∞ k=1 o w kp a pj = R 1 χ 1 j δ 1 p (δ 1α -iδ 2α ) , ( 23
)
where δ kp is the Kronecker delta, χ pj =

(V 2 +V I ) p (κ 1 j +κ 2 )(1-κ 1 j )+V p 2 (κ 1 j -κ 2 )(1+κ 1 j ) (V 2 +V I ) p (κ 1 j +κ 2 )(1+κ 1 j )+V p 2 (κ 1 j -κ 2 )(1-κ 1 j ) , κ 2 = p 2 / p 1
and κ 1 j = p I j / p 1 are the components of shear ratios for each region S I j (θ ) of the mesophase with the shear property of the matrix.

Effective properties for three-phase composites under nonuniform mesophase

The Eqs. ( 11)-( 13) are easily transformed applying Green's theorem to the area integrals. Replacing N (γ ) = u (γ ) /l, dy i = ldξ i and neglecting the line integral over the radial path in S I because its thickness h 1, the effective coefficients C * 55 , C * 54 , C * 45 and C * 44 are connected by the following relations:

C * 55 -iC * 54 = p - n j=1 p 1 -p I j V ⎛ ⎜ ⎝ 1 j
u (1)

dx 2 + i 1 j
u (1)

dx 1 ⎞ ⎟ ⎠ + p I j -p 2 V ⎛ ⎜ ⎝ 2 j u (2) dx 2 + i 2 j u (2) dx 1 ⎞ ⎟ ⎠ , ( 24
)
C * 45 -iC * 44 = -i p - n j=1 p 1 -p I j V ⎛ ⎜ ⎝ 1 j
u (1)

dx 2 + i 1 j
u (1)

dx 1 ⎞ ⎟ ⎠ + p I j -p 2 V ⎛ ⎜ ⎝ 2 j u (2) dx 2 + i 2 j u (2) dx 1 ⎞ ⎟ ⎠ , ( 25
)
where

p = p 1 V 1 + p 2 V 2 + n j=1 p I j V I j is the Voigt average, V I = n j=1 p I j V I j
is the volume fraction of the mesophase and V I j denotes the volume fraction of the region S I j . The volume fractions of matrix, fiber and the parallelogram periodic cell (V 1 , V 2 and V ) were defined previously. In [START_REF] Pobedrya | Mechanics of composite materials[END_REF], u (γ ) is the solution of the local problem 13 L, whereas in ( 25) u (γ ) is the solution of 23 L. Taking into account ( 18)-( 22), the analytical formulae for effective properties are obtained from [START_REF] Pobedrya | Mechanics of composite materials[END_REF] and [START_REF] Bakhvalov | Homogenisation: Averaging Processes in Periodic Media[END_REF] depending only on the unknown a 1 j as follows:

C * 55 -iC * 54 = p - n j=1 13 Q 1 j θ j+1 -θ j - n j=1 13 Q1 j e -i(θ j+1 +θ j ) sen θ j+1 -θ j -2 n j=1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ∞ p=3 13 Q pj e i( p-1) ( θ j+1 +θ j ) 2 p-1 sen ( p-1)(θ j+1 -θ j ) 2 + ∞ p=3 13 Q pj e -i( p+1) ( θ j+1 +θ j ) 2 p+1 sen ( p+1)(θ j+1 -θ j ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 26
)
C * 45 -iC * 44 = -i p - n j=1 23 Q 1 j θ j+1 -θ j - n j=1 23 Q1 j e -i(θ j+1 +θ j ) sen θ j+1 -θ j -2 n j=1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ∞ p=3 23 Q pj e i( p-1) ( θ j+1 +θ j ) 2 p-1 sen ( p-1)(θ j+1 -θ j ) 2 + ∞ p=3 23 Q pj e -i( p+1) ( θ j+1 +θ j ) 2 p+1 sen ( p+1)(θ j+1 -θ j ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 27
)
where

α3 Q pj = p 1 -p I j R 1 - p 2 -p I j κ 1 j R A pj χ pj + 1 ā j p -R 1 χ 1 j (δ 1α -iδ 2α ) δ 1 p 2V χ pj - p 2 -p I j 2 A 1 j V R 2 -R 2 1 R 1 R κ 2 -κ 1 j 2 R 2 (δ 1α -iδ 2α ) δ 1 p , A pj = κ 1 j + κ 2 R 2 p 1 + κ 1 j -κ 2 R 2 p 2R p 1 R p .
The coefficients ( 26) and ( 27) represent the effective axial coefficients of a reinforced three-phase composite with nonuniform mesophase S I . The analytical expressions of the effective properties ( 26) and ( 27) are functions of the properties and volume fractions of constituents of a composite with periodic cell w 1 , w 2 . As a limit case, the two-phase composite with perfect contact on interface reported in Guinovart-Díaz et al. [START_REF] Guinovart-Díaz | Influence of parallelogram cells in the axial behaviour of fibrous composite[END_REF] is obtained from ( 26) and ( 27) by taking the parameters κ 1 j = p I j / p 1 → ∞. The complete separation is obtained when κ 1 j = p I j / p 1 → 0.

3 Local problems α3 L (α = 1, 2) and effective properties for two-phase composites under nonuniform spring imperfect contact

The statement and the solution of the local problems as well as the effective coefficients for nonuniform imperfect spring model are calculated in an analogous way to the three-phase model reported in the previous section, and the principal details can be found in López-Realpozo et al (2014). Therefore, in this section we shall not repeat the same procedure and only some fundamental ideas will be explained, and the closed-form expressions of the effective coefficients will be written. Now, the family of boundary value problems depending on ε for the linear elastic heterogeneous media can be formulated in general form as follows:

∇ • C ε : ∇w ε = 0, x ∈ , ( 28
) [[T(w ε )]] = 0, x ∈ ε , (29) T (γ ) (w ε ) = (-1) γ +1 K[[w ε ]], x ∈ ε , (30) w ε (x) = Au 0 (x), x ∈ ∂ , (31)
where C ε (x) = C(y) is the S-periodic stiffness elastic tensor and y = x/ε is the local variable. Moreover, w ε = Au ε , A is a rotation matrix, the spring stiffness diagonal matrix is given by K, the displacement and the traction vectors are denoted by w and T, respectively, where

K = ⎛ ⎝ K n 0 0 0 K t 0 0 0 K s ⎞ ⎠ , w = ⎛ ⎝ w n w t w s ⎞ ⎠ , T = ⎛ ⎝ T n T t T s ⎞ ⎠ . ( 32
)
The corresponding tractions remain continuous [START_REF] Dumont | An asymptotic approach to the adhesion of thin stiff films[END_REF], while the effect of mechanical imperfection is incorporated through the mechanical displacements jumps across the interface (30). The spring constants K are the measures for the magnitude of the associated continuities, and they are a function of the position at the interface which is called proportional interface parameter [START_REF] Shodja | A piezoelectric medium containing a cylindrical inhomogeneity: Role of electric capacitors and mechanical imperfections[END_REF]. The vanishing value of Kn and Kt , Ks corresponds to pure debonding (normal perfect debonding), in-plane pure sliding and out-of-plane pure sliding, respectively. The status of the mechanical bonding is completely determined by appropriate values of these constants. For large enough values of the constants, the perfect bonding interface is achieved.; w n , w t , w s are the normal and the two tangential components of the rotated mechanical displacement vector w ε ; T n , T t , T s are the normal and tangential components of the traction vector T (T i = σ i j n j ), n is the outward unit normal on .

Here also we are interested to consider the antiplane case, and then, dimensionless local problems α3 L (α = 1, 2) related to ( 28)-(31) over the periodic cell S are written as follows: Laplace equation:

∇ 2 u (γ ) = 0, in S γ , ( 33
)
imperfect contact condition:

p γ (u (γ ) ,1 n 1 + u (γ ) ,2 n 2 ) + p γ n α p 1 = K s (θ ) R [[ u ]] , on , (34)
continuity condition for the stress:

p γ (u (γ ) ,1 n 1 + u (γ ) ,2 n 2 ) p 1 = -(1 -κ)n α , on , (35)
uniqueness condition:

u = 0, in S = S 1 ∪ S 2 , ( 36
)
p 11 -i p 21 = p 1 ⎡ ⎣ 1 - V 2 π R n-1 j=0 13 ā j+1 1 θ j+1 -θ j ⎤ ⎦ - p 1 V 2 π R n-1 j=0 ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∞ p=1 * e - i( p+1) ( θ j +θ j+1) 2 p+1 sen ( p+1)(θ j+1 -θ j ) 2 2 13 a j+1 p -(1-κ) δ 1 p R + 2 ∞ p=3 * 13 ā j+1 p e i( p-1) ( θ j+1 +θ j ) 2 p-1 sen ( p-1)(θ j+1 -θ j ) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ , ( 37
)
p 12 -i p 22 = -p 1 ⎡ ⎣ i + V 2 π R n-1 j=0 23 ā j+1 1 θ j+1 -θ j ⎤ ⎦ - p 1 V 2 π R n-1 j=0 ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∞ p=1 * e - i( p+1) ( θ j +θ j+1) 2 p+1 sen ( p+1)(θ j+1 -θ j ) 2 2 23 a j+1 p + i (1-κ) δ 1 p R + 2 ∞ p=3 * 23 ā j+1 p e i( p-1) ( θ j +θ j+1) 2 p-1 sen ( p-1)(θ j+1 -θ j ) 2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ . ( 38
)
The unknown constants a j p = α3 a j p are computed as solution of the following infinite system of algebraic equations for each arc:

ā j p + β j 1 R 2 H 1 δ 1 p ā j 1 + β j 1 R 2 H 2 δ 1 p a j 1 + β j p ∞ k=1 * w kp a j k = β j 1 δ 1 p (δ 1α -iδ 2α ) R, ( 39
)
where δ kp is the Kronecker delta and β j p = K j (1-κ)+κ p K j (1+κ)+κ p and K j is the tangential component of the imperfect parameter for each arc.

In the case where the area of mesophase S I is very thin, very good coincidence exists between the imperfect spring (37) and ( 38) and three-phase model ( 26) and [START_REF] Shodja | A piezoelectric medium containing a cylindrical inhomogeneity: Role of electric capacitors and mechanical imperfections[END_REF]. A condition between the interface and mesophase parameters is found. Thereby, the volume fraction of the interphase in the three-phase model is

V I = π(h 2 +2h R)/V ,
where h is the aforementioned thickness of the layer S I (see Fig. 2) and the relationship between the nondimensional debonding parameters K j and the interphase property κ 1 j can be introduced as in Hashin [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF], in the following form:

K j = R 1 κ 1 j / h. ( 40
)

Numerical algorithm based on a finite element method

In this part, the numerical method proposed to solve the problem (33)-( 36) is based on the finite element method. Even if this technique is quite standard, it is rapidly outlined here. For that purpose, it is first necessary to write the variational formulation of the interface problem (33)-( 36). First of all, it is convenient to make the change of variable in each phase ũ(γ

) (ξ ) = u (γ ) (ξ ) + ξ α , ( 41
) for α = 1, 2, ξ ∈ S. Then, ũ(γ ) is the solution of the problem ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∇ 2 ũ(γ ) = 0 in S, p γ p 1 ∂ ũ(γ ) ∂n = K S (θ ) R ũ in , p γ p 1 ∂ ũ(γ ) ∂n = 0 in , ũ = ξ α inS, ũ -ξ α S-periodic. (42)
Then, choosing a test function v, which can be discontinuous along the interface , multiplying the imperfect contact condition by the jump of this test function and integrating this relation among the interface, one obtains

v K S (θ ) R ũ ds - v p γ p 1 ∂ ũ(γ ) ∂n ds = 0. (43)
Then, on each sub-domain S γ , one can write

S γ p γ p 1 ∇v (γ ) • ∇ ũ(γ ) dξ + p γ p 1 ∂ ũ(γ ) ∂n v (γ ) ds = 0, (44)
and, adding the three last equalities (on , S 1 and S 2 ), a weak formulation of the problem can be written:

S γ p γ p 1 ∇ ũ • ∇vdξ + K S (θ ) R ũ v ds = 0, ũ = ξ α in S, ũ -ξ α S-periodic. (45)
Finally, using a standard finite element on each subdomain, and a "flat" finite element on , which has all its nodes on , the first ones related to S 1 and the other ones related to S 2 , it is possible to write a rigidity matrix of this problem that is invertible, with standard error estimates (see, for example, [START_REF] Nairn | Numerical implementation of imperfect interfaces[END_REF][START_REF] Dumont | An asymptotic approach to the adhesion of thin stiff films[END_REF] for more details). Due to finite element discretization, the integrals (see formulae ( 11)-( 13) for example) for the computation of the average stress and strains are replaced by sums over all elements contributions. In all the computations below, P2 triangular finite elements (conforming piecewise quadratic Lagrangian finite elements) are used. These elements are characterized by six nodes (one on each vertex and one on the middle of each edge). An example of meshed RVE is provided in Fig. 5 for fiber volume fraction 0.4 and fiber arrangement μ = 60 • .

Analysis of results

In order to illustrate the applicability of the theory developed above, several examples are chosen. As a validation of the model, we study some different particular cases of two-phase composites with hexagonal and square configurations of the periodic cell (Fig. 2). From now on, we shall use the classic short notation for the aforementioned effective elastic constants, i.e., p 11 = C * 55 , p 12 = C * 54 , p 21 = C * 45 and p 22 = C * 44 . The uniform imperfect elastic spring contact is obtained as a particular case of a stiffness matrix, taking

K 1 = K 2 = K 3 = • • • K j = K with 0 < K < ∞. Now,
some comparisons with other theoretical approaches are given in order to validate the present work.

Table 1 presents the normalized effective elastic modulus C * 44 /C [START_REF] Huang | Effects of void geometry on elastic properties of unidirectional fiber reinforced composites[END_REF] 44 , for perfect contact (K = 10 12 ) computed by the present model using AHM, the Finite Element Method proposed in this work and the results reported in Jiang et al. [START_REF] Jiang | A rigorous analytical method for doubly periodic cylindrical inclusions under longitudinal shear and its application[END_REF]. The material parameters used for this calculation are C [START_REF] Mansonwakeman | Composite processing and manufacturing: an overview[END_REF] 44 /C [START_REF] Huang | Effects of void geometry on elastic properties of unidirectional fiber reinforced composites[END_REF] 44 = 120 for Table 1 Effective elastic modulus obtained by Jiang et al. [START_REF] Jiang | A rigorous analytical method for doubly periodic cylindrical inclusions under longitudinal shear and its application[END_REF], AHM-uniform interface and FEM for composites with perfect contact, square and hexagonal periodic cell and ratio (C hexagonal and square periodic cell. Table 2 shows the behavior of the effective longitudinal shear elastic modulus C * 44 = C * 55 for the same arrangement of the periodic cell obtained by Jiang et al. [START_REF] Jiang | A rigorous analytical method for doubly periodic cylindrical inclusions under longitudinal shear and its application[END_REF], the present model by AHM for void volume fraction (K = 10 -12 ) and the Finite Element Method. The matrix shear modulus used is C

(2) 44 /C (1) 44 = 120) V 2 C * 44 /C (
(1) 44 = 30 GPa. In order to show the novelty of this work, we study now the nonuniform imperfect contact for two-phase composites with square cell. We consider the interface of the two-phase composites divided into two pieces = 1 ∪ 2 and different spring imperfect parameters K 1 and K 2 , which are assigned constants on 1 and 2 , respectively (see Fig. 6), according to the following distribution function

K s (θ ) = K 1 , 0 ≤ θ < θ 1 , K 2 , θ 1 ≤ θ < 2π.
Consequently, for illustrating the performance of the considered approach, a composite with hexagonal periodic cell and constituent properties C 3. Two different situations for soft interface are modeled with K 1 = 10 5 , K 2 = 0; K 1 = 0, K 2 = 10 5 using the AHM (spring and three-phase model) and FEM. Good agreement can be appreciated between the three methods where the relative errors are always smaller than 0.03 %. 3 Effective elastic moduli obtained by AHM and FEM approximations for nonuniform imperfect contact and hexagonal periodic cell (μ = 60 • ). Ratio between the constituent properties is C (2) 44 /C (1) 44 = 0.0002168

V 2 C * 44 /C
(1)

44 (K 1 = 10 5 , K 2 = 0) C * 44 /C (1)
44 (K 1 = 0, K 2 = 10 44 = 120. Good agreement between the two models is observed, where the difference is attained for low fiber volume fraction.

Figure 8a, b illustrates a comparison for the overall normalize properties C * 55 and C * 45 between three-phase and spring models, both for uniform and for nonuniform interface computed by FEM and AHM approaches. The composite has square periodic cell, and the interface for the nonuniform case is partitioned into two arcs with length π/2 and 3π/2, respectively. Hence, the three-phase composite for the nonuniform analysis has the normalized shear interphase property in each element S I 1 and S I 2 , κ 11 = p I 1 / p 1 and κ 12 = κ 2 = p 2 / p 1 , respectively. Moreover, the effective properties are considered as a function of κ 11 which varies in the range (0, 10 6 ), whereas κ 12 = κ 2 = p 2 / p 1 is constant in the whole range, and the contrast between the fiber and matrix is κ 2 = 120 (Fig. 8a) and κ 2 = 6 (Fig. 8b). Consequently, the spring model exhibits a similar pattern and the spring parameters K 1 , K 2 assigned to each arc are different. The parameters of imperfection used for the nonuniform analysis (FEM and AHM) are K j = κ 1 j p 1 R 1 / h. The composite with uniform interface means that the interphase S I is homogeneous. Accordingly, in the case of the uniform spring model, the imperfect constant is K 1 = K 2 = κ 11 p 1 R 1 / h. The fiber and interphase volume fractions are V 2 = 0.4 and V I = 0.001 in the computation. Notice that the three-phase and spring models coincide in the whole range of the contrast of the stiffness parameter κ 11 . Furthermore, when the parameter κ 11 → +∞, all models describe the effective properties for a perfect contact of a bi-phase fibrous composite. In the case κ 11 → 0, the difference between the uniform and nonuniform models is remarkable. In this case, the effective property is weaker in comparison with perfect contact (κ 11 → +∞) and it is even more weaker when the stiffness κ 2 = p 2 / p 1 is low. Hence, the spring uniform model as κ 11 → 0 reflects a total disbond between the fiber and matrix. These models predict the lower and upper bounds for the effective properties in the range of the κ 11 . The two approaches (FEM and AHM) are consistent in the numerical analysis. The differences between them for small value of κ 11 are due to the singularities that appear around the zone where the rigidity of the interface/interphase changes (θ = 0 • and θ = 90 • ) (these singularities appear for large value of rigidity). It is well known that the global axial behavior in fibrous composites with periodic square cell, isotropic components along the axis x 3 under perfect contact condition, is transversely isotropic; however, due to the nonuniform imperfection at the interface, the global behavior of this composite changes and it is a monoclinic symmetry so that the effective coefficient C * 45 is different from zero for the antiplane problem (see Fig. 8b).

A comparison between the three-phase and spring model, both for non-uniform case, is shown in

Conclusions

An approach using the asymptotic homogenization method for simulating the nonuniform imperfect bonding in elastic composite materials with oblique angle of the cell is proposed. We introduce different elastic mesophase constants for each partition of the interphase S I , between the matrix and inclusions, where a piecewise constant function p I (θ ) = p I j , θ j < θ < θ j+1 transmits a load from the matrix to the inclusion. The computation of the analytical expressions ( 26) and ( 27) of the shear moduli for three-phase model with nonuniform mesophase, derived by the asymptotic homogenization approach, is realized with good accuracy. The exactitudes of the solutions ( 26) and ( 27) were tested with other theoretical models (AHM spring model (37) and (38) and FEM), and good agreement was observed. The nonuniform imperfect interface and the angle of inclination of the cell play an important role in the global behavior of the composite. When the thickness of the interphase is very small (h → 0), the three-phase composites behave like two-phase composites. These models can be extended to more complex composites in order to predict the overall properties in the design and manufacturing of one-directional fibrous composite materials.

Fig. 1 Fig. 2

 12 Fig. 1 Heterogeneous medium with periodic cell

Fig. 4

 4 Fig. 4 Partition of the interface in n arcs for the spring model

Fig. 5

 5 Fig. 5 Meshed RVE for 0.4 fiber volume fraction and μ = 60 • fiber arrangement

  44 = 0.0002168 is taken. For the distribution of the imperfect parameters by arcs, we use θ 1 = π/2 in the calculations of Table

Fig. 6

 6 Fig. 6 Two-phase composite with interface divided by two arcs Table3Effective elastic moduli obtained by AHM and FEM approximations for nonuniform imperfect contact and hexagonal periodic cell (μ = 60 • ). Ratio between the constituent properties is C(2) 44 /C(1) 44 = 0.0002168

Fig. 7 Figure 7

 77 Fig. 7 Comparison of the axial effective properties computed by AHM and FEM where K 2 is fixed and the other parameter K 1 changes

Fig. 8

 8 Fig. 8 Different models are compared for the effective properties C * 55 and C * 45 in a composite with square periodic cell, uniform/nonuniform imperfect contacts and two different values of contrast between the fiber and matrix κ 2 = 120 (a) and κ 2 = 6 (b)

Fig. 9 .

 9 This figure shows the effective normalized properties C * 44 , C * 55 , C * 45 = C * 54 varying the length of the arcs of the

Fig. 9 Fig. 10

 910 Fig.[START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF] Behavior of the effective properties by AHM spring and three-phase models for a composite with square cell and nonuniform imperfect interface

Table 2

 2 Variations of effective longitudinal shear elastic modulus obtained by Jiang et al.[START_REF] Jiang | A rigorous analytical method for doubly periodic cylindrical inclusions under longitudinal shear and its application[END_REF], AHM and FEM for a composite with uniform imperfect contact with empty fibers and C

		1)
		44
		μ = 60 •			μ = 90 •
		Jiang et al.	AHM	FEM	Jiang et al.	AHM	FEM
	0.1	1.21815	1.21815	1,21815	1.21816	1.21816	1,21816
	0.2	1.48971	1.48971	1,48971	1.49000	1.49000	1,49000
	0.3	1.83711	1.83711	1,83711	1.83990	1.83990	1,83990
	0.4	2.29764	2.29764	2,29764	2.31343	2.31343	2,31343
	0.5	2.93931	2.93931	2,93931	3.00863	3.00863	3,00863
	0.6	3.90319	3.90319	3,90319	4.18637	4.18637	4,18638
				(1) 44 = 30 GPa
	V 2	C * 44 ≡ C * 55
		Hexagonal cell (μ = 60 • )		Square cell (μ = 90 • )
		Jiang et al.	AHM	FEM	Jiang et al.	AHM	FEM
	0.1	24.54550	24.54545	24.54546	24.54530	24.54530	24.54531
	0.2	20.00000	19.99996	20.00000	19.99590	19.99592	19.99592
	0.3	16.15330	16.15326	16.15326	16.12740	16.12740	16.12740
	0.4	12.85340	12.85336	12.85336	12.76050	12.76054	12.76054
	0.5	9.98430	9.98427	9.98427	9.73960	9.73964	9.73964
	0.6	7.45030	7.45029	7.45029	6.90960	6.90954	6.90954

Acknowledgments

The funding of CONACYT project number 129658 and Coordinación de la Investigación Cientifica (COIC-STIA-312-14, 313-14) de la UNAM is gratefully acknowledged. The authors gratefully acknowledge to the project SHICHAN, supported by FSP (Cooperation Scientifique Franco-Cubaine) Cuba 2011-26, PROJET No 29935XH and to the project Composite Materials from University of Havana. Thanks to Departamento de Matemáticas y Mecánica, IIMAS-UNAM, for its support and Ramiro Chávez Tovar and Ana Pérez Arteaga for computational assistance.

imperfection in the composite where the parameters of imperfections K 1 = 1, K = 50 are considered. The material parameters used in the computation are κ 2 = C (2) 44 /C (1) 44 = 120, the fiber volume fraction is V 2 = 0.3 and the interphase volume fraction is V I = 0.001. Notice that the effective normalized properties C * 44 and C * 55 are decreasing and different, except when the amplitudes of the arc are 90 • , 180 • , 270 • and 360 • . In Fig. 10, the effective axial properties of a composite with rhombic cell versus the volume fraction are studied. The angle of the periodic cell is μ = 75 • , and the contrast between the fiber and matrix used in the calculus is κ 2 = p 2 / p 1 = 120. Some different cases are presented via AHM: first, a composite with perfect contact conditions at the interface, denoted by black dotted line, is considered. The coefficients C * 44 and C * 55 possess the biggest rigidity because the fiber is much more rigid than the matrix. When in Fig. 10a, b the imperfect uniform parameter K = 50 (blue circles and dashed line) is included, the present models indicate that the composite begins to be weaker. The nonuniform interphase considered possesses two arcs of circumference with different length, with angles θ 1 = 90 • , θ 2 = 360 • -θ 1 and the corresponding imperfect parameters K 1 and K 2 , respectively (see Fig. 6). The coefficients C * 44 and C * 55 when K 1 < K 2 (black and red circles) are weaker than the uniform imperfection and reciprocally in the case K 2 < K 1 (black and red square) the fragility increases. The curve corresponding to the imperfect uniform parameter K = 0.1 (blue square and dashed line) approximates to the total separation of materials or porous fibers. One remarkable difference of the composite with rhombic periodic cell in comparison with hexagonal or square cell is that it exhibits monoclinic symmetry, i.e., C * 44 = C * 55 , C * 45 = 0 and C * 54 = 0. Moreover, in Fig. 10c, d it can be observed that the absolute value of C * 45 and C * 54 increases as the nonuniform interphase is considered. The total coincidence between three-phase and spring models can be obtained when the interphase parameter η = h/R 2 = 10 -5 [START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF]. For major values of η, some differences appear between the used models, because the thickness h is increased and the mesophase S I has more influence in the properties of the composite. Besides, the validation with FEM indicates the accuracy of the theoretical models calculated by AHM.


Raffaella Rizzoni
email: raffaella.rizzoni@unife.it

Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Imperfect interfaces as asymptotic models of thin curved elastic adhesive interphases

We obtain a limit model for a thin curved anisotropic interphase adherent to two elastic media. Our method is based on asymptotic expansions and energy minimization procedures. The model of perfect interface is obtained at the first order, while an imperfect interface model is obtained at the next order. The conditions of imperfect contact, given in a parallel orthogonal curvilinear coordinate system, involve the interphase material properties, the first order displacement and traction vectors, and their derivatives. An example of implementation of the imperfect interface condition is given for a composite sphere assemblage.

Introduction

During the last few years gluing techniques have become usual in engineering structural assembly and they are now frequent in aeronautics industry, where the use of composite material is necessary to lighten structures. These techniques have to be taken into account in structural modeling, in order to develop implementable predictive models in computational structural analysis softwares.

The thickness of the glue is usually very small when compared with the characteristic dimensions of the structure. A classical modeling approach is to consider the thickness as a small parameter and to study the limit problem when this parameter tends to zero. Traditionally, the stiffness of the glue is assumed to be smaller than those of the adherents. Many authors focused on the assumption of soft adhesively bonded joining, with fundamental contributions given in [START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF] and [START_REF] Klarbring | Asymptotic modelling of adhesive joints[END_REF] and in works addressing the assembly of elastic plates [START_REF] Geymonat | Analyse asymptotique du comportement en flexion de deux plaques collées[END_REF][START_REF] Zaittouni | Etude théorique et numérique du comportement d'un assemblage de plaques[END_REF], the limit behavior of soft thin interphases [START_REF] Licht | Comportement asymptotique d'une bande dissipative mince de faible rigidité[END_REF][START_REF] Licht | Comportement asymptotique d'une bande dissipative mince de faible rigidité[END_REF]Michaille, 1996, 1997;[START_REF] Ould-Khaoua | Thin layers in elasticity: a theoretical and numerical study[END_REF][START_REF] Ganghoffer | Modelling of the mechanical behaviour of joints bonded by a nonlinear incompressible elastic adhesive[END_REF][START_REF] Lebon | Numerical study of soft adhesively bonded joints in finite elasticity[END_REF]Lebon et al., , 2004;;[START_REF] Geymonat | Mathematical analysis of a bonded joint with a soft thin adhesive[END_REF][START_REF] Lenci | A limit model in the analysis of steel-concrete beams[END_REF][START_REF] Lebon | Asymptotic analysis of Mohr-Coulomb and Drucker-Prager soft thin layers[END_REF][START_REF] Lebon | Asymptotic study of soft thin layer: the non convex case[END_REF][START_REF] Cognard | Analysis of the nonlinear behavior of adhesives in bonded assemblies -comparison of TAST and Arcan tests[END_REF][START_REF] Schmidt | Modelling of adhesively bonded joints by an asymptotic method[END_REF][START_REF] Pelissou | Asymptotic modeling of quasi-brittle interfaces[END_REF][START_REF] Rekik | Homogenization methods for interface modeling in damaged masonry[END_REF][START_REF] Monchiet | Interfacial models in viscoplastic composites materials[END_REF], and, more recently, the dynamics of laminated beams [START_REF] Serpilli | Asymptotic modelling of the linear dynamics of laminated beams[END_REF].

On the other hand, the ratio between the stiffness of the glue and the stiffness of the adherents can be larger than 1/50, about 1/35 for two aluminum plates joined by an epoxy glue. Thus, it is interesting to consider the case in which the elastic moduli of the glue and the adherents are comparable. This case has been studied in [START_REF] Caillerie | The effect of a thin inclusion of high rigidity in an elastic body[END_REF] and [START_REF] Abdelmoula | Comportement asymptotique d'une interphase élastique mince[END_REF] and in some recent works focusing on interphase modeling in elasticity [START_REF] Benveniste | Imperfect soft and stiff interfaces in two-dimensional elasticity[END_REF][START_REF] Hashin | Thin interphase/imperfect interface in elasticity with application to coated fiber composites[END_REF][START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF][START_REF] Lebon | First order numerical analysis of linear thin layers[END_REF][START_REF] Serpilli | Limit models in the analysis of three different layered elastic strips[END_REF][START_REF] Benveniste | On two models of arbitrarily curved three-dimensional thin interphases in elasticity[END_REF]Lebon andRizzoni, 2010, 2011;[START_REF] Lebon | Asymptotic modelling of interface taking into account contact conditions: asymptotic expansions and numerical implementation[END_REF][START_REF] Rizzoni | Asymptotic analysis of an adhesive joint with mismatch strain[END_REF] and on structural interfaces [START_REF] Bigoni | Statics and dynamics of structural interfaces in elasticity[END_REF]Bertoldi et al., 2007a,b).

In particular, in [START_REF] Lebon | Asymptotic behavior of a hard thin linear elastic interphase: an energy approach[END_REF] a two-level model of a flat imperfect interface was derived using an energy asymptotic method. At the first level, the adhesive is replaced by a model of perfect interface, for which the stress and the displacement vectors are continuous. At the second level, the jumps in the displacements and in the stress vector along the interface are related to the interphase material properties and to the displacement and stress fields obtained at the first level.

In the present paper, we extend the results obtained in [START_REF] Lebon | Asymptotic behavior of a hard thin linear elastic interphase: an energy approach[END_REF] to the case of a curved, thin, linear elastic, anisotropic and homogeneous interphase joining two elastic adherents. In Section 2, we introduce the general notations and the three-dimensional equilibrium problem posed in a variational form. In Section 3, we obtain two different types of interface condition, by using a method based on asymptotic expansions and energy minimization procedures. These conditions extend to the curvilinear case the interface conditions obtained in [START_REF] Lebon | Asymptotic behavior of a hard thin linear elastic interphase: an energy approach[END_REF] for a straight interface.

As found in [START_REF] Lebon | Asymptotic behavior of a hard thin linear elastic interphase: an energy approach[END_REF], the first order term of the expansion yields the model of perfect interface. The second order term yields an imperfect interface model simulating the presence of the interphase. We interpret this term as a correction of the leading solution corresponding to the perfect interface model. The conditions of imperfect contact, obtained in a parallel orthogonal curvilinear coordinate system for a general anisotropic homogeneous interphase material, are specialized in Section 4 to the case of isotropic interphase. Section 5 is devoted to the implementation of the imperfect interface conditions and an example of a composite spheres assemblage is given in Section 6.

Formulation of the three-dimensional equilibrium problem

Let us consider a thin curved adhesive interphase B ε ∈ R 3 of constant thickness ε bonding two adherents occupying the regions ε + , ε -∈ R 3 . Two parallel interfaces, denoted S ε + and S ε -, separate ε + and ε -from B ε , respectively. Let S denote the middle surface between S ε + and S ε -. These interfaces are assumed to be perfect with the usual assumption that the displacement vector, u ε , and the stress vector, ε n, are continuous on S ε ± . Our goal is to substitute the thin interphase B ε with the surface S, geometric limit of the interphase as ε → 0 + , and to determine the imperfect interface conditions on the displacement and tractions across S which are equivalent to the three phase configuration with perfect interface conditions. The three regions are assumed to be anisotropic, homogeneous and linear elastic. We take a ± to denote the elasticity tensors of the materials occupying the regions ε ± , and b the elasticity tensor of the interphase material. The elasticity tensors are assumed to satisfy the following assumptions:

a ± ∈ L ∞ ( ε + ∪ ε -), b ∈ L ∞ (B ε ), ∃Á ± , Á > 0 : a ± (e) • (e) ≥ Á ± |e| 2 , b(e) • (e) ≥ Á|e| 2 , ∀ e : e = e T , (1)
and the standard assumptions of major and minor symmetries. After taking e to denote the strain tensor

e(u ε ) = 1 2 (∇u ε + (∇u ε ) T ), (2)
linear elasticity gives the Cauchy stress tensor ε as follows:

ε = a ± (e) in ε ± , (3)
ε = b(e) in B ε . (4)
A body force density f ∈ L 2 ( ε + ∪ ε -) 3 is assumed to be applied to ( ε + ∪ ε -) and a surface force density g to ε g ⊂ (∂ ε

+ ∪ ∂ ε -), with g ∈ L 2 ( ε g ) 3 .
Homogeneous boundary conditions are prescribed on ε u := (∂ ε

+ ∪ ∂ ε -) \ ε g : u ε = 0 on ε u . (5)
The equilibrium configurations of the three-phase composite body are the minimizers of the total energy

E ε (u) = ε ± 1 2 a ± (e(u ε )) • e(u ε ) -f • u ε dV - ε g g • u ε dA + B ε 1 2 b(e(u ε )) • e(u ε ) dV (6)
in the space of kinematically admissible displacements

V ε = {u ∈ H( ε ; R 3 ) : u = 0 on ε u }, (7)
where H( ε ; R 3 ) is the space of the vector-valued functions on the set

ε := ( ε + ∪ ε -∪ B ε ∪ S ε + ∪ S ε -)
, which are continuous and differentiable as many times as necessary. In view of the above regularity assumptions on a ± , b, f, and g, the existence of a unique minimizer u ε in V ε is ensured (Ciarlet, 1988, Theorem 6.3).

Thin interphase asymptotic analysis

In this section, we follow the asymptotic method based on the energy minimization and introduced in [START_REF] Lebon | Asymptotic behavior of a hard thin linear elastic interphase: an energy approach[END_REF]. The first step of the method consists in reformulating the equilibrium problem for a rescaled interphase domain made independent of ε via a change of variables. To do so, it is convenient to introduce a curvilinear orthogonal system ( 1 , 2 , 3 ) in ε , where 1 and 2 are two coordinate lines on the surface S and 3 is the coordinate line along the normal to the surface (Fig. 1). We assume that 3 takes the values 30ε/2, 30 , and 30 + ε/2, on S ε -, S, and S ε + , respectively. Let (ê 1 , ê2 , ê3 ) denote the local triad of the curvilinear coordinate system. The metric coefficients of this orthogonal parallel curvilinear coordinate system are denoted by h ε 1 , h ε 2 , h ε 3 , with h ε 3 = 1. The displacement gradient in the orthogonal parallel curvilinear system can be written in the form:

∇u ε = (v ε 1 |v ε 2 |u ε ,3 ) (8)
where a comma is used to denote partial differentiation with respect to 1 , 2 , 3 , the notation (a|b|c) indicates the matrix whose columns are the vectors a, b, c, and with

v ε ˛:= u ε ,˛( h ε ˛)-1 + Á ε ˛(u ε ), ˛= 1, 2, (9)
Á ε 1 := h ε 1,2 h ε 1 h ε 2 (ê 1 ⊗ ê2 -ê2 ⊗ ê1 ) + h ε 1,3 h ε 1 (ê 1 ⊗ ê3 -ê3 ⊗ ê1 ), ( 10
)
Á ε 2 := h ε 2,1 h ε 1 h ε 2 (ê 2 ⊗ ê1 -ê1 ⊗ ê2 ) + h ε 2,3 h ε 2 (ê 2 ⊗ ê3 -ê3 ⊗ ê2 ). ( 11
)
We now introduce the following change of variables:

( ˆ , 3 ) := p( ˆ , 3 ) := ( ˆ , 30 + ε -1 ( 3 -30 )) in B ε , ( 12
) ( ˆ , ˜ 3 ) = p( ˆ , 3 ) := ( ˆ , 3 ) ± ε 2 ∓ 1 2 ê3 in ε ± , (13)
where ˆ := ( 1 , 2 ) and ˆ := ( 1 , 2 ) are the surface coordinates. We take ± , B, S ± , g and u to denote the domains corresponding to ε ± , B ε , S ε ± , ε g and ε u , respectively, after the change of variables (Fig. 1). Let f := f • p-1 and g := g • p-1 denote the rescaled external forces. Let ũε ± = u ε • p-1 denote the displacement vector field from the adherents adjacent to the rescaled interphase and let u

ε = u ε • p -1
denote the displacement vector field from the rescaled interphase. In view of the condition of perfect interfaces between the adherents and the rescaled interphase, we have that ũε

± = u ε on S ± . Finally, let h ε i := h ε i • p -1 and hε i := h ε i • p-1 , i = 1, 2, 3
, denote the rescaled metric coefficients of the curvilinear coordinates system. Note that h 3 = 1 = h3 .

Within these assumptions and in view of relations ( 8)-( 11), the rescaled energy of the composite takes the following form:

E ε ( ũε ± , u ε ) : = ± 1 2 a ± (e( ũε ± )) • e( ũε ± ) -f • ũε ± hε 1 hε 2 dV ˜ - g (g • ũε ± ) hε s dA ˜ + B 1 2 1 ε B 33 (u ε ,3 ) • (u ε ,3 ) + 2B ˛3(v ε ˛) • (u ε ,3 ) + εB ˛ˇ(v ε ˛) • (v ε ˇ) h ε 1 h ε 2 dV , (14)
where summation convention for repeated indices is used, the indices ˛, ˇtake the values 1, 2, hε s is the Jacobian of the change of variables ( 12) and ( 13) in g , and B ij , i, j = 1, 2, 3 are the matrices whose components are defined by the relations

B ij hk := b hikj . ( 15
)
Note that the symmetry properties of b imply that (B ij ) T = B ji . The equilibrium problem of the rescaled three phase configuration can now be formulated as follows: find the pair ( ũε ± , u ε ) minimizing the energy (14) in the set of displacements

V = {( ũ± , u) ∈ H( ± ; R 3 ) × H(B; R 3 ) : ũ± = u on S ± , ũ± = 0 on u }. ( 16
)
The asymptotic method that we use in this paper is based upon the existence of the following expansions for the rescaled displacement vector fields ũε ± , u ε , for the rescaled metric coefficients, hε i , h ε i , i = 1, 2, 3, and for the Jacobian hε s :

ũε ± = ũ0 ± + ε ũ1 ± + ε 2 ũ2 ± + o(ε 2 ), (17)
u ε = u 0 + εu 1 + ε 2 u 2 + o(ε 2 ), ( 18
)
hε i = h0 i + ε h1 i + ε 2 h2 i + o(ε 2 ), i = 1, 2, 3, (19)
h ε i = h 0 i + εh 1 i + ε 2 h 2 i + o(ε 2 ), i = 1, 2, 3, (20)
hε s = h0 s + ε h1 s + ε 2 h2 s + o(ε 2 ). ( 21
)
Substituting these expansions into the rescaled energy ( 14), we obtain

E ε ( ũ± , u) = 1 ε E -1 (u 0 ) + E 0 ( ũ0 ± , u 0 , u 1 ) + εE 1 ( ũ0 ± , ũ1 ± , u 0 , u 1 , u 2 ) + ε 2 E 2 ( ũ0 ± , ũ1 ± , ũ2 ± , u 0 , u 1 , u 2 , u 3 ) + o(ε 2 ), (22)
where

E -1 (u 0 ) := B w -1 H 0 dV , ( 23
) E 0 ( ũ0 ± , u 0 , u 1 ) := ± w0 H0 dV ˜ - g w0 s h0 s dA ˜ + B (w -1 H 1 + w 0 H 0 )dV , ( 24
) E 1 ( ũ0 ± , ũ1 ± , u 0 , u 1 , u 2 ) := ± ( w0 H1 + w1 H0 )dV ˜ - g ( w0 s h1 s + w1 s h0 s )dA ˜ + B (w -1 H 2 + w 0 H 1 + w 1 H 0 )dV , ( 25
)
E 2 ( ũ0 ± , ũ1 ± , ũ2 ± , u 0 , u 1 , u 2 , u 3 ) := ± ( w0 H2 + w1 H1 + w2 H0 )dV ˜ - g ( w0 s h2 s + w1 s h1 s + w2 s h0 s )dA ˜ + B (w -1 H 3 + w 0 H 2 + w 1 H 1 + w 2 H 0 )dV , (26)
being w0 ( ũ0

± ) := 1 2 a ± (e( ũ0 ± )) • e( ũ0 ± ) -f • ũ0 ± , (27)
w1 ( ũ0

± , ũ1 ± ) := a ± (e( ũ0 ± )) • e( ũ1 ± ) -f • ũ1 ± , (28)
w2 ( ũ0

± , ũ1 ± , ũ2 ± ) := 1 2 a ± (e( ũ1 ± )) • e( ũ1 ± ) + a ± (e( ũ0 ± )) • e( ũ2 ± ) -f • ũ2 ± , (29)
wi s ( ũi

± ) := -g • ũi ± , i = 0, 1, 2, (30)
w -1 (u 0 ) := 1 2 B 33 (u 0 ,3 ) • (u 0 ,3 ) (31) w 0 (u 0 , u 1 ) := B 33 (u 0 ,3 ) • (u 1 ,3 ) + 1 2 B ˛3(v 0 ˛) • (u 0 ,3 ), (32)
w 1 (u 0 , u 1 , u 2 ) := B 33 (u 0 ,3 ) • (u 2 ,3 ) + 1 2 B 33 (u 1 ,3 ) • (u 1 ,3 ) + 1 2 B ˛3(v 0 ˛) • (u 1 ,3 ) + 1 2 B ˛3(v 1 ˛) • (u 0 ,3 ) + 1 2 B ˛ˇ(v 0 ˛) • (v 0 ˇ), (33)
w 2 (u 0 , u 1 , u 2 , u 3 ) := B 33 (u 0 ,3 ) • (u 3 ,3 ) + B 33 (u 1 ,3 ) • (u 2 ,3 ) + 1 2 B ˛3(v 0 ˛) • (u 2 ,3 ) + 1 2 B ˛3(v 1
˛) • (u 1 ,3 ) + 1 2 B ˛3(v 2
˛) • (u 0 ,3 ) + B ˛ˇ(v 0
˛) • (v 1 ˇ). (34)
In ( 23)-( 26) the following quantities were also defined:

H0 := h0 1 h0 2 , (35)
H1 := h0 1 h1 2 + h1 1 h0 2 , ( 36
) H2 := h0 1 h2 2 + h1 1 h2 2 + h2 1 h1 2 + h3 1 h0 2 , (37)
and analogous definitions were introduced for the quantities H i , i = 0, 1, 2.

In the next subsections, we obtain the conditions satisfied by the stationary points of the energies E -1 , E 0 , E 1 and E 2 . These conditions identify the relations between u 0 , u 1 , u 2 , . . . and the corresponding stress vectors arising from ũ0 ± , ũ1 ± , ũ2 ± , . . . at the interfaces S + , S -.

Stationary points of E -1

We minimize the energy E -1 in the class of displacements u 0 ∈ H(B; R 3 ). Since b is a positive definite tensor, the matrix B 33 is also positive definite and thus the energy E -1 is non negative. The stationarity of E -1 with respect to u 0 gives the condition

u 0 ,3 = 0, a.e. in B (38)
i.e., minimizers are independent of the coordinate 3 along the interphase thickness. Based on this result and on the initial assumption of perfect contact between the rescaled interphase and adherents, we obtain the following condition: ũ0 ˆ , 30 + 1 2

+ = ũ0 ˆ , 30 - 1 2 - , ˆ ∈ S. (39)
If we introduce the following definition of "jump" of a function f : + ∪ -→ R 3 across the rescaled interphase

[f]( ˆ ) := f ˆ , 30 + 1 2 + -f ˆ , 30 - 1 2 - , ˆ ∈ S, (40)
then condition (39) can be rephrased as follows:

[ ũ0 ] = 0 on S.

(41)

Stationary points of E 0

In view of (38), the energy E 0 simplifies as follows:

E 0 = ± 1 2 a ± (e( ũ0 ± )) • e( ũ0 ± ) -f • ũ0 ± H0 dV ˜ - g (g • ũ0 ± ) h0 s dA ˜ , (42)
and it becomes independent of u 0 and u 1 . We seek the energy minimizer of (42) in the class of displacements

V 0 = ( ũ± ) ∈ H( ± ; R 3 ) : ũ+ ˆ , 30 + 1 2 + = ũ- ˆ , 30 - 1 2 - , ˆ ∈ S, ũ± = 0 on u . (43)
Using standard arguments, we obtain that stationary points of the energy (42) satisfy the following equilibrium equations: div(a ± (e( ũ0 ± ))

+ f = 0 in ± , (44)
a ± (e( ũ0 ± ))n = g on g , (45)
a ± (e( ũ0

± ))n = 0 on ∂ ± \ g , (46)
[ ˜ 0 ê3 ] = 0 on S.

From the mechanical viewpoint, relations (41) and ( 47) correspond to conditions of perfect interface for the rescaled interphase modeling.

Stationary points of E 1

The relations ( 38), ( 44)-( 47) allow to simplify E 1 as follows:

E 1 = B 1 2 B 33 (u 1 ,3 ) • (u 1 ,3 ) + 1 2 B ˛3(v 0
˛) • (u 1 ,3 ) H 0 dV . ( 48
)
Stationary points u 1 of E 1 in H(B ; R 3 ) satisfy the condition

˜ 0 ê3 = B 33 (u 1 ,3 ) + B ˛3(v 0 ˛), (49)
with ˜ 0 ê3 the common value of the tractions on S ± . We integrate relation (49) along 3 taking into account that the fields at the order zero are independent of the coordinate 3 , and we use the continuity of the rescaled displacement vector fields across S ± , implying that v 0 ˛= ṽ0 ǫn S, to obtain

[ ũ1 ] = (B 33 ) -1 ( ˜ 0 ê3 -B ˛3(ṽ 0 ˛))
on S.

(50)

Stationary points of E 2

Using ( 38) and ( 44)-( 47), we simplify E 2 as follows:

E 2 = ± 1 2 a ± (e( ũ1 ± )) • e( ũ1 ± ) H0 dV ˜ + B (B 3˛( u 1 ,3 ) + B ˇ˛(v 0 ˇ)) • (v 1 ˛)H 0 dV . (51)
Noting the relation

v 1 ˛= u 1 ,h 0 ˛+ Á 0 ˛(u 1 ) - h 1 (h 0 ˛)2 u 0 ,˛+ Á 0 ˛(u 0 ), (52)
we observe that the second integral in ( 51) is given by

B 1 2 B 3˛( u 1 ,3 ) + B ˇ˛(v 0 ˇ) • u 1 ,h 0 ˛+ Á 0 ˛(u 1 ) H 0 dV . ( 53
)
up to (constant) terms in u 0 and its first derivatives. In view of (49) and under suitable regularity assumptions, the field u 1 admits the following representation

u 1 ( ˆ , 3 ) = [ ũ1 ] 3 + w( ũ1 )( ˆ ), (54)
where w(f

) := (1/2)(f( 30 + (1/2)) + ) + f( 30 -(1/2)) -
)) for a generic function f : + ∪ -→ R 3 . Integrating (54) along 3 simplifies the term linear in 3 and gives

E 2 = ± 1 2 a ± (e( ũ1 ± )) • e( ũ1 ± ) H0 dV ˜ + S (B 3˛( u 1 ,3 ) + B ˇ˛(v 0 ˇ)) • 1 h 0 ˛( w( ũ1 )) ,˛+ Á 0 ˛( w( ũ1 ))H 0 dV , (55)
up to in u 0 , which are considered constant in the minimization procedure. To "complete" the gradient of w( ũ1 ), we add to the second integral the term (( ˜ 0 ê3 ) • ( w( ũ1 )) ,3 ), which vanishes being w( ũ1 ) independent of 3 . We obtain

E 2 = ± 1 2 a ± (e( ũ1 ± )) • e( ũ1 ± ) H0 dV ˜ + S (M • ∇ w( ũ1 ))H 0 dV , (56)
where

M := (B 31 (u 1 ,3 ) + B ˇ1(v 0 ˇ)|B 32 (u 1 ,3 ) + B ˇ2(v 0 ˇ)| ˜ 0 ê3 ). ( 57
)
Finally, using standard arguments, we find that stationary points of ( 56) satisfy the following equilibrium equations:

div(a ± (e( ũ1 ± ))) = 0 in ± , (58)
a ± (e( ũ1

± ))n = 0 on ∂ ± \ S ± , (59)
[

˜ 1 ê3 ] = -div M on S, (60)
Mn = 0 on ∂S.

(61)

Relations ( 50) and ( 60) are imperfect interface conditions modeling the mechanical behavior of the anisotropic interphase at the first order of the asymptotic expansion. As a final remark, we note that condition (61) implies that the asymptotic expansions ( 17) and ( 18) do not hold in a neighborhood of ∂S (Lebon and Rizzoni, 2010, Section 3).

Imperfect interface conditions arising from an isotropic and homogeneous interphase

In this section the conditions ( 50) and ( 60) obtained in Section 3 are specialized for the case of an isotropic and homogeneous interphase with Lamé coefficients and . In view of the definition (15), the matrices B ij are given by

B ii = (2 + ) êi ⊗ êi + (ê j ⊗ êj + êk ⊗ êk ), i / = j / = k (62) B ij = êi ⊗ êj + êj ⊗ êi , i / = j. ( 63
)
Substituting these relations into (50) gives the following expressions for the jumps in the displacement components:

[ ũ1 1 ] = 1 ˜ 0 13 - 1 h 0 1 ũ0 3,1 + h 0 1,3 h 0 1 ũ0 1 , ( 64
) [ ũ1 2 ] = 1 ˜ 0 23 - 1 h 0 2 ũ0 3,2 + h 0 2,3 h 0 2 ũ0 2 , ( 65
) [ ũ1 3 ] = 1 (2 + ) ˜ 0 33 -(2 + ) 1 h 0 1 ũ0 1,1 + 1 h 0 2 ũ0 2,2 + h 0 1,2 h 0 1 h 0 2 ũ0 2 + h 0 2,1 h 0 1 h 0 2 ũ0 1 + h 0 1,3 h 0 1 + h 0 2,3 h 0 2 ũ0 3 . ( 66
)
Substituting ( 62) and ( 63) into (60) and using the expression of the divergence of a tensor given in [START_REF] Bolton | A simple notation for differential vector expressions in orthogonal curvilinear coordinates[END_REF], we have the following relations for the jumps of the stress components at order one:

[ ˜ 1 13 ] = - 1 h 0 1 M 11,1 + 1 h 0 2 M 12,2 + h 0 1,2 h 0 1 h 0 2 (M 12 + M 21 ) + h 0 2,1 h 0 1 h 0 2 (M 11 -M 22 ) + h 0 1,3 h 0 1 M 31 + h 0 1,3 h 0 1 + h 0 2,3 h 0 2 M 13 , (67)
[ ˜ 1 23 ] = - 1 h 0 1 M 21,1 + 1 h 0 2 M 22,2 + h 0 1,2 h 0 1 h 0 2 (M 22 -M 11 ) + h 0 2,1 h 0 1 h 0 2 (M 12 + M 21 ) + h 0 2,3 h 0 2 M 32 + h 0 1,3 h 0 1 + h 0 2,3 h 0 2 M 23 , (68)
[ ˜ 1 33 ] = -

1 h 0 1 M 31,1 + 1 h 0 2 M 32,2 + h 0 1,2 h 0 1 h 0 2 M 32 + h 0 2,1 h 0 1 h 0 2 M 31 + h 0 1,3 h 0 1 + h 0 2,3 h 0 2 M 33 - h 0 1,3 h 0 1 M 11 - h 0 2,3 h 0 2 M 22 , (69)
where

M 11 = (2 + ) ˜ 0 33 + 4 ( + ) (2 + ) 1 h 0 1 ũ0 1,1 + h 0 1,2 h 0 1 h 0 2 ũ0 2 + h 0 1,3 h 0 1 ũ0 3 + 2 (2 + ) 1 h 0 2 ũ0 2,2 + h 0 2,1 h 0 1 h 0 1 ũ0 2 + h 0 2,3 h 0 2 ũ0 3 , (70)
M 22 = (2 + ) ˜ 0 33 + 2 (2 + ) 1 h 0 1 ũ0 1,1 + h 0 1,2 h 0 1 h 0 2 ũ0 2 + h 0 1,3 h 0 1 ũ0 3 + 4 ( + ) (2 + ) 1 h 0 2 ũ0 2,2 + h 0 2,1 h 0 1 h 0 1 ũ0 2 + h 0 2,3 h 0 2 ũ0 3 , (71)
M 12 = M 21 = 1 h 0 2 ũ0 1,2 + 1 h 0 1 ũ0 2,1 - h 0 1,2 h 0 1 h 0 2 ũ0 1 - h 0 2,1 h 0 1 h 0 2 ũ0 1 , (72)
M i3 = M 3i = ˜ 0 i3 , i = 1, 2, 3. (73)

Implementation of the imperfect interface conditions

In this section we show how the conditions (41), ( 47), ( 50) and ( 60) can be used to obtain imperfect interface conditions for the original (unrescaled) problem. Consider a minimizer, u ε , of the energy (6), and a smooth extension of u ε to 0 ± , 0 g and 0 u , which are the domains obtained by taking the geometric limit of ε ± , ε g and ε u as ε → 0 + (Fig. 1c). Assume also the existence of an expansion in powers of ε also for u ε :

u ε = u 0 + εu 1 + ε 2 u 2 + o(ε 2 ). ( 74
)
In view of the change of variable (( 12) and ( 13)), we have

u ε ˆ , 30 ± ε 2 ± = ũε ± ˆ , 30 ± 1 2 ± , ˆ , ˆ ∈ S. (75)
Substituting the expansions ( 17) and ( 74) into (75), we obtain

u 0 ˆ , 30 ± ε 2 ± + εu 1 ˆ , 30 ± ε 2 ± + o(ε 2 ) = ũ0 ± ˆ , 30 ± 1 2 ± + ε ũ1 ± ˆ , 30 ± 1 2 ± + o(ε 2 ). ( 76
)
Now we assume that the terms on the left-hand side can be expanded in powers of ε

u 0 ( ˆ , ± 30 ) ± ε 1 2 u 0 ,3 ( ˆ , ± 30 ) + u 1 ( ˆ , ± 30 ) + o(ε 2 ) = ũ0 ± ˆ , 30 ± 1 2 ± + ε ũ1 ± ˆ , 30 ± 1 2 ± + o(ε 2 ), ( 77
)
and we identify the terms in ε to get

u 0 ( ˆ , ± 30 ) = ũ0 ± ˆ , 30 ± 1 2 ± , (78)
± 1 2 u 0 ,3 ( ˆ , ± 30 ) + u 1 ( ˆ , ± 30 ) = ũ1 ± ˆ , 30 ± 1 2 ± . (79)
Therefore, we have

[[u 0 ]] = [ ũ0 ] (80) [[u 1 ]] = [ ũ1 ] -w(u 0 ,3 ) (81) rr = 3IJ-4 ř3 ,  = = 3IJ+ 2 ř3 , (94)
where ˛and ˇare initially undetermined constants. These solutions apply inside the three phases, with appropriate values of the elastic moduli and with different values of the constants ˛, ˇ. Let ˛ε ± , ˇε ± , ˛ε, ˇε denote the six constants of the three-phase elastic problem. The condition ˇε -= 0 is necessary to avoid singularities at r = 0. The remaining five constants are determined by imposing the boundary condition u r (R e ) = ı and the continuity of the radial displacement u r and of the radial stress rr at r = R 0 ± ε 2 . We obtain the following conditions:

˛ε -R 0 - ε 2 = ˛ε R 0 - ε 2 + ˇε (R 0 -ε/2) 2 , (95)
3Ä -˛ε -= 3IJε -4 ˇε (R 0 -ε/2) 3 , ( 96
) ˛ε R 0 - ε 2 + ˇε (R 0 + ε/2) 2 = ˛ε + R 0 - ε 2 + ˇε + (R 0 + ε/2) 2 , (97)
3IJε -4 ˇε

(R 0 + ε/2) 3 = 3Ä + ˛ε + -4 ˇε + (R 0 + ε/2) 3 , ( 98
)
˛ε + R e + ˇε + R 2 = ı. ( 99
)
the perfect interface conditions calculated at order zero (cf. problem (P 0 )) are applied to the same problem. Let ˛0 ± , ˇ0 ± denote the four constants of this two-phase elastic problem. Again, the condition ˇ0 -= 0 is necessary to avoid singularities at r = 0. The remaining three constants are determined by imposing the boundary condition and the continuity of the radial displacement and of the radial stress at r = R 0 . We obtain the following conditions:

˛0 -R 0 = ˛0 + R 0 + ˇ0 + R 2 0 , (100)
3Ä -˛0 -= 3Ä + ˛0 + -4 + ˇ0 + R 3 0 , (101)
˛0 + R e + ˇ0 + R 2 e = ı. (102)
Lastly, we solve two-phase elastic problem with the imperfect interface conditions calculated at order one (cf. problem (P 1 )). In the spherical coordinate system, we have

h 1 := h = r sin Â, h 2 := h  = r, h 3 := h r = 1. (103)
Using ( 64)-( 73) and ( 103), the imperfect interface conditions (91) 4 and (91) 5 assume the form

[[u 1 r ]] = 3 (3Ä + 4 ) 0 rr (R 0 ) -2 (3Ä -2 ) (3Ä + 4 ) u 0 r (R 0 ) R 0 - 1 2 (u 0 r,r (R + 0 ) + u 0 r,r (R - 0 )), (104)
[[ 1 rr ]] = - 12 (3Ä + 4 ) 0 rr (R 0 ) R 0 + 36Ä (3Ä + 4 ) u 0 r (R 0 ) R 2 0 - 1 2 ( 0 rr,r (R + 0 ) + 0 rr,r (R - 0 )). ( 105
)
Let ˛1 ± , ˇ1 ± denote the four constants of the two-phase elastic problem with the imperfect interface conditions. We still find that the condition ˇ1 -= 0 avoids singularities at r = 0. The remaining three constants are determined by imposing ( 104) and ( 105), and the boundary condition at r = R e with ı = 0. We obtain the following conditions:

˛1 + R 0 + ˇ1 + R 2 0 -˛1 -R 0 = 9Ä -˛0 - (3Ä + 4 ) -2 (3Ä -2 ) (3Ä + 4 ) ˛0 -- 1 2 ˛0 + - 2ˇ0 + R 3 0 + ˛0 -, (106)
3Ä + ˛1 + -4 + ˇ1 + R 3 0 -3Ä -˛1 -= 36 (Ä -Ä -) (3Ä + 4 ) ˛0 - R 0 - 6 + ˇ0 + R 4 0 , (107)
˛1 + R e + ˇ1 + R 2 e = 0. (108)
Comparisons between the radial and circumferential stresses as functions of = log 10 Ä/Ä + on the external surface r = R e calculated from the exact three-phase solution, the perfect and the imperfect interface approximations are given in Figs. 2 and 3 for the two cases ε = 0.1R 0 and ε = 0.01R 0 . For drawing the plots, the following values for the material and geometrical constants were assumed:

-= 10 + , (109)
-= 0.20, + = 0.35, = 0.3, ( 110)

R e = 2R 0 , (111)
where ± are the Poisson coefficients of the materials "±" and is the Poisson coefficient of the interphase material. It should be noted that, for the values of ε considered in Figs. 2 and3, the imperfect interface model is more accurate than the perfect interface model for / + ranging from 10 -1 to 10 1 .

Conclusions

In this study, a model of a curved thin elastic interphase is formulated which is characterized by two-level interface conditions. At the first level, the interphase is replaced by a model of perfect interface, for which the traction and the displacement vector fields are continuous (cf. conditions (90) 4 and (90) 5 ). At the second level, the jumps in the displacements and in the stress vector along the interface are related to the interphase material properties and to the displacement and stress fields obtained at the first level (cf. conditions (91) 4 , (91) 5 ). The interface conditions are obtained in a system of parallel orthogonal curvilinear coordinates and are given in terms of the displacements and tractions belonging to the adherents and evaluated at both sides of the interphase. The main advantage of the interface model is that it makes possible to solve for the fields in the adherents without having to solve for the fields in the interphase. To illustrate its implementation, a spherical composite assemblage is studied and comparisons between the radial and circumferential stresses computed from the exact three-phase solution, the perfect and the imperfect interface approximations are given. A motivation for the study of interface models is the complexity of the finite element modeling due to the thinness of the thin layer, whose presence causes an ill-conditionment of the tangent matrices and affects the precision of the computations. The computational cost may be amplified by the differences between the material properties of the interphase and those of the adherents. The use of the imperfect interface model brings the intrinsic advantage of reducing the numerical analysis to a two-phase equilibrium problem. The imperfect interface laws obtained in this paper could also been advantageously used to evaluated the effect of coating interphases on the effective properties of composites. Implementation of the two-level model in a finite element software and a study the singularities at the edges of the interface will be addressed in a future work.

Fig. 1 .

 1 Fig. 1. (a) Original three phase reference configuration "adherent +"/interphase/"adherent -", (b) rescaled three phase reference configuration after the change of variables and (c) limit two phase reference configuration.

Fig. 2 .

 2 Fig. 2.Comparative plots for ε = R0/10 of the radial and circumferential stresses on the external surface calculated from the exact three-phase solution ( ε rr , ε ÂÂ ), the perfect ( 0 rr , 0 ÂÂ ) and the imperfect ( 0 rr + ε 1 rr , 0 ÂÂ + ε 1 ÂÂ ) interface approximations.

Fig. 3 .

 3 Fig. 3. Comparative plots for ε = R0/100 of the radial and circumferential stresses on the external surface calculated from the exact three-phase solution ( ε rr , ε ÂÂ ), the perfect ( 0 rr , 0 ÂÂ ) and the imperfect ( 0 rr + ε 1 rr , 0 ÂÂ + ε 1 ÂÂ ) interface approximations.

Acknowledgments

RR thanks the financial support of the Italian Ministry of Education, University and Research (MIUR) through PRIN project "Multi-scale modeling of materials and structures" (code 2010BFXRHS).

where [[f]]( ˆ ) := f( ˆ , + 30 ) -f( ˆ , - 30 ) is taken to denote the jump across S of a generic function f : 0 + ∪ 0 -→ R 3 and w(f)( ˆ ) := (1/2)(f( ˆ , + 30 ) + f( ˆ , - 30 )). Similar relations can be obtained for the tractions:

[

Substituting the conditions (41), ( 47), ( 50) and ( 60) into ( 80), ( 81), ( 84), (85) and using ( 78) and ( 82) we obtain the following interface conditions for the original limit problem:

with 0 ê3 the common value of the tractions on both sides of S. Using again the change of variables ( 12) and ( 13) in ( 44)-( 46), ( 58) and ( 59), we obtain the following two equilibrium problems posed only on the limit configurations of the adherents and each in which the interphase is treated as a surface:

(91)

Composite spheres assemblage

To illustrate the implementation of the imperfect interface conditions developed so far, we consider the boundary value problem of a composite spheres assemblage subject to uniform radial displacement of magnitude ı on the outer boundary. The problem will be solved in two ways: first, exactly, as a three-phase elastic problem where the phases are a sphere of material "+", a concentric interphase, and an inner spherical inclusion of material "-", and second, as a two-phase problem, where the phases are the two concentric spheres of materials "+" and "-" with in between imperfect interface conditions.

We take R e , R 0 + ε/2, to denote the exterior and the core radii of the outer sphere, and R 0ε/2 to denote the exterior radius of the inclusion. Let Ä + , + denote the elastic bulk and shear moduli of the outer spherical shell, Ä, the elastic bulk and shear moduli of the interphase, and Ä -, -the elastic bulk and shear moduli of the inclusion.

The general three-dimensional equations of homogeneous, isotropic linearized elasticity have the following general solution in terms of spherical coordinates r, Â, [START_REF] Lakes | Dramatically stiffer elastic composite materials due to a negative stiffness phase[END_REF]:


Ange Kongo-Konde
email: ange.kongo-konde@airbus.com

Iulian Rosu
email: rosu@lma.cnrs-mrs.fr

Frédéric Lebon
email: lebon@lam.cnrs-mrs.fr

Lucie Seguin

Olivier Brardo
email: olivier.brardo@airbus.com

Fabien Troude
email: florian.troude@airbus.com

Bernard Devesa
email: bernard.devesa@airbus.com

A Kongo Kondé

B Devésa

Thermomechanical couplings in aircraft tire rolling/sliding modeling

Keywords: Aircraft tire, Finite element, Friction law, Rubber, Steady State Rolling, Temperature

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Landing and take-off operations are the critical phases of a flight. Indeed, the high speed, the violent impact when aircraft touches the ground, the need of braking efficiently, the ground irregularities, the climatic conditions (crosswind, contaminations, …) as well as other parameters can make these maneuvers on ground very perilous.

Knowing at each time the forces in the contact patch could allow mastering the behaviour in the tire/ground interface. For aeronautical applications, the experiment on real tire, at moderated speeds, yields to knowledge of tire/ground behaviour in this reduced field. However, the high nonlinearities related to materials, geometry due to large deformations, coupled phenomena (contact with friction, thermal-mechanical coupling, wear) reduce the trust in extrapolation capabilities of these results at high speeds.

When we study tire yawed rolling, it is relevant to well estimate the Self Aligning Torque and the ratio, lateral force/ normal force. These respectively express the ability of tire to return in its steering direction when it is no longer steered, and the lateral grip potential of the rubber/ground couple.

Indeed, the ratio µ Y (between lateral force and vertical one) defines the equilibrium which exists between steering lock action and the trajectory of the vehicle. During µ Y upward phase, there is equilibrium, the more the driver or the pilot turns the wheel the more β increases. It corresponds to pseudo-sliding in the contact patch up to β max (around β=4 to 7° for automobile pneumatics, around β=6 to 10° for trucks and around 15 to 25° for aircraft, see Fig. 1).

Beyond these values, a loss of equilibrium sets up: the more the driver tries to increase the drift, the more the vehicle slides. This phase corresponds to pure sliding in the contact. The rubber is warming up, the temperature increases and implies decreasing of the friction coefficient µ. Thus, we observe the decreasing of µ Y and M Z as shown in Fig. 1. The specific values are normalized for analysing the observed phenomena. In skid, the result of friction forces is at the rear of the contact area center. This intrinsic property of rubber creates resistant torque which faces drift: this is the Self -Aligning Torque, M Z . It represents the pneumatic ability to return to his steering direction when the driver releases vehicle controls. It is maximal for slip angle β=6° and rapidly vanishes towards zero (around β=18° for aircraft tire), see Fig. 1.

Numerical modelling using Finite Element Analysis and the current computer capabilities open the way to produce finer and finer simulations of industrial problems formerly limited in terms of degrees of freedom. This new way allows also reducing the operating cost of experiments.

Aim of the study

The aim of this work is to highlight the effects of loading parameters (vertical load, velocity, inflating pressure) on the mechanical behaviour of tire in contact with the ground. In extreme conditions (crab landing due to crosswind for example), the high energy dissipation in the contact zone results in an increasing of temperature and wear. The local friction coefficient (µ local ) depends, a priori, highly on temperature, wear, sliding velocity and contact pressure. To emphasize the effects of these variations on the global mechanical behaviour of µ Y and M Z , we must consider a thermal-mechanical coupling. To do that, we propose a thermal model of tire in which the aim is to define the evolution of temperature during tire rolling and we also focus on the temperature diffusion in its thickness.

Finite Element Model for Tire Analysis FEA Model

Finite Element Analysis (FEA) is used nowadays to predict behaviour and performance of structures. Regarding our study, an accurate 3D tire model is essential. The fundamental requires of such modelling process are the accuracy and adequate information of the tire cross-section [START_REF] Olatunbosum | Finite element simulation of the tyre burst test[END_REF], particularly the internal features of plies and their modelling. The simulation process includes these following features:

• Tire geometry: Some image acquisition technologies developed for inspection operations in the tire manufacturing process are now being adapted by researchers for acquiring tire geometric data for relevant studies in the absence of any CAD tire geometry data [START_REF] Ghoreish | Finite Element Analysis of Steel-Belted Radial Tyre with Tread Pattern under Contact Load[END_REF]. This method provides the 2D section of tire including all tire inner features. So, the tire data are then stored in appropriate format and imported into ABAQUS software for further FE Modeling. Suitable points are then chosen as the nodes before including them into FE input file.

Fig. 2: 2D-section and 3D-partial models.

• Rubber and Reinforced fiber behaviours: The hyperelastic model of Mooney Rivlin [START_REF] Lahellec | Second-order estimate of the macroscopic behaviour of periodic hyperelastic composites: theory and experimental validation[END_REF] is largely used to model the rubber behaviour while orthotropic elastic model defines the reinforced fiber behaviour. The material properties considered in this study are derived from experimental tension test performed on samples cut in each different layers of the tire (see Fig. 2). We assume that the dependence of mechanical characteristics on temperature is neglected.

• Tire/rim contact: The contact of the tire to the rim is simplified with assumption that tire sticks to the rim. The nodes at the tire/rim contact are fixed at the rim by the tie option in ABAQUS.

• Coulomb friction at tire/ground contact: The tire/ground contact complicates the finite element model, since contact and friction problems are highly non-linear. Nonetheless, both effects cannot be ignored. The contact problem in FEM model is described by a deformable body (tire) to rigid body (ground) option and modeled by soft contact using penalty regularization. The Coulomb friction is modeled by a stiffness method.

Table 1: Mechanical assumptions.

The 3D finite element model is created in two steps. Firstly, the 2D mesh is created. Based on this mesh, 3D model is generated by revolution. Fig. 2 describes the FEM model. It shows the FE mesh of 2D model. It consists on 4-nodes elements and some 3-nodes required for the modeling of geometrical transition between the tread and sidewall and between sidewall and shoulder.

These 3-and 4-nodes elements are transformed, respectively, in 3D elements with 6-and 8-nodes.

The reinforced layers are modeled using in the same way. Fig. 2

Thermomechanical Coupling

In grip phenomenon, two mechanisms take place under the influence of relative sliding between rubber and ground: indentation and adhesion. In both cases, rubber viscoelastic characteristics play an important role. Some studies on rubber samples show significant influence of contact pressure, of the temperature and of the sliding velocity on the friction coefficient. Generally, this friction coefficient reaches a maximum value around 90°C before decreasing with the increasing of temperature [START_REF]Le pneu -L'adhérence[END_REF]. The temperature evolution in the contact patch is, thus, an important step which allows knowing the local friction coefficient µ local .

During a contact with friction between two bodies (Ω1 and Ω2), there is generation of heat in the contact due to energy dissipation. Generally, to solve a thermal problem leads to finding the solution of the heat equation (Eq. 1) with boundary conditions: Dirichlet and Neumann (convection, adiabatic and given temperature).

Q gradT k div dt dT c p + = ) . ( . ρ (1)
Where ρ is the mass density (kg/m 3 ), c p is the specific heat (Jkg -1 K -1 ), k is the thermal conductivity (Wm -1 K -1 ), T is the temperature (K), Q is the internal and external heat sources (Wm -3 ) and t is the time (s). In our study, the heat equation and material characteristics are only applied on the elastic body and the dissipated flux due to friction Q f can be written as a function of α. So, this flux is written on the contact boundary as:

T N local f u T Q . ). ( . λ µ α = (2)
With λ N contact pressure, T u tangential speed in the contact, µ local (T) the local friction coefficient depending on temperature, and α the share coefficient. We assume that material properties do not vary with the temperature. The friction coefficient only depends on the temperature.

Methods and Analysis approaches

Two approaches [START_REF] Kongo Kondé | Etude du comportement en roulement d'un pneu d'avion[END_REF] can be used to resolve this kind of problem. The first one is called "steady state analysis". This is a mixed Lagrangian/Eulerian approach in which, for an observer in the reference frame, the tire is a fixed set while materials flow through the refined stationary mesh.

The second one is called "Classical Dynamic analysis". This Lagrangian approach allows simulating real rolling. In this approach, for an observer in the reference frame, the tires rolls and, turn after turn, the elements touch and leave the contact area.

We have already discussed about the advantages and disadvantages of the two strategies [START_REF] Kongo Kondé | Etude du comportement en roulement d'un pneu d'avion[END_REF]. Numerically, both approaches are equivalents.

Static Analysis

The static analysis is performed in consecutive stages. In the first stage, the tire is modeled under inflation pressure using 2D axisymmetric model. Inflation analysis of axisymmetric model is followed by 3D static analysis of vertically loaded tire. The analysis consists of several steps. In this step, the 3D model is created, by rotation of an axisymmetric model around tire axis. In step one tire is inflated, while in step two ground surface gets gradually moved towards tire axis, until it reaches the goal value. In step three, equilibrium is reached, with constant force acting on surface. This step is needed as the basis for subsequent steady state rolling analysis or dynamic rolling analysis.

Steady State Analysis "SS"

In this analysis, the tire was modeled under steady state rolling conditions [START_REF] Nackenhorst | The ALE-formulation of bodies in rolling contact -Theoritical foundations and finite element approach[END_REF][START_REF] Korunovic | Finite Element Model for Steady-state Rolling Tire Analysis[END_REF]. This stage is restarted from the results of the second step (footprint analysis). A constant given ground velocity was assumed for the tire. The tire spinning velocity was changed to determine the tire free rolling condition and its dynamic loaded radius. In addition, the friction between tire tread surface and ground was taken into account using the simple linear Coulomb's law. The last analysis step is devoted to the application of the slip angle to the tire for a cornering situation. Slip angle is defined as the angle between the direction of travel and a normal plane to the axis of rotation. Analyses were performed using these values in conjunction with tire free rotational speed found in the previous step. The cornering is obtained by divided the rolling speed v 0 into v x =v 0 cosβ and v y =v 0 sinβ, where β is the slip angle.

Results

Static validation of numerical model

This static analysis is primarily performed because of model integrity checked, finding the deformed shape of the model. Fig. 3 shows the load-deflection curves from our tire model, manufacturer data and experimental results. The comparison reveals that there is a good matching between reference data (manufacturer data and experimental results) and computed values of tire model. So, what about the contact area?

We do not model the tread pattern features. So, the tire tread is smooth and without grooves, (see Fig. 2). The deletion of the tire features modifies the contact area and makes it greater. The increasing of contact area due to this observation is about 21% (see Table 2). Nevertheless, this value remains large. Further investigations have to be done to better match the contact area.

Simulation of tire rolling

Simulations are carried for different rolling velocities (v =0.6 v 0 to 6v 0 ), vertical loads (F Z = 0.15F 0 ; 0.39F 0 and 0.69F 0 ), nominal inflating pressure (P = P 0 ), slip angle (β = 16°), constant friction coefficient µ 0 and varying friction µ=µ 0 .f (T).

Loading parameters effects

• Rolling velocity effects: Fig. 4 shows the comparison between simulations performed at constant vertical load, pressure and friction coefficient by varying the rolling velocity. We observe that µ Y is independent of rolling velocity, but there exists a threshold velocity at which the magnitude of M Z decreases exponentially with the growing speed. The value of β max related to the maximum of M Z does not change. So, the high rolling velocities affect the magnitude of M Z .

Fig. 4: Rolling velocity effects on µ Y and M Z

• Vertical load effects: Fig. 5 shows the evolution of µ Y and M Z with the vertical load F Z .

The influence of vertical load on the mechanical tire response is important. We observe that µ Y is not sensitive to F Z while β max increases. We also observe a sensitive increasing of M Z magnitude and a shift of β max (Fig. 5) toward greater values while the vertical load increases. This observation was predictable. Indeed, the forces in the contact patch increase proportionately with the vertical load. De facto, β max corresponding to the slip angle for which contact nodes slide, increases. • Thermal -Mechanical Simulations: Thermal characteristics (conductivity, specific heat) are derived from [START_REF] Yuksel | A computation model to predict the thermomechanical behavior of automobile tires, Constitutive Models for Rubber III[END_REF]. The pure sliding (Fig. 6) is obtained by pulling the tire. We use the ABAQUS thermal -mechanical capabilities to performe the following simulations. The tire is pulled on ground along 4 meters at constant speed of v 0 , v 0 is the reference velocity. In this case, we observe that the temperature increases significantly when the tire is lengthways slided down even on a short distance. The rubber reaches a superficial temperature of 190°C (Fig. 6) close to critical temperature for which the rubber is deposited on concrete road [START_REF] Sakai | Thermal Engineering Analysis of Rubber Vulcanization and Tread Temperatures during Severe Sliding of a Tire[END_REF]. The running velocities considered are v 0 and 2v 0 at β=16°. The Fig. 7 shows the temperature distribution of the tire running at v 0 . In every instance, the rise in temperature is located on tire tread surface with the temperature gradient within the thickness (see Fig. 7 a). In Fig. 7 b, we observe that the warming is in the sliding direction. The temperature is higher from right side to left side on tire tread surface. The reinforced zone seems to be a thermal barrier. However, the thermalmechanical analysis remains too expensive. Fig. 8 shows the evolution of temperature of a tire tread node for the two velocities. This evolution is a step-like curve where the plateau represents the node temperature when it is out of the contact patch and the slope characterizes the increasing node temperature in the contact patch. When the node is out of the contact region, its temperature keeps constant because we do not consider the convective effects between tire and surrounding air. Normally, we should have a slight decrease due to convection which is a very complex phenomenon, and naturally depends on velocity k = f(v).

The main conclusion of this analysis is that the temperature values depend on the speed at which the tire is rolling. • Effects of considering a variable friction coefficient: The warm-up to due friction in the contact patch has the effect of modifying the local friction coefficient. To take these effects into account in the modelling of the tire behaviour, we consider the following friction law (Eq. 3) [START_REF] Awrejcewicz | Nonsmooth Dynamics of Contacting Thermoelastic Bodies[END_REF] in which friction coefficient depends on temperature.

        + - + = 1 ) ( ) ( 2 0 3 2 1 0 T T f f f T local µ µ (3)
We performed the above test rig (Fig. 10) inspired from Coulomb apparatus. Samples are cut in tread rubber and mounted on an aluminium block to ensure optimal contact between sample and the surface and to keep the sample at constant temperature.

Conclusion

The present research is a first step on the modelling of aircraft tire in which thermal effects are described via a friction law depending only on temperature. A Coulomb-like test rig allows the measuring of the static friction coefficient of rubber/concrete and rubber/asphalt. Indeed, taking into account these effects in the friction coefficient, we emphasize the decreasing of µ Y and the rapid vanishing of M Z towards zero beyond β max for which the full grip potential is reached as we did by considering, in first approximation, a Chichinadze-like law and by assuming that temperature map is known. In the further work, we will focus on the development of a cost-effective alternative thermo-mechanical approach (coupling steady state mechanical analysis and transient thermal analysis) and work on the automation of the proposed approach and improve the implementation of this algorithm in ABAQUS.

Fig. 1 :

 1 Fig. 1: Self -Aligning Torque and Lateral Friction coefficient for aircraft tire.

Fig. 3 :

 3 Fig. 3: Vertical deflection and contact area for β=16°C

Fig. 5 :

 5 Fig. 5: Load effects on µ Y and M Z

Fig. 6 :

 6 Fig. 6: Thermal result in pure sliding.

Fig. 7 :

 7 Fig. 7: Thermal result in cornering.

Fig. 8 :

 8 Fig. 8: Rise in temperature of a tire node at vo and 2vo

Fig. 9 :

 9 Fig. 9: Aluminium block and sample

Fig. 12 .

 12 Fig. 12. Comparison of mechanical response with constant and varying friction coefficient (Chichinadze) on µ Y and M z

Table 2 :

 2 Static validation of tire model.