Marco Gigliotti
email: gigliotti@emse.fr

Frédéric Jacquemin

Alain Vautrin

Internal Stresses Induced by Cyclical Environmental Conditions in Laminated Composite Plates

Internal stresses induced by cyclical environmental conditions in laminated composite plates

INTRODUCTION

Environmental conditions induce internal stresses in composite materials. At the mesoscopic scale, these stresses are mainly induced by the mismatch in the coefficients of hygroscopic expansion between adjacent plies and are usually calculated assuming uniform and constant hygrothermal fields.

Real hygrothermal conditions are of variable amplitude and, due to the diffusion processes, produce internal distributions which are far from being uniform. There are many applications, for example aeronautical structures subject to real flight conditions, in which the external hygrothermal conditions are approximately cyclical, the period of a cycle being the time of a flight. These conditions produce transient internal concentrations and transient internal stresses which deserve to be detailed [START_REF]Aeronautical Technologies for the Twenty-First Century[END_REF].

Calculation of internal stresses due to transient conditions has been performed for cylindrical and plate structures [START_REF] Benkeddad | On the Transient Hygroscopic Stresses in Laminated Composite Plates[END_REF][START_REF] Paul | Transient Hygrothermal Stresses in Laminated Cylinders[END_REF], the problem of a composite pipe subjected to cyclical hygrothermal conditions has been also solved [START_REF] Jacquemin | A Closed-form Solution for the Internal Stresses in Thick Composite Cylinders Induced by Cyclical Environmental Conditions[END_REF].

In the present paper both analytical and difference finite approaches are provided to determine the hygrothermal concentration fields inside a composite plate subject to cyclical conditions. Stresses are then calculated by using the classical lamination theory in a formulation adapted for taking into account non-uniform moisture concentrations [START_REF] Benkeddad | On the Transient Hygroscopic Stresses in Laminated Composite Plates[END_REF]. Residual curing stresses are also taken into account by considering the thermal differential between cure and room conditions. The uncoupled hygrothermoelastic model is presented in section 2.

The model is then used to simulate environmental induced stress in supersonic flight. The cycle consists of a mixed cycling of external temperature and moisture, conditioning due to maintenance is also taken into account. Discussion of results with emphasis on the complexity of the resulting state of stress is provided, some qualitative consideration on the effects of such stresses on the resulting damage, which has been experimentally observed [START_REF] Weitsman | Moisture in Composites: Sorption and Damage[END_REF], is also given.

MODEL APPROACH

The thermal and hygroscopic fields are in some sense coupled, the coefficient of moisture diffusivity depends on the temperature by an Arrhenius law, as it will be detailed later. Once the hygrothermal field is obtained, strains and stresses are calculated by assuming hygrothermoelastic behaviour. Attention is paid to plate structures which are sufficiently thin to enable us to formulate the following hypotheses:

-the thermal field is uniform along the entire structure. This hypothesis relies also on the fact that the time for thermal diffusion is much smaller than the time for moisture diffusion, -moisture diffusion is unidirectional along the thickness of the plate and follows the Fick's law [START_REF] Crank | The Mathematics of Diffusion" 2 nd Edition[END_REF], -stresses and strains are calculated according to the classical lamination theory [START_REF] Herakovitch | Mechanics of Fibrous Composites[END_REF].

HYGROTHERMAL FIELD

