Jean-Sébastien Affagard

Florent Mathieu

Jean-Mathieu Guimard

François Hild

Identification method for the mixed mode interlaminar behavior of a composite thermoset using displacement field measurements and load data

Keywords: A. Carbon Fiber Reinforced Plastics (CFRPs), B. Delamination, C. Finite element analysis (FEA), D. Mechanical testing

To analyze the delamination behavior of a thermoset composite, the parameters of a mixed mode cohesive model are identified by combining experimental displacement fields measured via global digital image correlation and load data in a single cost function written within a Bayesian framework using Finite Element Model Updating. A sensitivity analysis shows that the load data are of prime importance to enable the sought parameters to be tuned. Two delamination tests are analyzed in a single procedure that is first decoupled to find good initial guesses, and then fully coupled to get a set that is compatible with both experiments. The sensitivity of the calibrated parameters to their initial guess is analyzed. A low model error is found thereby validating the proposed approach (enabling to choose the best strategy for identification steps).

Introduction

Composite laminates are extensively used in aeronautics and aerospace industries [START_REF] Irving | Polymer Composites in te Aerospace Industry[END_REF][START_REF] Baker | Composite Materials for Aircraft Structures[END_REF] and thus challenge industry on its ability to accurately predict through numerical simulations and properly identify the behavior of specimens and structures up to failure. One of the critical degradation mechanisms is delamination between plies [START_REF] Davies | Standard Test Methods for Delamination Resistance of Composite Materials: Current Status[END_REF][START_REF] Allix | Mesomodeling of delamination: towards industrial applications[END_REF]. To characterize the interfacial properties, various experiments have been proposed. For instance, double cantilever beam (DCB) tests [START_REF]ASTM. D5528 -01e3 Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites[END_REF] are considered to evaluate the mode I properties, and crack lap shear (CLS) configurations [START_REF]ASTM. D5868 -01 Standard Test Method for Lap Shear Adhesion for Fiber Reinforced Plastic (FRP) Bonding[END_REF] for a priori mode II dominant propagations. Global fracture properties (e.g., critical energy release rates, stress intensity factors) are classically determined by using point data (e.g., applied load and remote displacement [START_REF]Instruction Générale de Contrôle pour détermination du G1c. Aérospatiale standard[END_REF][START_REF] Aitm | Determination of interlaminar fracture toughness energy[END_REF]). An alternative route consists of utilizing full-field measurements [START_REF] Ireman | Damage propagation in composite structural elements-coupon experiments and analyses[END_REF][START_REF] Abanto-Bueno | Investigation of crack growth in functionally graded materials using digital image correlation[END_REF][START_REF] Fedele | Identification of adhesive properties in GLARE assemblies by Digital Image Correlation[END_REF][START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF].

Cohesive zone models (CZMs) are also used to predict delamination in composite materials [START_REF] Cui | A combined stress-based and fracture-mechanics-based model for predicting delamination in composites[END_REF][START_REF] Allix | Damage analysis of interlaminar fracture specimens[END_REF][START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF][START_REF] Chaboche | On the interface debonding models[END_REF][START_REF] Reddy | Modeling the initiation and growth of delaminations in composite structures[END_REF]. Many of them are extensions of models initially applied to brittle materials [START_REF] Barenblatt | The Mathematical Theory of Equilibrium of Crack in Brittle Fracture[END_REF], quasibrittle materials [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF], or ductile materials [START_REF] Dugdale | Yielding of Steel Sheets Containing Slits[END_REF]. Compared with the previous (global) approaches, CZMs account for crack initiation and propagation. Full-field measurements techniques are very useful to identify material parameters of such models. Döll and Könczöl [START_REF] Döll | Micromechancis of fracture: optical interferometry of crack tip craze zone[END_REF] have shown how optical interferometry can be used to identify the strength of Dugdale model for PMMA. Abanto-Bueno and Lambros [START_REF] Abanto-Bueno | Experimental Determination of Cohesive Failure Properties of a Photodegradable Copolymer[END_REF] used a multicamera system to determine the cohesive properties of a copolymer via digital image correlation (DIC) analyses at two different scales. Fedele et al. [START_REF] Fedele | Identification of adhesive properties in GLARE assemblies by Digital Image Correlation[END_REF] calibrated the parameters of modified Xu and Needleman [START_REF] Xu | Void nucleation by inclusions debonding in a crystal matrix[END_REF] model for a glare assembly by updating finite element simulations (FEMU) with Q4-DIC (i.e., global DIC based upon finite element discretizations of the displacement field using 4-noded quadrangles [START_REF] Besnard | Finite-element'' displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]). Hong et al. [START_REF] Hong | Cohesive-zone laws for void growth -I. Experimental field projection of crack-tip crazing in glassy polymers[END_REF] used a field projection method to extract cohesive parameters from full-field measurements and applied the technique to analyze cracking of PMMA. Ferreira et al. [START_REF] Ferreira | On the analysis of notched concrete beams: From measurement with Digital Image Correlation to identification with Boundary Element Method of a cohesive model[END_REF] determined the traction profile in concrete when boundary element simulations are combined with Q4-DIC measurements. The authors have shown that Hillerborg's model is a very good approximation of the traction/separation law. Similarly, Shen

and Paulino [START_REF] Shen | Direct Extraction of Cohesive Fracture Properties from Digital Image Correlation: A Hybrid Inverse Technique[END_REF] tuned the traction profile with FEMU-U (meaning that FEMU uses only displacement data) for a ductile adhesive. Fuchs and Major [START_REF] Fuchs | Experimental determination of Cohesive Zone Models for Epoxy Composites[END_REF] utilized the J-integral to determine the cohesive stresses in glass/epoxy composites. Réthoré and Estevez [START_REF] Réthoré | Identification of a cohesive zone model from digital images at the micronscale[END_REF] updated X-FEM calculations to identify the parameters of a trapezoidal cohesive law applied to PMMA.

The present paper is a follow up on the identification of global fracture mechanics parameters (i.e., crack tip position, energy release rate, and stress intensity factors) and their associated uncertainties using experimental/numerical analyses of DCB and CLS tests [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]. With such a procedure, the experimental mode mixity could also be analyzed. More local analyses will be performed herein by considering a built-in cohesive zone model in the commercial finite element code Abaqus [START_REF] Simulia | Abaqus Analysis User's Manual[END_REF]. As proposed for the identification procedure of global parameters [START_REF] Fedele | Identification of adhesive properties in GLARE assemblies by Digital Image Correlation[END_REF][START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF], kinematic boundary conditions are prescribed in the finite element code in addition to the applied load, which is seldom considered. The paper is organized as follows. Section 2 describes the material, the specimen geometries, the mechanical tests and the full field measurements. The identification method updating the Finite Element (FE) model is introduced in Section 3. The sensitivity analysis is presented in Section 4 where the results will bring some elements on the procedure for parameter identification. Section 5 deals with the identified results obtained on the studied material. Last, Section 6 discusses the whole identification procedure and the results.

Set-up and measurements

In a previous study, the characterization of delamination properties was performed on a 0/0° interface configuration of a thermoset composite material T700/M21 (low grade made from prepregs cured with industrial quality process [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]). Three analyses were conducted to evaluate the mode I and II interlaminar fracture properties through DCB (Double Cantilever Beam) and CLS (Crack Lap Shear) tests. These tests were analyzed separately. The experiments were monitored by a Canon EOS 350D camera with a Sigma lens of 180 mm focal length. Displacement fields were measured using Q4-DIC [START_REF] Besnard | Finite-element'' displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF]. The specificity of this technique is to use a finite element discretization with 4-noded bilinear elements (Q4) and thus to ensure a continuity of the measured displacement field. The commercial Correli STC code was used [31]. In the present work, only the DCB test and one standard CLS experiment are considered to calibrate the cohesive properties thanks to a weighted Finite Element Model Updating (FEMU) approach. One of the main points is to use both tests at the same time to identify the parameters of the cohesive model in only one single step.

DCB configuration

In this configuration, the sample is loaded in an opening mode as presented in Fig. 1a. The dimensions of the sample are b = 20 mm, h 1 = 2.17 mm, h 2 = 1.85 mm, L = 250 mm, and t = 20 mm.

The loading in the x-direction is performed along the perpendicular direction of the plies with a tension/compression servohydraulic testing machine. Before the experiment is started, the sample was pre-cracked over a length a. In practice, it is performed manually by the operator while carefully clamping the specimen ahead of the desired location. An alternative is to perform a first propagation step with the testing machine; this step is then not taken into account in the subsequent investigation. In the present case the first route was chosen. The sample is then loaded under a displacement controlled procedure in which the force is measured and an image is acquired at each of the 28 steps of interest (Fig. 1b). Figures 1c and1d show the picture in the reference configuration and the last image in the deformed configuration of the ROI (Region of Interest, see Fig. 1a). The different images are used to measure the displacement field via Q4-DIC. In the present case, the edge size of each Q4 element is equal to 16 pixels (≈ 200 µm). The longitudinal (x-axis) and the transverse (y-axis) displacement fields between these two images are shown in Figs. 1e and1f, where the pixel size is equal to 12 µm.

