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

F Lachaud

J Morlier

S Schwartz

A Direct Method for the Assessment of Cohesive Zone Models for Thin Adhesive Layers Loaded in Mode I, Mode II, and Mixed-Mode I/II

Keywords: Adhesively bonded joint, cohesive zone model, macro-element, mode I, mode II, mixed-mode I/II, inverse method, direct method

In the context of increasing the strength-to-mass ratio of lightweight structures, the adhesively bonded joining technology appears to be an attractive solution. Nevertheless, the adhesive bonding method is important when the structural integrity of joints has to be ensured. In the literature, the cohesive zone models (CZMs) are shown to be able to predict both the static and fatigue strengths of adhesively bonded joints. The strength prediction is dependent on material laws and associated material parameters, characterizing the bondline behaviour mainly under pure mode I, mode II and mixed-mode I/II. The characterization methods are thus crucial. This paper aims at assessing the capabilities to identify the parameters of a particular CZM for both the inverse method, based on the energy balance associated with the path independent J-integral, and of a direct method described in this present work. The particular CZM has a classical shape based on the definition of a bilinear law for each of both pure modes, associated with pure mode interaction energy laws for initiation and propagation under mixed-mode I/II. The methodology used in this paper is based on a numerical test campaign only, involving the macro-element (ME) technique. A new approach for the fast formulation and implementation of ME modelling of two bonded beams is described.

Introduction

In the frame of structural design, the proper choice of joining technology is decisive for guaranteeing the integrity of the manufactured structure. The mechanical Rev. Adhesion Adhesives CC BY-NC -Creative Commons DOI: 10.7569/RAA.2018.097301 fastening, such as riveting or screwing, appears to be the reliable solution for the designers. Nevertheless, alone or in combination with the 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]. This higher level of mechanical performance allows for the design of lighter joints. In other words, the adhesive bonding offers the possibility to reduce the structural mass while ensuring the mechanical strength. The optimization of the strengthto-mass ratio is a challenge for several industrial sectors, such as aerospace, automotive, rail or naval transport industries. But, the reduction of structural mass makes sense only if the structural integrity is ensured. So to obtain the benefit from the adhesive bonding in view of mass reduction, it is required to be able to predict the strength of bonded joints.

The strength prediction consists in the comparison of computed strength criteria with allowable design values. The strength criteria could be based on theoretical, empirical, semi-empirical investigations and possibly including in-service feedback. The stress analysis allows for the computation of input data, necessary for the assessment of strength criteria. The experimental characterization allows then for the definition of allowable design values as well as of mechanical behavior to be used as input data for the mechanical analysis. As highlighted in [START_REF] Jumel | Instrumented end notched flexure -Crack propagation and process zone monitoring. Part I: Modelling and analysis[END_REF], the strength of a joining system at the macroscale depends on the experimental test specimen and procedure used. According to recent literature [START_REF] Khoramishad | A generalized damage model for constant amplitude fatigue loading of adhesively bonded joints[END_REF][START_REF] Khoramishad | Fatigue damage modeling of adhesively bonded joints under variable amplitude loading using cohesive zone model[END_REF][START_REF] Li | Mixed-mode cohesive zone models for fracture of an adhesively bonded polymer matrix composite[END_REF][START_REF] Da Silva | Advances in numerical modeling of adhesive joints[END_REF][START_REF] Campilho | Modelling adhesive joints with cohesive zone models: Effects of the cohesive law shape of the adhesive layer[END_REF] the cohesive zone model (CZM) appears to be one of the most suitable approaches that is able to model both the static and the fatigue behaviors of adhesive joints. A CZM offers thus the possibility to experimentally investigate both local and global mechanical behaviors. Based on Damage Mechanics and Fracture Mechanics, the CZM enables a diagnostic of the current state of the adhesive damage along the overlap. The damage associated with micro-cracks and/or voids coalescence results in a progressive degradation of the material stiffness before failure. Contrary to approaches based on Fracture Mechanics, the CZM does not require the hypothesis of the existence of an initial flaw. The localization of critical sites is an output of CZM-based approaches.

The characterization of CZM parameters makes use of a pre-cracked bonded overlap specimen loaded under various loadings, which are classically used in the frame of Fracture Mechanics.

Indeed, the end notched flexure (ENF) and double cantilever beam (DCB) test configurations are most frequently used for characterizing the CZM parameters of the adhesive employed in a thin adhesive layer loaded in pure mode I and pure mode II [START_REF] Da Silva | Advances in numerical modeling of adhesive joints[END_REF]. The characterization under mixed-mode I/II loading through mixed mode cantilever beam (MCB) and mixed mode bending (MMB) test configurations offers the possibility to subject the adhesive layer to a wide range of mixed-mode ratios while almost keeping the same experimental settings [START_REF] Hadavinia | The prediction of crack growth in bonded joints under cyclic-fatigue loading II. Analytical and finite element studies[END_REF][START_REF] Reeder | Mixed-mode bending method for delamination testing[END_REF][START_REF] Kenane | Mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites under fatigue loading[END_REF][START_REF] Högberg | Mixed-mode testing of adhesive layer[END_REF][START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF].

To characterize the CZM parameters, the inverse method based on the computation of the path independent J-integral [START_REF] Alfredsson | An experimental method to determine the complete stress-deformation relation for a structural adhesive layer loaded in shear[END_REF] along a closed contour of specifically CC BY-NC -Creative Commons 3 DOI: 10.7569/RAA.2018.097301 designed specimens has been suggested [START_REF] Högberg | Mixed-mode testing of adhesive layer[END_REF][START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF][START_REF] Alfredsson | An experimental method to determine the complete stress-deformation relation for a structural adhesive layer loaded in shear[END_REF][START_REF] Alfredsson | On the instantaneous energy release rate of the end notched flexure adhesive joint specimen[END_REF][START_REF] Anderson | The stress-elongation relation for an adhesive loaded in peel using equilibrium of energetic forces[END_REF][START_REF] Leffler | Shear behavior of adhesive layers[END_REF][START_REF] Cui | Constitutive law of adhesive layer measured with mixed-mode bending test[END_REF][START_REF] Cui | Delamination and debonding failure of laminated composite T-joints[END_REF]. The inverse method is based on experimental measurements of relative displacements at the crack tip of the bonded interface as a function of the reaction force. The use of the classical model of a beam on an elastic foundation allows for a simple computation of the J-integral and thus for an assessment of the CZM. Nevertheless, it suffers from some limitations mainly due to the hypotheses on which the J-integral is defined. A direct method has been previously suggested in [START_REF] Lélias | Mechanical behavior of adhesively bonded joints: Modeling, simulation and experimental characterization[END_REF]. Similarly to the inverse method, it makes use of both the experimental measurements of the displacement field at the bonded interface crack tip and of the model of beam on an elastic foundation. It allows for the characterization of CZM parameters assuming sufficiently accurate experimental measurements.

This paper aims at assessing the capabilities to identify the parameters of a CZM by the inverse method and the direct method. The CZM selected in this paper has a classical shape, based on the definition of a bilinear law for each of the two pure modes, associated with pure mode interaction energy laws for initiation and propagation under mixed-mode I/II. The methodology is based on a numerical test campaign performed only on the MCB test configuration. In other words, numerical tests are employed instead of experimental tests to evaluate both characterization methods. The numerical analysis uses the macro-element (ME) technique [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] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF], through a home-made computer code implemented on SCILAB. This technique has already been successfully assessed against the Finite Element (FE) analysis in a similar context [START_REF] Lélias | Mechanical behavior of adhesively bonded joints: Modeling, simulation and experimental characterization[END_REF][START_REF] Paroissien | A simplified stress analysis of bonded joints using 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]. Moreover, a new approach for an easy and fast implementation of the ME technique for the modelling of a bonded overlap is provided.

Numerical Test Campaign

Test Configuration

In the frame of the numerical test campaign presented in this paper, the MCB test configuration has been selected. It has been suggested by Högberg and Stigh [12,[START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF]. Similarly to the DCB test configuration, the loading consists in a pair of forces (termed F), being of the same magnitude but in opposite directions. Nevertheless, the action direction of the pair of forces is defined by an angle a, which allows for the adhesive layer to be subjected to pure mode I, pure mode II and mixed-mode I/II (see Figure 1). In this paper, contrary to the MCB test configuration described in [START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF], the two forces are not applied at the middle plane of the specimen but at neutral lines of each adherend. This MCB test configuration with a = p/2 corresponds to a DCB test configuration, while it leads with a = 0 to zero peel stress at the crack tip where the CZM parameters are measured.

The selected specimen design, including geometrical and material parameters, corresponds to the one described by Högberg and Stigh [START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF]. The crack length a=0 Rev. Adhesion Adhesives CC BY-NC -Creative Commons DOI: 10.7569/RAA.2018.097301 is then chosen (see Figure 1). The geometrical parameters are provided in Table 1 in conjunction with Figure 1.

The adherends are made of steel with a Young's modulus E = 200 GPa and a Poisson's ratio n = 0.3. The design is such that the adherends will remain in their linear elastic domain. The adhesive is assumed to have a classical bilinear damage evolution law following [START_REF] Allix | Damage mechanics of interfacial media: Basic aspects, identification and application to delamination, in Damage and Interfacial Debonding in Composites[END_REF], involving interaction energy laws for both initiation and propagation under mixed-mode:

G G G G G G G G I Ie n II IIe n I Ic n II IIc n       +       =       +       = 1 1         (1)
where n is a material parameter to be identified, G Ic and G IIc are the fracture energies in mode I and mode II, G Ie and G IIe are the elastic strain energies stored in mode I and mode II and G I and G II are related to the strain energy release rates in mode I and mode II, respectively. The fracture energies in mode I and mode II and the elastic stiffnesses under peel and shear, termed k I and k II respectively are the same as those used by Högberg and Stigh [START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF]. Nevertheless, the adhesive maximal peel and shear stresses, termed S max and T max , are different, to ensure a right energy dissipation during loading [START_REF] Turon | Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: Definition of interlaminar strengths and elastic stiffness[END_REF]. It is indicated that the law by Allix and Ladevèze [START_REF] Allix | Damage mechanics of interfacial media: Basic aspects, identification and application to delamination, in Damage and Interfacial Debonding in Composites[END_REF] already includes this condition. It is then chosen to keep the same maximal shear stress T max = 26 MPa, resulting in a Table 1 Geometrical parameters of the MCB specimen [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF].

a

Numerical Analysis

Macro-Element Technique Overview

The numerical analysis is performed using the macro-element (ME) technique for the modelling of bonded overlap [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] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF]. This technique has been developed for the simplified stress analysis of hybrid (bolted/bonded) joints [START_REF] Paroissien | Analytical two-dimensional model of a hybrid (bolted/bonded) single-lap joint[END_REF][START_REF] Paroissien | Simplified stress analysis of hybrid (bolted/bonded) joints[END_REF]. The ME technique is inspired by the Finite Element (FE) method and differs in the sense that the interpolation functions are not assumed, since they take the shape of the solutions of the governing differential equations system. A direct consequence is that only one ME is sufficient to mesh a complete bonded overlap in the frame of a linear stress analysis. The bonded overlap is then modelled by a four-node ME -also called bonded-beams element -the nodes of which are located at the extremities of the overlap on the neutral axes of adherends (see Figure 2). This ME involves 3 degrees of freedom per node or a total of twelve for a 1D-beam analysis.