We consider the problem of a composite plate in a Euclidean reference frame (x,y,z) subjected to cyclical hygrothermal conditions at its boundaries. By cyclical we mean that they vary in a periodic manner as a function of time. The plate can be considered as infinite along two directions, x,y for instance, so that the hygrothermal diffusion is essentially unidirectional along the z direction. As mentioned introducing the model, the thermal equilibrium is supposed instantaneous, the temperature T inside the plate is uniform, even if it may vary with time T=T(t). The moisture diffusion coefficient D depends on the temperature through the Arrhenius law:

        = ) ( exp ) ( t T B A t D (1)
where A and B are material constants. The moisture concentration c, which is the density of solute absorbed by the composite plate, is obtained by solving the following field equation:

0 , 0 ), , ( ) ( ), , ( > < < ∀ = t e z t z c t D t z c zz t (2)    = = = ∞ ) ( ) 0 , ( ) ( ) , ( ) , 0 ( z c z c t c t e c t c i (3)
where e is the thickness of the plate, ) (z c i is the initial concentration and ) (t c ∞ is the concentration at the lateral surfaces, which depends on the external relative humidity HR(t) by the law:

b t HR C t c ) ( ) ( ⋅ = ∞ (4)
C and b are constants. All the above functions have the same period τ.

The problem is relatively complex, as boundary conditions and diffusion coefficient depend explicitly on the time. As already noted in [START_REF] Carslaw | Conduction of Heat in Solids" 2 nd Edition[END_REF], a diffusion/conduction problem along a semi-infinite region with periodic boundary conditions is characterised by a typical dimension, indicated by λ, which measure the distance from the surface where a fluctuating regime exists. λ is given by:

τ π λ k 2 = ( 5
)
where k is a conduction/diffusion coefficient and τ is the period. λ can be exactly defined as the distance, taken from the surface, at which the amplitude of all oscillations is reduced by a factor of 0019 . 0 2 = -π e

. By analogy with this problem, a length e 0 can be defined as (see also [START_REF] Jacquemin | A Closed-form Solution for the Internal Stresses in Thick Composite Cylinders Induced by Cyclical Environmental Conditions[END_REF][START_REF] Gigliotti | Internal Stresses in Composite Laminates due to Cyclical Hygrothermal Loading. Application to Supersonic Flight[END_REF]):

∫ = τ π 0 0 ) ( 2 dt t D e (6)
which has the same meaning as λ, but it takes into account the effect of a varying diffusion coefficient. Some important considerations need to be addressed concerning e 0 :

-e 0 practically measures the distance from the surfaces in which hygrothermal oscillations exist. During the absorption phase, internal zones outside e 0 are in a transient regime. When "saturation" is reached, spatial regions inside e 0 keep fluctuating with period τ, while spatial regions outside e 0 reach a permanent (thus uniform) regime, -e 0 depends on the diffusion coefficient and the period of the cycle. An important conclusion is that a plate with 0 2e e < is in a fluctuating regime all along its thickness. An analytical solution of equation ( 2) with boundary conditions (3) can be obtained by employing a change of variable, in order to eliminate the explicit time dependence, and a recursive technique, in order to express the concentration as a function of the number of cycles. Details of such technique can be found in reference [START_REF] Jacquemin | A Closed-form Solution for the Internal Stresses in Thick Composite Cylinders Induced by Cyclical Environmental Conditions[END_REF] for laminated cylinders and in reference [START_REF] Gigliotti | Internal Stresses in Composite Laminates due to Cyclical Hygrothermal Loading. Application to Supersonic Flight[END_REF] for laminated plates. In the following equations the solution for a laminated plate is presented, by assuming that the diffusion coefficient and the moisture concentration are not discontinuous at the interfaces between adjacent plies. In other terms, concerning the diffusion problem, the plate behaves as if it was homogeneous. By indicating as N the number of cycles at which the internal regions of the plate reach a permanent state, the concentration field after (N-k) cycles is finally given by:

( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) 2 2 ' 2 2 ' 1 2 cos ' 1 2 1 2 2 1 2 cos 2 4 1 2 exp 1 2 2 2 ) ( ) , (
i ≤ ≤ - ∀                                 + + + -               +                                   - ∆ + - +       ∆ + = - ∫ ∑∑ ∞ + ∞ = - = ∞ π π π τ π π τ (7)
where:

( ) ( ) ( ) ( ) ( ) ∫ ∫ ∞ ∞ ∆ = = ∆ τ τ τ τ 0 0 1 , dt t c t D c dt t D (8) When ( ) 0 = z c i
, equation [START_REF] Herakovitch | Mechanics of Fibrous Composites[END_REF] becomes:

( ) ( ) ( ) ( ) ( ) ( ) ∑∑ ∞ = - = ∞ ∞         +                 - + - + - - = - 0 1 0 2 2 2 2 1 2 cos 1 2 exp 1 2 1 4 , n k i n e z n i N n n c c k N z c π β π β π 2 2 e z e ≤ ≤ - ∀ (9)
where:

( ) τ β ∆ = 2 e
(10) Equations ( 7) and ( 9) employ average concentration values at the boundaries (eq.8), thus they describe correctly the solution in the interior regions of the plate only, that is, outside e 0 . The finite difference method is employed besides the analytical one, to better represent the concentration distribution in regions close to the lateral surfaces.

STRESS FIELD

The stress field is obtained by employing the classical lamination theory, adapted for taking into account transient moisture concentration fields, such as those given by equations ( 7) and ( 9). In the following paragraph we present the main features of the model, see reference [START_REF] Gigliotti | Internal Stresses in Composite Laminates due to Cyclical Hygrothermal Loading. Application to Supersonic Flight[END_REF] for further details.

The main assumption of the model is to approximate the hygrothermal fields given by equations ( 7) and (9) by equivalent fields, which are linear over spatial sub-regions of a physical ply, called sub-plies. This is done so that the total strains, including the free expansion strains of hygrothermal nature, satisfy automatically the compatibility equations which otherwise are not satisfied, at least in the context of a plate theory. The employment of sub-plies has some shortcomings. As the number of sub-plies is finite, the concentration and the stress fields are approximate: the accuracy increases as the number of sub-plies increases, but the maximum number of sub-plies itself is limited by some physical consideration. In the present communication a macroscopic approach is proposed and, physically, a sub-ply must be able to represent the homogenised behaviour of the material, this imposes a limitation on its minimum dimensions.

EXAMPLE : APPLICATION TO SUPERSONIC FLIGHT

In order to illustrate the technique presented in section 2, we consider the hygrothermal cycle illustrated in figure 1. This cycle simulate the environmental conditioning encountered by skin aircraft structures during a flight. In the cycle some peculiar zones can be distinguished:

-point A, at room temperature and ground humidity conditions, which simulate ground environmental conditions -point B, at -55°C and 0% HR, which simulate the subsonic phase of the flight -point C, at 130°C and 0% HR, which simulate the supersonic phase of the flight point D, with the same conditions as in point A, which is located at the end of a flight An initial conditioning at 23°C and 50% HR for about 3 months is included in the calculations, in order to simulate a long term ground maintenance, which leads to a pseudohumid state.

A 4 mm thick [0 4 /90 4 ] s composite plate is considered, each ply "block" being 1 mm thick. Inside each ply block the spatial domain is divided into sub-domains 0.1 mm thick, such subdivision is sufficient to give a good accuracy of the calculated concentration fields. Figure 2 presents the concentration field at each point of cycle 1: point A for such cycle corresponds to the pseudo-humid state. For cycle 1 only the finite difference solution is employed.

Figure 3 presents the concentration field at each point of cycle 50. For such cycle, the analytical solution is fast to be estimated, as few terms in the series [START_REF] Herakovitch | Mechanics of Fibrous Composites[END_REF] are required, and such a solution is presented in figure 3 for comparison with the finite difference one. In both figures the symbol m(%) is employed, which indicates the percentage of solute density divided by the density of the reference dry solvent material.

Figures 2 and3 show that hygrothermal fluctuations are close to the lateral surfaces and extend over spatial regions of value e 0 , in this case equal to around 0.5 mm. The average effect of the supersonic cycles consists of drying the plate until a quasi-dry state, which is reached after around 300 cycles. The analytical solution, as already mentioned in section 2, respects average boundary conditions and is adequate for the internal regions only.

The initial mechanical state, before any conditioning, is not free of stress, residual stresses due to manufacturing are present in the structure. The stress free state is supposed to be achieved at a temperature of 210°C, which is around the final T g of a high-temperature composite system, the temperature differential is then given by ∆T= -187°C. Mechanical properties used for simulations are given in table 1, they are representative of a IM7/977-2 carbon fibre-epoxy ply and are given in the usual orthotropic material axes.

E 1 (GPa) E 2 (GPa) G 12 (GPa) ν 12 α 1 (ε o C -1 ) α 2 (ε o C -1 )

DISCUSSION

It can be noted from simulations that stresses are quite high after manufacturing. Residual curing stresses, which are presented as a reference in all figures, are uniform through plies, by the assumption of uniform ∆T, and take a value of around 42 MPa.

The internal stresses which arise after a long period of maintenance due to the absorption of moisture, point A in cycle 1, counteract only partially the residual curing stress and have an average value of 35 MPa, which is quite high when compared to the reference dry state.

Thermal fluctuations from -55°C to 130°C, points B and C respectively, extend over the whole thickness of the plate and exhibit an average amplitude of around 40 MPa, the mean stress is about 30 MPa. Low temperature conditions are quite critical, transverse stresses can reach, over the cycling, values as high as 60 MPa close to the lateral surfaces.

Hygroscopic fluctuations induce stress gradients over a region close to the lateral surfaces and of thickness e 0 and fluctuating stresses whose magnitude can be as high as 40 MPa.

Drying induced by supersonic cycles promotes a progressive resurgence of the residual curing stresses, the mean of all fluctuations increase with increasing cycles.

The calculated stresses, which values are quite high especially at low temperatures, may potentially damage the structure. In particular, the mixed hygrothermal cycling close to the lateral surface, whose extent is around 0.5mm (two or more physical plies), could be particularly aggressive.

CONCLUSIONS

A methodology for simulating hygrothermal fields and internal stresses induced by cyclical environmental conditions in laminated composite plates is presented. This methodology assumes uniform temperature distribution over the plate, employs the Fick's law for calculating the moisture concentration profiles and the classical lamination theory, adapted for transient conditions, for simulating the induced stresses.

Simulations are provided for a plate structure subjected to almost realistic supersonic flight conditions and give values of fluctuating stresses which can potentially induce damage in the structure without any external load.

The model is quite simple, based on straight and plausible hypotheses and is supported by analytical calculations. It can be adapted to constitutive laws different from the elastic ones, it can be generalised for taking into account the effect of damage.

The model is also intended as a tool for preliminary design and for conceiving experimental tests.

Fig. 1 :

 1 Fig. 1: Hygrothermal cycle

Figure 2 :

 2 Figure 2: Moisture concentration during cycle 1

Figure 3 :

 3 Figure 3: Moisture concentration during cycle 50

Figure 4 : Internal stresses xx σ during cycle 1 Figures 4

 414 Figure 4: Internal stresses xx σ during cycle 1 Figures 4 and 5 illustrate stresses xx σ and yy σ at points A, B, C and D of cycle 1. Figures 6and 7show the evolution of internal stresses with increasing number of cycles.

Figure 5 :

 5 Figure 5: Internal stresses yy σ during cycle 1

Figure 6 :

 6 Figure 6: Internal stresses xx σ at the start of cycles 1, 2, 50, 100, 300

Table 2 :

 2 Material properties used in the simulations (IM7/977-2)

	β 1 β 2

ACKWNOLEDGEMENTS

The present study was supported by the French Ministry of Research as part of the French Supersonic research program. All partners, ONERA, CCR EADS and Ecole Centrale Lyon are gratefully acknowledged.


Abbas Krayani

Frédéric Dufour

Gilles Pijaudier-Cabot

Elastoplastic nonlocal damage model for concrete and size effect analysis

Elastoplastic nonlocal damage model for concrete and size effect analysis A. Krayani, F. Dufour & G. Pijaudier-Cabot ERT R&DO, GeM, Ecole Centrale de Nantes, France ABSTRACT: This work deals with a combination of plasticity and nonlocal damage formulation for modelling concrete structure behaviour subjected to mixed mode failure. Plastic effect, driven by effective stress, is an hardening process and accounts for the development of irreversible strains while softening is controlled by damage to describe the degradation of the material stiffness. A regularization technique, based on the implicit gradient definition of the nonlocal strain tensor, is introduced to overcome the deficiencies induced by the softening damage relation i.e., spurious strain localization and dependance of the energy dissipation on mesh refinement. A 3D tensile bar localization benchmark is used to validate the effectiveness of the regularization technique. Moreover, simulation of reinforced concrete continuous deep beam is studied to illustrate the improvements achieved by coupled law compared to scalar damage models. Finally, size effect in mixed mode failure is investigated by mean of three asymmetric four-point bending tests on notched specimens.

INTRODUCTION

One of the specific behaviours of cementitious materials is that in uniaxial loading the traction strength is much smaller than compression strength. From that consideration, [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF] has developed a damage model which takes only the positive principal strains into account. This criterion is widely used for its simplicity and relative accuracy to represent mode I failures of concrete structures under monotone loading especially.

For some specific structures (thick and/or prestressed) the stress state may be locally bi-or tri-axial. As for any geomaterials the confinement has a positive effect on the structural strength and depending on the level of confinement the failure mode may change to mode II. For this reason, the modified von Mises criterion [START_REF] De Vree | Comparison of nonlocal approaches in continuum damage mechanics[END_REF] seems to be more attracive as the second invariant of the stress tensor is introduced into the yield criterion. However, this criterion is symetric (traction-compression) and does not represent the reality.

The idea followed by [START_REF] Jason | An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model[END_REF] is to couple plasticity, based on the effective stress, and local damage formulation in one single constitutive relation in order to represent a correct volumetric response of concrete, i.e. the confinement effect and the non symetric response of the material. The damage process, controlled by elastic strain , is described by the isotropic model developed by [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF]. The plastic process is described by means of a yield sur-face, inspired from [START_REF] Etse | Fracture energy formulation for inelastic behavior of plain concrete[END_REF] and modified by [START_REF] Crouch | Application of a stress return algorithm for elasto-plastic hardening-softening models with high yeild surface curvature[END_REF]. In this approach, damage controls softening, while plasticity controls harding, in compression especially.

This model generates ill-posed mathematical problems due to the loss of ellipticity of the governing differential equations. Numerically, the results suffer from pathological sensitivity to the size and the orientation of the finite element mesh and the total energy dissipated by the fracture process tends to zero upon mesh refinement. Several nonlocal constitutive laws, which can be regarded as a remedy of the pathalogical mesh dependence, can be found in the literature. Whether they are in an integral form [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] or in a gradient form [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], a salient characteristic of both types is the presence of a characteristic length in the constitutive relation.

The aim of this paper is to extend the previous approach to a non local model in order to simulate properly concrete structure behaviour with mixed mode failure. The regularization technique, based on the implicit gradient definitions of the nonlocal strain tensor, is adopted to overcome the mathematical problems due to the softening constitutive law (damage part), while the plastic part remains local.

In this contribution, the constitutive law of the model is first presented. Then, the validation and the performance of the model is evaluated in some structural cases.

MODEL FORMULATION

The effective stress approach has been chosen for coupling plasticity and damage effects. The total stress σ is then a function of the damage variable D and of an effective stress σ ′ that solves the equation of a plastic yield surface

σ = (1 -D)σ ′ (1)
2.1 Plasticity Elastoplastic constitutive law is defined by the following set of equations in which the effective stress σ ′ has been substituted to the applied stress σ :

ε σ ′ εp κ = ε e + ε p = Eε e = λm (σ ′ , κ) = λh (σ ′ , κ) (2)
where ε, ε e and ε p are the total, elastic and plastic strains respectively, E is the elastic stiffness tensor, m is the flow vector, κ is the set of internal variables and h is the plastic modulus. The plastic multiplier λ is given by the loading-unloading criterion (Kuhn-Tucker form)

F (σ ′ , κ) ≤ 0, λ ≥ 0, F (σ ′ , κ) λ = 0, (3)
where F (σ ′ , κ) is the yield function that defines the effective stress state. The plastic flow vector m is defined as (associated law)

m = ∂F ∂σ ′ (4)
The plastic yield surface, inspired from [START_REF] Etse | Fracture energy formulation for inelastic behavior of plain concrete[END_REF] and modified by [START_REF] Crouch | Application of a stress return algorithm for elasto-plastic hardening-softening models with high yeild surface curvature[END_REF], depends on four main functions ρ (effective stress invariant), k (hardening function), ρc (deviatoric parameter) and r (deviatoric shape function) (for more details, see [START_REF] Jason | An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model[END_REF]))