The presence of a crack is clearly visible in the image of the deformed configuration (Fig. 1d) and on the transverse displacement field (Fig. 1f). Figure 1g shows the mode I crack opening displacement for the first and the last analyzed picture. It illustrates the change of crack opening and the resulting compression, in the right part, where the displacement jump is negative. In the present case, the crack opening displacement is evaluated as the difference of mean displacement of the two layers of the sample.

CLS configuration

Figure 2a shows the loading for the CLS configuration. The dimensions of the sample are b = 10 mm, h 1 = 1.61 mm, h 2 = 1.57 mm, and L = 250 mm. It consists of a tensile test in a servohydraulic testing machine where the plies are aligned along the loading direction. In this configuration, a pre-crack is manually created in the specimen by the operator while carefully clamping the two arms of the specimen. Thirteen loading steps are analyzed (Fig. 2b) and the corresponding images are acquired.

Figures 2c and2d show the picture in the reference configuration and in the deformed configuration at the end of the analysis. In this present case the edge size of each Q4 element is equal to 32 pixels (≈ 200 µm). The physical size of one pixel is now 3.5 µm since a higher magnification was selected (i.e., to enable for a more local analysis of the cohesive zone model). In this configuration, a mixed mode condition is induced by the experiment itself [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF], which evolves during propagation. The presence of the crack is observed when analyzing the displacement fields. However, it is impossible to detect manually on the pictures since the mode I opening displacement is very small.

Identification procedure of cohesive parameters

3.1.

Mesh generation

The numerical model is implemented in the commercial finite element code Abaqus [START_REF] Réthoré | Identification of a cohesive zone model from digital images at the micronscale[END_REF]. A 2D Finite Element (FE) model is developed considering the plane stress assumption (Fig. 3a). Each part of the carbon-epoxy composite is meshed with CPS4 elements. The interface is modelled thanks to zero thickness 4-noded cohesive elements (COH2D4). The Dirichlet boundary conditions are based on the measured displacement field u meas from Q4-DIC analyses (Figs. 1 and2). Therefore, the whole experiment is not modeled but only the ROI monitored for DIC purposes [START_REF] Fedele | Identification of adhesive properties in GLARE assemblies by Digital Image Correlation[END_REF][START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]. In Ref. [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF], a minimum mesh refinement of two is advised. In the present case, a factor four on the refinement is selected to ensure an accurate description of the nonlinear stress shape behind the crack tip. It was checked that a finer mesh would lead to very similar results. In the numerical model, the prescribed displacements of these additional nodes are obtained with the Q4 interpolation functions of the measurement mesh. Since the Q4-DIC code Correli STC provides measured displacement fields using meshes made of Q4 elements, this interpolation is direct and does not add any additional error contrary to local DIC approaches.

Constitutive behavior

The behavior of the composite layers is assumed to be elastic. The elastic properties of each 0degree laminate are E 1 = 120GPa, E 2 = E 3 = 8.9GPa, G 12 = 5.3GPa, ν 13 = ν 12 = 0.33, ν 23 = 0.35 where 1 is the longitudinal direction of the ply, 2 is the transverse one, and 3 the out-of-plane direction. These values were obtained from previous identification campaigns where several tensile tests with loading and unloading sequences were performed on [0°], [±45°], [±67.5°] coupons tested on an electrohydraulic testing machine, in the same manner as Ref. [START_REF] Ladevèze | Damage modelling of the elementary ply for laminated composites[END_REF].

The cohesive zone model is characterized by a traction-separation law, which depends on the fracture energy and the strength of the interface (Fig. 3b). The area under the traction-displacement jump curves for a pure mode I or II is the respective fracture energy dependent on the critical crack

opening displacement, final i  i i c final i d t G      0 (1)
The nominal traction, t, consists of two components in 2D problems, namely, t n and t t that represent the normal and the shear tractions, respectively. The corresponding displacement jumps are δ n and δ t . Denoting by h 0 the original thickness of the cohesive element, the nominal strains are defined as

h i i 0    (2)
where i = n or t. The elastic behavior (Fig. 3b) of the interface is written as

                   t n tt nt nt nn t n K K K K t t   (3)
where the off-diagonal terms in the stiffness matrix (i.e., K nt ) are set to zero in order to express the uncoupled behavior between normal and shear components.