The main work is thus the formulation of the elementary stiffness matrix of the bonded-beams element. Indeed 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 distributions along the overlap of adhesives stresses, internal forces and displacements in the adherends.

Hypotheses and Formulation of the Elementary Stiffness Matrix

An approach for the formulation of the stiffness has already been described in detail in previous papers [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] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively 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 a mathematical software such as SCILAB for example. Compared with the early approach, the shape of solutions in terms of displacements and internal loads is not provided. Nevertheless, in the frame of nonlinear material analyses such as the one presented in this paper, the bonded overlap has to be meshed in order to locally update the material parameters within an iterative computation procedure (see section 2.2.3). As a result, the displacements and internal loads will be read directly at nodes. Moreover, the following description is useful for the derivation of the direct method. It is assumed that the thickness of the adhesive is constant along the length ∆ of the macro-element. Moreover, the adherends are simulated as linear elastic Euler-Bernoulli laminated beams.The general shape of the constitutive equations for the adherend j=1,2 provides the six first differential equations: where N j (M j ) is the normal force (bending moment) of the adherend j and u j (v j , q j ) is the longitudinal displacement (deflection, bending angle) of the adherend j. For the adherend j, A j represents the extensional stiffness, D j the bending stiffness and B j the coupling stiffness and D j =A j D j -B j B j ≠0. In this paper, the coupling stiffnesses remain equal to zero and A 1 =A 2 =Etb and D 1 =D 2 =Et 3 b/12. It is indicated that the Euler-Bernoulli kinematics is employed by Högberg and Stigh in [START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF].

The adhesive layer is simulated by an infinite number of elastic shear and transverse springs attached at both adherend interfaces. The adhesive shear stressdenoted T -and the adhesive peel stress -denoted S -are then given by:

S k v v T k u u h h

I II = - ( ) = -- - ( )      1 2 2 1 2 2 1 1 q q (3)
where h j is the half thickness of adherend j. In this paper, h 1 =h 2 =t/2. The classical local equilibrium from Goland and Reissner [START_REF] Goland | The stresses in cemented joints[END_REF] is used and provides the six last differential equations for j=1,2 (see Figure 3):

dN dx bT dV dx bS dM dx V bh T j j j j j j j = - ( ) = - ( ) + + =          + 1 1 0 1 ( 4
)
where V j is the shear force of the adherend j. Considering the local equilibrium equations Eqs. ( 4), the adhesive stresses are replaced by their expressions as functions of adherend displacements Eqs. (3). In conjunction with Eqs (2), it results in a system of twelve linear first-order ordinary differential equations:

Tbdx N 1 +dN 1 M 1 +dM 1 V 1 +dV 1 N 2 +dN 2 M 2 +dM 2 V 2 +dV 2 V 1 V 2 M 2 N 2 N 1 x M 1 y + Tbdx Sbdx
du dx D N B M du dx D N B M dv dx dv dx d d 1 1 1 1 1 1 1 2 2 2 1 2 2 2 1 1 2 2 1 = + = + = = ∆ ∆ ∆ ∆ q q q x x B N A M d dx B N A M dN dx bk u bk h b II II = + = + = + - 1 1 1 1 1 1 2 2 2 2 2 2 2 1 1 11 ∆ ∆ ∆ ∆ q q k k u bk h dN dx bk u bk h bk u bk h dV dx II II II II II II 2 2 2 2 1 11 2 2 2 1 + = - - + - q q q = = - = - + = + - bk v bk v dV dx bk v bk v dM dx bh k u bk h h b I I I I II II 1 2 2 1 2 1 1 1 1 1 1 q h h k u bk h h V dM dx bh k u bk h h bh k u b II II II II II 1 2 1 2 2 1 2 2 1 2 2 1 2 2 + - = + - + q q k k h h V II 2 2 2 2 q -                                    (5)
This system can be written as

dX dx AX =
where A is 12×12 matrix with real constant components and the unknown vector X is such that

t X=(u 1 u 2 v 1 v 2 q 1 q 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 [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] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF], such as:

CC BY-NC -Creative Commons

( ) -( ) ( ) ( ) -( ) -( ) ( ) ( ) -( ) - ( ) N N N N V V V V M M M 1 2 1 2 1 2 1 2 1 2 1 0 0 0 0 0 0 ∆ ∆ ∆ ∆ ∆ ∆ ∆ ( ) ( ) ( )                                     = M K u u BBe 2 1 2 0 0 ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )      u u v v v v 1 2 1 2 1 2 1 2 1 2 0 0 0 0 ∆ ∆ ∆ ∆ ∆ ∆ q q q q                                 (6)
The fundamental matrix of A, termed F A , is computed at x=0 and x=∆; using the SCILAB software, the associated command is "expm":

F F A A x expm A x expm A = ( )= ( ) = ( )= ( )    0 0 . . ∆ ∆ (7)
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 F A (x=0) and the second block of six lines and twelve rows comes from F A (x=∆), such that: ,,,,,; :

M x x U A i j A i ' , , , , , , ; : , = ( )=
j A i j x = = = = ( )            F ∆            (8)
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 (∆) q 1 (0) q 2 (0) q 1 (∆) q 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

(0) -N 2 (0) N 1 (∆) N 2 (∆) -V 1 (0) -V 2 (0) V 1 (∆) V 2 (∆) -M 1 (0) -Μ 2 (0) Μ 1 (∆) M 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 [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] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF]:

K BBe =N.M -1 .
Even if it is not the topic of this paper, it is obvious that this previous approach can be easily used to develop ME, under different local equilibrium equations (e.g. Hart-Smith [START_REF] Hart-Smith | Adhesive bonded single lap joints[END_REF], Luo and Tong [START_REF] Luo | Analytical solutions for nonlinear analysis of composite single-lap joints[END_REF]) or under different constitutive equations (e.g. Tsaï et al. [START_REF] Tsai | Improved theoretical solutions for adhesive lap joints[END_REF]) and/or including different number of layers of adhesives and adherends (e.g. double lap joint configuration).

Mesh and Boundary Conditions

The bonded overlap is regularly meshed with a parametrical number n_ME of bonded-beams elements. One extremity is clamped and the load is applied under displacement (termed H) at the other extremity [START_REF] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF] (see Figure 4). This choice of boundary conditions does not correspond to those experienced during experimental tests [START_REF] Högberg | Mixed-mode testing of adhesive layer[END_REF]. However, it is necessary to prescribe relevant degrees of freedom to run any analysis. In view of the application of the inverse method, it is mandatory that the adhesive layer does not deform at the joint extremity where the load is not applied. Clamping conditions avoid both peel and shear deformations. The inverse and direct methods are then applied by taking into account these boundary conditions (see section 3.1).

The results are not presented in this paper but a study on the influence of mesh size up to a maximal mesh density of twenty MEs per mm was performed under a pure linear elastic analysis under pure mode I (a = p/2). The conclusions are that (i) the original approach and the present approach for the formulation of the elementary stiffness matrix of ME provides exactly the same results, and (ii) the computed reaction as well as the adhesive peak stresses do not vary at all with the mesh density.

Nonlinear Computation Management

The use of a nonlinear adhesive material implies that the computation is nonlinear. A detailed description of the nonlinear algorithm used is provided in [START_REF] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF]. Only a brief overview is given here. The algorithm is based on the Newton-Raphson method and uses the secant stiffness matrix with an update at each iteration. In particular, the damage parameter is computed at each nodal abscissa according to the introduced adhesive material law. The norm of displacement jump (in mm) of interface l is defined by:

l d d = ( ) + ( ) v u 2 2 (9)
where d v (d u ) is the displacement jump of the interface (see Table 2) along the y-axis (x-axis). A mixity parameter b is defined by:

b d d q q = = -- - - u v u u h h v v 2 1 2 2 1 1 1 2 (10)
The mixity parameter is representative of the level of sliding displacement at the interface compared with the opening displacement at the interface. At each iteration, the mixity parameter b is updated. Under the current local mixity parameter, it is assumed that the material law is bilinear, such that the damage parameter d is:

d f e f = - ( ) - ( ) l l l l l l ( 11
)
where l e (l f ) is the displacement jump (in mm) of the interface at initiation (propagation). In order to compute l e (l f ), the interaction laws Eq. ( 1) are used while classically assuming that the projections on pure modes of the mixed mode evolution law under the current local mixity are bilinear (see Table 2): (

1 1 1 2 2 2 1 2 2
) ) + ( ) ( ) + ( )                      2
The damage parameter is computed only if d v is positive. Each ME is then updated with the damaged elastic stiffness taken as the maximal value of both damage parameters computed at each extremity of the ME.