F = ρ2 (σ ′ ) - k (σ ′ , k h ) ρc (σ ′ ) r 2 (σ ′ ) (5)
where k h is the hardening variable given as a function of confinement:

kh = √ 2 3 ε p :ε p ζ(σ ′ ) if k h < 1, kh = 0 if k h = 1, (6)
where ζ depends on the first normalized invariant

ζ = -A h + A 2 h -B h ξ + C h if ξ ≤ 0, ζ = -A h + A 2 h + C h if ξ > 0, (7)
A h ,B h ,C h are three model parameters. Note that equation ( 7) assumes that k h ranges between 0 and 1. When k h = 1, the yield surface becomes a limit surface with no hardening in pure elasto-plastic model. Combining equations ( 2) and ( 7) yields the following expression for the plastic modulus h:

h = √ 2 3 ∂F ∂σ ′ : ∂F ∂σ ′ ζ(σ ′ ) if k h ≤ 1, h = 0 if k h = 1, (8)
The plastic part of the constitutive relation contains 10 parameters in which four parameters remains at a fixed value: α = 0.5, γ = 0.99, A h = 7 × 10 -5 and k 0 = 0.1 (for more details, see [START_REF] Jason | An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model[END_REF]).

Only two parameters r c and r t will be given in the applications where the others are taken as p = 0.4,