In this work, a built-in CZM in Abaqus code is chosen even though the final goal is to use "bell-shape" models [START_REF] Allix | Mesomodeling of delamination: towards industrial applications[END_REF][START_REF] Allix | Damage analysis of interlaminar fracture specimens[END_REF][START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF]. This is not detrimental to the present work whose purpose is the feasibility of the methodology rather than the details of the chosen CZM. For pure mode I or II loadings, after the interfacial normal or shear tractions reach their respective interlaminar tensile or shear strengths, the stiffnesses are gradually reduced to zero (Fig. 3b) via a damage parameter [START_REF] Allix | Damage analysis of interlaminar fracture specimens[END_REF][START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF][START_REF] Chaboche | On the interface debonding models[END_REF][START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]. As a first proposal and among the few equations proposed within the FE code, the coupling of the fracture energy under mixed mode condition is defined in terms of the Benzeggagh-Kenane criterion [START_REF] Benzeggagh | Measurement of Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites With Mixed-Mode Bending Apparatus[END_REF] so that the total fracture energy, C G , is expressed as

            T t C n C t C n C G G G G G G ) ( (4)
where

t n T G G G  
and  is set to 1 (meaning that the 2D propagation space is a triangle, which is valid as a first approximation).

Inverse identification method

To identify the cohesive parameters, which are gathered in vector θ, a cost function  (θ) is based on the comparison of the measured and simulated displacement fields and resultant forces

                                  2 1 _ 2 _ _ _ _ _ _ _ _ 2 _ _ _ _ _ _ _ _ 2 _ _ _ _ _ _ _ _ 2 _ _ _ _ _ _ _ 2 2 2 2 ~                                                                                                 DCB step
CLS step CLS dof N N F F F F N N N N N N N N N N U U U U N N N N N N N N N N F F F F N N N N N N N N N N U U U U N N N N N N N N DCB DCB CLS CLS  θ θ θ (5)
where dof N is the number of kinematic degrees of freedom,

Boundary condition selection

It is worth noting that in all analyzed cases, experimental (i.e., measured) Dirichlet boundary conditions are applied to the numerical model. Consequently, any mode mixity induced by the stacking sequence (e.g., if not symmetric) or by any delamination plane not exactly located at midthickness (i.e., h1  h2) and/or loading conditions will be accounted for. These real conditions are often not taken into account in classical methods. In a previous work, it was shown that global mode mixities could be extracted by following such an approach [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF].

A first identification route consists of prescribing the displacement on the whole boundary of the analyzed ROI in the CLS experiment as for the determination of global parameters [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]. Since displacements are prescribed, they will give rise to equilibrium gaps [START_REF] Claire | Identification of damage fields using kinematic measurements[END_REF]. The equilibrium gap is separated into 8 contributions. Along the laminate (x-) direction the right and left reaction forces should be equal to the applied load (Fig. 3c-d). The right and left forces should be equal to 0 as well as the upper and lower ones. The identification leads to results with displacement residuals close to levels obtained for a global model [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF] (i.e., a root mean square difference between measured and computed displacements ≈ 0.3 pixel). For the equilibrium constraint, the load residual is of the order of 200 N on each boundary except on the boundaries where forces are transmitted (right and left boundaries along the laminate direction). On these boundaries, the cumulated residual is greater than 1,400 N. The same conclusion is drawn for the DCB test.

For the DCB experiment, it is possible to circumvent this issue by comparing the resultant of the cohesive load with the measured load (Fig. 3c). Conversely, for the CLS experiment, the upper and lower boundary conditions are stress-free in the simulations. Thus, the equilibrium gap is only built with the resultant load on the right side, which is equal to that of the lower left boundary (Fig. 3d).

Sensitivity analyses

To determine the feasibility of the identification, a sensitivity (S) analysis is performed on the whole set of cohesive parameter θ k [START_REF] Tarantola | Inverse Problems Theory[END_REF], for all the spatial domain (x) and for all relevant time steps (t) of the tests. This pre-analysis is carried out by means of the sensitivity matrix and allows for the best choice of identification strategy (e.g. two sets of parameters identified consecutively or only one set identified in one step). These parameters are the 6 cohesive properties and for the DCB test the crack initiation force. The latter is needed to account for the fact that the interface has been initially loaded at a level that is unknown. The sensitivity of the displacement field and resultant forces, collectively denoted by P, is written as

) , , ( ) , , ( ) , ( k k k t t t k     x P x P x S    (6)
and quantifies the effect of a variation θ k of each parameter θ k on the measured quantifies P. In the following analyses, a parameters variation of 1% is prescribed. The higher the sensitivity the easier the identification of the corresponding parameter provided the levels of the sensitivity fields can be made greater than the measurement resolutions [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization[END_REF].

The following analysis is performed with the initial values of the sought parameters, namely, K nn = 39 GPa/mm, K tt = 15 GPa/mm, σ n = 15, MPa, σ t = 37 MPa, G I = 350 J/m², and G II = 1500 J/m². These parameters are representative of values used at Airbus Group Innovations for such class of materials. Moreover, a pre-loading at 60 N is chosen for the DCB test. This force corresponds to a preloading performed before the measurement to localize the crack tip. Figure 4 shows the displacement sensitivity maps of the last loading step of the DCB experiment for each parameter. As expected, very small influences of the tangential parameters (i.e., K tt , σ t , G II ) are observed on the displacements (< 10 -5 pixel). When compared with the standard displacement resolutions they are significantly smaller, thereby proving that this experiment does not give any indication on the shear contribution of the CZM. On the contrary, the normal parameters lead to higher sensitivities, in particular the strength σ n (10 -4 pixel in the longitudinal direction, 10 -3 pixel in the transverse direction), but also the stiffness K nn and the toughness G I in the longitudinal direction (10 -4 pixel).

However, these levels are still lower than the measurement resolution (greater than 10 -2 pixel). The fact that all three normal parameters exhibit the same type of sensitivity is due to crack propagation that occurs for the last 16 steps of the experiment [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF].

Without any additional information, there is no hope of identifying the cohesive parameters with the considered approach (i.e., cost function only driven by displacement). The same analysis is carried out for the load sensitivities in the loading direction for which measured data are available.

Figure 5 shows that the sensitivity levels of the normal parameters are higher than the load resolution (of the order of 1N) in particular for the strength and toughness and to a lesser extent for the stiffness. This result explains why it is mandatory (at least for this configuration) to use a specific cost function that combines displacement and load measurements.

A similar analysis is performed for the CLS test. Figure 6 illustrates the displacement sensitivity maps of the last step of the CLS experiment, and Figure 7 the load sensitivity for all loading steps. Extremely low sensitivities are observed for both toughnesses. This result is to be expected since it is believed that propagation was very limited for the last recorded picture [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF] and the present FEMU process averages all steps. Conversely, the sensitivities to tangential and normal stiffnesses and strengths are more pronounced. This is due to the fact that the experiment was not in pure mode II loading [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]. There are strong correlations between normal and tangential parameters when displacements alone are analyzed. This is to be expected since no absolute stress (or load) scale is accounted for. As for the DCB test, the load sensitivities are very informative since their level is again higher than the load resolution for four out of 6 parameters (Figure 7). In both cases, the load information is crucial if cohesive parameters are sought.

In order to illustrate the global sensitivity to each parameter, the sensitivity matrix [M], which corresponds to the approximate Hessian of the nonlinear least squares formulation [START_REF]ASTM. D5528 -01e3 Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites[END_REF], is expressed as

   j t ij i M   S S  (7)
where   j  S ˆ is either the vector gathering one component of the nodal displacements over space and time, or all the load sensitivities over time for given parameter variations. The sensitivity matrices are shown for the DCB and CLS experiments in Figure 8. A high level in the sensitivity matrix represents a significant influence of the observed parameter, and a high level of off-diagonal terms suggests inter-dependencies between parameters as well.

For the DCB test, the displacement along the vertical direction shows that only the parameters K nn and σ n are sensitive and strongly correlated (i.e., nonzero off-diagonal terms). For the horizontal direction, the parameters K nn , K tt , σ n and σ t are the most sensitive. Moreover, the parameters K nn and σ n are strongly correlated in addition to K tt with σ t . Concerning the reaction forces, only the normal parameters (i.e., K nn , σ n , G I ) are linked and very influent (i.e., significantly higher than the load variance).

For the CLS experiment, the parameter K tt is the most sensitive in the horizontal direction. In the vertical direction, K nn and K tt are the most sensitive and not too (anti-)correlated. Concerning the reaction forces, all the parameters (i.e., normal and tangential) are influential. This result confirms again the benefit to use such a combined cost function, i.e. mixing displacement and load measurements.

Identification of cohesive parameters

This section presents the results of the parameter identification. Several protocols are considered based on the previous sensitivity analysis. First, using only the CLS experiment (i.e., a priori mode II dominant [START_REF]ASTM. D5868 -01 Standard Test Method for Lap Shear Adhesion for Fiber Reinforced Plastic (FRP) Bonding[END_REF]) the tangential behavior is analyzed with a high stiffness and strength in the normal direction (i.e., only the part related to the CLS test is considered in the cost function ( 5) with three unknown parameters). Then, a similar analysis is performed for the normal parameters using the DCB experiment (i.e., only the part related to the DCB test is considered in the cost function ( 5) with 3 unknown parameters). Once a first set of the six parameters is obtained, they will be used separately to update them by following the same approach as before, considering a sequential identification protocol. Last, with the two previously determined levels of cohesive parameters, the complete set of cohesive parameters is simultaneously identified by analyzing both tests with the coupled cost function  introduced in Equation (5).

First series of independent identifications

This first part aims to identify the tangential parameters using the CLS test. In order to identify them, the normal behavior is kept constant and chosen close to a purely elastic behavior (i.e., K nn = 500 MPa/mm, σ n = 250 MPa, G I = 400 J/m²) with a high level for initiation. The results are presented in Table 1 for different sets of bounds and initial values. Even though the three different sets of parameters are not identical, their level of variation remains small and that of the cost function is very close. The values of each parameter with the lowest cost function is considered for the subsequent identification steps, namely, K tt = 32 GPa/mm, σ t = 17 MPa, and G II = 1200 J/m².

The identified toughness G II is consistent with the level found by a global elastic identification when extrapolated to the failure load (1300 ± 50 J/m² [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]). Let us also note that the cost function has a level that is equal to about 7 times the standard resolutions associated with displacement and load measurements. This is an indication of a model error since if no model error had occurred, the cost function would be unitary. Conversely, the fact that its level does not depart too much from 1 is

an indication that the model should not be disqualified at this stage. Second, the normal parameters are identified using the DCB test. As in the previous case, the tangential behavior is assumed constant and purely elastic (i.e., K tt = 500 MPa/mm, σ t = 250 MPa, and G II = 1300 J/m²). The results gathered in Table 2 show that depending on the initial set of parameters the converged solution is not identical. However, all lie within the chosen bounds. In the present case, the three solutions are not as close in terms of overall level of the cost function as for the CLS experiment. This is presumably due to the existence of secondary minima in the cost function since its level at convergence is very close for the three sets of parameters. The set of values corresponding to the lowest cost function is chosen for subsequent analyses, namely, F = 58 N, K nn = 36 GPa/mm, σ n = 15 MPa, and G I = 310 J/m².

The identified toughness G I has a level lower than that identified with an elastic analysis (410 ± 50 J/m² [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]) yet with an identical order of magnitude. The level of the cost function is very close to the standard resolutions associated with displacement and load measurements. A very small model error is obtained in the present case, thereby validating it in the present analysis.

Second series of independent identifications

In this second loop, the tangential parameters are tuned by considering the values of the previously identified normal parameters, which still remain constant in the present minimization (F = 58 N, K nn = 36 GPa/mm, σ n = 15 MPa, and G I = 310 J/m²). The results are presented in Table 3. As in the previous cases, depending on the initial guess the converged solution is different. The identified tangential parameters are K tt = 12 GPa/mm, σ t = 39 MPa, G II = 1100 J/m².

The identified toughness G II is still close to the level found with a global elastic identification (1300 ± 50 J/m² [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]). Interestingly, the level of the cost function has been reduced by 20% (i.e., from 7.5 to 5.9) with this step. This decrease is in part due to the change of the tangential parameters so that the minimization procedure is initialized with another set of parameters. This result confirms what was expected from the sensitivity analysis, namely, that normal and tangential parameters had the same level of influence for the displacements. In particular, the stiffnesses K nn and K tt were correlated as well the strengths σ n and σ t . In order to identify the normal parameters, the tangential behavior is described with the previously identified parameters (i.e., K tt = 12 GPa/mm, σ t = 39 MPa, G II = 1100 J/m²). The results are presented in Table 4. The identified tangential parameters are F = 58 N, K nn = 38 GPa/mm, σ n = 15

MPa, and G I = 350 J/m². The identified parameters and the cost function levels are close to those of Table 2 except for the stiffness K nn . Although the sensitivity is not strongly correlated between K nn and the tangential parameters (i.e., K tt , σ t , G II ), the K tt parameter is slightly correlated with K nn and σ n , which may explain the difference of identified values. The levels of the cost function have been slightly reduced with this additional step. Being already close to 1, it was not expected to see significant changes when compared to those observed for the tangential parameters.

Combined identification

The last identification step consists of simultaneously identifying the normal and tangential cohesive parameters using both tests in a single analysis. The cost function is weighted accordingly. The results are gathered in Table 5. The level of the cost function is about five times the resolution limit, which validates the chosen model. The fact that this value lies between the two levels found by the two-step sequential identification procedure is expected thanks to the weighting introduced in Equation ( 5).

It is worth noting that the initial guesses analyzed in this part are all slightly different from the set found in Tables 3 and4. The reason for this choice is to get some information about the sensitivity to the initial guess. The two lower residuals lead to parameters that are close (Table 5). The force used to initialize the simulation of the DCB test is very stable with a level of 59 ± 0.1 N. Among the six cohesive parameters, the level of confidence is very different. For the toughness G I its level (350 ± 15 J/m 2 ) becomes closer to that found with a purely elastic identification (410 ± 50 J/m² [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]). In the experiment, crack propagation is observed making its uncertainty very small. This phenomenon also explains the small fluctuation of the normal strength σ n (13 ± 1.5 MPa).

Conversely, the identification of the normal stiffness is less stable even in the normal direction (K nn = 18 ± 3 GPa/mm).

The mean level G II (1250 ± 1 J/m²) is in very good agreement with a purely elastic identification when extrapolating the results for the failure load (1300 ± 50 J/m² [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]). The shear strength σ t (28 ± 6 MPa) is more difficult to capture. The relative fluctuation of the tangential stiffness (K tt = 14 ± 4 GPa/mm) is of the same order of magnitude as that observed in the normal direction.

In all analyzed cases, three different initial guesses were considered. The aim was to investigate the sensitivity of the CZM parameters to such initial values, and the level of the identification residual at the end of the optimization procedure. In general the level of the latter is very close, which indicates that overall calibration quality is similar. However, the values of some parameters are not necessarily close (i.e., K nn , K tt and  t ), which proves that they cannot be calibrated very precisely with the two considered experiments. Conversely, the evaluation of F DCB , G I , G II , and  n is more robust.

Interestingly, the two energy release rates are very well captured and consistent with analyses based upon linear elastic fracture mechanics [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]. The two stiffnesses are more difficult to evaluate quantitatively. They call for more accurate displacement measurements, but an order of magnitude is obtained, which is new for these types of tests. When the experiment allows stable propagation to occur, the strength of the CZM model can be estimated (i.e.,  n with the DCB experiment). For the CLS test, the propagation steps were very limited and did not yield a robust estimate of the strength  t .

Conclusion

In this paper, two delamination tests are analyzed to determine the parameters modeling the cohesive behavior in mode I and II. A combined cost function is introduced to study both tests at the same time with a weighting based upon a Bayesian framework to account for all sources of uncertainty. In the present setting, the displacement correlations were neglected and only the standard displacement resolution was considered. It could be extended to include the various correlations thanks to the covariance matrix associated with the measured degrees of freedom [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization[END_REF][START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF].

A sensitivity analysis is developed to illustrate the identification feasibility of the cohesive parameters. For the DCB test, the analysis shows a high sensitivity only for the normal behavior that illustrates a mode I loading. For the CLS test, which is a priori expected to be mode II-dominant for nominal boundary conditions, the sensitivity has the same level between the normal and tangential parameters, thereby indicating a mixed mode configuration, which is confirmed by previous analyses based upon full-field measurements [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF]. One of the key conclusions of this analysis is that without load data the cohesive parameters are not identifiable with the two experimental configurations studied herein.

The identification procedure has been conducted by following a three-step approach. The first step consists of identifying separately and sequentially the normal and tangential cohesive parameters by analyzing one test at a time. The other parameters are set to levels representative of a pure elastic interfacial behavior. Once a first set of parameters is obtained, it can be used to initialize a second step that still analyses the two tests independently. The identification residuals allow the chosen model not to be disqualified. Last, the two tests are studied simultaneously and the seven parameters are varied at the same time. For this last step, the identification residuals amount to about 5 times the standard resolution of displacement and load data. This low level validates the cohesive model used herein. It is also shown that all the parameters can be identified and with an error indication that can be less than 1 % in the most favorable case up to 30 % in the most difficult instance. This result was obtained for different initial guesses of the sought parameters.

It is worth noting that having set the 7 parameters free directly, the identification would have been unstable. This is one of the reasons of the developed 3-step strategy thanks to the sensitivity analysis. The latter proves that the initial guess is very important. This calls for simplified and fast techniques to get initial estimates faster than by performing finite element simulations, which are very demanding when dealing with cohesive zone models. Since the mode I and II toughnesses identified with the present procedure are in good agreement when compared to linear elastic fracture mechanics analyses [START_REF] Mathieu | Identification of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF], such type of approach is beneficial and leads to a more robust identification route. Systematic sensitivity analyses with quantitative estimations of uncertainties and errors are useful whatever the mode mixity of tests or configurations.

The present framework may also be used to identify initiation and propagation parameters under more complex configurations (e.g., interfaces between [±] plies) since it is very generic (i.e., it does not rely on closed-form solutions whose applicability is restricted to very simple cases). It has been shown that it is possible not only to extract global propagation parameters (e.g., G I and G II ) but also initiation values. Future analyses will also have to be extended to cases of rapid propagation to study delamination under dynamic loading conditions. One critical issue is related to the way the finite element simulations are driven with measured displacement fields that are corrupted by measurement uncertainties. Spatiotemporal DIC analyses may be used to tackle such difficult problems [START_REF] Besnard | Analysis of Image Series through Digital Image Correlation[END_REF].

F.

  sensor N the number of load sensors and step N the number of loading steps for each experiment.   U ~ and   F ~ are the vectors gathering the measured displacements and forces, respectively,   U and   F , which depend on the parameters θ, the calculated displacement and force vectors,u and the standard displacement and force resolutions, respectively. In the present case, weighted residuals are considered that stem from a Bayesian foundation to account for the standard resolution levels associated with Q4-DIC and the load sensors. The underlying hypothesis is that noise is Gaussian and white. The quadratic form is the argument of the exponential in the Gaussian probability distribution. The additivity of the functionals is the counterpart of the statistical independence of the two measurements (i.e. load and displacements), which implies that probabilities are multiplied. The chosen normalization by the noise amplitude in all functionals guarantees that no additional prefactors are to be considered. If only noise is present (i.e., no model error), the expectation of  (θ) is equal to 1[START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization[END_REF][START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF][START_REF] Bertin | Integrated Digital Image Correlation applied to elasto-plastic identification in a biaxial experiment[END_REF].The standard displacement resolution pixel 04 0.  u  in both cases (or  0.5 µm for the DCB experiment, and  0.15 µm for the CLS experiment), whereas the standard load resolution is The identification of θ is performed through the minimization of the cost function through Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm[START_REF]Matlab. Matlab, the Language of Technical Computing, version R2015a: the MathWorks[END_REF].

FiguresFigure 1 :

 1 Figures

  Figure 1: (a) DCB sample with the position of the Region Of Interest (ROI) used to measure

Figure 2 :

 2 Figure 2: (a) CLS sample with the position of the Region Of Interest (ROI) used to measure

Figure 3 :

 3 Figure 3: (a) Finite element mesh in which the cohesive nodes are marked in red.

Figure 4 :Figure 5 :

 45 Figure 4: Displacement (expressed in pixel) sensitivity maps of the DCB test for the last loading step

Figure 6 :

 6 Figure 6: Displacement (expressed in pixel) sensitivity maps in the CLS test for the last loading step in

Figure 7 :

 7 Figure 7: Reaction force sensitivity as a function of the loading step in the CLS experiment for a

Figure 8 :

 8 Figure 8: Sensitivity matrix for the DCB experiment (left) for the displacement field in the horizontal

Table 1 :

 1 Tangential cohesive parameters identified when analyzing the CLS experiment for different bounds and initial values. The normal parameters are set to high values

		K tt (GPa/mm)	15 -70	10 -60	15 -80
	Bounds	σ t (MPa)	10 -60	5 -70	15 -70
		G II (J/m²)	700 -1800 700 -1800 700 -1800
		K tt (GPa/mm)	35	20	30
	Initial values	σ t (MPa)	18	23	32.5
		G II (J/m²)	1100	1600	1200
		K tt (GPa/mm)	34	32	47
	Identified values	σ t (MPa)	19	17	27
		G II (J/m²)	1100	1200	890
	Cost function	 (-)	7.55	7.53	7.57

Table 2 :

 2 Normal parameters identified when analyzing the DCB experiment for different bounds and initial values. The tangential parameters are set to high values

		F DCB (N)	56 -62	55 -63	57 -63
	Bounds	K nn (GPa/mm) σ n (MPa)	5 -30 5 -25	5 -45 7 -35	7 -30 5 -25
		G I (J/m²)	250 -500 250 -600 250 -500
		F DCB (N)	59,5	59	60
	Initial values	K nn (GPa/mm) σ nn (MPa)	25 15	15 22	8.5 13
		G I (J/m²)	350	350	310
		F DCB (N)	59	58	60
	Identified values	K nn (GPa/mm) σ nn (MPa)	28 13	37 15	21 13
		G I (J/m²)	350	310	360
	Cost function	 (-)	1.34	1.25	1.26
	5.2.

Table 3 :

 3 Tangential parameters identified when analyzing the CLS experiment for different bounds and initial values. The normal parameters are set as the mean values of the previous identification

			step
		K tt (GPa/mm)	5 -70	5 -80	5 -50
	Bounds	σ t (MPa)	5 -60	5 -70	5 -50
		G II (J/m²)	700 -1800 700 -1800 700 -1800
		K tt (GPa/mm)	35	30	32
	Initial values	σ t (MPa)	18	33	17
		G II (J/m²)	1100	1200	1200
		K tt (GPa/mm)	12	12	11
	Identified values	σ t (MPa)	39	39	39
		G II (J/m²)	1100	1100	1000
	Cost function	 (-)	5.924	5.925	5.926

Table 4 :

 4 Normal parameters identified when analyzing the DCB experiment for different bounds and initial values. The tangential parameters are set as the mean values of the previous identification

			step
		F DCB (N)	55 -63	57 -63	57 -63
	Bounds	K nn (GPa/mm) σ n (MPa)	1 -45 7 -35	7 -30 5 -25	7 -40 5 -25
		G I (J/m²)	250 -600 250 -500 250 -500
		F DCB (N)	59	60	60
	Initial values	K nn (GPa/mm) σ nn (MPa)	15 22	8.5 13	21 13
		G I (J/m²)	350	310	360
		F DCB (N)	58	59	60
	Identified values	K nn (GPa/mm) σ nn (MPa)	38 15	20 14	20 14
		G I (J/m²)	350	370	370
	Cost function	 (-)	1.19	1.21	1.20
	5.3.

Table 5 :

 5 Cohesive parameters identified when analyzing both experiments for different bounds and

		initial values
		F DCB (N)	57 -62	57 -62	57 -62
		K nn (GPa/mm)	10 -40	10 -40	10 -40
		K tt (GPa/mm)	9 -30	5 -30	5 -30
	Bounds	σ n (MPa)	5 -25	5 -25	5 -25
		σ t (MPa)	15 -45	15 -45	15 -45
		G I (J/m²)	250 -450	250 -450	250 -450
		G II (J/m²)	900 -1800 600 -1800 600 -1800
		F DCB (N)	59	58	60
		K nn (GPa/mm)	20	15	20
		K tt (GPa/mm)	15	20	12
	Initial values	σ n (MPa)	15	20	14
		σ t (MPa)	35	25	39
		G I (J/m²)	350	300	370
		G II (J/m²)	1300	1100	1100
		F DCB (N)	59	60	60
		K nn (GPa/mm)	23	16	21
		K tt (GPa/mm)	12	19	11
	Identified values	σ n (MPa)	12	15	12
		σ t (MPa)	33	22	34
		G I (J/m²)	340	350	360
		G II (J/m²)	1250	1250	1250
	Cost function	 (-)	4.92	5.10	4.91

Acknowledgement

The support of this research by "Agence Nationale de la Recherche" is gratefully acknowledged (VULCOMP Phase 2 project, Grant No. ANR-12-RMNP-0018)


A NEW MODELLING OF THE MULLINS EFFECT AND THE VISCOELASTICITY OF ELASTOMERS BASED ON A PHYSICAL APPROACH

Keywords: Natural Rubber, hyperelastic, Mullins effect, hysteresis, physical approach

 in which macroscopic constitutive equations are build in regards with the physics of polymeric chains. In this context of physical considerations, the aim of the present paper is to study independently each phenomenon involved in rubber-like materials and to assemble them in a global constitutive equation.

First, the hyperelastic behaviour of elastomers is modelled by the physical approach of [START_REF] References Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF], widely known as the eight-chains model. This model accurately reproduces the large strains elastic behaviour of elastomers under different types of deformation. Second, the hysteretic time dependent behaviour is approached by the model developed by [START_REF] Bergström | Constitutive modelling of the large strain timedependent behavior of elastomers[END_REF] that considers the separation of the network in two phases: an elastic equilibrium network and a viscoelastic network that captures the nonlinear rate-dependent deviation from equilibrium. This model is quite simple and successfully reproduces the rate-dependent hysteretic properties of elastomers. Last, as shown in the bibliography, the Mullins stress-softening effect can be considered as a damage phenomenon which only depends on the maximum stretch attained during the deformation history [START_REF] Govindgee | Mullins' effect and the strain amplitude dependence of the storage modulus[END_REF]. In the present approach, the physical theory of [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] based on an alteration of the polymeric network is adopted. This theory was introduced in the eight-chains hyperelastic model and successfully simulates the decrease of the material stiffness between the first and the second loading curves under cyclic loading. As these three models are based on the physics of the polymeric network, they are gathered in a new efficient constitutive equation. This model is able to reproduce simultaneously the Mullins effect and the time-dependent hysteretic behaviour of elastomers. Finally, the constitutive parameters of this new model are identified by fitting experimental data.

INTRODUCTION

The fatigue of elastomers has become a key point in automotive industry. Indeed, the prediction of the fatigue behaviour of automotive parts, such as engine mounts or bump stoppers, is a necessary prerequisite to the development of new products. In this context, the numerical modelling, for example the finite element method, takes an important role and the development of efficient constitutive equations for elastomers is of prime necessity.

Elastomers are known to exhibit a highly non linear elastic behaviour under static load; a rate dependent or viscoelastic behaviour including hysteresis under cyclic loading; and a socalled Mullins effect which can be described as a stress-softening phenomenon observed during the first cycles under fatigue loading conditions. Another phenomenon, the Payne effect, involves the thixotropic behaviour under dynamic loading [START_REF]Thixotropic behavior of rubber under dynamic loading histories: Equipments and theory[END_REF]. Our study is devoted to the modelling of the behaviour of a carbon black filled natural rubber. It is well known that the addition of filler particles leads to an increase of the material stiffness and to an improvement of the ultimate properties of the material (ultimate stress and maximum extensibility). Moreover, it has to be noted that the Mullins effect and the hysteretic nature of the rubber is amplified by the presence of fillers. Considering cyclic experiments, the Mullins effect is measured through the area between the first and the second loading curves and the hysteretic behaviour is characterised by the area between loading and unloading curves.

The aim of the present paper is the development of a rubber-like constitutive equation that is able to simulate both these two phenomena. A physical approach has been adopted for the development of the model. In section 2, classical hyperelastic constitutive equations are first recalled. Then, the models adopted for the Mullins effect and the hysteresis are independently presented. Finally, a combination of these two sub-models is developed. Section 3 presents the identification of the previous model. The comparison between experimental data and theoretical results is performed for both uniaxial tensile and pure shear tests. Section 4 is devoted to concluding remarks.

CONSTITUTIVE EQUATION

HYPERELASTIC MODEL

Rubber like materials are classically considered isotropic, incompressible and hyperelastic. They are characterised by a strain energy function W. Thus, the definition of a hyperelastic model consists in defining the form to W. Two different points of view can be adopted to define W: phenomenological or physical. In the former case, many authors proposed models based on strain invariants. [START_REF] Mooney | A theory of large elastic deformation[END_REF] developed a linear model of the strain invariants which was generalised by [START_REF] Rivlin | Large elastic deformations of isotropic materials -VII. Experiments on the deformation of rubber[END_REF]. Others forms of W were considered by [START_REF] Gent | Forms of the stored (strain) energy function for vulcanized rubber[END_REF] and [START_REF] Hart-Smith | Elasticity parameters for finite deformations of rubber-like materials[END_REF]. Principal stretches based model were proposed by [START_REF] Valanis | The strain-energy function of a hyperelastic material in terms of the extension ratios[END_REF] and [START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubberlike solids[END_REF]. All these constitutive equations are based on both experimental data and mathematical developments. In the later case, the definition of W is based on the physics of the polymeric chains network. [START_REF] Treloar | The elasticity of a network of long chain molecules (1 and 2)[END_REF] used Gaussian statistics applied to the network in order to describe the macroscopic behaviour of rubber-like materials. These physical considerations lead to the neo-Hookean constitutive equation. This model is not dedicated to large deformations. In order to overcome the limitations of the neo-Hookean approach, authors used non-Gaussian statistics to describe molecular chains deformations. [START_REF] Kuhn | Beziehunger zwichen elastischen konstanten und dehnungsdoppelbrechung hochelastischer stoffe[END_REF] described the stretching limit of the chains. Their developments are the basis of the chains models like the full network model of [START_REF] Wu | On improved 3-D non gaussian network models for rubber elasticity[END_REF], the three chains model of [START_REF] James | Theory of the elastic properties of rubber[END_REF], the four chains model of [START_REF] Flory | Network structure and the elactic properties of vulcanized rubber[END_REF] and the eight chains model of [START_REF] References Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF]. The network configuration corresponding with the eight-chains model is presented in Figure 1. In our study, we consider the eight chains model because of its simplicity and its good capability of describing the strain-hardening of the material. The corresponding expression of W is given:

      + = β β β sinh ln . Nl r N C W chain r (1) Where [ ]        = + + = = - Nl r l N r nk C chain chain r / ) ( 3 1 1 2 / 1 2 3 2 2 2 1 L β λ λ λ θ (2)
In which L is the Langevin function defined by

β β β 1 ) coth( ) ( - = L .

MULLINS EFFECT

It is well known that elastomers exhibit a reduction of stiffness after the first loading cycle of a fatigue experiment. This phenomenon was observed, first, by Bouasse ad Carrière (1903), studied in details by [START_REF] Mullins | Softening of rubber by deformation[END_REF] and is widely known as the Mullins effect. Under quasistatic condition, the viscoelastic behaviour of elastomers does not influence the material response and the behaviour can be considered as purely elastic and submitted to the Mullins effect. The corresponding stress-strain curves are presented in Figure 2. The undamaged (or virgin) material is first stretched to the extension ratio λ I and the stress follows the path I in Figure 2(b). Then the unloading from λ I to 0 follows the path I'. The second loading from 0 to λ II >λ I first follows the path I' until λ=λ I then it follows the path II. The second unloading from stretch ratio λ II to 0 follows the path II' which is different than path I'. At a given stretch, the stress on II' is lower than the stress on I'. Repeating this process, the loading path corresponding to the increase of stretch from 0 to λ II is the path that joins II' and the part III of the virgin curve. Finally, the corresponding unloading follows the path III'. These observations illustrate Mullins' remarks, as accommodation occurs only from strain lower than the maximum strain obtained in the material history, and when strain reaches the maximum strain ever attained, the behaviour becomes the one of an undamaged material.

Figure 2. Mullins effect under quasi-static conditions

Different approaches have been proposed to model the Mullins effect, and there is no unanimous explanation of the physical causes of it. [START_REF] Mullins | Theoretical model for the elastic behavior of filledreinforced vulcanizated rubbers[END_REF] proposed a model in which the rubber could exist in either one of two possible phases, so-called hard and soft phases. They consider that the rubber is initially in its hard phase, and that it would degrade into soft phase as it deforms. A damage parameter is introduced to measure the transition between the two phases. This parameter is related to the strain amplification in the remaining hard phase of the rubber. [START_REF] Mullins | Softening of rubber by deformation[END_REF] suggested that this stress-softening is due to disentanglements of the network chains produced by the breakdown of filler particles/ rubber interactions. The same idea was suggested by [START_REF] Bueche | Molecular basis of the Mullins effect[END_REF][START_REF] Bueche | Mullins effect and rubber-filled interaction[END_REF] that developed a probabilistic model of contact between chains and fillers.

More recently, [START_REF] Chagnon | Experimental identification and rheological modeling of the Mullins' effect for carbon black-filled rubber[END_REF] analysed the Mullins effect from an experimental point of view. They studied the stiffness decrease under deformation during tensile tests in order to determine a damage evolution law. Following this experimental work, [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] considers that the Mullins effect is mainly due to the breakage of links between rubber chains. During stretching, they observed the evolution of chains length between two junction points and concluded that the number of monomers between two junctions points increases. As a consequence, the chain density decreases. The evolution of these two parameters, said the number of monomers per chains and the chain density, are introduced in the eight-chains model. Constitutive parameters of the original eight chains model are replaced by functions of the maximum strain and a phenomenological law is proposed to simulate the increasing number of monomers between two junction points:

( ) max 1 0 exp chain N N N λ = (3)
The conservation of the number of monomers per unit of volume implies that the two parameters of the eight-chains model must be related by:

Cte C N r = (4)
For details on this model, the reader can refer to [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF].

HYSTERESIS

Under cyclic loading conditions, elastomers reveal their viscoelastic nature. Two different viscous phenomena can be identified. The first one is the hysteretic behaviour: after accommodation to the Mullins effect, loading and unloading curves do not follow the same path, as shown in Figure 3. The second phenomena is the well-known long-term relaxation.

Both phenomena are due to viscoelasticity but we can say that they do not take place at the same time scale. Consequently, hysteresis and long-term relaxation can not be modelled by the same method. In the present work, our study is restricted to the hysteretic behaviour. In the last few years, several studies attempted to model the hysteretic nature of rubber-like material both experimentally and theoretically [START_REF]A constitutive model for carbon black filled rubber: experimental investigation and mathematical representation[END_REF][START_REF] Bergström | Constitutive modelling of the large strain timedependent behavior of elastomers[END_REF][START_REF] Banks | Modeling of nonlinear hysteresis in elastomers under uniaxial tension[END_REF][START_REF] Miehe | Superimposed finite elastic-viscoelastic-plastoelastic stress response with damage in filled rubbery materials: experiments, modeling and algorithmic implementation[END_REF]. These models are based on the same approach: the superposition of a rate-dependent response on an equilibrium elastic response. They only differ in the details of the stress-strain relationships.

Here, we adopt the model developed by [START_REF] Bergström | Constitutive modelling of the large strain timedependent behavior of elastomers[END_REF]. As mentioned above, their model is based on the decomposition of the elastomeric behaviour into two responses, an elastic one and a rate-dependent one. The elastic part of the behaviour is described by the eight chains model, presented above. The non-equilibrium part reproduces the chains slippage phenomenon which ensures that the material behaviour is similar under loading and unloading conditions. This slippage of the chains is modelled by a classical tube approach. This constitutive equation is schematised in Figure 4.

More precisely, the model is governed by three major equations; the two eight-chains models stress-strain relationships (presented above) and the slippage law given by:

[ ] m B C Bp chain B C τ λ γ 2 1 1 - = • (5)
where B C and m are three parameters to be determined.

Elastic Part

Hyperelastic behaviour Eight chains model

Rate-dependent Part

Viscous damping Tube approach

Hyperelastic behaviour Eight chains model Figure 4. The Bergström and Boyce model.

PRESENTATION OF THE FINAL CONSTITUTIVE EQUATION

In order to consider both Mullins effect and hysteresis simultaneously, it is necessary to combine the different models described in the previous paragraphs. The key point is that these models are based on the same physical approach. As we considered that the Mullins effect is the stiffness decrease between the first and the second loadings curves, it only affects the material behaviour for loading paths. The elastic part of Bergström and Boyce model corresponds to loading paths and the rate-dependent part to unloading paths, so that our modelling of the Mullins effect, detailed in Marckmann et al. (2000), is only introduced in the elastic part of the Bergström and Boyce constitutive equation. Finally, this final assembled model can be represented schematically as shown in Figure 5. Note that it depends on eight parameters.

Elasticity with Mullins effect

Hyperelastic behaviour with Mullins effect Marckmann et al. model

Rate-dependent Part

Viscous damping Tube approach

Hyperelastic behaviour Eight chains model

PARAMETERS IDENTIFICATION AND SIMULATION

Uniaxial tension and pure shear cyclic experiments were conducted on a carbon black filled natural rubber. These tests were performed at different maximum strains. A global identification of the parameters, said using uniaxial tensile and pure shear data, leads to the values of the material parameters presented in Table 1.

Parameters Values

Mullins effect model

Table 1. Material parameters

To verify these results, the experiments are simulated. Comparisons between numerical results and experimental data are presented in Figure 6 for uniaxial tensile tests and in Figure 7 for pure shear tests. The numerical results show that our model is able to simulate both Mullins effect and hysteresis. Nevertheless, the quantitative results are not as good as expected. In our model, the Mullins effect highly dominates the hysteretic effect of stress-strain curves. For a given loading stretch, the first cycle is well reproduced by the model. But the second cycle exhibits a too small amount of hysteresis, said the cycle is flattened. This difficulty can be explained by the fact that the first cycle includes both the Mullins effect and the hysteresis, so that the stiffness decreases between loading and unloading paths is well reproduced. For the second cycle, the material is accommodated to the Mullins effect and the hysteresis is underestimated. This might be a consequence of the fitting methodology that was performed globally. It might be more appropriate to identify various effects separately.

Stretch in tensile direction

Finally, note that the relaxation behaviour is not taken into account in the present model: the stiffness decreases after each successive loop due to viscoelasticity is not considered. In order to improve this approach, it is necessary to add long-term time dependent phenomena. Indeed, the relaxation highly influences the stress-strain curve during the first loading cycles and can be neglected for next cycles. Addition of relaxation parameters in the model would reduce the influence of the Mullins effect that is overestimated for the first cycles. Moreover, the shape and the size of stress-strain loops would be reproduced in a more precise manner.

CONCLUDING REMARKS

The different phenomena that constitute the behaviour of elastomers are classically studied separately. In order to model the global behaviour of these materials, it is necessary to combine these separated approaches. A major problem induced by the addition of models is the large number of material parameters to be determined, so that these parameters can not be simply calculated. That is the reason why the choice of efficient models with few numbers of parameters for each phenomenon is highly important. Moreover, in order to be assemble these sub-models into a single global constitutive equations, the approaches adopted must be compatible.

The present study is based on a physical approach that focus on the chains behaviour under loading conditions. The hyperelastic basis of our model is the eight-chains constitutive equation. The Mullins effect modelling considers the breakage of secondary links between elastomeric chains. And finally, the hysteretic part of the behaviour is based on the slippage of chains. Consequently, all these approaches can be mixed together in order to define a global model.

This final model successfully describes the hysteresis and Mullins effect of elastomers. It can be improved by adding the other characteristic phenomena observed in elastomeric materials, the most important being the long-term viscoelasticity.

Figure 1 .

 1 Figure 1. Spatial representation of the eight chains model.

  Figure 3. Hystereric behaviour

  τ and Bp chain λ are respectively the stress and the chain length of the rate-dependent part of the model, and 1 C , 2

Figure 5 .

 5 Figure 5. Final constitutive equation including large strains, hysteresis and Mullins effect.

Figure 7 .

 7 Figure 6. Tensile test simulation using the global model.


M Aït-Bachir

E Verron

Energy release rate of small cracks under finite multiaxial straining

The energy release rate of small cracks governs fatigue crack nucleation. A method is presented here to efficiently and accurately evaluate the energy release rate of such cracks, arbitrarily oriented, under general conditions of finite multiaxial loading. As a motivation, the dependence on crack length is then investigated. It is demonstrated that the energy release rate of small cracks is proportional to the crack length and that the proportionality factor is a function of the far-field parameters only. An attempt is then made to search for a general expression of this proportionality factor under simple loading conditions.

INTRODUCTION

Studies carried out in the recent years have brought to light the physical phenomena governing fatigue life of rubber [START_REF] Cam | Mechanism of fatigue crack growth in carbon black filled natural rubber[END_REF][START_REF] Gorju | Fatigue life of rubber components: 3d damage evolution from x-ray computed microtomography. In Constitutive Models for Rubber V. Mars[END_REF]. It turns out to be driven by the growth of small cavities transforming then into small cracks and propagating throughout the material up to a critical size that leads to a major loss of material properties. The propagation of these small cracks represents the main stage of fatigue life.

Regarding small crack growth prediction two distinct approaches are usually considered. The crack initiation approach is based on the evaluation of the mechanical fields of a crack-free material in order to study how a small flaw would propagate when subjected to these mechanical conditions; this approach lies thus within the framework of continuum mechanics. Some predictors have already been developed to predict crack initiation (Mars 2002;[START_REF] Verron | Definition of new predictor for multiaxial fatigue crack nucleation in rubber[END_REF]. To the contrary, the crack propagation approach studies the propagation of an existing small crack embedded in the material [START_REF] Gent | Cut growth and fatigue of rubbers. i-the relationship between cut growth and fatigue[END_REF][START_REF] Lake | The mechanical fatigue limit for rubber[END_REF].

Actually, in some simple cases it is possible to reconcile these two approaches. Indeed, considering a small crack of length c under plane stress uniaxial tension, Rivlin and Thomas have been able to factorize the energy release rate T [START_REF] Rivlin | Rupture of rubber[END_REF]:

T = 2kWc (1)
where W stands for the strain energy density of the far-field region and k is a factor that depends only on the far-field loading conditions. This result has later been extended to both pure shear and equibiaxial tension [START_REF] Yeoh | Relation between crack surface displacements and strain energy release rate in thin rubber sheets[END_REF]. Thus, it turns out that for a small crack subjected to simple loading conditions, the energy release rate is proportional to the crack length and the proportionality factor is a function of the far-field parameters only.

In this paper we aim at developping a method to easily estimate the energy relase rate of small cracks under arbitrary plane stress loading conditions and then at investigating the factorization of the energy release rate with respect to the crack length.

METHODS

A simple tool to model a small crack under arbitrary loading conditions

A method is needed first to model a small crack, arbitrarily oriented, under general conditions of far-field multiaxial loading.

Present definition of a small crack

In order to define the concept of small crack we introduce first the notion of "boundary" as a circular region drawn around the crack, on which we compare solutions for the "body" and the "crack neighbourhood" (see Figure 1):

• "body" refers to the solution obtained for a deformed body without considering the effects of a crack. As we move away form the edge, the gradients induced by edge effects get smaller and imply an upper limit on the size of the "boundary". • "crack neighbourhood" refers to the solution obtained for a crack embedded in an infinite medium, with specified far-field loading conditions. As we move away from the crack, the gradients induced by the crack get smaller and imply a lower limit on the size of the "boundary".

A crack is then said to be small when one can draw a boundary around it such that the solution on the boundary, in both the "body" and in the "crack neighbourhood", is constant to within a given tolerance.

Note here that this definition is very general because edge effect is considered. However, in the present study, the body is made up of an isotropic hyperelastic material and is subjected to multiaxial homogeneous loading conditions; thus gradients of the mechanical fields in the body are nul. The body-crack neighbourhood boundary is then chosen at a location where gradients tend to zero.

Modelling boundary displacements under plane stress multiaxial straining

The small crack assumption permits us to impose the far-field state of stress directly via the displacement at the boundary. Indeed, except in a very small region around the crack, the stress and strain fields are homogeneous and are the same as in the crackfree model for which the analytical expressions of strain and stress fields along with the strain energy density are completely known.

The original model consists of a small throughcrack with orientation θ (in the undeformed configuration) under a far-field equibiaxial straining (λ 1 , λ 2 ) in the frame of reference ( , ) e e 1 2

(see Figure 2). Throughout the deformation the crack rotates and ends up with an orientation ψ with respect to its orientation in the undeformed configuration. The problem with this method is that a new geometry is required whenever the crack orientation θ is to be changed. However, from the perspective of the crack, changing the crack orientation boils down to changing the far-field loading (or equivalently the far-field straining) (see Figure 3). Thus, by performing first a change of basis and then a pull-back to maintain the crack orientation fixed throughout the deformation, we can express the associated far-field strain in the crack-based frame and capture the same range of conditions with a single FE model. Moreover, no restriction is made on the geometry of the far-field boundary which requires working out a very general expression of the displacement to be imposed to a node at the boundary. In order to conveniently study the effect of biaxiality, far-field loading parameters (λ 1 , λ 2 ) are replaced by (λ, B) with λ 2 = λ and λ 1 = λ B . Thus B = log λ 2 /log λ 1 quantifies the biaxiality and can be called the biaxiality factor (Mars 2002). The relationship between the undeformed coordinates (X, Y) and the deformed coordinates (X′, Y′) of a point P in the crack-based frame ( , )

′ ′ e e 1 2 (see Figure 2) under the far-field loading conditions described above is given by:

′ ′ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ X Y R R R X Y T T B ( ) ( ) ( ) ψ θ λ λ θ 0 0 (2)
where

R( ) cos sin sin cos θ θ θ θ θ = - ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ (3) and R( ) cos sin sin cos ψ ψ ψ ψ ψ = - ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ (4)
It should be pointed out that the implementation of the previous calculation (Eq. 2) and further computations carried out in this study were all obtained with the finite element software Abaqus in which the displacement of a node at the far-field boundary was imposed via a DISP subroutine.

Basic validations

First we have validated the calculations described in 2.1.2 by verifying that a crack-free model gives homogeneous loading for all B, λ and θ.

We then focused our attention on the "small crack" requirement. Indeed, far from the crack, it is desired to get a homogeneous stress field. This condition can never be met rigorously because of the finite dimension of the model and the presence of the crack. However, on using the method presented in section 2.1.1 a criterion was set to check accuracy of the approximation. The latter was considered acceptable when the relative error on the far-field stress when compared to the crack free model was less than 1%. From the experience a crack length 40 times smaller than the dimension of the model fits well into this criterion.

Note here that an additional satisfactory way to validate the calculations is to compare the energy release rate computed by Abaqus (J-integral) with the fracture mechanics solutions and from some studies carried out in finite strain [START_REF] Yeoh | Relation between crack surface displacements and strain energy release rate in thin rubber sheets[END_REF].

Energy release rates of small cracks and their factorization

It is first desired to develop a method that enables us to capture the influence of the presence of a small crack on the variation of the mechanical fields in the crack neighbourhood and that reveals how that variation makes up for the energy release rate. Using this method, the factorization of the latter with respect to the crack length is then investigated. Finally, on the basis of the numerical results obtained via finite element analysis, we examine the dependence of the energy release rate of small cracks on the farfield parameters B, λ and θ.

J-integral

The J-integral represents a way to calculate the energy release rate. It was first introduced within the framework of planar small strain as a contour path integral around the crack tip [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF]). Thereafter, it was extended to planar finite strain. Indeed, on considering a contour Γ surrounding the crack tip and leaning on both faces of the crack in the undeformed configuration, the J-integral writes:

J q NdS = ⋅ ∑ ∫ Γ (5)
where ∑ stands for the Eshelby stress tensor, q → is the crack direction vector and N → is the outward normal vector to the surface element dS in the under formed configuration.

We remind the reader that, for an elastic material, the J-integral is path-independent which permits the arbitrary choice of the contour surrounding the crack tip to calculate the energy release rate [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF]. For the sake of simplicity, the far-field Eshelby stress tensor is denoted ∑ ∞ throughout the rest of this paper.

Let us focus now our attention on the contour. A rectangular contour of characteristic dimension R has been chosen (see Figure 4). On choosing R such that R/c → ∞ the evaluation gets highly simplified (see Table 1). Indeed, Γ A , Γ B and Γ C lie then in the far-field region wherein the Eshelby stress tensor is uniform and is equal to ∑ ∞ .

Figure 3. Original and new models after performing the appropriate change of basis (undeformed configuration).

It follows that

J l d l R R R = - ( ) → + ∞ ∞ -∑ ∑ ∫ lim ( ) 11 11 (6)
Using the symmetry with respect to the crack line, we deduce that:

J ld l = - ( ) ∞ + ∞ ∑ ∑ ∫ 2 11 11 0 ( ) (7)
Note here that the expression above is not an approximation but the exact formula of the energy release rate. However, because the model has finite dimension, the upper bound for the integral is necessarily < + ∞. Thus, we can only compute an approximation of the energy release rate:

J ld l L ≅ - ( ) ∞ ∑ ∑ ∫ 2 11 11 0 ( ) ( 8
)
where L is the maximum distance to the crack face in the ′ e 2 -direction.

Proportionality of J with respect to the crack length

Let us consider two different small cracks: one with length cand the other one with length kc where k is a real number. Both are subjected to the same farfield loading conditions. Using Eq. ( 7), the energy release rate of a crack of length cwrites:

J c l dl c ( ) ( ) = - ( ) ∑ ∑ ∫ ∞ + ∞ 2 11 11 0 (9)
while the energy release rate of a crack of length kc writes:

J kc l dl kc ( ) ( ) = - ( ) ∑ ∑ ∫ ∞ + ∞ 2 11 11 0 (10)
where ∑ 11 c (respectively ∑ 11 kc ) denotes the Eshelby stress associated with the crack of length c (respectively of length kc).

By using the substitution l = kl′, we obtain that dl = kdl′. It follows that:

J kc k kl dl kc ( ) ( ) = - ′ ( ) ′ ∑ ∑ ∫ ∞ + ∞ 2 11 11 0 (11)
Let us focus now on the effect of multiplying the crack length by a factor k on the transformation of the mechanical fields with respect to the original problem of a crack of length c. We recall that the latter consists in a small crack embedded in a medium that has infinite dimension; thus changing the crack length by a factor k boils down to performing a homothetic transformation with a scale factor k and where the homothetic center is the center of the crack. And as the segment Γ + goes by the center of the crack, the mechanical fields along its length obey the homothetic transformation with respect to the crack length. Hence we have:

( ) ( ) kl l c kc = ∑ ∑ 11 on + 11 Γ (12)
and we deduce that:

J kc k kl dl kc ( ) ( ) = - ( ) ∑ ∑ ∫ ∞ + ∞ 2 11 11 0 (13) = - ( ) ∑ ∑ ∫ ∞ + ∞ 2 11 11 0 k k l d l c ( ) (14) = kJ c ( ) (15)
Note here that no limits were placed on the particular plane stress loading condition. Thus, the argument we have just made is completely general and is valid for all λ, B and θ.

This result implies that it is sufficient to work out the energy release rates for one crack size, and Hence we have:

N q Nd ⋅ ∑ ∫ Γ Γ Γ + -′ e 1 -∑ ∫ 11 0 ( ) l dl R Γ A ′ e 2 ∑ ∞ 12 R Γ B ′ e 1 2 11 ∑ ∞ R Γ C -′ e 2 -∑ ∞ 12 R Γ - -′ e 1 - ∑ - ∫ 11 0 ( ) l dl R
that the energy release rate of all other crack sizes are then determined.

NUMERICAL RESULTS

The behaviour of the isotropic hyperelastic material was modeled by an incompressible neo-Hookean model.

Variation of the energy release rate versus crack length

In order to validate the proportionality law demon strated previously (see Eq. ( 15)), the dependence of the energy release rate on the crack length has been investigated for various far-field loading conditions (B ranging from -0.5 to 1, λ ranging from 1 to 5 and θ ranging from 0° to 90°). Indeed, for each and every far-field state of stress, the energy release rate has been computed twice for two different crack lengths: c and 2c. All the results obtained so far show that, for Bincreasing from -0.5 to 1 with an increment of 0.25, for λ increasing from 1.1 to 5.0 with an increment of 0.1 and θ increasing from 0° to 90° with an increment of 15°, we always verify:

J B c J B c ( , , , ) ( , , , ) λ θ λ θ 2 2 = (16)

Variation of the energy release rate versus crack orientation

For various biaxiality factors (B ranging from -0.5 to 1), the influence of λ on the variation of the energy release rate versus the crack orientation has been investigated. As we can observe (see Figure 5 and Figure 6), for B in [-0.5, 1] the energy release rate strictly decreases as the crack orientation increases in both small and Figure 5. Energy release rate under simple loading conditions (uniaxial tension ---, pure shear ---and equibiaxial tension -) in small strain (λ = 1.1). finite strain: the maximum is reached at θ = 0° and the minimum is reached at θ = 90°. However, under equibiaxial tension (B = 1), the energy release rate is independent from the crack orientation in both small and finite strain and thus only depends on λ.

DISCUSSION: FAR-FIELD PARAMETER OF THE SCALE LAW

The proportionality of the energy release rate with respect to the crack length (section 2.3.2) permits the following factorization:

J (λ, B, θ, c) = f (λ, B, θ) c ( 17
)
where f is a function of the far-field parameters only.

Comparison with linear elastic fracture mechanics (LEFM)

From the solution of LEFM for a rubberlike material [START_REF] Mars | Heuristic approach for approximating energy release rates of small cracks under finite strain, multiaxial loading[END_REF]), f can be factorized into:

f B f B ( , , ) ( , ) (cos sin ) λ θ λ θ β θ = + 0 2 2 2 (18)
where

β = + + 2 1 2 B B ( 19
)
and f 0 is a function of λ and B only. f 0 (λ, B) is basically the value of f(λ, B, θ) at θ = 0 °. Thus, regarding the original problem (see Figure 2), f 0 depends only on the far-field state of stress.

Note here that all the energy release rates computed in small strain for B ranging from -0.5 to 1 and θ ranging from 0° to 90° completely match this theoretical factorization.

Factorization of the energy release rate in finite strain under simple loading cases

Under equibiaxial tension, because the energy release rate is independent from the crack orientation, the previous factorization of f remains obviously valid. Under uniaxial tension, we can draw the same conclusion from the numerical results. Indeed, on comparing the factorization with the energy release rate values computed for λ ranging from 1.1 to 5, we always have:

f (λ, -0.5, θ) = f (λ, -0.5, 0°) cos 2 θ (20)

However, under pure shear and actually for all B in [-0.5, 1.0], a simple comparison with numerical results clearly shows that the previous factorization (Eq. 18) cannot be extended to finite strain.

As for now, all the arguments that we made in this subsection are based only on finite element results but they have the benefit to guide us towards the simplification of the expression of f under simple loading cases. From the remarks above we deduce that: J (λ, -0.5, θ, c) = (f 0 (λ, -0.5) cos 2 θ) c (21)

and

J(λ, 1, θ, c) = f 0 (λ, 1) c (22)
wherein the general expression of f 0 still needs to be determined.

CONCLUSIONS

A new method for modelling small cracks under arbitrary loading states and finite straining has been presented. The latter was successfully tested against the most known cases (uniaxial tension, equibiaxial tension and pure shear).

We have also investigated how the energy release rate is balanced by the distribution of configurational stresses. It has then been proved that the energy release rate of a small crack (suitably defined) always follows a linear scale law with respect to crack size, regardless of loading state.

Finally some progress has been made towards a general-purpose expression for energy release rate under arbitrary loading.

Figure 1 .

 1 Figure 1. Body (whitish grey), crack neighbourhood (dark grey) and boundary (grey).

Figure 2 .

 2 Figure 2. Small through-crack under finite biaxial straining with arbitrary orientation.

Figure 4 .

 4 Figure 4. Contour for the evaluation of the J-integral.

Figure 6 .

 6 Figure6. Energy release rate under simple loading conditions (uniaxial tension ---, pure shear ---and equibiaxial tension -) in large strain (λ = 2.5).

Table 1 .

 1 Contributions to the energy release rate.

	Segment