Finally, the displacement is linearly applied as a function of the numerical time. All the numerical test results presented in this paper are obtained from a simulation run involving one hundred time steps, with a constant numerical time step d t .

CC BY-NC -Creative Commons DOI: 10.7569/RAA.2018.097301

Characterization Method

The Inverse Method

The inverse method is based on the energy balance associated with the computation of the path independent J-integral [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF] on a closed contour G:

J Wdy T dU dx ds = - ∫ ( 13
)
where W is the strain energy density, T n = s is the traction vector, σ is the stress tensor, U is the displacement vector, n is the normal unit vector directed outward to the path G, and (x,y) is the specified two-dimensional coordinate system. From the fundamental work by Fraisse and Schmit [START_REF] Fraisse | Use of J-integral as fracture parameter in simplified analysis of bonded joints[END_REF] it is shown that the J-integral parameter can be computed from stress analysis based on a model of beam on an elastic foundation as:

J T d S d u v u v u u v v u v d d d d d d d d d d , , ,
( )= ( ) + ( ) ∫ ∫ 0 0 (14)
In the frame of the inverse method:

(i) the adhesive peel stress is obtained from experimental tests under pure mode I loading as [START_REF] Anderson | The stress-elongation relation for an adhesive loaded in peel using equilibrium of energetic forces[END_REF]:

S J v u v v d d d d ( )= ∂ ( ) ∂ , (15)
(ii) the adhesive shear stress is obtained from experimental tests under pure mode II loading as [START_REF] Alfredsson | An experimental method to determine the complete stress-deformation relation for a structural adhesive layer loaded in shear[END_REF]:

T J u u v u d d d d ( )= ∂ ( ) ∂ , (16)
(iii) the adhesive peel and shear stresses are obtained from experimental tests under mixed-mode I/II loading as [START_REF] Högberg | Mixed-mode testing of adhesive layer[END_REF][START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF]:

S J v u v v d d d d ( )= ∂ ( ) ∂ , ( 17
)
T J u u v u d d d d ( )= ∂ ( ) ∂ , (18)
An advantage of this method is that it offers the possibility to monitor the evolution of the adhesive stress at the crack tip from the measurements of macroscopic quantities possibly measurable from experimental test fixtures, such as the CC BY-NC -Creative Commons 13 DOI: 10.7569/RAA.2018.097301 applied load (in N) or the evolution of displacement jump (in mm) at the crack tip. The related expressions can be found in [START_REF] Högberg | Mixed-mode testing of adhesive layer[END_REF][START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF][START_REF] Alfredsson | An experimental method to determine the complete stress-deformation relation for a structural adhesive layer loaded in shear[END_REF][START_REF] Alfredsson | On the instantaneous energy release rate of the end notched flexure adhesive joint specimen[END_REF][START_REF] Anderson | The stress-elongation relation for an adhesive loaded in peel using equilibrium of energetic forces[END_REF][START_REF] Leffler | Shear behavior of adhesive layers[END_REF][START_REF] Cui | Constitutive law of adhesive layer measured with mixed-mode bending test[END_REF][START_REF] Cui | Delamination and debonding failure of laminated composite T-joints[END_REF]. In order to address the constraints inherent to the numerical analysis (see section 2.2.3) the following approximate expression for the J-integral parameter, deduced from the work by Fraisse and Schmit [START_REF] Fraisse | Use of J-integral as fracture parameter in simplified analysis of bonded joints[END_REF] after the relevant adaptations to the Euler-Bernoulli framework (the same as ME) is used:

J b N L A M L D b N L A M L D u v d d , ( )= ( ) + ( )         + ( ) + ( )     1 2 1 2 1 2 1 2 2 2 2 2      - ( )+ ( ) ( ) + ( )+ ( )+ ( ) -( ) ( )      1 2 0 0 2 0 0 2 0 0 1 2 2 1 2 2 1 b N N A M M h N N               + ( ) ( )+ ( ) ( )     - ( ) ( ) 2 1 1 2 2 1 1 8 1 1 0 0 D b V L L V L L b V q q q + + ( ) ( )     V 2 2 0 0 q (19)
This expression is valid when A 1 =A 2 =A, D 1 =D 2 =D, B 1 =B 2 =0 and t 1 =t 2 =t=2h. Even if the previous expression can be computed using the output from the stress analysis, it is nevertheless suitable for the analysis of experimental tests, which needs the application data reduction scheme [START_REF] Sarrado | An experimental data reduction method for the mixed mode bending test based on the J-integral approach[END_REF][START_REF] Azari | Fracture load predictions and measurements for toughened epoxy adhesive joints[END_REF]. In addition, the required differentiation of the J-integral parameter to obtain the adhesive peel and shear stresses at x=L is taken as the ratio of the difference between two consecutive computed J over the difference between two consecutive computed displacement jump d u or d v :

S t J t J t t t T t J t J t v t v v t u t u d d d d d d d d ( ) ( ) = ( ) -- ( ) ( ) - - ( ) ( ) ( ) = ( ) -- ( ) t t t u t ( ) - - ( )        d d (20)
This is the reason why the chosen numerical time step is elevated (d t = 100).

The Direct Method

This method is presented in [START_REF] Lélias | Mechanical behavior of adhesively bonded joints: Modeling, simulation and experimental characterization[END_REF]. It is based on the measurement around the crack tip of the displacement of the neural axis according to the x-axis and the y-axis. Contrary to the inverse method, no spatial integration of equilibrium equations is required. It is assumed that the coupling stiffness remains equal to zero B 1 =B 2 =0;

CC BY-NC -Creative Commons DOI: 10.7569/RAA.2018.097301 the expressions for the adhesive shear stress including the coupling stiffnesses are provided in Appendix A.

In the case of pure mode I loading, the adhesive shear stress vanishes so that the local equilibrium of adherends can be reduced to the following set of differential equations for j=1,2:

+ =          + 0 1 0 1 ( 21
)
As a result:

S b d M dx j j = - ( ) 1 1 2 2 (22)
Using the constitutive relationship, the adhesive peel stress can be expressed as:

S D b d w dx j j j = - ( ) 1 4 4 (23)
In the case of pure mode II loading, the adhesive peel stress vanishes so that the local equilibrium of adherends can be reduced to the following set of differential equations for j = 1,2:

dN dx bT dV dx dM dx V bh T j j j j j j = - ( ) = + + =          1 0 0 (24)
As a result:

T b dN dx j j = - ( ) 1 1 (25)
Using the constitutive relationships, the adhesive shear stress can be expressed as: In the case of mixed-mode I/II loading, the local equilibrium of adherends is given by Eq. ( 4). The following expressions for the adhesive peel and shear stresses are obtained:

1 2 2 ( 28
)
In the frame of this paper, the stress analysis provides the nodal displacements only, due to the approach chosen for the formulation of the elementary stiffness matrix. The consecutive differentiation of nodal displacements at the crack tip requires assessing the adhesive peel and shear stresses which are computed as the ratio of the difference between the displacements at the node located at x=L and the node at x=L-L/n_ME over L-L/n_ME.

Test Campaign

The assessment of both inverse and direct methods is performed through five different loading conditions, corresponding to five different values of the angle α (see Figure 4): i. pure mode I with a = p/2; ii. mixed-mode I/II with a = p/4; iii. mixed-mode I/II with a = p/8; iv. mixed-mode I/II with a = p/16; v. pure mode II with a = 0.

For each of the five loading conditions, four virtual tests are performed and post-processed. Each of the virtual tests is associated with a mesh refinement or a density of mesh per mm: based on the constitutive behavior of the adhesive material. The results denoted as computed are related to direct output data of the numerical analysis as the adhesive stress distribution (using Eq. ( 3)), the reaction force or the J-integral parameter. Finally, results obtained from the inverse method and from the direct method are provided on the basis of computed results as described in section 3.1 and section 3.2 respectively. For each loading condition, the results are given up to the propagation beginning at x=L (it means d=1 at x=L).

Results and Discussion

Pure Mode I and Pure Mode II

As a function of the displacement jump at x=L and up to the propagation beginning at x=L (it means d=1), the computed reaction forces F under pure mode I and pure mode II are provided for the four mesh densities in Figure 5 and Figure 6. It appears that the mesh density has a significant influence on the accuracy of predictions under pure mode I: a maximal relative difference of -11% in the reaction force is obtained for a mesh density of eight MEs per mm. The reaction force tends to stabilize with increase in mesh density. On the contrary, this influence appears as very limited under pure mode II, for the range of mesh density selected; a maximal relative difference of -1.1% in the reaction force is obtained for a mesh density of eight MEs per mm. However, the difference in the computed reaction forces with one ME per mm and eight MEs per mm remains similar for the pure mode I and pure mode II: 82 N and 84 N respectively. But the reaction force in pure mode II is about ten times higher than in pure mode I. As a result, the relative difference in pure mode II appears as much reduced. The peel (shear) stress distributions as a function of the displacement jump at x=L predicted by the inverse method and by the direct method are provided in Figure 7 (Figure 8) for a mesh density of eight MEs per mm. Moreover, they are compared with the theoretical and the computed ones. Under pure mode I, the four peel stress distributions appear as superimposed. Similarly, under pure mode II, the four shear stress distributions appear as very close. Only the inverse method overestimates the maximal stress, before a slight underestimation at higher jump displacements. The relative differences in the maximal peel (shear) stresses obtained with the inverse and direct methods from the theoretical peel (shear) stress peaks are given in Figure 9 (Figure 10) as a function of the mesh density. Under pure mode I, the increase of mesh density tends to reduce the difference from the theoretical value for both methods and to reduce the rate of this evolution. This behaviour is observed for the direct method under pure mode II too. For the inverse method, if the increase in the mesh density reduces the evolution rate, the relative difference from the theoretical value tends to increase by up to 6.7%, which is not considered as significant.