B h = 2 × 10 -2 , C h = 2 × 10 -6 and A = -0.5.

Damage model

Damage growth, controlled by elastic strain tensor, is governed by the following loading function:

g(ε e , k d ) = ε eq (ε e ) -k d (9)
where ε eq is the equivalent strain. k d initially equals the damage threshold ε D 0 and during the damage process equals the largest ever reached value of ε eq . The evolution of damage is governed by the standard loading/unloading condition:

g(ε e , k d ) ≤ 0, kd ≥ 0, kd g(ε e , k d ) = 0 (10)
For quasi-brittle materials such as concrete, the damage propagates much more easily under tension than under compression. Taking that into account, one possible choice is the so-called modified von Mises definition according to de Vree et al. ( 1995)

ε eq = k -1 2k (1 -2ν) I 1 (11) + 1 2k (k -1) 2 (1 -2ν) 2 I 2 1 + 6k (1 + ν) 2 J 2
where I 1 is the first invariant of strain tensor, J 2 is the second deviatoric invariant of the strain tensor and k is the ratio between uniaxial compressive and tensile strengths. This model uses only one single damage function

D = 1 - 1 -A ε eq + A exp (B (ε eq -ε D 0 )) (12)
where A andB are two constant of the model. In this contribution, Mazars' definition of equivalent strain is used

