Eric Paroissien
email: eric.paroissien@isae-supaero.fr

Lucas F M Da Silva

Frederic Lachaud

Lucas Fm Da Silva

Frédéric Lachaud

Simplified stress analysis of functionally graded single-lap joints subjected to combined thermal and mechanical loads

Keywords: functionally graded adhesive, single-lap bonded joint, stress analysis, macroelement, thermoelasticity, dissimilar adherend

published or not. The documents may come

Simplified stress analysis of functionally graded single-lap joints subjected to combined thermal and mechanical loads

NOMENCLATURE AND UNITS

𝜒 𝛼 adherend thermal unbalance parameter (-) (X,Y,Z) element reference system of axes

Introduction

In the frame of structural design, the proper choice of joining technology is decisive for the integrity of the manufactured structure. Mechanical fastening, such as riveting or screwing, appears to be a reliable solution for the designers. Nevertheless, alone or in combination with mechanical fastening, the adhesive bonding technology may offer significantly improved mechanical performance in terms of stiffness, static strength and fatigue strength [START_REF] Hart-Smith | Design methodology for bonded-bolted composite joints[END_REF][START_REF] Kelly | Quasi-static strength and fatigue life of hybrid (bonded/bolted) composite single-lap joints[END_REF][START_REF] Da Silva | Handbook of Adhesion Technology (2 volumes)[END_REF].

Indeed, unlike the discrete load transfer of mechanical fasteners, the load transfer between structural bonded components is continuous all along the overlap. This higher level of mechanical performance allows for lighter joints. In other words, adhesive bonding offers the possibility to reduce the structural mass while ensuring the mechanical strength. The optimization of the strength-to-weight ratio is a challenge for several industrial sectors, such as aerospace, automotive, rail or naval transport industries.

Nevertheless, stress gradients at both overlap ends appear in bonded joints, due to the relative deformation of the adhesive layer with regards to the adherends. It leads to a load transfer restricted on a small length at the overlap ends. In order to increase the load capability of bonded joints, the reduction of adhesive peak stresses is wanted. The specimen design for the thick adherend shear test [START_REF] Astm D | Standard Test Method for Thick-Adherend Metal Lap-Shear Joints for the Stress-strain Behavior of adhesives by Tension Loading[END_REF] leads to both a homogenization of the adhesive shear stress and a drastic reduction the adhesive peel stress, all the more when care is taken to reduce the edge effects [START_REF] Cognard | Influence of adhesive thickness on the behaviour of bonded assemblies under shear loadings using a modified TAST fixture[END_REF]. Another approach is to make the material and/or geometrical properties of the adherends and/or the adhesive layer vary along the overlap. Several design solutions have been published [START_REF] Da Silva | Handbook of Adhesion Technology (2 volumes)[END_REF]. For example, a solution is the tapering of adherends at overlap ends, which allows for a progressive increase of the neutral line lag and a reduction of adhesive peel stress [START_REF] Hart-Smith | Adhesive-bonded scarf and stepped-lap joints[END_REF][START_REF] Oterkus | Nonlinear analysis of bonded composite joints[END_REF]. A more local solution is the rounding of adherend corner associated with adhesive spew fillets [START_REF] Zhao | Single lap joints with rounded adherend corners: Stress and strain analysis[END_REF][START_REF] Zhao | Single lap joints with rounded adherend corners: Experimental results and strength prediction[END_REF]. The mixed adhesive solution which is a rough version of a graded joint consists in the use of various different adhesives along the overlap to increase the joint strength [START_REF] Raphael | Variable-adhesive bonded joints[END_REF][START_REF] Da Silva | Joint strength optimization by the mixed-adhesive technique[END_REF][START_REF] Breto | Finite Element Analysis of Functionally Graded Bond-Lines for Metal/Composite Joints[END_REF][START_REF] Machado | Influence of low and high temperature on mixed adhesive joints under quasi-static and impact conditions[END_REF]. In recent past years, functionally graded adhesive (FGA) have been more

and more considered [START_REF] Durodola | Functionally graded adhesive joints -A review and prospects[END_REF][START_REF] Kawasaki | Functionally Graded Adhesive Joints Bonded by Honeymoon Adhesion Using Two Types of Second Generation Acrylic Adhesives of Two Components[END_REF]. FGA joints involve a continuous variation of the adhesive properties along the overlap allowing for the homogenization of the stress distribution and load transfer. When dissimilar adherends have to be bonded, the adhesive stress distribution is asymmetrical, so that one of the overlap ends is overstressed. Moreover, this overstressing is magnified under thermal loads due to the mismatch in coefficient of thermal expansion (CTE) of adherends. The capability of a local graduation of the adhesive stiffness is a promising solution to optimize the strength of multimaterial joints under severe loads, such as combined thermal and mechanical. This situation occurs very often in multi-material structures found in the transport industry. That is why the development of dedicated stress analyses to predict the stress distribution is fundamental. The Finite Element (FE) method is able to address the stress analysis of FGA joints [START_REF] Breto | Finite Element Analysis of Functionally Graded Bond-Lines for Metal/Composite Joints[END_REF][START_REF] Durodola | Functionally graded adhesive joints -A review and prospects[END_REF]. Nevertheless, since analyses based on FE models are computationally costly, it would be profitable both to restrict them to refined analyses and to develop simplified approaches, enabling extensive parametric studies and optimization processes. Moreover, numerous simplified stress analyses of bonded joints are available and provide accurate predictions [START_REF] Van Ingen | Stress analysis of adhesively bonded single lap joints[END_REF][START_REF] Tsaï | An evaluation of analytical and numerical solutions to the single-lap joint[END_REF][START_REF] Da Silva | Analytical models of adhesively bonded joints-Part I: Literature survey[END_REF]. In 2014, Carbas et al. published a first analytical approach for 1D-bar stress analysis of FGA joints [START_REF] Carbas | Modelling of functionally graded adhesive joints[END_REF]. This stress analysis is based on the shear-lag approach by Volkersen [START_REF] Volkersen | Die Nietkraftverteilung in Zugbeanspruchten Nietverbindungen mit konstanten Laschenquerschnitten[END_REF] associated with a resolution scheme making use of Taylor expansion in power series (TEPS) to solve the involved differential equations. This stress analysis is restricted to half of the overlap length of balanced joints with a linear graduation of the adhesive shear modulus. Stein et al. presented a 1D-bar analysis using TEPS resolution able to address unbalanced bonded joints under any adhesive properties graduations [START_REF] Stein | Stress solution for functionally graded adhesive joints[END_REF][START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF]. This analysis is called by the authors Improved Shear Lag Model (ISLM). Moreover, Stein et al. provided a sandwich-type analysis using TEPS resolution, taking into accounts both in-plane and out-of-plane load, termed General Model (GM). The sandwich-type analysis concept comes from the analysis methodology by Goland and Reissner [START_REF] Goland | The stresses in cemented joints[END_REF] who provided the first closed-form solution for the adhesive stress distribution for simply supported balanced joint made of adherends undergoing cylindrically bending. Goland and Reissner took into account the geometrical non linearity due to the lag of neutral line to assess the bending moment at both overlap ends through a bending moment factor. This methodology was then employed by other researchers to improve the initial model [START_REF] Hart-Smith | Adhesive bonded single lap joints[END_REF][START_REF] Williams | Stresses in adhesive between dissimilar adherends[END_REF][START_REF] Bigwood | Elastic analysis and engineering design formulae for bonded joints[END_REF][START_REF] Oplinger | A Layered Beam Theory for Single-Lap Joints[END_REF][START_REF] Tsaï | Improved Theoretical Solutions for Adhesive Lap Joints Int[END_REF][START_REF] Högberg | Mechanical behavior of single-layer adhesive joints -An Integrated approach[END_REF][START_REF] Nemes | Modeling of cylindrical adhesively bonded joints[END_REF][START_REF] Luo | Fully-coupled nonlinear analysis of single lap adhesive joints[END_REF][START_REF] Weißgraeber | A general sandwich-type model for adhesive joints with composite adherends[END_REF] leading to various forms of the bending moment factor [START_REF] De L'armée | Bending moment calculation for single lap joints[END_REF]. In 2017, Stapleton et al. used a joint element (JE) for the stress analysis of FGA joints under various geometrical configurations, including in-plane and out-of-plane load as well as non-linear material behavior [START_REF] Stapleton | Design of functionally graded joints using a polyurethane-based adhesive with varying amounts of acrylate[END_REF]. A JE is a 4-nodes brick element allowing for the modelling of two bonded adherends [START_REF] Stapleton | Design of functionally graded joints using a polyurethane-based adhesive with varying amounts of acrylate[END_REF][START_REF] Stapleton | The analysis of adhesively bonded advanced composite joints using joint finite elements[END_REF]. Over a similar period of time, the first and third authors of the present papers and co-workers have been working on the development of the macro-element (ME)