It can be concluded firstly that the adhesive stress evolution computed as a function of displacement jump at x=L, referred to as the adhesive constitutive evolution, corresponds to the theoretical ones under pure mode I and pure mode II. Secondly, the direct method provides the theoretical adhesive constitutive evolution under pure mode I and pure mode II. Thirdly, if the inverse method provides the theoretical adhesive constitutive evolution under pure mode I, it provides under pure mode II a relevant adhesive constitutive evolution along the linear range and along a large part of the softening part of the constitutive behavior with possibly a slight overestimation of the maximal stresses. Nevertheless, from Figure 10, it is shown that the inverse method is able to predict the theoretical maximal stress, by selecting a coarser mesh density.

Mixed-Mode I/II

In this section, quadratic initiation and propagation criteria are used; this means n=2 in Eq. ( 1). The computed reaction forces F for the three mixed-mode loading conditions defined by a = p/4, a = p/8 and a = p/16 are provided as functions of the norm of the displacement jump of both interfaces at x = L for the four mesh densities in Figure 11. Similarly to the pure mode II case, the influence of the mesh density on the computed reaction force appears as very limited, for the range selected, due to the increase of the level of the reaction force (see section 4.1). The reaction force computed for a mesh density of one ME per mm is different from the one for a mesh density of eight MEs per mm, i.e., -2.4%, -1.8% and -1.2% for α=π/4, α=π/8 and a = p/16, respectively.

The adhesive peel and shear stress distributions as a function of the displacement jump at x = L predicted by the inverse method and by the direct method are provided in Figure 12 and Figure 13 respectively, for a mesh density of eight It is shown that the predictions by the direct method and the inverse method are different. Of course, the direct method predicts the computed results, since the direct method and the numerical simulations based on ME use the same hypotheses. The accuracy of direct method predictions is related to the post-processing approach and thus to the mesh refinement (see Figure 14). The adhesive constitutive evolutions in peel and in shear under mixed-mode identified by the inverse method are significantly different from the computed ones. The reason is that the way the shear and peel stresses are obtained in the inverse method under mixedmode are not mathematically true in general. Indeed, Eq. ( 14) does not lead to Eq. ( 17) and Eq. ( 18) in general. The particular CZM used in this paper to model the behaviour of the adhesive layer in the bonded joint is thus a case for which it is not true. When it is true, the adhesive peel and shear stresses are derivable from a potential energy [START_REF] Goustanios | Path dependence of truss-like mixed mode cohesive laws[END_REF]. For the particular CZM used in this paper, the adhesive peel and shear stresses appear as not derivable from a potential energy so that the work of separation is path-dependent.

Moreover, the constitutive adhesive evolutions in peel and shear computed and predicted by the direct method are given for a mesh density of eight MEs per mm and for the three loading conditions in Figure 15 and Figure 16 respectively. Under the mixed-mode loading conditions a = 11 p/4 and a = p/8, the adhesive constitutive evolution in peel remains very close to the pure mode I (see Figure 15 more than a small difference of -1.7%. from the maximal theoretical peel stress. On the contrary, the variation in mixed-mode loading conditions is clearly shown in the adhesive constitutive evolution in shear (see Figure 16). It can thus be deduced that loading conditions close to the pure mode II are relevant choices for the experimental characterization of the mechanical behavior of the adhesive layer under mixedmode loading conditions, if the MCB test configuration is chosen. The actual level of mixity at x=L has been investigated. A local mixity angle y relevant to the actual description of loading conditions (angle α) is defined as:

y p b = - - 2 1 tan (29)
As already shown in [START_REF] Högberg | Specimen proposals for mixed-mode testing of adhesive layer[END_REF], the local mixity angle at x=L differs from the loading condition a for the crack length selected (see Figure 17). Moreover, as a function of the increase of the damage parameter at x=L, the local mixity angle at x=L decreases to reach a value close to a. This evolution of the local mixity angle should lead to the shapes of the peel and shear stresses (see Figure 15 and Figure 16) which deviate from those of pure mode ones.

Assessment of Mixed-Mode I/II Criterion

There are several shapes of initiation and propagation [10-11, 21, 36]. An approach could consist in choosing a shape for the initiation and propagation and then to best fit the related parameters to the assessed peel and shear stresses. The energy criterion for initiation and propagation in this paper uses only one material parameter, the exponent n (see Eq. ( 1)). A comparison of the reaction force, adhesive peel and shear stress as a function of jump displacement is then provided between n=1 and n=2 in Figure 18 the one with n=2. As a result, a small difference is thus shown on the macroscopic measurement due to the choice of the mixed-mode model. Qualitatively, the shape of the adhesive constitutive evolution in peel and shear is similar for n=1 and n=2, although a slight variation in the slopes after the initiation is observed. Nevertheless, a clear difference appears quantitatively in terms of level of stress and jump displacement. A relative difference of -10% (-17%) for the case with n=2 is obtained for the maximal peel (shear) stress, respectively. As a result, it could be thought that models for the initiation and propagation could be assessed through the use of both the direct method and the numerical analysis. Finally, since the inverse method cannot always be used under mixed mode loading, it could be a priori. As a result, a preliminary and iterative step including experimental and numerical tests has to be conducted to design the test specimen to demonstrate the ability to obtain relevant experimental data. The use of the ME technique could help in reducing the computation time for numerical tests to virtually explore the design possibilities associated with the loading conditions geometries and adhesive constitutive behaviour laws. The implementation of an ME using the Timoshenko model instead of Euler-Bernoulli model could improve the relevance of the assessment of constitutive laws, especially when using thicker adherends.

Nomenclature and Units

Figure 1

 1 Figure 1 MCB test configuration considered in this paper and the definitions of geometrical parameters.

Figure 2

 2 Figure 2Modelling of a bonded overlap by a bonded-beams element[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] Lélias | An extended semi-analytical formulation for fast and reliable mode I/II stress analysis of adhesively bonded joints[END_REF].

Figure 3

 3 Figure 3 Free-body diagram of infinitesimal elements of the adherend 1 (top) and adherend 2 (bottom).

9 DOI: 10

 910 .7569/RAA.2018.097301

-

Figure 4

 4 Figure 4 Applied displacement H and fixed displacement applied to an MCB test configuration.

  CC BY-NC -Creative Commons15 DOI: 10.7569/RAA.2018.097301

  i. one ME per mm; ii. two MEs per mm; iii. four MEs per mm; iv. eight MEs per mm. The results of the test campaign are presented and discussed in section 4. The results denoted as theoretical are related to the input data provided to the model Rev. Adhesion Adhesives CC BY-NC -Creative Commons DOI: 10.7569/RAA.2018.097301

Figure 5

 5 Figure 5 Reaction force as a function of the displacement jump at x=L for the four mesh densities under mode I.

Figure 6

 6 Figure 6 Reaction force as a function of the displacement jump at x=L for the four mesh densities under mode II.

Figure 7

 7 Figure 7 Peel stress at x=L as a function of the displacement jump. Comparison of distributions among the theory, computation, the direct method and the inverse method for a mesh density of 8 MEs per mm.

Figure 8 Figure 9

 89 Figure 8 Shear stress at x=L as a function of the displacement jump. Comparison of distributions among the theory, computation, the direct method and the inverse method for a mesh density of 8 MEs per mm.

Figure 10

 10 Figure 10Relative difference in % in the maximal shear stress obtained with the inverse and direct methods from the theoretical shear stress peak as a function of mesh density.

Figure 11 Figure 12

 1112 Figure 11 Reaction force as a function of the norm of displacement jump at x=L for the four mesh densities under mixed-mode loading conditions a = p/4, a = p/8 and a = p/16.

Figure 13

 13 Figure 13 Shear stress at x = L as a function of the displacement jump under a loading condition a = p/16. Comparison of distributions among the theoretical pure mode II, computation, the direct method and the inverse method for a mesh density of 8 MEs per mm.

[ 8 Figure 14 Figure 15

 81415 Figure 14 Peel stress computed from the direct method at x = L as a function of the displacement jump under a loading condition a = p/16 for 1 ME per mm and 8 MEs per mm.

to Figure 20 . 8 aFigure 16

 20816 Figure 16 Shear stress at x=L as a function of the displacement jump for the three mixedmode loading conditions.

Figure 17 Figure 18

 1718 Figure 17 Local mixity angle at x = L as a function of the damage parameter at x = L.

Figure 19

 19 Figure 19 Computed peel stress at x = L as a function of the displacement jump under mixed-mode loading condition a = p/16 for n = 1 and n = 2 with a mesh density of eight MEs per mm.

Figure 20

 20 Figure 20 Computed shear stress at x = L as a function of the displacement jump under mixed-mode loading condition a = p/16 for n = 1 and n = 2 with a mesh density of eight MEs per mm.

1 ) 2 ) 4 )

 124 Under pure mode II, the adhesive peel stress vanishes, so that: Under pure mixed-mode I/II, the adhesive peel and shear stresses are given by: Rev. Adhesion Adhesives CC BY-NC -Creative Commons DOI: 10.7569/RAA.2018.097301

  d u displacement jump (mm) of the interface along the x-axis d ue displacement jump (mm) of the interface along the x-axis at initiation d uf displacement jump (mm) of the interface along the x-axis at propagation d v displacement jump (mm) of the interface along the y-axis d ve displacement jump (mm) of the interface along the y-axis at initiation d vf displacement jump (mm) of the interface along the y-axis at propagation l norm of displacement (mm) of the interface l e norm of displacement jump (mm) of the interface at initiation l f norm of displacement jump (mm) of the interface at propagation n adherend Poisson's ratio q j bending angle (rad) of the adherend j around the z direction y local mixity angle (rad) (x,y,z) system of axes CZM cohesive zone model DCB double cantilever beam ENF end notched flexure FE Finite Element MCB mixed mode cantilever beam ME macro-element

in mm b in mm e in mm t in mm L in mm

  MPa, instead of 20 MPa. The choice of keeping S max to its original value instead of T max does not change qualitatively the results provided in his paper. The material parameters of the adhesive layers are given in Table2.

	0	4	0.2	8	100
		maximal peel stress S max = 36.6
					DOI: 10.7569/RAA.2018.097301
				CC BY-NC -Creative Commons	5