ε eq = 3 i=1 ε e i + 2 (13)
where ε e i + are the positive principal elastic strains. The damage variable D is determined as a linear combination of two damage variables D t and D c , that represent tensile damage and compressive damage respectively, by the help of two coefficients α t and α c computed from the elastic strain tensor and depending on the type of effective stress state:

D = α t D t + α c D c (14) D t,c = 1 - 1 -A t,c ε eq + A t,c exp (B t,c (ε eq -ε D 0 )) (15) α t,c = 3 i=1 ε t,c i ε e i + ε eq (16)
where A t , B t , A c , B c are four model parameters. ε t and ε c are the principal values of strains associated to the positive and negative effective stresses respectively.

The combination of plasticity and damage exhibits strain-softening and all the inherent difficulties attached to this specific material behavior, i.e., spurious strain localization and dependence of the energy dissipation on mesh refinement. To avoid this undesirable behaviour, the model may be regularized considering a nonlocal approach. In the present paper, a regularization technique, similar to that presented by [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], is chosen. It is based on implicit gradient definition of the nonlocal strain tensor which is calculated for each component of the elastic strain tensor as:

ε e ij = εe ij -c∇ 2 εe ij ( 17
)
where ∇ 2 is the Laplacian operator. The parameter c is of the dimension length squared and characterizes the non local interaction. This equation must be completed by the homogeneous Neumann boundary conditions