technique for the simplified stress analysis of bonded, bolted and hybrid (bonded/bolted) joints [START_REF] Paroissien | Contribution aux Assemblages Hybrides (Boulonnés/Collés) -Application aux Jonctions Aéronautiques[END_REF][START_REF] Paroissien | Hybrid (bolted/bonded) joints applied to aeronautic parts: Analytical one-dimensional models of a single-lap joint[END_REF][START_REF] Paroissien | Analytical two-dimensional model of a hybrid (bolted/bonded) single-lap joint[END_REF][START_REF] Paroissien | A simplified stress analysis of bonded joints using macro-elements[END_REF][START_REF] Paroissien | Elasto-Plastic Analysis of Bonded Joints with Macro-Elements[END_REF][START_REF] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF][START_REF] Lélias | Mechanical behavior of adhesively bonded joints: Modeling, simulation and experimental characterization[END_REF][START_REF] Paroissien | Simplified Stress Analysis of Hybrid (Bolted/Bonded) Joints[END_REF]. Dedicated 4-nodes Bonded-bars (BBa) and Bonded-beams (BBe) have been formulated. As for the JE model, only one BBa or BBe, depending on the chosen kinematics, is sufficient to be representative for an entire bonded overlap in the frame of a linear elastic analysis (see Figure 1). When the geometrical or material properties of the adherends or the adhesive layer vary along the overlap a mesh is necessary along the overlap length direction only. The ME technique is inspired by the FE method and differs in the sense that the interpolation functions are not assumed. Indeed, they take the shape of solutions of the governing ordinary differential equations (ODEs) system, coming from the constitutive equations of the adhesive and adherends and from the local equilibrium equations, related to the simplifying hypotheses. The main work is thus the formulation of the elementary stiffness matrix of the ME. Once the stiffness matrix of the complete structure is assembled from the elementary matrices and the boundary conditions are applied, the minimization of the potential energy provides the solution, in terms of adhesive stress distributions along the overlap, internal forces and displacements in the adherends. The ME technique can be regarded as mathematical procedure allowing for the resolution of the system of ODE, under a less restricted application field of simplifying hypotheses, in terms of geometry, material behaviours, kinematics, boundary conditions and loads.

Stress analyses of bonded joints under thermal loads can be found in the literature linked to the aerospace [START_REF] Hart-Smith | Adhesive-Bonded Double-Lap Joints[END_REF] or to the emergence of the industry of electronic packaging [START_REF] Zeyfang | Stresses and strains in a plate bonded to a substrate: Semiconductor devices[END_REF][START_REF] Chen | Thermal stresses in Bonded Joints[END_REF][START_REF] Suhir | Stresses in bi-material thermostat[END_REF][START_REF] Marques | Adhesive joints for low and high temperature use: An overview[END_REF].

However, to the best knowledge of authors, there is not any published stress analyses of FGA joints under thermal load. The present paper aims at presenting simplified stress analysis of FGA single-lap joints under combined mechanical and thermal loads. Under the 1D-bar kinematics, the resolution scheme by TEPS and by ME is used. As the 1D-bar TEPS and ME analyses provide the same predictions, the resolution scheme by ME under the 1D-beam kinematics is employed only. Indeed, the ME technique offers the possibility to extend the application field to analyses involving nonlinear material behaviors, various geometries and various applied boundary conditions [START_REF] Paroissien | A simplified stress analysis of bonded joints using macro-elements[END_REF][START_REF] Paroissien | Elasto-Plastic Analysis of Bonded Joints with Macro-Elements[END_REF][START_REF] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF][START_REF] Lélias | Mechanical behavior of adhesively bonded joints: Modeling, simulation and experimental characterization[END_REF][START_REF] Paroissien | Simplified Stress Analysis of Hybrid (Bolted/Bonded) Joints[END_REF]. The developed stress analyses are then assessed against reference stress analyses for bonded and FGA joints on several test cases. The three stress analyses presented need dedicated computer codes, which are provided as supplementary materials with the present papers. These codes run on the MATLAB commercial software. Moreover, for the comfort of readers, this paper provides the useful mathematical steps, even if some elements have eventually been published elsewhere. The following hypotheses are taken (i) the adherends are linear elastic materials simulated as bars, (ii) the adhesive layer is simulated by an infinite number of linear elastic shear springs linking both adherends, and (iii) the shape of graduation of the adhesive layer shear modulus is considered. As a result, it is supposed that all the adhesive stress components vanish except the in-plane shear. The case of a single-lap joint subjected to combined mechanical and thermal loads is considered. The geometrical parametrization is provided in Figure 2. The subscript 1 (2) refers to the upper (lower) adherend. The origin of the global reference system is taken at the centre of the overlap, with the x-axis along the overlap length direction, the only axis according to which displacements u are possible. The joint is submitted to an uniform variation of temperature  T , to a tensile force F at one extremity and is fixed at the other one. The stress analysis is conducted in force but could similarly be made in tensile flow

F/b.

Governing equations

The local equilibrium of both adherends (see Figure 3) provides the following equations:

𝑑𝑁 𝑗 𝑏𝑑𝑥 = (-1) 𝑗 𝑇(𝑥), 𝑗 = 1,2 (1)
where b is the overlap width, N j the normal force in the adherend j and T the adhesive shear stress. The total strain is equal to the mechanical strain plus the thermal strain such as:

𝑑𝑢 𝑗 𝑑𝑥 = 𝑁 𝑗 𝐴 𝑗 + 𝛼 𝑗 Δ 𝑇 , 𝑗 = 1,2 (2)
where  j is the coefficient of thermal expansion of the adherend j. A j is the membrane stiffness of the adherend j, given by: where e j is the thickness of the adherend j and E j the Young's modulus of the adherend. The displacement u j (x) is the normal displacement of points located at the abscissa x on the neutral line of adherend j (see Figure 2).

𝐴 𝑗 = 𝐸 𝑗 𝑏𝑒 𝑗 (3)
The constitutive equation for the adhesive layer is provided by:

𝑇 = 𝐺 𝑎 𝑢 2 -𝑢 1 𝑒 𝑎 = 𝑘 𝐼𝐼 Δ𝑢 (4)
with:

Δ𝑢 = 𝑢 2 -𝑢 1 (5)
where e a is the adhesive thickness, G a the adhesive shear modulus and k II =G a /e a the adhesive shear relative stiffness. u is the differential displacement of the adherend interface. The stress analyses presented use k II andu, so that they can be directly applied when the thickness of the adhesive layer varies along the overlap.

TEPS resolution

The approach using TEPS resolution scheme is firstly used. The differentiation of Eq. (2) with respect to x provides:

𝑑 2 𝑢 𝑗 𝑑𝑥 2 = 1 𝐴 𝑗 𝑑𝑁 𝑗 𝑑𝑥 , 𝑗 = 1,2 (6)
By using the local equilibrium equation Eq. ( 1) and the adhesive constitutive equation Eq. ( 4), it comes:

𝑑 2 𝑢 𝑗 𝑑𝑥 2 = 1 𝐴 𝑗 (-1) 𝑗 𝑏𝑇(𝑥) = 𝑏 𝐴 𝑗 (-1) 𝑗 𝑘 𝐼𝐼 Δ𝑢, 𝑗 = 1,2 (7)
As a result, a second order differential equation in the slipping displacement (relative horizontal displacement of the interface) can be written:

𝑑²Δ𝑢 𝑑𝑥² -𝜂 ̃2𝑘 𝐼𝐼 Δ𝑢 = 0 (8)
with:

𝜂 ̃2 = 1 𝑒 1 𝐸 1 + 1 𝑒 2 𝐸 2 = 1+𝜒 𝐴 𝐴′ 2 (9)
𝜒 𝐴 = 𝐴′ 2 𝐴′ 1 = 𝑒 2 𝐸 2 𝑒 1 𝐸 1 = 𝐴 2 𝐴 1 (10)
A' j is the membrane stiffness of the adherend j per unit of width. 𝜒 𝐴 is representative for the stiffness unbalance of the joint. The differential equation Eq. ( 8) is relevant to the one obtained by Stein et al. for the ISLMwhich does not consider any thermal loadbut written in slipping displacement instead of shear strain [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF]. A solution is then searched for any x included between -c and c under the shape of TESP:

Δ𝑢(𝑥) = ∑ 𝑢 𝑛 𝑥 𝑛 ∞ 𝑛=0 (11)
For the series terms to have the same unit as the function approximated, the following variable change is made in the present analysis:

ζ = 𝑥 𝑐 (12)
As result, the solution is searched for any X included between -1 and 1 under the shape:

Δ𝑢(𝑥) = ∑ 𝑢 𝑛 (𝑐ζ) 𝑛 ∞ 𝑛=0 = ∑ 𝑢 𝑛 𝑐 𝑛 ζ 𝑛 ∞ 𝑛=0 = ∑ 𝑈 𝑛 ζ 𝑛 ∞ 𝑛=0 (13)
with:

∀𝑛, 𝑈 𝑛 = 𝑐 𝑛 𝑢 𝑛 (14)
The m th derivative of u is then assessed as follows

𝑑 𝑚 Δ𝑢 𝑑𝑥 𝑚 = 1 𝑐 𝑚 𝑑 𝑚 Δ𝑢 𝑑ζ 2 = 1 𝑐 𝑚 ∑ ∏ (𝑛 + 𝑖) 𝑚 𝑖=1 𝑈 𝑛+2 ζ 𝑛 ∞ 𝑛=0 (15)
The graduation of adhesive properties is then described under the shape of a TESP:

𝑘 𝐼𝐼 (ζ) = ∑ 𝐾 𝑛 ζ 𝑛 ∞ 𝑛=0 = ∑ 𝑘 𝑛 𝑥 𝑛 ∞ 𝑛=0 = 𝑘 𝐼𝐼 (ζ) (16)
with:

∀𝑛, 𝐾 𝑛 = 𝑐 𝑛 𝑘 𝑛 (17)
The expressions for u and k II are then replaced in the second order differential equation Eq.

(8) leading to:

∑ (𝑛 + 1)(𝑛 + 2)𝑈 𝑛+2 ζ 𝑛 ∞ 𝑛=0 -𝑐 2 𝜂 ̃2 ∑ 𝑈 𝑛 𝑋 𝑛 ∞ 𝑛=0 ∑ 𝐾 𝑙 ζ 𝑙 ∞ 𝑙=0 = 0 (18) ∀𝑛, 𝑈 𝑛+2 = 𝑐 2 𝜂 ̃2 (𝑛+1)(𝑛+2) ∑ 𝑈 𝑙 𝐾 𝑛-𝑙 𝑛 𝑙=0 (19)
This recursion relationship is of course similar to the one obtained by Stein et al. [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF]. Two boundary conditions have to be applied to be able to determine the terms of the series. In this paper, the values of the normal forces in the lower adherend are used, which reads:

𝑁 2 (𝑥 = -𝑐) = 0 = 𝑁 2 (𝑋 = -1) (20)
𝑁 2 (𝑥 = 𝑐) = 𝐹 = 𝑁 2 (𝑋 = 1) (21)
A relationship between the sipping displacement and the normal force in the lower adherend is then established as follows. From Eq. ( 2) and since the sum of the normal force of the upper adherend and that of the lower adherend is equal to F at any position along the overlap, the slipping displacement is written such as:

𝑑Δ𝑢 𝑑𝑥 = (1 + 𝜒 𝐴 ) 𝑁 2 𝐴 2 -𝜒 𝐴 𝐹 𝐴 2 + (1 - 1 𝜒 𝛼 ) 𝛼 2 Δ 𝑇 (22)
where the adherend thermal unbalance is characterized by:

𝜒 𝛼 = 𝛼 2 𝛼 1 (23)
As a result, from the application of boundary conditions in Eq. ( 20) and Eq. ( 21), both last required equations are obtained.

∑ (𝑛 + 1)𝑈 𝑛+1 ∞ 𝑛=0 = 𝑐 𝐹 𝐴 2 + 𝑐 (1 - 1 𝜒 𝛼 ) 𝛼 2 Δ 𝑇 (24)
∑ (𝑛 + 1)(-1) 𝑛 𝑈 𝑛+1 ∞ 𝑛=0 = -𝑐𝜒 𝐴 𝐹 𝐴 2 + 𝑐 (1 - 1 𝜒 𝛼 ) 𝛼 2 Δ 𝑇 (25)
These last equations allow for the introduction of the thermal load in the model. The solution is finally obtained by truncation of the series at an order n_max. Equations Eq. ( 19), Eq. ( 24)

and Eq. ( 25) allow for the writing of a linear system, the size of which is (n_max+1)², which can be solved using a mathematical programming software such as MATLAB.

ME resolution

Firslty, the elementary stiffness matrix of a BBa element is formulated. The length of the BBa is , on which the material and geometrical properties are supposed constant. The element reference system of axis is denoted (X,Y,Z). The elementary stiffness matrix of the BBa, termed K BBa , describes the interaction between the four nodal forces and the force nodal displacements (see Figure 4), such as:

( -𝑁 1 (0) -𝑁 2 (0) 𝑁 1 (Δ) 𝑁 2 (Δ) ) = 𝐾 𝐵𝐵𝑎 ( 𝑢 1 (0) 𝑢 2 (0) 𝑢 1 (Δ) 𝑢 2 (Δ) ) ⟺ 𝐹 𝑒 = 𝐾 𝐵𝐵𝑎 𝑈 𝑒 ( 26
)
where F e (U e ) is the nodal force (displacement) vector of the BBa element. In the frame of the 1D-bar analysis, the closed-form expressions for each component of K BBa can be obtained [START_REF] Paroissien | Contribution aux Assemblages Hybrides (Boulonnés/Collés) -Application aux Jonctions Aéronautiques[END_REF][START_REF] Paroissien | Hybrid (bolted/bonded) joints applied to aeronautic parts: Analytical one-dimensional models of a single-lap joint[END_REF]. Even if the mathematical description has already published under another shape, it provided in Appendix A.

As expected, the obtained stiffness matrix is the same as the one obtained without considering the thermal load. The method to take into account a linear variation of shear stress in the adherend thickness following [START_REF] Tsaï | Improved Theoretical Solutions for Adhesive Lap Joints Int[END_REF] is described in Appendix D.

Since the material properties of the adhesive vary along the overlap, the approach using the ME technique consists in regularly meshing the overlap with BBa elements with n_ME BBa elements (see Figure 5). Each BBa element has a length =L/n_ME, on which the material to be approximated by a stepped function. In this work, it is supposed that the adhesive shear relative stiffness for the n th BBa is equal to its value at the centre of the BBa. The elementary stiffness matrices of the each BBa can then be assessed. The bars outside the overlap are simulated as bar elements associated with the elementary stiffness matrix K bar,j for adherend j given in Appendix A. From the elementary stiffness matrices, the stiffness matrix of the joint, termed K s , is then built. The boundary conditions are applied using the classical FE rules. In particular, the thermal load is replaced by an equivalent mechanical one, such that the equivalent nodal force vector provides the same nodal displacement as the thermal load. The potential energy leading to the classical linear system F s =K s U s , the size of which is (2n_ME+4)², where F s is the nodal force of the structure and U s the nodal vector.

Figure 5. ME model of the FGA single lap joint.

Under 1D-beam kinematics 2.1.1 Hypotheses

The model is based on the following hypotheses: (i) the thickness of the adhesive layer is constant along the overlap, (ii) the adherends are simulated by linear elastic Euler-Bernoulli laminated beams and (iii) the adhesive layer is simulated by an infinite number of elastic shear and transverse springs linking both adherends. Similarly to the 1D-bar analysis, the case of a single-lap joint subjected to combined mechanical and thermal loads is considered, for which the geometrical parametrization is provided in Figure 6. The joint is simply supported at both extremities and submitted to a uniform variation of temperature  T and to a tensile force F. It is indicated that any boundary conditions could be applied when the ME technique is applied, even if simply supported is chosen in this paper. The stress analysis is conducted in force but could similarly be made in tensile flow F/b.

Governing Equation

The local equilibrium selected for the formulation of the presented BBe element is related to the one used by Luo and Tong [START_REF] Luo | Fully-coupled nonlinear analysis of single lap adhesive joints[END_REF] and allows for a coupling between the in-plane and out-ofplane load. Moreover, the formulation presented can be easily modified to correspond to the Goland and Reissner [START_REF] Goland | The stresses in cemented joints[END_REF] or to the Hart-Smith [START_REF] Hart-Smith | Adhesive bonded single lap joints[END_REF] local equilibrium.

The local equilibrium of both adherends (see Figure 7) provides the six following equations:

𝑑𝑁 𝑗 𝑑𝑋 = (-1) 𝑗 cos 𝜃 𝑗 𝑏𝑇, 𝑗 = 1,2 (27)
𝑑𝑉 𝑗 𝑑𝑋 = (-1) 𝑗+1 𝑏𝑆 + (-1) 𝑗 sin 𝜃 𝑗 𝑏𝑇, 𝑗 = 1,2

𝑑𝑀 𝑗 𝑑𝑋 + 𝑉 𝑗 + cos 𝜃 𝑗 𝑏 (ℎ 𝑗 + 𝑒 𝑎 2 ) 𝑇 -sin 𝜃 𝑗 𝑁 𝑗 = 0, 𝑗 = 1,2 (28)
ℎ 𝑗 = 𝑒 𝑗 2 , 𝑗 = 1,2 (30)
where V j is the shear force in the adherend j, M j the bending moment in the adherend j,  j the bending angle in the adherend j and S is the adhesive peel stress. The constitutive equations can then be written as:

𝑁 𝑗 = 𝐴 𝑗 𝑑𝑢 𝑗 𝑑𝑋 -𝐵 𝑗 𝑑𝜃 𝑗 𝑑𝑋 -𝑁 𝑗 Δ 𝑇 , 𝑗 = 1,2 (31)
𝑀 𝑗 = -𝐵 𝑗 𝑑𝑢 𝑗 𝑑𝑋 + 𝐷 𝑗 𝑑𝜃 𝑗 𝑑𝑋 + 𝑀 𝑗 Δ 𝑇 , 𝑗 = 1,2 (32)
𝜃 𝑗 = 𝑑𝑣 𝑗 𝑑𝑋 (33)
where (see Appendix A) A j is the membrane stiffness of adherend j, B j the coupling membrane-bending stiffness of adherend j, D j the bending stiffness of adherend j, 𝑁 𝑗 Δ 𝑇 the

X Y + N 1 +dN 1 M 1 +dM 1 V 1 +dV 1 V 2 +dV 2 N 2 +dN 2 M 2 +dM 2 M 1 V 1 V 2 M 2 N 1 N 2 S.bdx T.bdx T.bdx  1  1  2  1 e a /2
thermal normal force in the adherend j and 𝑀 𝑗 Δ 𝑇 is the thermal bending moment in the adherend j. In the case of a lay-up characterized by a mirror symmetry, B j =0 and 𝑀 𝑗 Δ 𝑇 = 0.

The constitutive equations for the adhesive layer are provided by:

𝑆 = 𝐸 𝑎 𝑒 𝑎 [𝑣 1 -𝑣 2 ] = 𝑘 𝐼 Δ𝑣 (34)
𝑇 = 𝐺 𝑎 𝑒 𝑎 [𝑢 2 -ℎ 2 𝜃 2 -(𝑢 1 + ℎ 1 𝜃 1 )] = 𝑘 𝐼𝐼 Δ𝑢 (35)
with:

Δ𝑢 = 𝑢 2 -𝑢 1 -ℎ 2 𝜃 2 -ℎ 1 𝜃 1 ( 36
)
Δ𝑣 = 𝑣 1 -𝑣 2 (37)
where E a is the adhesive peel modulus and k I =E a /e a the adhesive peel relative stiffness. v is representative of the opening displacement of the adherend interface. Contrary to the 1D-bar analysis, the presented analysis cannot be directly applied when the thickness of the adhesive layer varies along the overlap. Indeed, the variation of the thickness induces a lag of the neutral axis, which has to be taken into account because of the deflection. Nevertheless, it is indicated that the variation of the adhesive thickness could be easily taken into account when using the ME technique.