Table 2

 2 Adhesive material parameters.

	Stress		Mode I	Stress		Mode II
	S max	G Ie	G Ic	T max	G IIe	G IIc
	kI			kII
	δ ve		δ vf	δ ue		δ uf
			Displacement jump			Displacement jump

G Ie in N/mm G Ie in N/mm G IIe in N/mm G IIe in N/mm

	0.76	3.128 E-2	3.464 E-2	2.30

d ve in mm d vf in mm d ue in mm d uf in mm

	1.71 E-3	4.15 E-2	7.28 E-3	1.77 E-1

k I in GPa/mm S max in MPa k II in GPa/mm T max in MPa

	21.4 DOI: 10.7569/RAA.2018.097301 36.6	3.57 CC BY-NC -Creative Commons

  strain energy release rate (energy per unit area: mJ or N/mm) in peel G II strain energy release rate (energy per unit area: mJ or N/mm)in shear G Ic adhesive fracture energy (energy per unit area: mJ or N/mm)in peel G Ie adhesive elastic strain energy stored (energy per unit area: mJ or N/mm) in peel G IIc adhesive fracture energy (energy per unit area: mJ or N/mm)in shear G IIe adhesive elastic strain energy stored (energy unit area: mJ or N/mm)in

	A j	extensional stiffness (N) of adherend j
	B j D j	extensional and bending coupling stiffness (N.mm) of adherend j bending stiffness (N.mm 2 ) of adherend j
	E	adherend Young's modulus (MPa)
	F	magnitude of applied force (N)
	G I
		shear
	H	magnitude of applied displacement (mm)
	J	J-integral parameter
	K BBe	elementary stiffness matrix of a bonded-beam element
	L	length (mm) of bonded overlap
	M j	bending moment (N.mm) in adherend j around the z direction
	N j	normal force (N) in adherend j in the x direction
	S	adhesive peel stress (MPa)
	S max	maximal adhesive peel stress (MPa)
	T	adhesive shear stress (MPa)
	T max	maximal adhesive shear stress (MPa)
	V j	shear force (N) in adherend j in the y direction
	a	crack length (mm)
	b	width (mm) of the adherends
	d	damage parameter
	e	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	power usd in the adhesive material law
	n_ME number of macro-elements
	t	adherend thickness (mm)
	u j	displacement (mm) of adherend j in the x direction
	v j	displacement (mm) of adherend j in the y direction
	D	overlap length (mm) of a macro-element
	D j	characteristic parameter of adherend j in N 2 .mm 2
	a	angle (rad) used for the definition of the load application in MCB test
	b	mixed-mode parameter
	d t	numerical time step (s)

CC BY-NC -Creative Commons

29 DOI: 10.7569/RAA.2018.097301

Acknowledgement

The author affiliated to Sogeti High Tech gratefully acknowledges the engineers and the managers involved in the development of JoSAT (Joint Stress Analysis Tool), which is an internal research program. Both first authors indicate that this paper is a contribution to the initiative TACCOS; TACCOS stands for "Toulouse Adhésion Cohésion COllage Structural, de la Chimie à la Mécanique" which means "Toulouse Adhesion Cohesion Structural Adhesive Bonding, from Chemistry to Mechanics" in English.

CC BY-NC -Creative Commons

DOI: 10.7569/RAA.2018.097301 thought that the assessment provided by the direct method could be restricted to a particular combination of geometrical parameters, material parameters and boundary conditions under consideration.

Summary and Conclusions

In order to identify the parameters of CZM for thin adhesive layers, the inverse method and the direct method have been tested through a numerical test campaign. The numerical test campaign is based on the ME technique. A new approach for a fast formulation of the elementary stiffness matrix of bonded-beams element is provided in this paper. It can be concluded firstly that both inverse method and direct method are able to predict the adhesive constitutive evolution in both pure modes I and II. Secondly, the direct method provides accurate predictions under mixed-mode loading conditions. Thirdly, the inverse method fails in the assessment of the adhesive shear and peel stresses as a function of displacement jump under mixed-mode loading, with the particular CZM (but having a classical shape) used in this paper to model the adhesive layer in the bonded overlap. The use of the inverse method should then be restricted to the assessment of the adhesive constitutive evolution under pure modes only. Considering the J-integral framework, other restrictions should be considered when using the inverse method such as loading and unloading scheme or material time-dependent characteristics. Although the direct method appears as attractive, it could be difficult to apply in practice. Indeed, third-order and fourth-order differentiations of measured displacements are required, so dedicated methods to process the recorded signals should be employed. Under pure mode I and pure mode II, the order of differentiation could be reduced if the strain field on the external skin of adherends is recorded [START_REF] Jumel | Instrumented end notched flexure -Crack propagation and process zone monitoring. Part I: Modelling and analysis[END_REF]. From [START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF], it is clear that the experimental measurement of values such as applied force as a function of displacement jump is not sufficient for an accurate assessment of debonding problems. Nevertheless, assuming that the practical means are available to post-process the experimental test results for a reliable assessment of successive differentiations of measured displacements, the direct method would offer the possibility to investigate the CZM for thin adhesive layers. In particular, the experimental characterization of the shapes of constitutive adhesive evolutions under pure mode I, pure mode II and mixed-mode I/II could lead to validate the CZM available in the literature or to develop other ones by changing some of the underlying hypotheses. For this objective, the shape of the constitutive adhesive stress is of highest interest. That is why the authors of this paper are currently investigating this topic. The use of a novel experimental test technique [START_REF] Costa | Static assessment of the mixed-mode behavior of three epoxy adhesives[END_REF] allowing for the application of a wide range of mixed-mode ratios could be considered in conjunction with the direct method. It is indicated that the use of a series of loading and unloading schemes could help in a refined understanding of the mechanical behaviour of thin adhesive layers. Finally, when an adhesive layer needs to be characterized, its material behaviour is unknown

Appendix A

The expressions for the adhesive peel and shear stresses including the coupling stiffnesses are obtained as explained in section 3.2 by using the constitutive equations of adherends (Eq. ( 2)) and the local equilibrium equations of adherends (Eq. ( 4)).

Under pure mode I, the adhesive shear stress vanishes, so that:


S Moindjie

J L Autran

D Munteanu

G Gasiot

P Roche

Multi-Poisson Process Analysis of Real-Time Soft-Error Rate Measurements in Bulk 65nm and 40nm SRAMs

Altitude and underground real-time soft error rate (SER) measurements on SRAM circuits have been analyzed in terms of independent multi-Poisson processes describing the occurrence of single events as a function of bit flip multiplicity. Applied for both neutron-induced and alpha particle-induced SERs, this detailed analysis highlights the respective contributions of atmospheric radiation and alpha contamination to multiple cell upset mechanisms.

It also offers a simple way to predict by simulation the radiation response of a given technology for any terrestrial position, as illustrated here for bulk 65nm and 40nm SRAMs.

Introduction

Real-time soft error rate (RTSER) technique is an experimental reliability method to determine the soft error sensitivity of a given component, circuit or system [START_REF] Ziegler | SER -History, Trends and Challenges[END_REF][START_REF] Autran | Soft Errors: from particles to circuits[END_REF][START_REF] Cannon | SRAM SER in 90, 130 and 180 nm Bulk and SOI Technologies[END_REF][START_REF] Heijmen | Altitude and Underground Real-Time SER Tests of Embedded SRAM[END_REF][START_REF] Autran | Soft-errors induced by terrestrial neutrons and natural alpha-particle emitters in advanced memory circuits at ground level[END_REF]. The methodology consists in a direct observation of a large number of devices working in parallel under standard operating conditions and exposed to ambient background radiation at ground level. From the detection and counting of single events upsets occurring during the experiment, one can estimate the soft error rate (SER) of the devices under test. In case of SRAM circuits, the result can be given in terms of "event" SER, "total bit flip" SER, "multiple cell upset" (MCU) SER, and can be numerically expressed in FIT/Mbit (one FIT equals one failure per 10 9 hours of operation) [START_REF] Autran | Introductory Invited Paper -Real-time soft-error rate measurements: A review[END_REF][START_REF] Autran | A Review of Real-Time Soft-Error Rate Measurements in Electronic Circuits[END_REF].

For modern CMOS circuits at terrestrial level, SER is the result of a weighting of two main components due to: i) the atmospheric particle flux and particularly the high-energy neutrons interacting with IC materials; ii) the traces of radioactive elements present in the circuit materials, mainly U or Th contaminations at sub-ppb concentrations. These two SER components are respectively labeled neutron-SER (n-SER) and alpha-SER (a-SER). Neutron-SER is proportional to the intensity of the high-energy neutron flux at test location whereas alpha-SER is fixed and linked to the concentrations of radioactive impurities. The measured SER at a given test location is generally expressed as [START_REF] Autran | Introductory Invited Paper -Real-time soft-error rate measurements: A review[END_REF][START_REF]Measurement and Reporting of Alpha Particles and Terrestrial Cosmic Ray-Induced Soft Errors in Semiconductor Devices[END_REF]:

SER test location = AF × n-SER NYC +α-SER (1)
where AF is the so-called acceleration factor [START_REF]Measurement and Reporting of Alpha Particles and Terrestrial Cosmic Ray-Induced Soft Errors in Semiconductor Devices[END_REF] which define the relative neutron flux at test location with respect to sea-level (reference New York City, NYC) and n-SER NYC is the value of the neutron-SER normalized for NYC conditions. From Eq. (1), it is obvious that a-SER can be directly measured during an underground experiment (AF = 0) and that measurements at two different locations is sufficient for determining the two unknown quantities n-SER NYC and a-SER.