∇ε e .n = 0 (18)
with n the outward unit normal to the physical boundary. Now the local equivalent strain in equation ( 13) is replaced by its nonlocal counterpart:

ε eq = 3 i=1 ε e i + 2 (19)
3 STRUCTURAL APPLICATION

Bar under uniaxial tension

The validation of the regularization technique introduced to the damage part of the model is illustrated by means of a tension bar with a defect in the middle.

The bar length is 1 m and the cross section is 2 × 2 cm 2 . In simple tension tests, when the yield surface is reached, the effective response of the plastic damage constitutive law presented in this paper is perfectly plastic and no further damage can evolve due to a constant elastic strain. For this reason, the plastic process in this example is described by the Von Mises' model with linear hardening and no limit surface. The following parameters are used in this analysis: E = 3.3 × 10 10 Pa, ν = 0.2, A t = 0.9, B t = 10000, ε D 0 = 10.5 × 10 -5 , σ y = 3 × 10 6 Pa and E t = 1.1 × 10 10 Pa, where σ y the yield stress and E t the hardening modulus are the marerial properties for plastic part. Fig. 1 shows the damage distribution over the bar close to failure for a non local parameter c equal to 0.002 m 2 . The properties of the plastic gradient damage model are assessed by carrying out the analysis with several values of the gradient parameter c and a fixed mesh of 80 pentagon elements. Fig. 2 and Fig. 3 show the load-displacement response and the profile of damage for several values of the internal length parameter c. As desired, both the ductility in the loaddisplacement response and the width of the final damage profile increase with the internal length while for c= 0 the response of the regularized version coincide with that of the local one where the damage is localized in the center.