ME resolution

The resolution scheme follows the same part as for the 1D-bar analysis (see section 2.1.4).

The single-lap joint is meshed in BBe elements along the overlap and beam elements for the parts outside the overlap. The boundary conditions, the mechanical and thermal loads are then applied. Similarly to the 1D-bar analysis, the thermal load is applied under the shape of an equivalent nodal force vector given by:

𝐹 𝑡ℎ = ( -𝑁 1 Δ 𝑇 -𝑁 2 Δ 𝑇 𝑁 1 Δ 𝑇 𝑁 2 Δ 𝑇 0 0 0 0 𝑀 1 Δ 𝑇 𝑀 2 Δ 𝑇 -𝑀 1 Δ 𝑇 -𝑀 2 Δ 𝑇 ) (38)
Contrary to the 1D-bar analysis, it is not possible to simply obtain closed-form expressions for the components of the stiffness matrix of the BBe element. An approach for the formulation of the stiffness matrix of BBe element under Goland and Reissner equilibrium has already been described in detail in previous papers [START_REF] Paroissien | Contribution aux Assemblages Hybrides (Boulonnés/Collés) -Application aux Jonctions Aéronautiques[END_REF][START_REF] Paroissien | Hybrid (bolted/bonded) joints applied to aeronautic parts: Analytical one-dimensional models of a single-lap joint[END_REF][START_REF] Paroissien | Analytical two-dimensional model of a hybrid (bolted/bonded) single-lap joint[END_REF][START_REF] Paroissien | A simplified stress analysis of bonded joints using macro-elements[END_REF][START_REF] Paroissien | Elasto-Plastic Analysis of Bonded Joints with Macro-Elements[END_REF][START_REF] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF][START_REF] Lélias | Mechanical behavior of adhesively bonded joints: Modeling, simulation and experimental characterization[END_REF][START_REF] Paroissien | Simplified Stress Analysis of Hybrid (Bolted/Bonded) Joints[END_REF]. Nevertheless, this approach could be long to set up. In this paper, a new approach is provided for a fast and easy implementation within mathematical software such as MATLAB for example. The present formulation ME has never been published. The element reference system (X,Y,Z) of axes is considered. Following Luo and Tong approach [START_REF] Luo | Fully-coupled nonlinear analysis of single lap adhesive joints[END_REF], a first approximation is made. The bending angle is supposed very small. The six local equilibrium equations become then:

𝑑𝑁 𝑗 𝑑𝑋 = (-1) 𝑗 𝑏𝑇, 𝑗 = 1,2 (39)
𝑑𝑉 𝑗 𝑑𝑋 = (-1) 𝑗+1 𝑏𝑆 + (-1) 𝑗 𝜃 𝑗 𝑏𝑇, 𝑗 = 1,2 (40)
𝑑𝑀 𝑗 𝑑𝑋 + 𝑉 𝑗 + 𝑏 (ℎ 𝑗 + 𝑒 𝑎 2 ) 𝑇 -𝜃 𝑗 𝑁 𝑗 = 0, 𝑗 = 1,2 (41)
A second approximation is made. It consists in neglecting the product of the adhesive shear stress with the bending angle 𝑇𝜃 𝑗 ≪ 1, 𝑗 = 1,2. The six local equilibrium equations become then:

𝑑𝑁 𝑗 𝑑𝑋 = (-1) 𝑗 𝑏𝑇, 𝑗 = 1,2 (42)
𝑑𝑉 𝑗 𝑑𝑋 = (-1) 𝑗+1 𝑏𝑆, 𝑗 = 1,2 (43)
𝑑𝑀 𝑗 𝑑𝑋 + 𝑉 𝑗 + 𝑏 (ℎ 𝑗 + 𝑒 𝑎 2 ) 𝑇 -𝜃 𝑗 𝑁 𝑗 = 0, 𝑗 = 1,2 (44)
Compared to the local equilibrium by Hart-Smith [START_REF] Hart-Smith | Adhesive bonded single lap joints[END_REF] only the bending moment is modified, involving a coupling between normal forces and bending moment. The following quotation is introduced for any functions f:

( 𝑓 + 𝑓 - ) = ( 1 1 1 -1 ) ( 𝑓 1 𝑓 2 ) (45)
A third and last approximation is made which reads

𝑁 - 2 (𝜃 1 ∓ 𝜃 2 ) ≪ 1.
Under this approximation and taking into account that the sum of normal forces at any abscissa is equal to the applied force F, a system of twelve first order linear ODEs is obtained:

𝑑𝑢 + 𝑑𝑋 = 1 2 ( 𝐷 1 Δ 1 + 𝐷 2 Δ 2
) 𝑁 + +

) 𝑁 -+

) 𝑀 + +

) 𝑀 -

𝑑𝑣 + 𝑑𝑋 = 𝜃 + (46)
𝑑𝜃 + 𝑑𝑋 = 1 2 ( 𝐵 1 Δ 1 + 𝐵 2 Δ 2 ) 𝑁 + + 1 2 ( 𝐵 1 Δ 1 - 𝐵 2 Δ 2 ) 𝑁 -+ 1 2 ( 𝐴 1 Δ 1 + 𝐴 2 Δ 2 ) 𝑀 + + 1 2 ( 𝐴 1 Δ 1 - 𝐴 2 Δ 2 ) 𝑀 - ( (47)
) 𝑑𝑢 - 𝑑𝑋 = 1 2 ( 𝐷 1 Δ 1 - 𝐷 2 Δ 2 ) 𝑁 + + 1 2 ( 𝐷 1 Δ 1 + 𝐷 2 Δ 2 ) 𝑁 -+ 1 2 ( 𝐵 1 Δ 1 - 𝐵 2 Δ 2 ) 𝑀 + + 1 2 ( 𝐵 1 Δ 1 + 𝐵 2 Δ 2 ) 𝑀 - ( 48
) 𝑑𝑣 - 𝑑𝑋 = 𝜃 - ( 49
) 𝑑𝜃 - 𝑑𝑋 = 1 2 ( 𝐵 1 Δ 1 - 𝐵 2 Δ 2 ) 𝑁 + + 1 2 ( 𝐵 1 Δ 1 + 𝐵 2 Δ 2 ) 𝑁 -+ 1 2 ( 𝐴 1 Δ 1 - 𝐴 2 Δ 2 ) 𝑀 + + 1 2 ( 𝐴 1 Δ 1 + 𝐴 2 Δ 2 ) 𝑀 - 50
where  j =A j D j -B j B j ≠0. By letting F=0, the stress analysis of the sandwich by Hart-Smith is deduced [START_REF] Hart-Smith | Adhesive bonded single lap joints[END_REF]. In addition, by letting e a =0, it corresponds to the one by Goland and Reissner [START_REF] Goland | The stresses in cemented joints[END_REF]. This system can be written as 𝑑𝑉 𝑑𝑋 = 𝐴𝑉 where A is 12x12 matrix with real constant components and the unknown vector

V is such that t V=(u 1 u 2 v 1 v 2  1  2 N 1 N 2 V 1 V 2 M 1 M 2 ).
But the elementary stiffness matrix corresponds to the relationship between the vector of nodal forces and the vector of nodal displacements, such as:

( -𝑁 1 (0) -𝑁 2 (0) 𝑁 1 (Δ) 𝑁 2 (Δ) -𝑉 1 (0) -𝑉 2 (0) 𝑉 1 (Δ) 𝑉 2 (Δ) -𝑀 1 (0) -𝑀 2 (0) 𝑀 1 (Δ) 𝑀 2 (Δ) ) = 𝐾 𝐵𝐵𝑒 ( 𝑢 1 (0) 𝑢 2 (0) 𝑢 1 (Δ) 𝑢 2 (Δ) 𝑣 1 (0) 𝑣 2 (0) 𝑣 1 (Δ) 𝑣 2 (Δ) 𝜃 1 (0) 𝜃 2 (0) 𝜃 1 (Δ) 𝜃 2 (Δ)) (58)
The fundamental matrix of A, termed  A , is computed at X=0 and X=; using the MATLAB software, the associated command is "expm":

{ Φ 𝐴 (𝑋 = 0) = 𝑒𝑥𝑝𝑚(𝐴. 0) Φ 𝐴 (𝑋 = Δ) = 𝑒𝑥𝑝𝑚(𝐴. Δ) (59)
From these two 12*12 matrices, two matrices M' and N' are extracted. M' (N') is composed of the lines related to the nodal displacements (forces). For each, a first block of six lines and twelve rows comes from  A (X=0) and the second block of six lines and twelve rows comes from  A (X=), such that:

{ 𝑀′ = Φ 𝑈 (0, Δ) = ( [Φ 𝐴 (𝑋=0)] 𝑖=1,
where i (j) indicates the line (row) number. As K BBe is defined according to ([u 1 (0) u 2 (0) u 1 ()

u 2 () v 1 (0) v 2 (0) v 1 () v 2 ()  1 (0)  2 (0)  1 ()  2 ()]
), a simple rearrangement of the order of lines of M' is performed to produce the matrix M. Similarly, the matrix N' is subjected to the same operation. In a similar way, the terms related to nodal forces at X=0 are multiplied by -1

to follow the arrangement ([-N 1 (0) -N 2 (0) N 1 () N 2 () -V 1 (0) -V 2 (0) V 1 () V 2 () - 1 (0) - 2 (0)  1 ()  2 ()]
). It leads to the definition of the matrix N. The elementary stiffness matrix K BBe is equal to the product of N and the inverse of M: K BBe =N.M -1 .The stiffness matrix of the beam element under a local equilibrium coupling the in-plane and out-of-plane load is described in Appendix C. As for the 1D-bar analysis, the minimization of the potential energy leads to a linear system the size of which is (6n_ME+12)².

Even if it is not the topic of this paper, it is obvious that this previous approach can be easily used to develop ME including different number of layers of adhesives and adherends (e.g.

double lap joint configuration), various beam models (e.g. Timoshenko beam model, see Appendix D) or taking into account for a linear variation of shear stress in the adherend thickness following [START_REF] Tsaï | Improved Theoretical Solutions for Adhesive Lap Joints Int[END_REF] (see Appendix D).

Comparative study

Overview

A comparative study of the ISLM by Stein et al. [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF], the present 1D-bar TEPS, 1D-bar ME and 1D-beam ME analysis is presented in this section, starting with a convergence study. This study is performed against three test cases (TCs):

(i) TC#1: a balanced joint configuration under a pure mechanical load;

(ii) TC#2: an unbalanced joint configuration under a pure thermal load;

(iii) TC#3: an unbalanced joint configuration under combined mechanical and thermal loads.

The joint configurations are almost inspired from to those found in [START_REF] Carbas | Modelling of functionally graded adhesive joints[END_REF][START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF]. The mechanical load is F=5 kN while the thermal load is  T =+50°K. The balanced joint configuration is made of two steel adherends. The unbalanced joint configuration has the same geometry as the balanced one, but the lower adherend is made in aluminum instead of steel. The geometrical and mechanical parameters are given in Table 1 and Table 2 respectively in accordance to

where E a,max (E a,min ) is the maximal (minimal) adhesive peel modulus in the graduation and G a,max (G a,min ) is the maximal (minimal) adhesive shear modulus in the graduation. In this work, the ratio between the maximal (minimal) adhesive peel modulus and the maximal (minimal) adhesive shear modulus through is constant and equal to 2(1+ a ), where  a is the adhesive Poisson's ratio. In this work, the adhesive peel modulus is then represented by the adhesive Young's modulus. The adhesive properties are then summarized in Table 3.

Convergence study

The convergence study is performed on the TC#1 and the TC#2 (FGA balanced joint under pure mechanical load and pure thermal load). The resolution scheme based on TEPS needs to truncation order (n_max) while the one based on the ME technique needs a mesh with n_ME BBa or BBe. A convergence study is then performed by recording the maximal adhesive stresses as a function n_max and n_ME.

The maximal adhesive shear stress (T max ) is provided in Figure 8 as function of the order of series truncation (n_max) for both TC#1 and TC#2 as predicted by the ISLM and the TEPS analysis. It is shown that T max tends to a finite value (12.56 MPa for TC#1 and 11.11 MPa for TC#2) for an order of series truncations lower than n_max=20. For the case of n_max=100

with the TEPS analysis, the ratio between each series term U n with the sum of series termswhich is equal to u(x=c)is provided in Figure 9 for both TC#1 and TC#2, illustrating the fast convergence of series. Moreover, the TEPS analysis provides a maximal adhesive shear stress relatively different of 1.68E-4 % from the one provided by the ISLM (for TC#1).

Figure 8. Maximal adhesive shear stress as function the order of series truncation for the ISLM by Stein et al. [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF] and the 1D-bar TEPS analysis for both TC#1 and TC#2. The relative difference in the maximal adhesive provided the 1D-bar ME analysis from the one by 1D-bar TEPS analysis as function of the number of MEs for both TC#1 and TC#2 is provided in Figure 10. It is shown that T max provided by the 1D-bar ME analysis tends to the one provided by the TEPS analysis when the number of MEs is increased. For n_ME=1000, the relative difference is 0.16% (0.22%) for TC#1 (TC#2). As a result, the TEPS resolution scheme is less costly in terms of CPU time than the ME one, since convergence is obtained at a lower size of the linear system to be inverted. This behavior is related to the meshing strategy associated with the graduation of adhesive properties. It is thought that the number of MEs could be reduced by adapting the length of each ME according to the current gradient of adhesive properties for example. However, the mesh optimization is not the topic of this paper. In Figure 11, the maximal adhesive shear stress provided by the 1D-bar and 1D-beam analysis as function the order of the number of MEs for both TC#1 and TC#2 is provided. As expected from Figure 10, it is shown that T max tends to a finite value.

Elements of validation

The ME technique is a particular resolution scheme allowing for the system of ODEs coming from simplifying hypotheses on which various modelssuch as Volkersen, Goland and Reissner, Hart-Smith, Luo and Tongare based. It was shown in that, for bonded joints with constant adhesive properties under mechanical or thermal loads, the predictions from the models using the ME resolution scheme provide the same results as those provided by the related reference models [START_REF] Paroissien | Contribution aux Assemblages Hybrides (Boulonnés/Collés) -Application aux Jonctions Aéronautiques[END_REF][START_REF] Paroissien | Hybrid (bolted/bonded) joints applied to aeronautic parts: Analytical one-dimensional models of a single-lap joint[END_REF][START_REF] Paroissien | Analytical two-dimensional model of a hybrid (bolted/bonded) single-lap joint[END_REF][START_REF] Paroissien | A simplified stress analysis of bonded joints using macro-elements[END_REF][START_REF] Paroissien | Elasto-Plastic Analysis of Bonded Joints with Macro-Elements[END_REF]. In other words, the same hypotheses lead to the same results. Moreover, it was shown that the predictions from the ME analysis are in close agreements with those from FE models built on bar or beam element linked by peel and/or shear springs, under mechanical and/or thermal loads, involving the update of adhesive properties for each ME to take into account for nonlinear adhesive material behaviors [START_REF] Paroissien | A simplified stress analysis of bonded joints using macro-elements[END_REF][START_REF] Paroissien | Elasto-Plastic Analysis of Bonded Joints with Macro-Elements[END_REF].

These FE models were developed to be the most representative for the ME analysis in order to validate the codes. As a result, the ME resolution scheme provide validated predictions relevant to the simplifying hypotheses. Similarly, the TEPS resolution scheme allows for the resolution of differential equations related to the simplifying equations. It was validated and assessed in the case of FGA single-lap joints under mechanical load by Stein et al. [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF].

In addition, the stress distributions at constant maximal and minimal adhesive properties are then provided from the use of ME models in the following sections. An order of truncation equal to 100 and a number of MEs equal to 500 is chosen in the following sections.

Considering the balanced joint configuration in the frame of the 1D-bar analysis, the adhesive shear stress distributions along the overlap are provided in Figure 12. The predictions considering homogeneous shear modulus G a,min and G a,max , as well as those from ISLM by Stein et al. [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF], the present 1D-bar TEPS and ME analyses are included. It appears that the predictions from the ISLM, the TEPS analysis and 1D-bar ME analysis provide the same predictions. Considering the balanced joint configuration in the frame of the 1D-beam analysis, the adhesive shear and peel stress distributions along the overlap are provided in Figure 13 and

Test cases

1D-bar analyses

The adhesive shear stress distributions along the overlap are provided in Figure 15 and Figure 16 for the 1D-bar analyses for TC#2 and TC#3 (unbalanced joint configuration under combined mechanical and thermal loads), respectively. The ISLM cannot then be applied. It is shown that the predictions by the 1D-bar TEPS and ME analyses are superimposed. For each case, the graduation of adhesive properties allows to reduce the peak stresses below those obtained from the case at constant minimal shear modulus. However, the reduction is less pronounced for the unbalanced cases (TC#2 and TC#3). For the TC#1, the reduction in adhesive shear peak stress at x=c of the FGA joint is -21.7% from the bonded joints with a constant shear modulus G a,min , while it is -14.5% (-15.3%) for TC#2 (TC#3).

1D-beam analyses

In the frame of 1D-beam analysis, the adhesive shear and peel stress distributions along the overlap are provided in Figure 17 to Figure 18, and in Figure 19 to Figure 20, for TC#2 and TC#3 respectively. The predictions come from the 1D-beam ME analysis only since the GM cannot be applied for these TCs. The adhesive shear stress distribution along the overlap of FGA joints from the 1D-bar analysis differs significantly from the one from 1D-beam analysis. As for the 1D-bar analysis, the reduction of adhesive peak stresses is shown while the stress distribution of FGA joints appears close to the one of the bonded joints with minimal and constant adhesive modulus. For TC#1, the reduction in adhesive shear stress at x=c for the FGA joint is -13.6% from the bonded joints with a constant minimal Young and shear modulus, while it is -10.0% (-9.61%) for TC#2 (TC#3). Similarly, in terms of peel stress at x=c, the reduction is -4.85% (-2.59%) for TC#1 (TC#3).

Reduction of adhesive peak stresses

This section aims at illustrating how the use of the 1D-beam ME analyses could help in the design of adhesive graduation to reduce the adhesive peak stresses, for the unbalanced joint configuration subjected to pure thermal load and combined mechanical and thermal loads in particular. According to Hart-Smith [START_REF] Hart-Smith | Design methodology for bonded-bolted composite joints[END_REF][START_REF] Hart-Smith | Adhesive bonded single lap joints[END_REF], the adhesive peel stress could lead to an anticipated failure of single-lap bonded joint, whereas the potential of shear deformation is not reached. For the unbalanced joint configuration, under a pure thermal load, it is shown that the level of adhesive peel stresses remain very low (see Figure 18), while the adhesive shear stress are symmetrical in absolute value (see Figure 17). Under combined mechanical and thermal loads, the level of adhesive peel stress is significantly higher due to the introduction of the mechanical load inducing a bending moment at both overlap ends (see Figure 20). Moreover, even if the adhesive peel stress distribution is asymmetrical due to the

The increase of the parameter p allow for the enlargement of the overlap length at higher modulus, while increasing the graduation slope at both overlap ends where the adhesive stress gradients are higher. The shape of various graduations is illustrated in Figure 21. stress is the highest, the reductions obtained become -5.22%, -8.42% and -10.6% with p=2, p=3 and p=4 respectively. The choice of adhesive graduation law associated higher power order allows then for a lag of the adhesive shear peak stress in direction of the center of the overlap.

Figure 22. Adhesive shear stress distribution along the overlap for the unbalanced joint configuration under a pure thermal load, for various adhesive graduations.

- Under a combined mechanical and thermal load, the adhesive shear and peel stress distributions along the overlap are provided in Figure 24 to Figure 25, respectively. As for the pure thermal load case, the adhesive peak shear stresses are not located at both overlap ends for p=2, p=3 and p=4. Relatively to the adhesive shear peak stresses with p=1 a reduction of -4.80%, -5.39% and -3.18% is obtained with p=2, p=3 and p=4 respectively. Relatively to the adhesive stress at overlap end, where the level of adhesive peel stress is the highest, the reductions obtained become -4.80%, -7.76% and -9.77% with p=2, p=3 and p=4 respectively.

The adhesive peak peel stresses are located at both overlap ends; the maximal peak is located at x=c. It is shown that with p=2, p=3 and p=4 respectively, both peaks at overlap ends are reduced. Relatively to the adhesive peel peak stresses at x=c with p=1, a reduction of -2.87%, -4.91% and -6.45% is obtained with p=2, p=3 and p=4 respectively  the present 1D-bar TEPS analysis restricted to a pure mechanical load provide the same results as the ISLM by Stein et al [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF];

 the present 1D-bar TEPS and the 1D-bar ME analysis provide the same results;

 the use of TEPS resolution scheme provides converged results at lower CPU cost than the ME resolution scheme;

 the graduation of the adhesive properties allows for the reduction of adhesive peak stresses;

 the present 1D-beam ME analysis restricted to a pure mechanical load provide similar results as the GM by Stein et al [START_REF] Stein | Analytical models for functionally graded adhesive joints: A comparative study[END_REF];

 the reduction of the adhesive shear peak stresses is found less pronounced when the 1Dbeam analysis is used instead of the 1D-bar analysis;

 the reduction of the adhesive peak stresses in less pronounced for an unbalanced joint than for a balance joint.

 higher level of reduction can be obtained by modifying the graduation law.

A dedicated validation campaign based on FE modelling should be undertaken in order to assess the relevance of the simplifying hypotheses and the performance of the resolution scheme for the stress analysis of FGA joints. In particular, the free stress state at both overlap ends cannot be captured with the simplifying hypotheses used. Besides, optimization processes could be used to optimize the graduation of adhesive properties as function of the adhesive stress distribution to minimize the adhesive peak stresses. Finally, in order to increase the strength of FGA single-lap joints, an idea could be to graduate the properties of both the adhesive and adherends. For example, the reduction of adhesive peel stresses at both overlap ends could be obtain by tapering the adherend edge, while increasing the ratio between the overlap length and the adherend thickness [START_REF] Hart-Smith | Design methodology for bonded-bolted composite joints[END_REF]. The elementary stiffness matrix is finally computed from the matrix M e and N e : The elementary stiffness matrix of the bar element, simulating the adherend j outside the overlap is:

𝐹
𝐾
𝐾 𝑏𝑎𝑟,𝑗 = 𝐴 𝑗 ( 1 -1 - 1 1
) , 𝑗 = 1,2 (A-20)

Appendix B

This appendix provides the derivation of the constitutive equations of laminated beams used in the 1D-beam analysis, in the (X,Y i ,Z) reference local axis of the adherend, the height origin of which is taken on the neutral axis. The normal force and the bending moment are written such as:

𝑁 𝑖 (𝑋) = ∫

  A j extensional stiffness (N) of adherend j B j extensional and bending coupling stiffness (N.mm) of adherend j D j bending stiffness (N.mm 2 ) of adherend j E a adhesive peel modulus (MPa) E a,min adhesive shear modulus (MPa) E a,max adhesive shear modulus (MPa) E j adherend Young's modulus (MPa) of adherend j F magnitude of applied force (N) F e element nodal force vector F e,therm element nodal force vector equivalent to thermal load G a adhesive shear modulus (MPa) G a,max maximal adhesive shear modulus (MPa) G a,min minimal adhesive shear modulus (MPa) K BBa elementary stiffness matrix of a bonded-bars element K BBe elementary stiffness matrix of a bonded-beams element K bar,j elementary stiffness matrix of a bar for the adherend j L length (mm) of bonded overlap M e element matrix linking the element nodal displacement to the constant integration vector M j bending moment (N.mm) in adherend j around the z direction M e element matrix linking the element nodal force to the constant integration vector 𝑀 𝑗 Δ 𝑇 thermal bending moment (N.mm) in adherend j around the z direction N j normal force (N) in adherend j in the x direction 𝑁 𝑗 Δ 𝑇 thermal normal force (N) in adherend j in the x direction S adhesive peel stress (MPa) T adhesive shear stress (MPa) T max maximal adhesive shear stress (MPa) U e element nodal displacement vector V j shear force (N) in adherend j in the y direction b width (mm) of the adherends c half-length (mm) of bonded overlap e a thickness (mm) of the adhesive layer h j half thickness (mm) of adherend j k I adhesive elastic stiffness (MPa/mm) in peel k II adhesive elastic stiffness (MPa/mm) in shear n_max order of truncation n_ME number of macro-elements p power of the graduation law u j displacement (mm) of adherend j in the x direction v j displacement (mm) of adherend j in the y direction  overlap length (mm) of a macro-element   T variation of temperature (K) u slipping displacement (mm)  j characteristic parameter (N 2 .mm 2 ) of adherend j   j coefficient of thermal expansion (K -1 ) of adherend j  j bending angle (rad) of the adherend j around the z direction 𝜒 𝐴 adherend stiffness unbalance parameter (-)

Figure 1 .

 1 Figure 1. Modelling of a bonded overlap by a macro-element.

2 .

 2 Description of simplified stress analyses of FGA single-lap joints 2.1. Under 1D-bar kinematics 2.1.1 Hypotheses

Figure 2 .

 2 Figure 2. Geometrical parametrization of the single-lap joint, boundary conditions and applied loads (1D-bar analysis).

Figure 3 .

 3 Figure 3. Free body diagram of infinitesimal pieces included between x and x+dx of both

F

Figure 4 .

 4 Figure 4. Free body diagram of infinitesimal pieces included between x and x+dx of both

  constant. The actual graduation of the adhesive shear relative stiffness has then

Figure 6 .

 6 Figure 6. Geometrical parametrization of the single-lap joint, boundary condition and applied loads (1D-beam analysis).

F

Figure 7 .

 7 Figure 7. Free body diagram of infinitesimal pieces included between x and x+dx of both

  𝑏. 𝑢 -+ 𝑘 𝐼𝐼 𝑏(ℎ + + 𝑒 𝑎 )𝜃 + + 𝑘 𝐼𝐼 𝑏ℎ -𝜃 - + + 𝑒 𝑎 )ℎ -) 𝜃 + + (

Figure 2 and

 2 Figure2and Figure6. A parabolic graduation of adhesive properties is assumed such as:

Figure 9 .

 9 Figure 9. Ratio between each series term and the sum of series terms for both TC#1 and

Figure 10 .

 10 Figure 10. Relative difference in % in the maximal adhesive provided the 1D-bar ME analysis

Figure 11 .

 11 Figure 11. Maximal adhesive shear stress provided by the 1D-bar and 1D-beam analysis as

Figure 12 .

 12 Figure 12. Adhesive shear stress distribution along the overlap associated with TC#1 for the

Figure 14 ,

 14 Figure14, respectively. The predictions considering homogeneous shear and peel modulus

Figure 13 .Figure 14 .

 1314 Figure 13. Adhesive shear stress distribution along the overlap associated with TC#1 for the

Figure 15 .

 15 Figure 15. Adhesive shear stress distribution along the overlap associated with TC#2 for the

Figure 16 .

 16 Figure 16. Adhesive shear stress distribution along the overlap associated with TC#3 for the

Figure 17 .

 17 Figure 17. Adhesive shear stress distribution along the overlap associated with TC#2 for the

Figure 18 .

 18 Figure 18. Adhesive shear stress distribution along the overlap associated with TC#2 for the

Figure 19 .

 19 Figure 19. Adhesive peel stress distribution along the overlap associated with TC#3 for the

Figure 20 .

 20 Figure 20. Adhesive peel stress distribution along the overlap associated with TC#3 for the

  ME peel stress at Ea,min and Ga,min 1D-beam ME peel stress at Ea,max and Ga,max [TC#3] 1D-beam ME peel stress graded unbalance of the joint, the adhesive peel peaks stresses located at both overlap ends are significant. As a result, for both previous load cases, a symmetrical adhesive graduation is kept in order to try to reduce the adhesive peak stresses. It is assumed to follow a symmetrical power law parametrized by p, with p=1,2,3,4, such as: 𝐸 𝑎 (𝑥) = 𝐸 𝑎,𝑚𝑎𝑥 (𝑥) -(𝐸 𝑎,𝑚𝑎𝑥 (𝑥) -𝐸 𝑎,𝑚𝑖𝑛 (𝑥)) ( 𝑥) = 𝐺 𝑎,𝑚𝑎𝑥 (𝑥) -(𝐺 𝑎,𝑚𝑎𝑥 (𝑥) -𝐺 𝑎,𝑚𝑖𝑛 (𝑥)) (

Figure 21 .

 21 Figure 21. Adhesive peel modulus along the overlap as a function of p.

Figure 23 .

 23 Figure 23. Adhesive peel stress distribution along the overlap for the unbalanced joint

Figure 24 .

 24 Figure 24. Adhesive shear stress distribution along the overlap for the unbalanced joint

Figure 25 .

 25 Figure 25. Adhesive peel stress distribution along the overlap for the unbalanced joint

Table 1

 1

		. Geometrical parameters of joint configurations
	b (mm)	e a (mm)	e 1 =e 2 (mm)	L (mm)	l 1 =l 2 (mm)
	25	0.2	2	25	75

Table 2 .

 2 Material parameters of adherends.

		Coefficient of thermal expansion (K -1 ) Young's modulus (GPa)
	steel	12E-6	210
	aluminum	24E-6	70

Table 3 .

 3 Adhesive material properties.

	E a,max (MPa)	E a,min (MPa)	 a

  𝜎 𝑖 𝑏𝑑𝑌 𝑖 the adherend i n i is the number of layers and h pi is the final height of the p i th layer.Moreover, the orthotopic behavior provides𝜎 𝑖 𝑝 𝑖 = 𝑄 𝑖 𝑝 𝑖 (𝜀 𝑖 𝑝 𝑖 -𝛼 𝑖 𝑝 𝑖 Δ 𝑇 ), 𝑖 = 1,2 (B-3)where, in the adherend I, 𝑄 𝑖 𝑝 𝑖 is the matrix of reduced stiffness in the p i th layer.As a result, the normal force and the bending moment are given by:

	+ℎ 𝑖 -ℎ 𝑖 𝑀 𝑖 (𝑋) = ∫ -𝑌 𝑖 𝜎 𝑖 𝑏𝑑𝑌 𝑖 = 𝑏 ∑ +ℎ 𝑖 -ℎ 𝑖 = -𝑏 ∑ ∫ ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 𝑛 𝑖 𝑝 𝑖 =1 𝑛 𝑖 𝑝 𝑖 =1 ∫ 𝑄 𝑖 𝑝 𝑖 (𝜀 𝑖 𝑝 𝑖 -𝛼 𝑖 𝑝 𝑖 Δ 𝑇 )𝑑𝑌 𝑖 𝜎 𝑖 𝑝 𝑖 𝑑𝑌 𝑖 ∫ 𝜎 𝑖 , 𝑖 = 1,2 𝑝 𝑖 𝑌 𝑖 𝑑𝑌 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 , 𝑖 = 1,2 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 𝑛 𝑖 𝑝 𝑖 =1 , 𝑖 = 1,2 𝑀 𝑖 (𝑋) = -𝑏 ∑ ∫ 𝑄 𝑖 𝑝 𝑖 (𝜀 𝑖 𝑝 𝑖 -𝛼 𝑖 𝑝 𝑖 Δ 𝑇 )𝑌 𝑖 𝑑𝑌 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 𝑛 𝑖 𝑝 𝑖 =1 𝑖 = 1,2 which finally leads to: 𝑁 𝑖 (𝑥) = ∑ 𝑄 𝑖 𝑝 𝑖 [∫ 𝑑𝑦 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 ] 𝑛 𝑖 𝑝 𝑖 =1 𝑑𝑢 𝑖 𝑑𝑥 -𝑏 ∑ 𝑄 𝑖 𝑝 𝑖 [∫ 𝑦 𝑖 𝑑𝑦 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 ] 𝑛 𝑖 𝑝 𝑖 =1 𝑑𝜃 𝑖 𝑑𝑥 -𝑏 ∑ (B-6) 𝑀 𝑖 (𝑥) = 𝑏 ∑ 𝑄 𝑖 𝑝 𝑖 [∫ 𝑦 𝑖 𝑑𝑦 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 ] 𝑛 𝑖 𝑝 𝑖 =1 𝑑𝑢 𝑖 𝑑𝑥 + 𝑏 ∑ 𝑄 𝑖 𝑝 𝑖 [∫ 𝑦 𝑖 2 𝑑𝑦 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 ] 𝑛 𝑖 𝑝 𝑖 =1 𝑑𝜃 𝑖 𝑑𝑥 + 𝑛 𝑖 𝑝 𝑖 =1 ∑ 𝑄 𝑖 𝑝 𝑖 𝛼 𝑖 𝑝 𝑖 [∫ 𝑦 𝑖 𝑑𝑦 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 ] Δ 𝑇 𝑛 𝑖 𝑝 𝑖 =1 (B-7) The parameters involving in the constitutive equations Eq. (31) to Eq. (33) are thus defined (B-1) (B-2) (B-4) (B-5) 𝑄 𝑖 𝑝 𝑖 𝛼 𝑖 𝑝 𝑖 [∫ 𝑑𝑦 𝑖 ℎ 𝑝 𝑖 ℎ 𝑝 𝑖 -1 ] Δ 𝑇 such as for i=1,2 𝐴 𝑖 = 𝑏 ∑ 𝑄 𝑖 𝑝 𝑖 (ℎ 𝑝 𝑖 -ℎ 𝑝 𝑖 -1 ) 𝑛 𝑖 𝑝 𝑖 =1 (B-8) 𝐵 𝑖 = 𝑏 2 ∑ 𝑄 𝑖 𝑝 𝑖 (ℎ 𝑝 𝑖 2 -ℎ 𝑝 𝑖 -1 2 ) 𝑛 𝑖 𝑝 𝑖 =1 (B-9) 𝐷 𝑖 = 𝑏 3 ∑ 𝑄 𝑖 𝑝 𝑖 (ℎ 𝑝 𝑖 3 -ℎ 𝑝 𝑖 -1 3 ) 𝑛 𝑖 𝑝 𝑖 =1 (B-10) 𝑁 𝑖 𝑇 ̅ = 𝑏 ∑ 𝑄 𝑖 𝑝 𝑖 𝛼 𝑖 𝑝 𝑖 (ℎ 𝑝 𝑖 -ℎ 𝑝 𝑖 -1 ) 𝑛 𝑖 𝑝 𝑖 =1 Δ 𝑇 (B-11) 𝑀 𝑖 𝑇 ̅ = 𝑏 2 ∑ 𝑄 𝑖 𝑝 𝑖 𝛼 𝑖 𝑝 𝑖 (ℎ 𝑝 𝑖 2 -ℎ 𝑝 𝑖 -1 2 ) 𝑛 𝑖 𝑝 𝑖 =1 Δ 𝑇 (B-12) where, in 𝑁 𝑖 (𝑋) = 𝑏 ∑ Appendix C

Acknowledgement

This work has not received any specific grant.

Appendix A

This appendix presents the mathematical description of the elementary stiffness matrix of the BBa element. Equation Eq. ( 7) can be explicitly written such as a system of a coupled second order ODE at constant coefficients:

This system is solved such as:

where c 1 , c 2 , c 3 and c 4 are integration constants. The boundary conditions at both extremities of the BBa element, in terms of displacements, lead to the expressions for the integration constants as functions of nodal displacements u i , u j , u k and u l (see Figure 4):

It can then be written under this shape:

With:

The normal forces are then computed from Eq. ( 2), Eq. (A-2) and Eq. (A-3):

The nodal normal forces are then deduced as function of the integration constants:

with:

In equation Eq. (A-15) the equivalent nodal force vector to the thermal load is appearing:

This appendix provides a brief description of the formulation of the beam element used in the 1D-beam analysis. Under the approximation of small bending angle, the local equilibrium equations of adherend i=1,2 outside the overlap are given by:

It corresponds to those obtained along the overlap when the adhesive stresses vanish. But, the normal force is equal to F at any X. The system of six first order linear ODEs to be solved is then found under the following shape:

The resolution is performed using the exponential matrix as described in section 2.2.3.

Appendix D

In [START_REF] Paroissien | A simplified stress analysis of bonded joints using macro-elements[END_REF][START_REF] Paroissien | Elasto-Plastic Analysis of Bonded Joints with Macro-Elements[END_REF], a path to take into account for a linear variation of the shear stress in the adherend thickness following [START_REF] Tsaï | Improved Theoretical Solutions for Adhesive Lap Joints Int[END_REF] in the formulation of MEs is described and reminded here. In the 1Dbar analysis, it is sufficient to modify the adhesive shear relative stiffness such as:

with:

where G 1 (G 2 ) is the shear modulus of the adherend 1 (2).

In the 1D-beam analysis, the constitutive equations of adherends to consider are: where H i is the shear stiffness taken into account the shear correction factor.


F Dufour

G Pijaudier-Cabot

NUMERICAL MODELING OF MATERIAL FLOW WITH AN INTERNAL LENGTH

Keywords: Cosserat model, Bingham rheology, fresh concrete flow

In this paper we propose a new numerical method to study material flow with an internal length. The finite element method with Lagrangian interpolation points (FEMLIP) is used for its capability to model large deformation processes for any type of materials. It is widely accepted in the litterature that fresh concrete rheology is very well estimated by Bigham model in the classical continuum approach. Nevertheless in most cases the form size and the distance between reinforced bars are about few times greater than aggregate size itself. This means that the microstructure has a non negligible effect on the flow (aggregate rotation independent of material rotation + boundary layer effect) and that it must be taken into account to get a good estimation of the filling capability of the material. We propose to use a Cosserat continuum with a Bingham law to tackle this work.

INTRODUCTION

The objectives of this work are to model heterogeneous two-phase fluid like fresh concrete flow during filling form processes with a particular concern about aggregate segregation and/or locking effect around reinforced bars. The aggregate mean size is only few times smaller than the distance between the bars or between the bars and the mould. Therefore, microstructure must be taken into account to get a good estimation of the concrete flow and boundary conditions must be addressed carefully.

We already have modelled [START_REF] Dufour | Développements de la méthode des éléments finis avec des points d'intégration lagrangiens : Application à la géomécanique[END_REF] concrete flow either with an homogeneous or an heterogeneous approach within the classical continuum framework. In the first approach, informations about material microstructure are clearly missing (e.g. material locking between bars cannot be modelled) whereas the second approach is very much demanding on computer ressources and cannot be used for large scale applications.

Therefore, we propose in this work to combine the Bingham model with a Cosserat continuum in order to enhance microstructure effects into the model without having to use very fine and numerically costly meshes. We wish to avoid as much as possible heterogeneous modelling for large scale problems because of memory and computational time limits. Along these lines, heterogeneous modelling will be used only for validations purposes.

NUMERICAL METHOD

To tackle these objectives, we use a finite element method with Lagrangian interpolation points (FEMLIP) implemented in our ellipsis code [START_REF] Moresi | Mantle convection models with viscoelastic/brittle lithosphere: Numerical methodology and plate tectonic modeling[END_REF] where the problem domain is represented by an Eulerian mesh with an embedded set of Lagrangian integration points or particles. Unknowns (i.e. velocity field) are computed at the mesh nodes and the Lagrangian particles carry history variables during the large deformation process including material properties. This method is ideally suited to model fluid-like behavior of continuum solids like two-phase fluid. From the material properties point of view, unlike in a pure Eulerian method, material interfaces do not show convection during computations as fictitious lines representing interfaces are preserved during flow between two adjacent particles with different properties. From the computational point of view, interfaces are taken into account in one finite element for which rheological properties are averaged for integration. The consequence is that the modelling of heterogeneous material needs very fine meshes with the FEMLIP same as for any numerical method based on integration over volumes.

RHEOLOGICAL LAWS

Concrete, in its fresh state, can be considered as a homogeneous fluid, provided that a certain degree of flow can be achieved and that the mix is quasi-homogeneous with regularly dispersed aggregate particles (in space and size). Modelling concrete rheology is a very difficult task as every concrete behaves differently depending on their composition and mixing method. Therefore it is impossible to write a relation between the stress and the strain rate based on every single parameter at the microscale, quantities and nature of each admixtures in particular, which are nowadays more and more numerous. Here, the concrete rheology is thought at a mesoscale and it is written as a relation between the shear stress τ and the shear strain rate γ. A complete review of rheological models for cement based materials has been done by [START_REF] Yahia | Applicability of rheological models to high-performance grouts containing supplementary cementitious materials and viscosity enhancing admixture[END_REF]. We discuss here only the three simplest ones, namely Newton's, Bingham's and Herschel-Bulkley's models.

The simplest model (one parameter) assumes that concrete is a Newtonian fluid:

τ = η γ (1)
where η is the viscosity but [START_REF] Tatersall | The rationale of a two-point workability test[END_REF] has shown that a single parameter is not enough to model concrete flow. At least two parameters are needed to model concrete with a slump lower than 30 cm. The first one needs to be related to the static behaviour (slump test) and the second one to the dynamic fluid response (rheometer test). Hence, a stress threshold has been added to the Newton's model to give the Bingham's model:

τ = τ 0 + η γ if τ ≥ τ 0 (2) γ = 0 if τ ≤ τ 0
where τ 0 is the yield stress and η is the plastic viscosity.

When experimental data are fitted for many different concretes and for self-compacting concrete especially, it often yields a negative value of the yield stress although very small. For this reason de Larrard et al. (1998) prefered to use a Herschel-Bulkley's model (3 parameters):

τ = τ 0 + a γb (3)
Ferraris and de Larrard (1998) have found a mean value for b of 1.53 for concrete without superplasticizer and 1.36 for others. However the Bingham's model (e.g. b = 1) is most widely used [START_REF] Tsong | Flow behaviour of high strength high-performance concrete[END_REF] because it needs less experimental tests to be calibrated. For the same reason, we have chosen to use the Bingham's model in this study.

One example for each rheological model is plotted in Figure 1. For an easier numerical implementation, the Bingham's model needs to be smoothed [START_REF] Papanastasiou | Flows of materials with yield[END_REF] in order to avoid the corner at τ = τ 0 on the material response and the use of two equations under and above the shear stress threshold. This yields:

η = η 0 + τ 0 γ 1 -e -m γ (4)
where m is the stress growth rate parameter. The higher m, the more accurate the approximation. We found out, for the range of strain rates reached during our numerical tests, that a value m=5000 is high enough as we observed no differences with results obtained with a value of m=10000. [START_REF] Dufour | Numerical modelling of concrete flow: homogeneous approach[END_REF]) have already shown the ability of this rheological law together with the FEMLIP to model concrete flow during slump test (Fig. 2) and Lbox filling (3).

The use of Bingham's law in for a classical continuum needs two different types of experimental tests. The first one, the slump test, to identify the yield stress and the second one, a rheometer test, to identify the plastic viscosity. For the rheometer test we have chosen to measure the flow time of fresh concrete in a funnel. The reason is that it can be easily modelled by ellipsis as the flow is nearly 2D.

ENHANCED MODEL

In a 2D conventional continuum the deviatoric part of an isotropic material is characterised by a shear viscosity η. For a Cosserat continuum we also have a Cosserat shear viscosity η c and a bending viscosity M . The constitutive relation for a generalised Newtonian fluid can be written in the usual pseudo-vector form:

τ = ΛD (5)
where τ is the deviatoric stress vector and the matrix Λ is expressed as:

        η η 0 0 0 0 η 0 0 0 0 η + η c η -η c 0 0 η + η c 0 0 symm. M 0 M         (6)
It means that when the Bingham's law for a Cosserat continuum is implemented two additional parameters are required. [START_REF] Choi | Distinct elements vs Cosserat theory: A comparison for the case of an infinite shearlayer[END_REF] have proposed relations between bending viscosity and viscosity η and internal length d such as:

M = 2ηd (7)
If we assume that the internal length of such a material is the granular size, then only one more parameter needs to be identified. One example of behaviour of this Cosserat continuum is studied for Newtonian material. We consider the discharge of a Cosserat material from a model silo (Fig. 4). In this case we are interested in the influence of the internal length parameter on the discharge velocity (Fig. 5) of a Cosserat fluid as described by ( 5 In Fig. 5 we plot the volume flowing out the reservoir versus time and for different values of the ratio α = d/a where a is the silo aperture and d is the internal length of the material. The flow rate is almost identical for all values of α (Fig. 5-d and5-e) larger than the aperture (1.0) which corresponds, in the elastic case, to a situation in which no flow can occur. Note that in purely viscous materials static equilibrium states do not exist, ensuring that stable arches do not form. For 0 ≤ α ≤ 1.0, the smaller the internal length the faster the outflow. Indeed [START_REF] Dufour | A particle-in-cell formulation for large deformation in cosserat continua[END_REF] have observed the same effects in the case of an infinite shearlayer. The first effect is a boundary layer effect along the boundaries. Here is a possibility to force the material rotation to zero for non spherical aggregates which can slide without rotating. The second effect is that for the same body force, a Cosserat fluid will flow slower than a classical one because energy is stored within the material due to bending stiffness. Frictional contacts between aggregates are modeled by Cosserat couple stresses.

BINGHAM'S FLOW WITHIN A COSSERAT CONTINUUM

To get the microstructure description of a Cosserat model together with the Bingham rheology we propose to model Bingham's flow within a Cosserat continuum. In the Cosserat constitutive matrix the viscosity η can be identified to the Bingham's viscosity, η c and M needs to be identified. The first term is a constant fraction of the Bingham's viscosity and the second one is calculated as a function of the internal length taken equal to the granular size D 50 . The coupling of the Bingham's law with a Cosserat model allows the use of coarser mesh with a homogeneous fluid approach compared to the modelling of heterogeneous fluid. The microstructure effect is brought into equations instead of the geometry. This coupling will be used in the future whenever boundary layer effects are of importance in the concrete flow for instance in the Lbox problem or in order to study concrete flow in the casting of densely reinforced strutural elements.

  FIG. 1. Different rheological models for fluid-like materials. (a) Newton (η = 650 Pa.s), (b) Bingham (τ 0 = 500 Pa, η = 500 Pa.s) and (c) Herschel-Bulkley (τ 0 = 400 Pa, a = 500 Pa.s b and b = 1.2).

  FIG. 2. Different snapshots of slump test on SCC at time 0 s. (a), 0.36 s. (b), 1.87 s. (c) and 7.24 s. (d).

  FIG. 3. Different snapshots of the Lbox test on SCC at time 0 s. (a), 0.45 s. (b), 0.85 s. (c) and 7.63 s. (d).

  FIG. 4. Snapshots of the model. (a) Initial conditions, (b) t = 4.75 10 -4 sec., (c) t = 1.42 10 -3 sec., (d) t = 2.85 10 -3 sec. and (e) t = 3.8 10 -3 sec.. The plot below each picture corresponds to the vertical velocity profile throughout the silo aperture.

FIG. 5 .

 5 FIG. 5. Outgoing volume for (a) α=0, (b) α=1/3, (c) α=2/3, (d) α=5/3 and (e) α=10/3