In this study, altitude (on the ASTEP platform [START_REF] Autran | Altitude SEE Test European Platform (ASTEP) and First Results in CMOS 130nm SRAM[END_REF][START_REF] Autran | ASTEP (2005-2015): Ten Years of Soft Error and Atmospheric Radiation Characterization on the Plateau de Bure[END_REF]), underground (at LSM, Modane [START_REF] Autran | Altitude and Underground Real-Time SER Testing of SRAMs Manufactured in CMOS Bulk 130, 65 and 40 nm[END_REF][START_REF] Autran | Combined Altitude and Underground Real-Time SER Characterization of CMOS Technologies on the ASTEP-LSM Platform[END_REF]) and sea level (at IM2NP laboratory, Marseille) real-time soft error rate (SER) measurements on bulk 65nm SRAM circuits have been analyzed in terms of independent multi-Poisson processes describing the occurrence of single events as a function of bit flip multiplicity. This approach has been applied for estimating the event rates induced by atmospheric neutrons and alpha-particle emitters as a function of event multiplicity. We used these extracted results to predict by simulation the radiation response of the 65nm SRAMs for a third test location at sea-level (Marseille, France), illustrating the advantage of this method not only to evaluate the global SER value but also to accurately estimate the occurrence of MCU events as a function of their event multiplicity. Additional results

concerning altitude and sea level measurements on bulk 40nm SRAMs extend the proposed analysis to the specific case of very large multiplicity events.

Theory

The occurrence of soft errors in electronics is a continuous time stochastic process. The principle of a RTSER experiment is to count the number of random events N(t) that have occurred up to some point t in time over a large population of circuits (the number of events is therefore small relative to the number of memory elements considered in the test). In such a real-time approach, soft errors are not permanent because they are eliminated when new data is written after detection: a RTSER test is thus equivalent to a life-test with replacement, i.e. in which a failing device is replaced with a new device immediately upon failure detection.

Assuming that soft errors are also random in space (within the circuits), independent of each other and occur with a fixed event rate l, the counting process {N(t), t ≥ 0} in a finite interval of time t consequently obeys the Poisson(lt) distribution [START_REF] Kingman | Poisson Processes[END_REF]:

P(N (t) = n) = e -λt (λt) n n! , n = 0, 1, 2 (2)
The maximum likelihood estimate for l is simply the number of detected events N divided by the experiment duration T and eventually normalized with respect to the number of units on test. Expressing this later quantity in Mbit for a memory test and time in MBit×h, l can be directly estimated from:

l = N/ Σ ( 3
)
where Σ is the number of MBit×h cumulated at time T given by:

Σ = MEM (t)dt 0 T ∫ (4)
where MEM(t) represents the number of Mbit under test at time t. This quantity may be variable in case of device automatic disconnection by the system due to abnormal current consumption or other anomaly detected at device, daugthercard or mothercard levels.

The counting process, previously introduced, can be applied to "single events" (whatever the type of these events), or more specifically to single bit upsets (SBU) or to multiple cell upsets (MCU) defined by a given event multiplicity (noted MCU(q), q = 2, 3,.. corresponding to the number of bitcells simultaneously upset during a same particle interaction event). Let N 1 (t), N 2 (t),…. N i (t), … N M (t) be M independent Poisson processes with event rates l 1 , l 2 ,…, l i ,… l M respectively corresponding to the counting of SBU, MCU(2), …, MCU(i),… MCU(M) events. The total event counting process N(t) = S N i (t) will be the merging of all these Poisson processes and will be also a Poisson process with rate l = l 1 + l 2 +…+ l i … + l M , as illustrated in Fig. 1.

Performing such a multi-Poisson process decomposition of the global "single event" process occurring during a RTSER experiment will allow us to reveal the respective contributions of atmospheric neutrons and alpha-particle emitters in the MCU SER response of the technology under test, as illustrated in the following.

Real-time experiments

For the purpose of this study, we considered raw data of several RTSER experiments performed on: i)

Results and discussion

Bulk 65nm SRAMs

From raw data (bit flip distributions) shown in Figs. 2, we first extracted in Table 2 the event rate per hour and per megabit of memory (Eq. ( 3)) of single events as a function of the event multiplicity (an event of unit multiplicity corresponding to a single bit upset).

From Eq. (1) also valid in terms of event rate, one can consider that for a given test location:

λ i = AF × n-λ i +α-λi (5)
where l i , n-l i and a-l i are respectively the total, neutron and alpha event rates for events of multiplicity i. i=1

corresponds to a single bit upset (SBU), i=2 to a double cell upset, etc.…).

Considering Eq. ( 5) for both ASTEP (AF = 6.5 averaged over the duration of the experiment) and LSM (AF = 0 when atmospheric cosmic rays are totally screened) experiments, it is then possible to extract these n-l i and a-l i values for the tested technology. This extraction is reported in Table 3 and illustrated in Fig. 4. These results

show that, for this particular bulk 65nm technology, soft errors are dominated at sea level by the contribution of alpha-particle radioactivity present in IC materials up to events of multiplicity equal to 6. In particular, the alpha event rate is 6 times larger than neutron event for SBU, 7 times larger for MCU(2), 5 times larger for MCU [START_REF] Cannon | SRAM SER in 90, 130 and 180 nm Bulk and SOI Technologies[END_REF] and 4 times larger for MCU [START_REF] Autran | Soft-errors induced by terrestrial neutrons and natural alpha-particle emitters in advanced memory circuits at ground level[END_REF] and MCU [START_REF] Autran | Introductory Invited Paper -Real-time soft-error rate measurements: A review[END_REF]. Neutrons and alphas have finally similar event rate for MCU [START_REF] Heijmen | Altitude and Underground Real-Time SER Tests of Embedded SRAM[END_REF] whereas neutrons are found to dominate for larger MCUs.

In order to verify the validity of such extractions, we use on more time Eq. ( 5) to estimate the alpha and neutron event rates for Marseille location from data of Table 2 extracted by the combination of ASTEP and LSM results. Taking AF = 0.9 for Marseille, Table 4 shows the results of the so-called "predicted" event rates using this procedure. These values are compared to the experimental (i.e. measured) event rates and are used to "predict" the number of events vs. multiplicity for a fixed duration of 18,697 h corresponding to the measurements. Table 4 shows that these predicted values are close to the measured ones (within 10% of uncertainty for the number of detected events), demonstrating here the interest of such a multi-Poisson process analysis to calculate not only the global event rate of the circuits but also to quantify (ideally to predict) the respective contributions of MCU events as a function of event multiplicity.

Bulk 40nm SRAMs

A similar analysis has been applied to raw data related to the bulk 40nm SRAM experiments. In this case, only data related to altitude and sea-level experiments are available (Fig. 5). Table 5 summarizes the number of detected events and the estimated event rates as a function of event multiplicity. Compared to the 65nm technology, 40nm SRAMs exhibit larger multiplicity events, up to 16 and 22 cells in size, respectively detected during altitude and sea level measurements. As previously explained [START_REF] Autran | Real-Time Soft-Error Testing of 40nm SRAMs[END_REF], these large multiplicity (≥5) cell upsets always involve adjacent cells preferentially distributed in columns in the physical memory array; they are the consequence of charge diffusion-collection and bipolar amplification that propagate into the silicon well directly impacted by the ionizing particle at the origin of the observed MCU.

The neutron and alpha event rates extracted from previous data are also reported in Table 5 and are shown in Fig. 5. Due to an evident lack of statistics for large MCU events that become more and more rare when their size increases, the extraction of neutron and alpha event rates has been only performed for MCU below or equal to 12 in size. Figs. 4 and 5 clearly shows the evolution of the event rates between the two technologies as a function of the event multiplicity. First, one can observe a global decrease of the alpha event rate for the 40nm with respect to the 65nm; this reduction is especially important for SBUs (reduction factor ~2) and significant up to MCU(4).

Secondly, one can note a slight increase of the neutron sensitivity concerning both SBU and all MCU events (up to MCU( 7)) for the 40nm technology with respect to the 65nm one. This tendency, previously reported for the global neutron and alpha SER values [START_REF] Autran | Altitude and Underground Real-Time SER Testing of SRAMs Manufactured in CMOS Bulk 130, 65 and 40 nm[END_REF], can be here quantified for a given event multiplicity, which is important in terms of event prediction and occurrence process simulation.

Conclusion

In summary, we presented in this study a detailed analysis of real-time soft error rate measurements describing the occurrence of single and multiple cell upset events as independent Poisson processes running in parallel, one for each even multiplicity (i.e. event size expressed in number of cells simultaneously upset by a single particle).

Combining a minimum of two distinct experiments conducted in different terrestrial environments, this approach allowed us to estimate the rates of all the events induced by both neutron-induced and alpha particle emitters as a function of their multiplicity. The proposed method has been successfully applied to quantify the respective contributions of neutrons and alpha contamination to multiple cell upset mechanisms for two SRAM technologies manufactured in CMOS bulk 65nm and 40nm, giving new lighting on the evolution of the soft error rate with integrated circuit integration. as extracted in Table 5.

Figure and table captions

Figure 5.

Neutron and alpha event rates for the bulk 40nm SRAM technology as a function of event multiplicity as extracted in Table 5.

Table 1. Summary of the real-time soft-error rate (RTSER) experiments considered in this study.

Table 2. Event rate extraction from RTSER data for the altitude, underground and sea level experiments on bulk 65nm SRAMs.

Table 3. Extracted values of neutron and alpha event rates as a function of event multiplicity for the 65nm SRAM technology (from data of Table 2 using Eq. ( 5) applied for both altitude and underground conditions).

Table 4.

Comparison between extracted (Table 2) and predicted event rates (from Table 3 using Eq.( 5) applied for Marseille conditions) for the sea-level experiment.

Table 5. Event rate extraction from RTSER data for the altitude and sea level experiments on bulk 40nm

SRAMs. Extracted values of neutron and alpha event rates (using Eq. ( 5) considering altitude and sea level data) as a function of event multiplicity are also reported.

Process 1, N 1 (t), SBU

Number of Events

Event rate ( × 10 -7 h -1 Mbit -1 )

Number of Events