Shear failure mode: local approach

To demonstrate the abality and the performance of our constitutive law to predict correctly the failure mode in structural applications, the experimental behavior of reinforced concrete Continuous Deep Beam (reported by [START_REF] Asin | Behavior of statically interminate deep beams[END_REF]) has been simulated. The specimen, whith a minimum amount of stirups, has been selected for numerical simulation. However, the vertical shear reinforcement was not taken into account to simplify the study. The material parameters for the concrete are: E = 3.47 × 10 10 Pa, ν = 0.2, A t = 1, A c = 2.75, B t = 13000, B c = 2470, ε D 0 = 1 × 10 -4 , r c = 120 × 10 6 Pa and r t = 11 × 10 6 Pa. Moreover, a comparative numerical investigation has also been performed with two isotropic damage models in order to show advantages achieved by the introduction of plasticity into the damage formulation. The first one using the classical damage formulation defined by [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF] with the following parameters: E = 3.47 × 10 10 Pa, ν = 0.2, The geometry and the loading setup are presented in Fig. 4 where the thickness is 0.15 m. Because of symmetry, only (the right) half of the specimen has been meshed with 3048 eight node cubic elements. With these geometry and parameters, the simulations are performed with the local approaches, since, using the nonlocal model presented in this paper would require a mesh small enough to take into account non local interaction and would reach computational lim-its due to the memory size of the calculation. A satisfactory agreement of the experimental and numerical results is observed. obtained by the damage models, where only flexural cracks are formed. Fig. 7 compares the load deflection curves of the plastic damage and damage models. As expected, the peaks for the damage models are underestimated as the failure is in mode I. This strong difference is due to the fact that for Mazars' criteria is written only as a function of positive elastic strain (Eq. ( 13)) which means that it is only due to microcracks in mode I, while for the modified von Mises definition, the strength under biaxial tensile-compressive stress for concrete is overestimated which leads to overestimation of the shear resistance of the material in this model (see e.g. [START_REF] Patzák | Adaptivity resolution of the localized damage in quasi-brittle materials[END_REF]).

A t = 1, A c = 2.75, B t = 8000, B c = 1768, ε D 0 = 1 × 10 -4 ,

Size effect in mixed mode tests