Event rate ( × 10 -

Figures 2 and 3

 3 Figures 2 and 3 shows the cumulative bit flips distributions versus test duration respectively obtained for the three 65nm experiments and the two 40nm experiments. Test has been conducted under nominal operating conditions: V DD = 1.2V for 65nm SRAMs (1.1V for 40nm SRAMs), room temperature, standard checkerboard test pattern. These results will serve as raw data for the multi-Poisson process analysis conducted in the next section.

Figure 1 .

 1 Figure 1. Merging M Poisson processes corresponding to the occurrence of single events of different cell upset multiplicities.

Figure 2 .

 2 Figure 2. Cumulated distributions of bit flips versus time for the altitude, underground and sea level experiments on bulk 65 nm SRAMs.

Figure 3 .

 3 Figure 3. Cumulated distributions of bit flips versus time for the altitude and sea level experiments on bulk 40 nm SRAMs.

Figure 4 .

 4 Figure 4. Neutron and alpha event rates for the bulk 65nm SRAM technology as a function of event multiplicity

  3,226Mbit of bulk 65nm Single-Port SRAMS between 2008 and 2017; ii) 7,168 Mbit of bulk 40nm Single-Port SRAMS between 2011 and 2017. Details about these experiments can be found in Refs.[START_REF] Autran | Introductory Invited Paper -Real-time soft-error rate measurements: A review[END_REF][START_REF] Autran | ASTEP (2005-2015): Ten Years of Soft Error and Atmospheric Radiation Characterization on the Plateau de Bure[END_REF][START_REF] Autran | Altitude and Underground Real-Time SER Testing of SRAMs Manufactured in CMOS Bulk 130, 65 and 40 nm[END_REF][START_REF] Autran | Altitude and Underground Real-Time SER Characterization of CMOS 65nm SRAM[END_REF][START_REF] Autran | Real-Time Soft-Error Testing of 40nm SRAMs[END_REF][START_REF] Autran | Real-time soft error rate measurements on bulk 40nm SRAM memories: a five-year dual-site experiment[END_REF]. Table

I summarizes these experiments conducted on the ASTEP platform (2,255 m of altitude), at the underground laboratory of Modane (LSM) and at sea-level in Marseille (IM2NP laboratory). 65nm altitude and underground experiments were conducted in parallel with two identical setups; the setup installed in Modane was transported in Marseille in 2014 for sea level measurements. The 40nm experiments were successively performed on ASTEP and in Marseille using the same setup. Both 65nm and 40nm experiments are always running in 2017 in Marseille.

)

  Insufficient event sta>s>cs

	Location ASTEP multiplicity Event	Bulk 65nm SRAMs Acceler. factor AF Effective cumulated hours of tests 01/21/2008 Start date Stop date 05/07/2009 6.5 11,278 Events ( × 10 -7 h -1 Mbit -1 ) Events 1 n-l 1 = 0.955 Number of Event rate Number of multiplicity ( × 10 -7 h -1 Mbit -1 ) Altitude experiment Underground experiment Bulk 40nm SRAMs Stop date factor AF 3/11/2011 6/10/2014 6.0 ( × 10 -7 h -1 Mbit -1 ) Events a-l 1 = 5.92 28,494 ( × 10 -7 h -1 Mbit -1 ) Event rate Number of Event rate ( × 10 -7 h -1 Mbit -1 ) hours of tests Sea level experiment Event Neutron event rate Alpha event rate cumulated Start date Acceler. Effective Bulk 65nm SRAMs Bulk 65nm SRAMs Bulk 65nm SRAMs -
	LSM 1 Marseille 2 3	04/11/2008 44 01/25/2015 03/13/2015 Ongoing * 7 5	0 l 1 = 12.1 2 0.9 l 2 = 1.93 3 l 3 = 1.38 4	57,058 109 n-l 2 = 0.096 18,697 24 n-l 3 = 0.112 12 n-l 4 = 0.112	_ l 1 = 5.92 a-l 2 = 1.30 _ 10/16/2014 Ongoing * 0.9 l 2 = 1.30 a-l 3 = 0.652 l 3 = 0.652 a-l 4 = 0.109	36 9 5	_ 21,676	l 1 = 5.97 l 2 = 1.49 l 3 = 0.829
	4	3	l 4 = 0.827 5	2 n-l 5 = 0.026	l 4 = 0.109 a-l 5 = 0.109	1		l 4 = 0.166
	5	1	l 5 = 0.276 6	2 n-l 6 = 0.026	l 5 = 0.109 a-l 6 = 0.109	1		l 5 = 0.166
	6	1	l 6 = 0.276	2	l 6 = 0.109	0		l 6 = 0
	7	1	l 7 = 0.276	0	l 7 = 0	0		l 7 = 0
				Moindjie et al. Table 2

Sea level measurements Event multi- plicity Number of Events Event rate ( × 10 -7 h -1 Mbit -1 )

			Bulk 40nm SRAMs
	Altitude experiment		Sea level experiment	Extracted event rates
	Event
	multiplicity	Measured Predicted Measured Predicted
	1	36	41	5.97	6.77
	2	9	8	1.49	1.39
	3	5	5	0.829	0.751
	4	1	1	0.166	0.207
	5	1	1	0.166	0.132
	6	0	1	0	0.132
	7	0	0	0	0.038
	Total events	52	57

* Results reported on 04/30/2017 Moindjie et al. Table1

n-l 7 = 0.004 a-l 7 = 0Moindjie et al. Table3


Guénolé Lallement
email: guenole.lallement@st.com

Fady Abouzeid

Martin Cochet

Jean-Marc Daveau

Philippe Roche

Jean-Luc Autran

A 2.7pJ/cycle 16MHz SoC with 4.3nW Power-Off ARM Cortex-M0+ Core in 28nm FD-SOI

This work presents a System-on-Chip designed for Energy-Harvested applications. It embeds an ARM Cortex -M0+ microcontroller, 4 KB RAM, 4 KB ROM, an ultra-low power frequency synthesizer, a custom power switch, and a Power Management Unit enabling Active and Sleep modes. The system fabricated in 28 nm FD-SOI technology achieves 2.7 pJ/cycle at 16 MHz during active mode, and the core consumes 4.3 nW during deep sleep mode. The system operates at a fixed voltage of 0.5 V, and can switch from Active and Sleep/Deep Sleep modes, adjusting the frequency from 16 MHz to 8 MHz or 32 kHz in one cycle upon energy availability. By combining frequency/power switching with extra Reverse Body-Biasing the system power consumption is reduced by 53% and 98% in respectively sleep and deep sleep modes.

I. INTRODUCTION

Pulsed by the wide variety of low cost, battery-operated and connected applications, Internet-of-Things is considered to be the next decade market with 30.7 billion devices expected in 2020 [START_REF] Lucero | IoT platforms : enabling the Internet of Things[END_REF]. Those actual systems' power consumption are in the milliwatt range, hardly sustainable by integrated energy harvesting solutions for autonomous sensor nodes or outdoor monitors [START_REF] Bahai | Ultra-low Energy systems: Analog to information[END_REF]. With the advent of new smart objects with perpetual operations, designers and manufacturers must go beyond the Koomey's law which conjectures that the amount of power required for a given computing load will fall by a factor of 2 every 1.5 years. [START_REF] Koomey | Implications of Historical Trends in the Electrical Eficiency of Computing[END_REF] By enabling efficient low-cost systems, this work aims to push mass market microcontrollers (MCU) down to the microwatts power consumption. The proposed design developed in 28 nm Fully Depleted Silicon-On-Insulator (FD-SOI) technology is based on optimized standard cells and memories in conjunction with design methodology improved for 0.5 V operation with strong tolerance to process variability. The integrated Power Management Unit offers an adaptive power reduction related to the core activity, which helps to achieve 2.7 pJ/cycle ACTIVE and 4.3 nW DEEP SLEEP mode power consumption. This confirms the System-on-Chip's compatibility with power budgets generated by energy-harvesting devices.

In Section II the chosen system architecture and implementation process is presented. The silicon measurements are then presented in Section III. Lastly, an adaptive power reduction analysis is provided based on the core activity in Sec. IV.

II. SYSTEM ARCHITECTURE AND IMPLEMENTATION A. System overview

As presented in Fig. 1, the System-on-Chip (SoC) designed is composed of three Ultra-Low Voltage (ULV) power domains:

• An Always-On (AO) integrating a custom Power Management Unit (PMU), a one-cycle switching frequency synthesizer (DDSS), and event triggered peripherals;

B. Frequency synthesizer operation

The proposed clock generator (Fig. 2) is a custom alternative to PLL which trades-off phase-locking for lower power and area [START_REF] Cochet | On-Chip 28nm FD-SOI digital frequency generator for SoC clocking[END_REF]. It operates on a free-running unregulated internal 32phases oscillator which period T RO is compared with a known 32 kHz reference via a simple counter. That stage output, W, is then used as the input command of a fractional phase-selection divider, providing output periods with 1/32 nd increments. The oscillator is custom designed, while the rest of the circuit is automatically placed and routed from standard cells. The total clock generator area is 981 µm 2 . The MEP was achieved at 0.5 V/16 MHz operating condition, using regular threshold-voltage and 34nm to 46nm gate-length standard cell libraries (cf. Fig. 3). Regular-threshold transistors were also selected with regards to the use of Reverse-Body Biasing (RBB) for static power consumption reduction in sleep states. The power switch was strictly designed to minimize the static current consumption. The Cortex -M0+ was modeled using a current load in the range of 1 µA to 10 µA, corresponding to the simulated current profile of the MCU for a given set of programs. The targeted on-resistance (R ON ) was defined with a 10 mV maximum voltage drop target, which induces a limited power impact margin of 2%-4%, that is R ON ∈ [1 kΩ; 10 kΩ]. The switch design is based on distributed LVT pMOS transistors with gate tied to the body, offering fast swing recovery and super shut-off performances (Fig. 4). The restoration time was characterized to remove the need for a feedback controller with acknowledgment, replaced by a hard-coded wait operation. These specifications ensure the minimum power consumption overhead due to the controller area and activity power savings, and the inherent switching time when the SoC goes into DEEP SLEEP. (cf. Sec.IV).

MEP

The memories used for ROM and RAM are based on 6transistors bitcells with ultra-low voltage capabilities, enabling operation at 0.5 V and retention at lower supply values. The full system was implemented with ultra-low voltage corners Regarding the yield, the cumulative percentage of working parts from 0.45 V to 0.5 V is reported in Fig. 5b, as extracted from the wafer. These silicon results show that the reference processor can be supplied down to 0.48 V, within a variability window of less than 2%.

B. Power evaluation

The power contribution of each parts of the system was measured separately at 0.5 V/16 MHz, plotted in Fig. 6 and the figures associated reported in the first row of Tab. I. The total energy consumption when the system is running is 2.67 pJ/cycle. The total power consumption due to the static current is 1.5 µW, dominated by the M0+ core.

IV. CORTEX M0+ WITH ADAPTIVE POWER REDUCTION A. Available low power modes

The ultra-low power application oriented Cortex -M0+ processor is a 32-bit ARMv6 architecture with a two-stage pipeline for improved response time and efficiency. Thus, it achieves lower power/higher performances when compared to the previous M0 core [START_REF] Yiu | The Definitive Guide to ARM R Cortex R -M0 and Cortex R -M0+ Processors[END_REF].

Architecturally, two sleep modes are proposed: normal SLEEP and DEEP SLEEP. In our implementation, the SLEEP behavior is defined by the clock gating of the M0+ core whereas the DEEP SLEEP enables the use of the power switch on the Power-Off domain. Based on these sleep mode features, a PMU was designed and merged with a Wake-up Interrupt Controller (WIC) as additional hardware level to enable mode transitions. With this architecture the SoC user has the opportunity -depending on their application requirements -to clock gate or power down the microcontroller and wake up on certain hardware events.

B. Resulting power savings

Table I demonstrates the leverage offered by the M0+ power modes combined with the PMU to reduce power consumption of the system. All measurements were performed at 0.5 V/16 MHz. Switching from ACTIVE to SLEEP power mode leads to reduce the core power consumption from 15.1 µW to 1.96 µW, reaching its lowest value of 4.3 nW in DEEP SLEEP. In low power modes, the fastest clock frequency is no longer necessary. Consequently by switching in one cycle the DDSS clock from 16 MHz to 8 MHz the power consumption of the Always-On (SoC) is reduced by 48% (13%). For further power saving, the DDSS is disabled in DEEP SLEEP and the 32 kHz reference clock is used. As reported in Fig. 7 this last technique decreases the AO/SoC power consumption by 95%.

Subsequently, trimming the SoC frequency offers room for Reverse Body-Biasing (RBB) usage leading to a SoC Fig. 8: Reduction of the leakage power using RBB.

By combining together, all these power saving approaches, we induced a power breakdown from ACTIVE to SLEEP with frequency scaling & RBB and DEEP SLEEP with power gating, frequency scaling & RBB of respectively 53% and 98%. This leads to a total power consumption of the SoC of 0.7 µW in DEEP SLEEP mode (cf. Fig 9). Lastly, the power modes selected and the associated actions are summarized in Table II. The time required to switch the processor between states results in longer response latency and consume a fair amount of energy. The most energy-efficient SoC inactive state depends of these power and timing penalties. Hence, to select the lowest power mode for a given timeout, it is necessary to evaluate the energy overhead due to these transitions. Fig. 10 gives an insight into power modes switching and the associated figures. Using these performance measurements, the designer can select the best mode according to the SoC activity. The IDLE state is given as the worst power saving mode. It ideally starts in 0 clock cycles and corresponds to an active polling mode of the MCU.

The switching time from ACTIVE to SLEEP or DEEP SLEEP and the associated energy are defined by design. They have been validated using RTL and Prime Time power simulations. This time results in a minimum number of cycle and energy reported in Fig. 10 using the red and orange dots, respectively for SLEEP and DEEP SLEEP modes. The slopes are given by the energy/cycle associated to each mode. Finally, the blue and purple dots define the time mode switching is beneficial in terms of energy.

Hence, for energy driven applications, the gains of entering the power-off modes is benchmarked: the SLEEP state can directly be selected. Indeed only 0.6 (≤1) cycle of overhead is reported. However, for sleep time over 55 cycles the DEEP SLEEP should be triggered because it leads to the lowest SoC energy mode.

In the case of duty cycled core operations, the sum of the times spent in ACTIVE mode T ON and sleep mode T OF F is constant. Knowing T ON through a simple counter gives us T OF F and therefore the best sleep mode can be chosen.

D. Detailed comparison with the state of the art

The performances of the system is compared with the latest state of the art 32-bit microcontrollers in Table III. This work combines a low power consumption SoC with excellent energy/cycle associated with a PMU allowing efficient dynamic power mode selection.

V. CONCLUSION This research presents a System-on-Chip fabricated in 28nm FD-SOI, optimized for energy-harvested applications and offering mass market MCU with very high efficiency: 2.7 pJ/cycle at 0.5 V/16 MHz in ACTIVE mode, and 4.3 nW MCU in DEEP SLEEP mode. A frequency synthesizer enabling one-cycle frequency scaling, a dedicated power switch, and a PMU offering efficient mode switching lead to reduce the SoC power consumption down to 0.7 µW. The system power signatures fit in the energy-harvesting ∼100 µW power budget, demonstrating autonomous system capabilities [START_REF] Bahai | Ultra-low Energy systems: Analog to information[END_REF].

Fig. 1 :

 1 Fig. 1: 28 nm FD-SOI SoC implementation and corresponding block diagram. Active area is 0.073 mm 2 .

Fig. 2 :

 2 Fig. 2: Block diagram the Direct Digital Sampling and Synthesis (DDSS) clock generator [4]. C. Low energy design considerations The Cortex -M0+ minimum active energy operating point (MEP) was extracted using synthesis loops. Various operating conditions and standard cell libraries were explored for V dd ∈ [0.45V ; 0.9V ], F ∈ [1 MHz; 750 MHz], threshold-voltage V th ∈ [Regular; Low], and gate-length L g ∈ [30 nm; 46 nm].The MEP was achieved at 0.5 V/16 MHz operating condition, using regular threshold-voltage and 34nm to 46nm gate-length standard cell libraries (cf. Fig.3). Regular-threshold transistors were also selected with regards to the use of Reverse-Body Biasing (RBB) for static power consumption reduction in sleep states.

FFig. 3 :

 3 Fig. 3: Subset of Cortex -M0+ Synthesis loops results for V dd ∈ [0.45V ; 0.9V ]. Using limited range of gate-lengths and Regular-Threshold Voltage (RVT) standard cell libraries leads to the optimal MEP, with an operating voltage of 0.5 V and frequency of 16 MHz.

Fig. 4 :Fig. 5 :

 45 Fig.4: Distributed power switches using LVT-transistor design and simulated performances at 0.5 V. Total area is 90 µm 2 . centered at 0.5 V, with adjusted clock derating, setup and hold uncertainties. The tolerance to variability was improved by checking the hold timing violations at low voltages, fast/slow process and high/low temperatures.III. SILICON MEASUREMENTSA. Benchmarking and measured performancesThe system was fabricated in 28nm FD-SOI technology, packaged and measured using a custom development board and a Kintex-7 FPGA from Xilinx. A C-code program running the Dhrystone benchmark, SLEEP and DEEP SLEEP operations is used as software for yield measurements on 80 dies at wafer level, and performances measurements on 5 packaged dies, at 25 • C. The speed performance obtained are given in Fig.5afor V dd ∈ [0.45 V; 0.65 V]. The mean clock frequency goes from 10 MHz at 0.47 V to 150 MHz at 0.65 V with the targeted 16 MHz at 0.5 V. The F max associated to a RBB bias voltage of 500 mV is also plotted. Its utilization in SLEEP and DEEP SLEEP modes is scheduled for static power consumption reduction (Sec. IV).

Fig. 6 :

 6 Fig. 6: Pie chart of the energy/cycle (left) and the leakage power (right) SoC's repartition.

Fig. 9 :

 9 Fig. 9: Total power breakdown of the SoC.

Fig. 10 :

 10 Fig. 10: Evolution of the energy consumption in pJ for selected power state.

  A retention domain combining 4 KB of RAM and 4 KB of ROM in the form of 2 SRAM banks. The PMU was designed as a Finite State Machine to enable ACTIVE and Standby Modes (SLEEP and DEEP SLEEP) at very low power cost. The peripherals included are an UART controller, a SPI interface, 10 GPIO ports, a Wakeup Interrupt Controller and the ARM Debug Access Port. They are memory-mapped to the AHB bus through the APB.

• A Power-Off domain that includes an ARM Cortex -M0+ core and a custom power switch; •

TABLE I :

 I ENERGY/CYCLE AND POWER BREAKDOWN

	State	M0+	Always-ON	DDSS	Memory	Units
	ACTIVE	0.94 15.1	0.39 6.23	0.89 14.5	0.44 6.96	pJ/c. µW
	IDLE 1	0.86 13.7	0.38 6.23	0.89 14.5	0.25 4.02	pJ/c. µW
	SLEEP	0.12 1.96	0.22 3.59	0.89 14.5	0.03 0.53	pJ/c. µW
	DEEP SLEEP	2.7e -4 4.3 nW	0.21 3.49	0.89 14.5	0.03 0.54	pJ/c. µW

1 IDLE given as the worst waiting option corresponding to active polling.

TABLE II :

 II SUMMARY OF POWER MODES AVAILABLE

	C. Adaptive mode selection according to time spent in a mode
	State option	ACTIVE	SLEEP +Freq. scal. +RBB	DEEP SLEEP +Freq. scal. +RBB
	AO clock	16 MHz	8 MHz	32 kHz
	DDSS	ON	ON	OFF
	M0+ clock	16 MHz	Gated	Gated
	Power switch	ON	ON	OFF
	Body bias	0 V	RBB 0.5 V	RBB 0.5 V

ACKNOWLEDGEMENTS

The authors would like to thank their colleagues David Bonciani and Janit Kumar, Manohara Mr and Amit Patel for wafer testing and testing board design respectively.