Quasibrittle materials like concrete are charactarized by gradual softening in a fracture process zone (FPZ) that is independent and not negligeable compared to the structure size D. Hence, the response of geometrically similar specimens is not proportional and yield a material size effect [START_REF] Bažant | Size effect in blunt fracture: concrete, rock, metal[END_REF]. It can be explained by stress redistribution and the associated energy release due to the developement of large cracks or a large FPZ prior to failure. A complete explanation of the size effect can be found in the textbook by [START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasi-Brittle Materials[END_REF]. From the modelling viewpoint, the size effect can be described only with failure models that contain internal length which is related to the characteristic length of the material, i.e. the width of the FPZ [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure[END_REF]. The aim of this section is to highlight the existence of the size effect on the mixed mode failure as it is predicted by the plastic gradient enhanced damage formulation presented in this paper. Asymmetric four-point bending tests proposed by [START_REF] Gálvez | Mixed mode fracture of concrete under proportional and nonproportional loading[END_REF] In order to avoid some possible mesh bias and ensure that the global energy dissipation in failure process is captured correctly, the size of an element in the process zone should be kept constant and small enough compared to the characteristic length. c n is a constant depending on geometry which plays no role in size effect analysis. The ultimate loads P u for the three sizes and its corresponding nominal strength are listed in the Table 1. The numerical results are interpreted with the help of Bazant's size effect law [START_REF] Bažant | Size effect in blunt fracture: concrete, rock, metal[END_REF].

σ N = Bf ′ t (1 + D/D 0 ) -1/2 (21)
where B is a dimentionless geometry-dependent parameter, D 0 is a characteristic size and f ′ t is the tensile strength of the material. Bf slope -1/2 where these two lines intersect at the abscissa D = D 0 . Obviously, it can be seen that there is a pronounced size effect that is in a good agreement with Bazant's size effect relation. The larger the specimen, the lower the nominal strength. It is also observed that the plasticity plays more important role for the smaller specimen (the hardening variable k h is more important as the specimen is smaller). On the other hand, for the larger specimen the developement of damage is dominant and its behavior is closer to linear elastic fracture mechanics and more different from a limit analysis solution.

CONCLUSION

In this work, the hardening plasticity is combined with a gradient enhanced damage formulation for modelling the structure behaviour with mixed mode failure. The example of tensile bar shows the effectivity of the regularization technique to overcome the mathematical problems due to the strain softening behaviour (damage part) where the damage distribution zone is controlled by a material heterogeneity parameter. Furthermore, the simulations of the Continuous Deep Beam show the improvements achieved by the coupled model compared to the damage models. The model predicts rather accurately the failure mode (e.g. mode II) and the ultimate load in complex structures.

Finally, the numerical simulations of three asymmetric four-point bending tests on notched specimens of different sizes show a dependence of the nominal strength on the specimen size. The larger the specimen, the lower the nominal strength. Thus, the nonlocal approach exhibits a size effect in mixed mode failure which is in a good agreement with Bazant's size effect relation.

Figure 1 :Figure 2 :

 12 Figure 1: Development of damage in a tensile bar with a defect in the middle.

Figure 3 :

 3 Figure 3: Damage profiles along the beam for different values of the internal length scale.

  and the other one employing the modified von Mises definition of the equivalent strain according to de Vree et al. (1995) with the following parameters: E = 3.47 × 10 10 Pa, ν = 0.2, A t = 0.9, B t = 8000 and ε D 0 = 1 × 10 -4 . The reinforcement bars are modelled with a von Mises plasticity model with E = 2.1 × 10 11 Pa, ν = 0.3, σ Y = 5 × 10 8 Pa (yield stress), Et = 4.5 × 10 10 Pa (tangent modulus).

Figure 4 :

 4 Figure 4: Geometry and loading setup of CDB.

  Fig. 5 provides a comparison of the experimental crack pattern with the numerical damage profile of the plastic damage model. As seen in the experiments, flexural damage band (mode I) starts first at midspan (stage A). Further load increase leads to the development of the flexural damage band and a propagation of damage (mode I & II) continues over the support (stage B), then a shear damage band (pure mode II) formed suddenly (stage C). Finally, failure occurs at the load introduction zone near the loading column (stage D).The ultimate failure load obtained experimentally was 1180 kN , the corresponding numerical value is 1285 kN , which is overestimated by 11%. Thus, the model can reproduce correctly the mode of failure and predict rather well the ultimate load.

Figure 5 :Figure 6 :

 56 Figure 5: Developement of crack pattern in CDB; Experiment versus simulation with plastic damage model.

Figure 7 :

 7 Figure 7: Load-deflection curves. Comparison between plastic damage model, Mazars' model and modified von Mises' model.

Figure 8 :

 8 Figure 8: Geometry and loading setup for asymmetic four- point bending tests.

Figure 9 :

 9 Figure 9: Profile of damage for small size specimen.

Fig. 10

 10 Fig. 10 shows the load deflection responses for the three different sizes. The nominal strength is obtained with the formula σ N = c n P u bD (20)

Figure 10 :

 10 Figure 10: Load-displacement curves for small size specimen (D=75 mm), medium size specimen (D=150 mm) and large size specimen (D=300 mm).

′tFigure 11 :

 11 Figure 11: Size effect on nominal strength.

Table 1 :

 1 Numerical results for three different sizes.

	Size	Depth Ultimate load Nominal strength
		D	P u	σ N
		(mm)	(kN)	(MPa)
	small	75	8.213	2.19
	medium	150	14.67	1.956
	large	300	24.57	1.638