Masoud Rostami

Vladimir Zeitlin
email: zeitlin@lmd.ens.fr

Eastward-Moving Convection-Enhanced Modons in Shallow Water in the Equatorial Tangent Plane

We report a discovery of steady long-living slowly eastward moving large-scale coherent twin cyclones, the equatorial modons, in the shallow water model in the equatorial beta-plane, the archetype model of the ocean and atmosphere dynamics in tropics. We start by constructing analytical asymptotic modon solutions in the non-divergent velocity approximation, and then show by simulations with a high-resolution numerical scheme that such configurations evolve into steady dipolar solutions of the full model. In the atmospheric context, the modons persist in the presence of moist convection, being accompanied and enhanced by specific patterns of water-vapour condensation.

As is known, the large-scale atmosphere and ocean dynamics in tropics is largely determined by equatorial waves. Indeed, superposition of atmospheric Outgoing Longwave Radiation (OLR) data onto the spectrum of equatorial waves gives a very good correspondence [1], except for a branch corresponding to slow eastward-moving motions associated with the Madden-Julian Oscillation (MJO) [2]. The Kelvin waves, the only species of the equatorial waves that can produce steady large-scale eastward-moving structures, have much faster velocity than MJO [3]. This fact motivates a search for coherent dynamical structures with slow eastward propagation which, thus, should be necessarily nonlinear. There do exist steady eastward-moving structures in a fluid layer on the mid-latitude f -plane, the modons which are exact dipolar solutions of the quasi-geostrophic (QG) equations [4]. Although the classical modons were found analytically in the QG approximation, it was shown first by a computer-assisted analysis [5] and then by direct numerical simulations [6] that corresponding solutions exist also in the parent rotating shallow water (RSW) model on the f -plane.

The QG modon solutions were extended to the full sphere [7,8] in the framework of equivalent barotropic model, which is structurally close to QG. The pioneering paper [9] advanced an idea that the MJO could be related to such a modon in spherical geometry and showed that, like in the midlatitude tangent plane [6], the equivalent barotropic modon persists in the RSW on the sphere. Yet, the archetype model for understanding dynamics in the tropics is RSW in the equatorial beta-plane, e.g. [10]. We, thus, look for eastward-moving modons in this model. However, it encounters an obstacle from the very beginning, because there is no consistent QG approximation at the equator, e.g. [START_REF] James | Introduction to Circulating Atmospheres[END_REF], so the known modon solutions can not be borrowed. Our main observation is that there is a dynamical regime in RSW, called long-wave approximation in oceanography [START_REF] Gill | Atmosphere -Ocean Dynamics[END_REF], which corresponds to small pressure variations, and gives, to the leading order, equations which do allow for modon solutions.

In the atmospheric context this regime is the RSW analog of the non-divergent equatorial balance model [START_REF] Charney | [END_REF]. The relevance of this model to large-scale tropical motions in the atmosphere is supported by scale and data analyses [14,15]. We construct the modon solutions in this regime and, following [6], use these asymptotic solutions to initialize high-resolution numerical simulations with the full RSW model, and show that coherent dipolar steady eastward-moving structures do arise and persist. We then show that, in the atmospheric context, inclusion of moisture, with condensation and evaporation, enhances the modons, without disrupting long-time coherence, and produces specific convective patterns.

We start with RSW equations in the equatorial beta-plane with no dissipation:

     ∂ t v + v • ∇v + β y ẑ ∧ v + g∇h = 0 , ∂ t h + ∇ • (vh) = 0 , (1)
where ∇ = (∂ x , ∂ y ), v = (u, v), u is zonal and v-meridional components of velocity, h is geopotential height (thickness), β is the meridional gradient of the Coriolis parameter, and ẑ is a unit vertical vector. The interpretation of the model in the oceanic context is direct, while in the atmospheric context it should be understood as vertically averaged primitive equations in pseudo-height pressure coordinates [START_REF] Zeitlin | Geophysical Fluid Dynamics: Understanding (almost) everything with rotating shallow water models[END_REF], and can be extended to include water vapor with condensation and related latent heat release [START_REF] Bouchut | [END_REF], as well as surface evaporation (see below).

We introduce a pressure perturbation parameter λ: h = H(1 + λη), where H is unperturbed thickness, and fix spatial, velocity and time-scales: (x, y) ∼ L, (u, v) ∼ V, t ∼ L/V . Under hypothesis that λ → 0, and gHλ/V 2 = O(1), hence V << √ gH, i.e. the characteristic velocity is much smaller than the phase velocity of the Kelvin waves c = √ gH, the nondimensional equations take the form:

∂ t v + v • ∇v + β y ẑ ∧ v + ∇η = 0 , (2)
λ(∂ t η + v • ∇η) + (1 + λη) ∇ • v = 0, (3)
where β = βL 2 /V , and v = v 0 + λv 1 + .... In the leading order in λ (3) gives ∇ • v 0 = 0, the motion is non-divergent, and u 0 = -∂ y ψ, v 0 = ∂ x ψ. Cross-differentiation of the zonal and meridional momentum equations results in the equation for the stream-function:

∇ 2 ψ t + J (ψ, ∇ 2 ψ) + βψ x = 0, (4)
where J denotes the Jacobian. The modon solutions are built following [4]. They are are obtained under hypothesis of steady motion with constant zonal velocity U , by supposing a linear relationship between the absolute vorticity and stream-function in co-moving frame, which gives inhomogeneous Helmholtz equation. It is solved by separation of variables in polar coordinates in terms of Bessel functions, first in the outer domain under condition of decay, and then in the inner domain, and matching the inner and outer solutions across a circle of a given radius a in the plane. The solution has the form: where J 1 and K 1 are ordinary and modified Bessel functions of order one, respectively, p is real, p 2 = β/U , so U > 0, and the motion is eastward. For each pair (a, p) there exist a series of eigenvalues k arising from matching conditions, of which the lowest corresponds to a dipole, to be called the asymptotic modon. The next eigenvalue gives a quadrupolar solution called "shielded modon", and so on. By construction, this solution can be centred anywhere in the equatorial beta-plane, yet it is only when its center is at the equator that it results in a cyclonic pair. The streamlines of the asymptotic modon are presented in Fig. 1.

           ψ ext = -U a K 1 (pa) K 1 (pr) sin θ, r > a, ψ int = U p 2 k 2 J 1 (ka) J 1 (kr) -r k 2 (1 + U + U k 2 ) sin θ, r < a, (5)
First-order corrections u 1 , v 1 can be found in the next order of the asymptotic expansion.

There is no guarantee that such expansion converges, but numerical results presented below suggest that there is indeed an exact solution of (1) corresponding to this asymptotic one.

We initialized the numerical simulations in the RSW model with the velocity field corresponding to the asymptotic modon solution (5) with U √ gH, and flat pressure. Simulations were performed with high-resolution well-balanced finite-volume numerical scheme [START_REF] Bouchut | Nonlinear Dynamics of Rotating Shallow Water: Methods and Advances[END_REF] in a rectangular domain with sponges at the boundaries, which allow to mostly evacuate short inertia-gravity waves. No explicit dissipation was added. The domain was chosen to be symmetric with respect to the equator and wide enough 9L d ×6.5L d , where L d = (gH) 1/4 /β 1/2 is the equatorial deformation radius, in order to minimize the influence of boundaries onto the modons. We varied spatial resolution from 600 × 600 to 1200 × 1200, to check numerical convergence. The natural units of length and time in the numerical scheme are L d and 1/βL d , diminishes with increasing resolution (several % loss at high resolution), not shown. The difference between the asymptotic and "exact" modons is clear from the comparison of their relative vorticity in Figure 5, which shows that the "exact" modon is more compact and has a larger peak vorticity. Although we did not explore in detail the space of parameters of the modon solutions, we varied the parameter a, the modon's radius. We observed that while for a < L d equatorial modons always emerge, for a > L d , with the same velocity scale, the initial asymptotic modon transforms into a packet of equatorial Rossby waves moving westward. We also checked the robustness of the equatorial modon to the details of initialization, by either removing the outer field of the asymptotic modon, or misaligning its axis. In all such simulations the eastward-moving equatorial modon was always emerging after leaving some "debris" and/or wobbling at initial stages. This indicates that equatorial modons are attracting solutions.

After having established existence of equatorial modons in the adiabatic environment, we switched on the effects of moisture within the diabatic atmospheric moist-convective RSW (mcRSW) [START_REF] Bouchut | [END_REF]. The equations of the model read:

           ∂ t v + v • ∇v + β y ẑ ∧ v = -g∇h, ∂ t h + ∇• (vh) = -γP, ∂ t Q + ∇•(Qv) = -P + E. ( 6
)
Here Q is a bulk amount of water vapor in the air column, γ is a parameter depending on the underlying stratification, P is the condensation sink, and E is the evaporation source of moisture, which are parameterized as follows:

P = Q -Q s τ H(Q -Q s ), E = α|v|(Q s -Q)H(Q s -Q). (7)
Q s is a saturation moisture threshold, τ is relaxation time, α is a parameter regulating evaporation, and H(...) denotes the Heaviside function. As in previous test simulations [START_REF] Bouchut | [END_REF],

we take a uniform initial moisture distribution Q i close to saturation:

Q i = Q s -0.01 with
Q s = 0.9, a short relaxation time τ equal to several time-steps of the numerical scheme, γ = 1, and α = O(10 -1 ). We performed exactly the same simulations as in the adiabatic case. They showed that the equatorial modons keep their coherence in the moist-convective environment. Moreover, as is known, and confirmed in mcRSW [19,20], condensation leads to intensification of the cyclonic vortices. Correspondingly, the potential vorticity (PV)

anomalies become stronger, the size smaller, and the phase speed of the convectively-coupled modon becomes higher than that of its adiabatic counterpart. A comparison of evolution of the same modon in adiabatic and moist-convective environments is presented in Fig. 6.

A typical pattern of the associated moisture field, and characteristic condensation patterns at the front and at the rear of the modon are shown in Fig. 7. Such patterns persist all along the simulation. We should stress that, although we initialized the simulations with non-divergent velocity field of the asymptotic modon, the "exact" modon develops and maintains a characteristic convergence/divergence pattern (not shown), which redistributes moisture in its core, and prompts condensation in specific zones. Notice that as condensation enhances cyclonic vorticity [19], the moist-convective modon intensifies, and moves faster than the dry one.

To conclude, we established existence of long-living, slow eastward-moving dipolar coherent structures, the modons, in the shallow-water dynamics in the equatorial beta-plane. In cannot make a firm link to MJO at this stage, some features of the equatorial modon do resemble those of dipolar structures associated with the MJO [21]. Work is in progress in this direction, in particular including improved parametrisation of moist processes [20] and vertical structure, which has specific features in MJO, cf. [22] and references therein. Independently of a possible link to MJO, the very existence of the eastward-moving equatorial modons, with their characteristic velocity pattern, slow eastward propagation, and relatively weak signature in pressure should be kept in mind in data analyses and simulations with "big" atmospheric and oceanic models. It is known [23] that long-wave pressure anomalies in RSW at the equator produce westward-propagating dipolar Rossby-wave packets and eastward-propagating Kelvin waves, in accordance with the Gill's scenario [24]. The genesis of eastward-moving dipoles is an important topic, to be addressed elsewhere. Our analysis above suggests that comparable zonal and meridional scales of the initial disturbance, and relative smallness of pressure anomaly are necessary for that.

Let us comment on the relation of equatorial modons we discovered to those known, or hypothesized, in the literature. First, while the relation of the spherical modons to the standard QG ones on the midlatitude tangent plane was established in the limit of small spherical modon size [7], this proof, which uses the structural resemblance of equivalent barotropic equations on the sphere and QG equations, both incorporating a finite barotropic deformation radius, is not directly transposable to the equatorial beta-plane case. So the proof that our modons is a limiting case of the spherical ones is pending. Second, westwardpropagating localized solutions, which were also called modons, were constructed in [25] as a generalization of westward-moving equatorial Rossby-wave solitons, and confirmed by direct numerical simulations in [23]. The long-wave scaling used in [25], [23] differs from ours, with the main difference residing in the disparity of zonal and meridional scales. A possibility of existence of eastward-propagating modons, if some ad hoc "heuristic" terms were added to the systematically derived asymptotic equations, was also evoked in [25]. Question remains whether our modons could be related to this conjecture, although the difference in scalings seems to prevent this. We should add, for completeness, that meridionally translating zonally-symmetric modon-like solutions in the vertical plane at the equator were derived in [26] in the context of the atmospheric equatorial boundary layer.

1 FIG. 1 .

 11 FIG. 1. Normalised streamlines and velocity field of an asymptotic modon in stationary (left ) and co-moving (right) frames. Dashed circle: separatrix of radius a.

FIG. 2 .Fig. 4 .FIG. 4 .

 244 FIG. 2. Initial adjustment of the asymptotic modon with non-dimensional a = 0.5, U = 0.15, as seen in the thickness and velocity fields. Time in units of 1/βL d .

FIG. 5 .

 5 FIG. 5. Relative vorticity of the asymptotic (left) vs "exact" (right) modons. "Exact" modon is more compact and more intense.

10 FIG. 6 .

 106 FIG. 6. Steady eastward propagation of the equatorial modon, as seen in PV field in moistconvective (left) and adiabatic (right) environments, with the same initial conditions. t = 0, 10, 20, 30, 40 [1/(βL d )], from top to bottom. (a = 0.5, U = 0.2, α = 0.05).

FIG. 7 .

 7 FIG. 7. (a = 0.5, U = 0.2, α = 0.05). Distribution of the water vapor (left) and superposition of pressure (black) and condensation (red) isopleths (right) at t = 30 [1/(βL d )].

ACKNOWLEDGMENTS

We are grateful to J.-I. Yano for his constructive remarks and useful discussions.


J P Crété

P Longère

J M Cadou

X-FEM ASSISTED SIMULATION OF DUCTILE DAMAGE INDUCED 2D CRACK PROPAGATION

Keywords: X-FEM, viscoplasticity, damage, computational mechanics

The present work is devoted to the numerical simulation of crack propagation in engineering materials whose failure results from void initiation, growth and coalescence. The behavior of the plate material is described via a Gurson type model accounting for the combined effects of strain hardening, thermal softening, viscoplasticity and void growth induced damage. The eXtended Finite Element Method has been retained to describe the kinematic consequences of the crack propagation across the mesh. The crack is assumed to propagate as soon as the stored energy around the crack tip reaches a critical value. The related crack length is estimated from the crack velocity considering the current time increment. The constitutive model and the extended finite elements were both implemented in the engineering FE computation code Abaqus as user subroutines. The numerical simulation of a notched plate under tension loading has been conducted. While making some simplifications, the present work reproduces numerically the 2D propagation of a crack resulting from void growth induced damage.

INTRODUCTION

The present work is devoted to the numerical simulation of crack propagation in engineering materials whose failure results from ductile damage The fracture of ductile materials is known to result from void initiation, growth and coalescence. During the process of void growth induced damage, the bulk material is subject to a progressive loss of its overall properties and to the appearance, in addition to the isochoric plastic deformation due to dislocation glide in the matrix material, of an inelastic dilatancy due to void growth. In order to describe these effects, the behavior of the material at stake is described via a modified version of the micromechanics based Gurson model, see Gurson ( [4]), Tvergaard and Needleman ( [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]) and Longère et al ( [START_REF] Longère | Dragon: Modelling the Transition between Dense Metal and Damaged (Micro-Porous) Metal Viscoplasticity[END_REF]), accounting for the combined effects of strain hardening, thermal softening, viscoplasticity and void growth induced damage. The numerical treatment of crack propagation in structures is not trivial. The first attempts to achieve it mainly consisted in meshing very finely the crack front to reproduce as accurately as possible the stress-displacement fields at the crack tip. This method implies a very large number of meshes and may be prohibitive when applied to engineering complex structures. Very recently, methods consisting in enriching the finite element kinematics to account for the aforementioned singular fields at the crack tip have emerged. One of these methods, namely the eXtended Finite Element method (X-FEM) [START_REF] Moes | Belytschko: A finite element method for crack growth without remeshing[END_REF], has been here retained to describe the kinematic consequences of the crack propagation across the mesh. In the present work, the crack is assumed to propagate as soon as the stored energy around the crack tip reaches a critical value. The related crack length is estimated from the crack velocity considering the current time increment. To estimate the ability of the adopted approach, the numerical simulation of a notched plate under tension loading has been conducted employing the engineering FE computation code Abaqus [START_REF] Hibbitt | ABAQUS Users Manual, Version 6.9[END_REF]. While making some simplifications, the present work reproduces numerically the 2D propagation of a crack resulting from void growth induced damage. The constitutive model and the extended finite elements were both implemented in the computation code Abaqus as user subroutines. The constitutive model of the material considered in the present study is outlined in Sect.2. The finite element enrichment technique adopted in the present work is described in Sect.3. The combination of the constitutive model with the finite enrichment method is detailed in Sect.4. The application of the approach to the numerical simulation of the tension loading of a notched plate is shown and commented in Sect.5.

MATERIAL CONSTITUTIVE MODELLING

Constitutive equations

The failure of the material considered in the present work is supposed to result from void initiation, growth and coalescence. To reproduce the process of ductile damage at stake, the material behavior is described via a Gurson type model, extended by Tvergaard and Needleman.

The GTN damage-plastic potential may be written as :

Φ GT N = ( σ eq σ y ) 2 + 2q 1 f cosh(- 3 2 q 2 p m σ y ) -(1 + q 3 f 2 ) = 0 (1)
where σ eq represents the equivalent stress, p m the pressure, σ y the rate dependent yield stress, f the volume fraction of voids, and where (q 1 , q 2 , q 3 ) are material constants. Motivated by physical considerations, Longère et al [START_REF] Longère | Dragon: Modelling the Transition between Dense Metal and Damaged (Micro-Porous) Metal Viscoplasticity[END_REF] recently introduced in GTN model a mean stress provoking a shift of the yield locus towards negative stress triaxialities allowing to describe damage growth under shear dominated loading. The modified damage-plastic potential is accordingly proposed in the form :

Φ GT N = ( σ eq σ y ) 2 + 2q 1 f cosh(- 3 2 q 2 p m + p r σ y ) -(1 + q 3 f 2 ) = 0 (2)
p r = b ln(q 1 f ) (3)
with b a positive constant value. The yield stress σ y accounts for the combined effects of strain hardening, thermal softening and viscoplasticity :

σ y (κ, κ, T, ...) = σ y (κ, T ) + σ vp ( κ, T, ...) (4)
where κ represent the matrix plastic strain, κ the matrix plastic strain rate and T absolute temperature. In the rate independent yield stress σ y (κ, T ), the isotropic strain hardening is described via a Voce type law and the thermal softening via a power law :

σ y = (R 0 + R ∞ [1 -exp(-kκ)] β ) 1 - T T ref m (5)
where (R 0 , R ∞ , k, β) are isotropic hardening related constants and (T ref , m) thermal softening related constants. The strain rate induced overstress reproduces the potential tension/compression asymmetry ( [START_REF] Graff | Strain localization phenomena associated with static and dynamic strain ageing in notched specimens: experiments and finite element simulations[END_REF]) :

σ vp = Y [ κ exp( V a p m k b T )] 1 n (6)
where (Y ,n) are viscosity related constants and (V a ,k B ) behaviour asymmetry related constants, with V a = V h β 3 where V h is a constant and β Burgers vector magnitude (β=2.5 Å), and with k B Boltzmann constant (k B =1.3804.10 -23J/K). Assuming the normality rule, the evolution laws are expressed as:

d ∼ p = Λ ∂Φ GT N ∂σ ∼ = Λ( ∂Φ GT N ∂σ eq n ∼ - 1 3
∂Φ GT N ∂p m δ ∼ ) = ˙ pD n ∼ + 1 3 ˙ pM δ ∼ ( 7
)
n ∼ = 3 2 s ∼ σ eq (8)
where the distortional and dilatational parts, namely ˙ pD and ˙ pM , respectively, of the inelastic strain rate d ∼ p are given by

˙ pD = Λ ∂Φ GT N ∂σ eq = 2Λ σeq σ y (9)
˙ pM = -Λ ∂Φ GT N ∂p m = 3q 1 q 2 f Λ sinh[-3 2 q 2 (p m + pr )] σ y (10)
with σeq = σeq σy , pm = pm σy and pr = pr σy . The evolution law of the isotropic hardening variable κ is deduced from the equality of the macroscopic plastic work rate with the microscopic one, see Gurson ( [4]) :

κ = σ eq ˙ pD -p m ˙ pM (1 -f )σ y (11)
According to the aforementioned considerations, adiabatic heating under dynamic evolution is evaluated from

ρC Ṫ = σ eq ˙ pD -p m ˙ pM -r κ (12)
The porosity rate ḟ is decomposed into a contribution due to growth of existing defects, namely ḟg , and a contribution due to the formation of new defects, namely ḟn :

ḟ = ḟg + ḟn (13) ḟg = (1 -f )Trd ∼ p = (1 -f ) ˙ pM (14) ḟg (0) = f 0 ḟn = B σy (15) ḟn (0) = 0 B = f sup p σ c Φ p-1 Io exp(-Φ p Io ) (16)
where p is a constant (p=2) and where f sup represents the upper bound of the nucleated secondary void volume fraction. Φ Io is the hole nucleation criterion, see Longère & al. ( [6]).

Numerical procedure

The material behavior outlined in the previous subsection is implemented as user material (UMAT) in the engineering finite element computation code Abaqus. The numerical integration is achieved using the classical return mapping procedure combined with the Newton-Raphson solving algorithm, see Aravas ( [10]) for further details. Time increments are small and the tangent operator is elastic. Dividing Eq.( 9) by Eq.( 10), or inversely, allows for eliminating the viscoplastic multiplier and defining from the incremental viewpoint :

Ξ ∆ = ∆ pD ∂Φ GT N ∂p m + ∆ pM ∂Φ GT N ∂σ eq = 0 (17)
The numerical integration consists thus in solving the following system of equations

Φ GT N (σ eq , p m ; H α ) = 0 (18) Ξ ∆ = 0 (19)
p m = p e m + K∆ pM (20)
σ eq = σ e eq -3µ∆ pD (21)

∆H = h(∆ pD , ∆ pM , σ eq , p m , H) (22)
where σ e eq and p e m represent the trial, equivalent stress and pressure respectively, and ∆H the system of complementary laws Eqs.( 9)-( 10) and Eqs.( 35)-( 16) written in the incremental form. The constants K and µ represent the bulk and shear moduli respectively. Adiabatic conditions are assumed for strain rates greater than 1s -1 .

FE ENRICHMENT

Principle

The extended finite element method (X-FEM) [START_REF] Moes | Belytschko: A finite element method for crack growth without remeshing[END_REF] consists in embedding strong discontinuities, induced by e.g. cracks, inside finite elements while enriching the regular displacement of the affected finite elements with singular functions attempting to reproduce the kinematic consequences of the discontinuities at stake. The current displacement field u(x, t) is generally expressed as :

u(x, t) = i∈I u i (t)N i (x) + j∈J b j (t)N j (x)H(x) + k∈K N k (x) 4 l=1 c l k (t)F l (x) (23)
where u i (t) represents the regular nodal displacement of node i, b j (t) the discontinuous nodal displacement magnitude of node j belonging to a crack-crossed finite element, and c l k (t) the singular nodal displacement magnitude of node k belonging to a crack tip-containing finite element. N i (x) represent the shape functions, F l (x) the four singular functions at the crack tip, and H(x) the generalized Heaviside funtion with H(x) = sign(d Γ (x)) where d Γ (x) denotes a signed distance function from the crack (H(x) takes the value +1 for the nodes located above the crack and -1 for the nodes located below the crack). The total set of nodes is denoted as I, the set of nodes belonging to crack-crossed finite elements as J, and the set of nodes belonging to crack tip-containing finite elements as K, see 1. In the case of a two-dimensional (2D) elastic problem, the singular functions F l (x) are expressed in a cylindrical frame (r, θ) as:

F l=1,2,3,4 (r, θ) = [ √ rsin θ 2 , √ rcos θ 2 , √ rsinθsin θ 2 , √ rsinθcos θ 2 ] (24)

Adopted approach

Regarding the singular displacement field at the crack tip

The X-FE method is particularly efficient when applied to 2D problems involving elastic-brittle materials for which the singular functions F l (x) can be determined analytically, see Eq.( 24). Dealing with elastic-(visco)plastic materials makes the resolution much more difficult because the aforementioned singular functions are generally not a priori known. In certain cases, the functions may be approximated by analytical studies considering simplified material behaviors or by numerical simulations. The fact that the material of the present study is subject to strain hardening, thermal softening, viscoplasticity and ductile damage, see Sect.2, increases significantly the complexity of the singular functions in question. An alternative way consists in using the singular functions defined in the elastic context so as to represent the presence of crack in the element containing the crack tip. Using only the Heaviside function without other enrichment functions, the crack grows element by element, i.e. the crack tip is always located on an edge of an element, so the position of the crack tip is approximate and the field the near crack tip can be perturbed. Fig. 2 represents the evolution of the equivalent stress which is averaged over a patch near the crack tip, see [START_REF] Haboussa | Combescure : X-FEM analysis of the effects of holes or other cracks on dynamic crack propagations[END_REF], for an elastic-plastic notched plate under tensile loading. According to Fig. 2, using the sole Heaviside function allows to describe the field near the crack tip with lightly finest mesh compared with X-FEM using singular functions. Moreover neglecting the singular functions allows to reduce drastically the required number of degrees of freedom (dof), 48 dof with the singular functions vs. 8 dof with the sole Heaviside function.

In absence of knowledge of the singular functions for the complex material behavior considered in the present work, and according to the aforementioned remarks, we are here using the reduced current displacement field X-FEM with singular functions (mesh size = 0,049mm) X-FEM without singular functions (mesh size = 0,033mm)

u(x, t) = i∈I u i (t)N i (x) + j∈J b j (t)N j (x)H(x) (25)

Regarding the numerical integration

There exist several techniques devoted to capture the crack propagation. The original method consists in subdividing the finite element in which the crack is currently propagating into adaptive sub-triangles, see Fig. 3. This however implies re-projecting the state variables values defined at the integration points of the original finite element onto the integration points of the sub-finite elements. The procedure which is easy when dealing with linear evolutions (as it is the case for elastic behavior) is far from being trivial when dealing with non-linear evolutions (it is the case for the present damage-plastic behavior). For this reason, the technique consisting in increasing significantly the number of integration points [START_REF] Elguedj | Combescure: Appropriate extended functions for X-FEM simulation of plastic fracture mechanics[END_REF] of the original finite element has been preferred in the present case. We are here using 64 integration points.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • a) First increment b) Second increment

Concerning the nodes belonging to cut and uncut finite elements

To ensure that the enrichment vanishes in all elements not cut by the crack, see [START_REF] Zi | Belytschko: New crack-tip elements for XFEM and applications to cohesive cracks[END_REF], the current displacement field u(x, t) takes thus the form

u(x, t) = i∈I u i (t)N i (x) + j∈J b j (t)N j (x)(H(x) -H(j)) (26)
where H(j) corresponds to the global discontinuous step function at node j.

Coupling X-FEM and ductile damage

The physical process of crack formation being very complex, in particular the transition between the stages of void coalescence and crack germination, we are here tentatively assuming that the transition in question is instantaneous. To avoid too strong mesh size dependence, we are adopting a damage-crack transition criterion based on the stored energy, averaging over a patch at the crack tip. The patch method has notably been used by Haboussa & al.( [8]) involving a mean stress based criterion.

Crack propagation criterion

The instantaneous material state is supposed to be well described via the Helmholtz free energy ω( e ∼ , κ, T ) where e ∼ represents the elastic strain tensor. The Helmholtz free energy is decomposed into a reversible and a stored part, namely ω r ( e ∼ ) and ω s (κ, T ) respectively:

ω( e ∼ , κ, T ) = ω r ( e ∼ ) + ω s (κ, T ) (27)
According to the constitutive model outlined in Sect.2, the two contributions to the Helmholtz free energy are expressed as

ω r ( e ∼ ) = 1 2 e ∼ : C ∼ ∼ : e ∼ ( 28
)
ω s (κ, T ) = h(κ)g(T )
where h(κ) represents the stored energy of cold work and g(T ) the thermal softening function. The thermodynamic forces are written as

σ ∼ = ∂ω ∂ e ∼ = dω r d e ∼ = C ∼ ∼ : e ∼ (29) r = ∂ω ∂κ = ∂ω s ∂κ = h (κ)g(T )
where r represents the matrix isotropic hardening force. According to Eq.( 5), the expressions of the strain hardening and thermal softening functions are given by

σ y = R 0 g(T ) + r(κ, T ) (30) h (κ) = R ∞ [1 -exp(-kκ)] β (31) g(T ) = 1 - T T ref m (32)
From the incremental viewpoint, the value of the stored energy at current step (n+1) is deduced from its value at previous step (n) via

ω n+1 s = ω n s + ∆ω s ( 33
)
∆ω s = ∂ω s ∂κ ∆κ + ∂ω s ∂T ∆T = h (κ)g(T )∆κ + h(κ)g (T )∆T
This quantity is then averaged over a patch containing a finite number of elements namely, p (according to Eq.( 35)), located at the crack tip and in the crack growth direction according to

W patch = 1 A p i=1 ω i s A i (34)
A = p i=1 A i
where A represents the patch area. Keeping constant the total patch area for configurations involving different mesh sizes allows for reducing the mesh size dependence of the numerical results. The principle is illustrated in Fig. 4. The crack propagates if

F (W patch ) = 1 - W c W patch > 0 (35)

Angle of propagation θ c

Patch for crack growth criteria

Crack propagation speed

The crack growth rate is estimated using the expression developed by [START_REF] Kanninen | Advanced Fracture Mechanics[END_REF] and adapted to the present approach :

ȧ = C r F (W patch ) = C r (1 - W c W patch ) (36)
where ȧ represents the crack growth rate and C r the Rayleigh wave celerity in the material. The value of the critical value W c depends on the material and on the patch radius. From the numerical viewpoint, the crack advance increment is deduced from

∆a = C r (1 - W c W patch )∆t (37)

Crack propagation orientation

Identifying the crack propagation orientation constitutes another challenge in the numerical treatment of crack propagation, in particular when dealing with elastic-(visco)plastic materials. Fracture mechanics based criteria, see e.g. Haboussa et al.( [8]), as well as bifurcation analyses, see e.g. Hill [START_REF] Hill | Acceleration waves in solids[END_REF] or Rice [START_REF] Rice | The Localization of Plastic Deformation[END_REF], are used depending on the application. In a first approach, we used the former works with some modifications. The presence of damage near the crack tip involves a stress drop making the use of the stress field unsuitable from the physics point of view. Therefore an approach based on plastic strain field near the crack tip seems to be more efficient. In a first step, we define a local plastic strain tensor averaged in a patch which is perpendicular to the crack. The crack propagation direction is tentatively deduced from

θ c = 2arctan( 1 4

APPLICATION TO MODE I CRACK PROPAGATION IN A NOTCHED PLATE MADE OF A DUCTILE MATERIAL

To evaluate the prediction ability of our approach, we are considering the crack propagation in a 2D (plane strain) notched plate made of a ductile material, namely the mild steel in Longère et al [START_REF] Longère | Dragon: Modelling the Transition between Dense Metal and Damaged (Micro-Porous) Metal Viscoplasticity[END_REF] study, submitted to Mode I loading, see Fig. 5. Finite elements are 4 node-quadrilateral with 64 integration points, enriched according to the formulation outlined in Sect.3. The material behavior is described via a user material subroutine integrating the rate equations given in Sect.2. For reason of confidentiality, some constant values are not given. The nodes of the plate upper side are submitted to a constant velocity (1m/s) whereas those of the lower side are constrained.

Tension without crack propagation

In a first step, we study the notched plate without crack propagation. The evolution of the energy W patch storaged in a patch (see Eq.35) with a radius of 6.25mm and a angle of 0 • for three different mesh size (the mesh 3 is the smallest mesh size and the mesh 1 the most coarse) is drawn in Fig. 6. According to Fig. 6 one can observe a convergence of the curves which denotes a weak influence of the mesh size on the crack propagation criterion proposed in Sect.4. In Fig. 7, the porosity map arround the crack tip at 0.001s is presented. According to Fig. 7, the porosity is concentrated near the crack tip, provoking the large mesh dependency observed when a failure porosity is used to represent the crack propagation [START_REF] Hambli | Comparison between Lemaitre and Gurson damage models in crack growth simulation during blanking process[END_REF].

The method presented in Sect.4 which allows to propagate a crack, needs two parameters (W c the crack growth criterion and the radius of the patch). They have to be defined from experimental data. In a first approach, to evaluate the efficiency of the model we use arbitrary values. The crack growth quantity W c is chosen to be equal to 4.10 6 J.m -3 and the radius value of the patch used to definine the crack propagation direction is 6.25mm.

Tension with crack propagation

The history of the upper side reaction force for configurations with and without crack propagation are drawn in Fig. 8. According to Fig. 8, one can clearly observe that the crack propagation provokes a drop in the plate response, as expected. The evolution of the averaged stored energy W patch is presented in Fig. 8. According to Fig. 8, Wpatch increases with the increasing displacement (or time), reaches W c , rises above W c , then drops suddenly to a value below W c , then increases, reaches W c and repeats the aforementionned scenario. The sequential evolution of Wpatch may be explained by the need for the crack to cut completely the element. Refining the mesh size would reduce the magnitude of the fluctuations of Wpatch around Wc and smooth the evolution.

The map porosity in the plate at 4ms is drawn in Fig. 9. According to Fig. 9, the porosity is concentrated along the crack and as expected the crack growth follows an angle of 0 approxi-mately. One can observe numerical artefacts consisting in a concentration of porosity in the two corners of the plate which are due to boundary conditions.

CONCLUDING REMARKS

• This work is devoted to the numerical simulation of crack propagation in engineering materials whose failure results from void initiation, growth and coalescence.

• Our numerical approach coupled a Gurson material ( [4], [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]) behavior type with extended finite element method ( [START_REF] Moes | Belytschko: A finite element method for crack growth without remeshing[END_REF]) to reproduce the crack propagation in a 2D notched structure submitted to Mode I loading.

• The crack growth criterion is based on the energy stored in a patch located near the crack tip. This criterion based on a patch allows to attenuate the mesh size dependence observed when dealing with softening materials. In a first approach the crack propagation direction is deduced from the plastic strain field near the crack tip. The crack length is defined using the definition introduced in reference ( [START_REF] Kanninen | Advanced Fracture Mechanics[END_REF]) in adequation with our approach.

• Experimental tests in modes I and II would be very useful to define more accurately the quantities used in our approach.

[1-3] The Authors would like to aknowledge the financial support of DGA (French Ministry of

Figure 1 :

 1 Figure 1: Principle of designation of enriched nodes

Figure 2 :

 2 Figure 2: History of the average equivalent stress with and without singular functions. Case of a notched plate submitted to Mode I loading

Figure 3 :

 3 Figure 3: Gauss Points evolution in finite elements during a crack growth

Figure 4 :

 4 Figure 4: Shape of the patch use for evaluation of the crack growth criteria

Figure 5 :

 5 Figure 5: Numerical test of a notched plate under tension loading

Figure 6 :Figure 7 :

 67 Figure 6: Comparison of the evolution of Wpatch during the test for an angle of 0 • for a patch with a radius of 6.25mm

Figure 8 :

 8 Figure 8: a) Evolution of reaction force with and without propagation and b) Evolution of Wpatch during crack growth

Figure 9 :

 9 Figure 9: Porosity field at time 0.004s

  [

	p		p	2
	22 p	-sign( p 12 ) 8 +	22 p	)	(38)
	12		12
	U(t)
		0.1m
	0.1m
	crack (0.02m)

Defence) as well as the fruitful discussions with Prof. M. Salaun from Toulouse University, ISAE/ICA, France.


Kaili Xie

Clement De Loubens

Frédéric Dubreuil

Deniz Z Gunes

Marc Jaeger

Marc Leonetti

Clément De Loubens

Marc Léonetti

Interfacial rheological properties of self-assembling biopolymer microcapsules

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Microcapsules are drops surrounded by thin elastic shells that have the ability to protect and deliver substances in biological or industrial processes. They have a great potential for the development of new multi-functional materials as paints 1 but also to deliver nutrients, bacteria, DNA or drugs in the human body via the blood microcirculation 2,3 or the digestive tract [4][5][6] . Controlling the mechanical properties of microcapsules is of prime importance to provide their stability in various physical and chemical a Aix-Marseille Université, CNRS, Centrale Marseille, M2P2 UMR 7340, 13451, Marseille, France. b Aix-Marseille Université, CNRS, Centrale Marseille, IRPHE UMR 7342, 13384, Marseille, France. E-mail: leonetti@irphe.univ-mrs.fr c Université Grenoble Alpes, LRP, F-38000 Grenoble, France. d CNRS, LRP, F-38000 Grenoble, France. E-mail: clement.de-loubens@univ-grenoblealpes.fr e Université Grenoble Alpes, CERMAV, F-38000 Grenoble, France. f CNRS, CERMAV, F-38000 Grenoble, France. g Nestlé Research Center, Vers-chez-les-Blanc, PO Box 44, 1000 Lausanne 26, Switzerland. † Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here]. See DOI: 10.1039/b000000x/ environnement and to release encapsulated molecules 7 .

In this regard, layer-by-layer assembly of polyelectrolytes on dissolvable colloids is considered as a gold standard method, since it allows to control accurately the thickness of the membrane [8][9][10] and so its mechanical properties 11,12 . But the drawback of this multi-steps method is its length when a large number of layers has to be deposited to get high mechanical resistance. Recently, facile routes of fabrication of biopolymer microcapsules have been proposed that takes advantage of the self-assembling of polylectrolytes at water and oil interfaces [13][14][15][16] . As example, polymer -surfactant complex has been used to encapsulate Vitamine E in oil droplets 17 . A biocompatible solution has been obtained by the complexation of chitosan, a water-soluble selfhealing biopolymer 18 , with a counter-charged oil-soluble ionic surfactant, in order to form membranes at the interfaces of microdroplets 18 . Chitosan [START_REF] Peter | Chitin and Chitosan in Fungi[END_REF] acting as a kind of polysaccharide carrying positively charged groups (Figure 1) presents good solubility in water under acid condition (e.g., pH<4) with a great potential in drug delivery due to its biocompatibility and mucoadhesive porperties 6 . The surfactant phosphatidic fatty acid is a type of food molecule well-known in chocolate manufacture, as well as a component of cell membranes 18 , referred to as 'PFacidYN', exhibiting negative charges (Figure 1), which is insoluble in water but well soluble in vegetable oils. The membrane obtained with these complexing molecules has a multi-layer aspect with a disordered structure and its thickness grows by a diffusion-controlled mechanism 18 , hence can be tuned by varying concentrations of complexing molecules and complexation time. However the robustness of this method in terms of stability and control of the mechanical properties of the microcapsules is unknown while it is the main condition for its scale-up. Indeed the resulting mechanical properties of these kind of self-assembling biopolymer microcapsules have not been characterized in a large range of deformations and conditions of fabrication.

Mechanical properties of thin membrane of microcapsules are described generally by bi-dimensional non-linear viscoelastic constitutive models [START_REF] Barthes-Biesel | [END_REF][21][22] , as for fluid interfaces 23,[START_REF] Brenner | Interfacial transport processes and rheology[END_REF] . Rheology of complex plane fluid-fluid interfaces can be characterized by specific tools that respond to small surface stresses in isolating them from bulk stresses as the double-wall ring geometry rheometer [START_REF] Fuller | [END_REF] . In contrast, the linear elastic properties of solid membrane of microcapsules can be inferred from direct experiments on microcapsules 26 , for example osmotic pressure tests [27][28][29] , atomic force microscopy (AFM) 11,30 or the analysis of wrinkles on deflated microcapsules 31 . Whereas advanced tools of interfacial rheology tend to dissociate the mechanical response of fluid-fluid interfaces from the fluid mechanics of the bulk phases [START_REF] Fuller | [END_REF] , Chang & Olbricht 32,33 have demonstrated the strength of hydrodynamics methods to deform millimetric capsules with simple flows and to determine surface elastic modulus. More recently, these methods have been extended to micrometric capsules [34][35][36][37][38] . By combination with advanced numerical simulations 21,22,39 , the nonlinear elastic constitutive laws in regime of large deformation 40 and the viscosity of the membrane 41 have been determined for microcapsules made from a biopolymer (HSA) that was physically cross-linked to form a thin membrane in hydrogel. It was shown that this kind of microcapsules can sustained up to 180% of stretching without plastic deformations or break-up. If microfluidic hydrodynamic stretching has been proved relevant and convenient to determine the mechanical properties of strongly cohesive membranes such as made with HSA, the method has never been applied to weakly cohesive membranes such as multilayers of polyelectrolytes. Mechanical properties of such microcapsules have mainly been characterized in small deformation regime by osmotic pressure tests, micromanipulation 12 or by applying stress / strain on plane surfaces 16,18 showing viscoelastic and self-healing properties of this kind of self-assembling membranes. Here we will demonstrate the advantage of hydrodynamic methods, namely the elongation flow, to characterize as well the linear as the non-linear interfacial rheological properties of polyelectrolytes chitosan / PFacidYN microcapsules. Moreover, we will show that the mechanical behaviour of this kind of polyelectrolyte microcapsules differ greatly from those of cross-linked biopolymer membranes.

The objective of this paper is to characterize the interfacial rheological properties of self-assembling biopolymer microcapsules that were obtained in one single facile step. The aim is to question their fabrication robustness and provide new insights of their non-linear rheological behavior. Firstly, suspensions of chitosan / PFacidYN microcapsules 18 were produced using microfluidics droplets generators with various concentrations of complexing molecules and size of droplets. We paid attention also to control the time of complexation of the two electrolytes at water / oil interface and the stability of the microcapsules in regard to osmotic shocks outside the fabrication system. Secondly, the linear and non-linear interfacial rheological properties of these microcapsules were characterized by stretching them under controlled hydrodynamic stresses in a milli-fluidic extensional flow chamber 38,42 . We show that surface shear elastic modulus depends on the composition and the complexation time. Moreover, the microfluidic fabrication method allowed us to get suspensions of monodisperse microcapsules with uniform mechanical properties. Finally chitosan / PFacidYN microcapsules show remarkable non-linear rheological behaviors as yield stress (i.e. solid / liquid transition) and break-up that can be tuned also via the formulation and are related to the thickness of the membrane measured by Atomic Force Microscopy (AFM). We conclude on the opportunities offered by theses results to 'engineer' in a facile route the mechanical properties of self-assembling biopolymer microcapsules for specific applications.

Materials and methods

Chitosan solutions

Chitosan (CAS number 9012-76-4, Sigma-Aldrich) with medium molecular weight and 75-85% deacetylation (Figure 1) was used for the fabrication of polyelectrolytes microcapsules. Chitosan powder was dissolved in Millipore water (resistivity > 18 MΩ cm) by adjusting the pH with hydrochloric acid (1 mol/L) at 3.0. The chitosan solution was then filtered to remove undissolved particles with syringe-filters Minisart (pore size of 5.0 µm). The density was measured with the densimeter DMA 4500 M (Anton Paar) and adjusted to 1.09 ± 0.01 g/mL by addition of pure glycerol at 22 • C (98% of purity, CAS number 56-81-5, VWR) to avoid sedimentation or creaming in silicon oil Wacker AP1000 that was kindly offered by IMCD group. All Chitosan solutions used in this investigation were stored at 4 • C. Before microcapsules fabrication, the pH value of the mixture of chitosan and glycerol solution was checked and eventually re-adjusted at 3.00 ± 0.02. Table 1 in supplementary materials summarizes the different solutions of chitosan / glycerol used for microcapsules fabrication. Their viscosities range from 13.9 to 794 mPa•s at 22 • C.

Solutions of phosphatidic fatty acids in oil

The surfactant used to precipitate capsule membrane is ammonium phosphatidic fatty acid (PFacidYN, Figure 1). It is comprised in the commercial lecithin YN (Palsgaard 4448, food-grade, E442) that was provided in kind by Palsgaard. In mass, the phosphatidic acids are 55% w/w, neutral triglycerides 40% and ammonium salts 5%, see Minifie [START_REF] Minifie | Cocoa and Confectionary: Science and Technology[END_REF] and Gunes et al. 18 for details. PFacidYN was dissolved in rapeseed oil (CAS number 8002-13-9, Sigma-Aldrich) at different concentrations (Table 1 in supplementary materials). Undissolved particles were removed by centrifugation. The viscosity of these solutions was about 65-70 mPa•s according to the concentration of PFacidYN.

Fabrication of monodisperse microcapsules

Monodisperse chitosan / PFacidYN microcapsules were fabricated in a home-made microfludic T-junction (Figure 1) that produced water-in-oil droplets [START_REF] Glawdel | [END_REF] . The chip was manufactured by a 3-D printer (Formlabs Inc., Form 1+) and made of a hydrophobic resin (Clear, FLGPCL03, Formlabs). A round glass capillary (from CM Scientific Ltd.) of 200 µm internal diameter and 65 µm thick was inserted into the side channel branch where the aqueous phase was pumped into. The intersection area in the flow direction of the continuous phase (rapeseed oil) was narrowed to strengthen the shear forces and to generate aqueous droplets. The two immiscible phases were driven by pressure and flow rate controllers (Fluigent). Consequently, we were able to produce monodisperse suspensions of spherical microcapsules: the standard deviation of the size distribution was within 4% (Figure 1). Various microcapsules radii have been obtained by adjusting the flow rates of both phases with the same set-up: from 80 to 200 µm.

To obtain uniform suspensions of microcapsules in term of interfacial rheological properties, it is necessary to control the time of membrane formation due to the complexation of chitosan with PFacidYN. Firstly, microcapsules were pre-formed in the microfluidic chip with the oil phase containing a very low concentration of PFacidYN (< 0.01% w/w) to stabilise the droplet and avoid coalescence, whereas the aqueous phase contained the desired concentration of chitosan /glycerol mixture. These capsules were collected for 5 minutes. During this time, the membrane formation was very limited as the concentration of PFacidYN was low. Secondly, the microcapsules were transferred into an oil phase with 0.066 to 2.75 % w/w PFacidYN. In this way, the time let for the complexation of chitosan with PFacidYN, t comp , can be controlled easily and varied from 2 to 25 min. Thirdly the reaction was stopped by diluting the oil phase with a large quantity of mineral oil (CAS number 8042-47-5, Sigma Aldrich) and cyclohexane (CAS number 110-82-7, assay 99.5%, Sigma Aldrich) with a ratio 4:1 v/v. Microcapsules were let to sediment within a few minutes and the supernatant organic phase was removed. The microcapsules were dispersed in silicone oil Wacker AP1000 to be streched in elongation flow. Density of the internal fluid (chitosan / glycerol mixture) was similar to the density of the silicon oil to avoid sedimentation or creaming. The viscosity of silicone oil AP1000 was 1.18 Pa•s at 22 • C for shear rates up to 500 s -1 .

Atomic Force Microscopy

The membrane thickness of dried microcapsules was measured by Atomic Force Microscopy (AFM). Firstly, microcapsules deposited on a silicon wafer were left to dry for a week at 35 • C. Microcapsules were so collapsed on themselves during this process (Figure 2-left). The topology of the surface of collapsed microcapsules was done by AFM. AFM measurements were carried out with a Dimension Icon (Brucker) and a Pico Plus (Molecular Imaging) equipped with silicon probes having a nominal spring constant of 0.2 N.m -1 and tip radius below 10 nm (All in One AFM probes from Budget Sensors) using standard procedure for contact mode in fluids. Data were treated and analysed using Gwyddion software. Topographic variations were obtained after baseline correction via profile measurements on various microcapsules 11 . The membrane thickness was estimated by measuring the height of the collapsed capsule flat regions (Figure 2) and ranged from 40 to 900 nm according to the formulation used to fabricate the microcapsules.

Elongation flow chamber

The interfacial rheological properties of microcapsules were determined by stretching individual microcapsules in an elongation flow chamber. 38,40,42 Planar elongation flow was generated in a cross-like channel made of PMMA with 1x1 mm 2 square-sections (Figure 3). The fluid was injected into the chamber by a glass syringe mounted on a home-made syringe pump driven by a PI actuator M235-52S. We paid attention to avoid any deformable The rate of elongation in the center of the cross was measured by Particle Tracking Velocimetry (PTV), see de Loubens et al. 38 for details. A suspension of 1.6 mg/mL spherical ploystyrene particles (from Bangs Lab) with average diameter of 10 µm was seeded in silicone oil AP1000. The rate of elongation was found to be uniform in a region of interest (ROI) of 250 µm around the stagnation point where the flow is zero and unstable. On the center-plane (z = 0), the velocity field v is given by v = ε(xxyy) where ε is the rate of elongation (in s -1 ), x and y the cartesian axis with origin the stagnation point (Figure 3). All capsule deformation analysis should be conducted in the ROI. Within the range of flow rate 1 to 400 µL/s, a linear variation of the elongation rate was found. In addition, the dependence of ε on the depth of the channel (z direction) was characterized. As expected, the flow had a Poiseuille-like profile. For each measurement, the z position of the microcapsule was measured to determine accurately the elongation rate using the law ε = f (z) that was determined experimentally 38 . To prevent any effect of the channel walls, only capsules close to the center (i.e. distance from the center below ± 100 µm) were studied. For radii of 200 µm, the capsule sees only a variations of 2% of the rate of elongation in these conditions. The projection of a deformed capsule on the xy plane was assumed to be an ellipse with the major and minor axis length L and S, respectively. In fact, this is the first shape mode which is excited by the elongation flow whatever the nature of the soft particle at the center. We defined the Taylor parameter by

D = L -S L + S (1)
For low values of elongation rate (ε < 30 s -1 ), capsules were trapped one by one at the stagnation point of the elongation flow. Two containers at the ambient pressure were plugged to the outlets of the channel, which allowed to stabilize the microcapsule at the stagnation point. Then we suddenly applied the flow and measured the deformation as function of the time (Figure 3, ε = 9.98 s -1 ). The steady-state deformation D ∞ was extracted from the plateau of the curve. Alternatively, for high values of rates of elongation (ε > 30 s -1 ), the microcapsule was sucked up in one branch of the channels up to 5 mm from the stagnation point and the flow was applied. The capsule came so from one of the two inlet branches of the channel, then passed by the stagnation point and went back in a perpendicular branch (outlet). In this way, it was possible to obtain a saturated dynamics and large deformations of the microcapsules (Figure 3, ε = 41.72 and 68.69 s -1 ) before the capsule moved out of the ROI. We checked that the two ways were equivalent for moderate rates of elongation 40 . The rate of elongation ε was varied from 0.5 to 300 s -1 .

The hydrodynamic stress σ (in Pa) applied on the microcapsule in the center of the cross-like channel is

σ = η ε ( 2
)
where η is the viscosity of the external phase, i.e. the silicone oil AP1000. For microcapsules with thin membranes, their interfacial rheological properties can be modelled by surface constitutive laws without bending resistance. Various constitutive laws 21,39 describe the elastic deformation of microcapsules in small and large deformations: generalized Hooke law, Skalak law, neo-Hookean law. In the limit of small deformations, all these laws reduced to the linear Hooke law with the same surface shear elastic modulus G s (in N/m) that is related to the stable deforma- Each symbol corresponds to a microcapsule. The surface shear elastic modulus G s is determined from the slope of the curve for each microcapsule (dotted lines, Equation 3). White symbols: 0.

G s = 25 6 σ R D ∞ (3)
Note that G s can be correlated with the three dimensional Young modulus E, the Poisson coefficient ν and the thickness of membrane h by

G s = Eh 2(1 + ν) = Eh 3 (4)
if the membrane is assumed to be made of an uniform incompressible elastic material (ν=0.5). The determination of the surface viscous modulus of the membrane is out of the scope of this paper. This parameter is difficult to determine by experiments on microcapsules as we need to track the movement of the membrane in shear flow and deduce the viscous modulus from inverse numerical simulations 22,41 .

Each microcapsule was tested for at least six different values of elongation rates covering a range with at least a factor two. Nineteen different formulations of microcapsules were tested (Table 2 in supplementary materials) with 4-7 microcapsules by batch. The radius was also varied from 80 to 200 µm. Finally, more than 120 microcapsules in total were analyzed in this paper.

Results and discussion

Elastic properties of microcapsules

Our experimental set-up allowed us to apply different hydrodynamic stress σ on a given microcapsule and to analyse the evolution of its deformation. The steady-state Taylor parameter D ∞ varied linearly with increasing of stress σ within the small deformation regime (Figure 4), as expected by asymptotic theories 45 ( D ∞ < 0.1). In this regime, the shape of the microcapsule recovered to initially spherical status after the flow stopped. Different microcapsules with the same formulation exhibited similar behavior, which means that capsules fabricated with our present method have the same elastic properties of membrane for a given formulation. The slope of the curve D ∞ versus σ provided us a constitutive rheological parameter, the surface shear elastic modulus G s , that has been determined for each capsule (Equation 3). According to the formulation, G s ranged from 25 mN/m to 1.02 N/m (Table 2 in supplementary materials). Each value of G s was the mean of 4 to 7 microcapsules. We estimated that the measurement uncertainty was within ±6% with our set-up. For a given formulation, the standard deviation was also of this order. Our microfluidic procedure of fabrication and collection of suspensions of microcapsules allowed us to obtain not only monodisperse microcapsules (Figure 1), but also uniform elastic properties of the membrane of microcapsules. It was also noteworthy that elastic properties of chitosan / PFacidYN microcapsules were independent of their radius. For example, human serum albumin (HSA) microcapsules fabricated from batch emulsification of two immiscible phases are polydisperse in size [START_REF] Edwards-Levy | Serum Albumin: Structure, Functions and Health Impact[END_REF] and present variations of their elasticity of several orders of magnitude 38,47 (Figure 5). For microcapsules prepared with 0.25%w/v chitosan solution, 0.36%w/w PFacidYN and a complexation time of 2 minutes, the mean G s was about 0.11 N/m (mean of 60 microcapsules) with a standard deviation of ± 0.02 N/m for radii ranging from 80 to 200 µm. On the contrary, a factor 3.5 in size gave a factor 2000 in G s for HSA microcapsules. These strong differences can be explained by the kinetics of membrane formation. For chitosan / PFacidYN capsules, the thickness of the membrane formation scaled with the time at power 1/2 on time scale of several hours in agreement with a growth by diffusion-controlled self-assembly of surfactant and/or chitosan chains 18,48 , while most of HSA molecules have been adsorbed at the water / oil interface and have reacted in 30 min 49 . Consequently, the maximum concentration of HSA molecules at the interface should be equal to R[HSA]/3 in first approximation.

The maximum surface concentration in HSA molecules thus increased both with the capsule size and with the initial bulk HSA concentration 38,47 . This result enables us to fabricate microcapsules with uniform interfacial rheological properties that do not depend on the size. A first way to control the interfacial rheological properties of the membrane of chitosan / PFacidYN microcapsules is to vary the time of complexation between both molecules. For fixed concentration of chitosan (0.25%w/v) and PFacidYN (0.36%w/w), the surface shear elastic modulus G s increased linearly with the complexation time t comp , G s 0.04t comp (Figure 6). For a complexation time of 2 and 25 minutes, G s was 0.111 ± 0.006 N/m and 1.023 ± 0.049 N/m, respectively, either a factor 10. This result is correlated to the fact that the membrane thickness increases also with the square-root of the time 18 .

Another way to control the interfacial rheological properties of these microcapsules is also to change the total concentration in complexing molecules. Similarly to the complexation time, G s increased linearly with the concentration of surfactant PFacidYN (Figure 7). For 0.25%w/v chitosan and 2 minutes complexation of membrane, G s ranged from 0.025 ± 0.001 to 0.842 ± 0.058 N/m when the concentration of PFacidYN increased from 0.066 to 2.75%w/w. In this case, the linear relationship between G s and the concentration of PFacidYN was G s 0.31[PFacidY N]. For lower concentration of surfactant, microcapsules coalesced easily which were prevented to explore in this study. This linear evolution of G s was correlated to an increase of the membrane thickness that has been measured by AFM. Measuring the topographies of collapsed microcapsules (Figure 2) allowed us to determine the magnitudes of the membrane thickness for different concentrations of PFacidYN. We found that the thickness of the membrane increased from about 40-50 nm up to 700-950 nm for concentrations of PFacidYN of 0.066 and 2.75% w/w respectively. At intermediate concentrations of surfactant (0.66 -1.32% w/w PFacidYN), the membrane thickness was about 100-170 nm. The uncertainty is relatively large as this kind of measurement is technically difficult due to poor adhesion of microcapsules on the substrate and the difficulty to determine the number of layers for each step. The uncertainty on the number of layers can be reduced if a region is scanned at the interface between the microcapsule and the substrate where there is more chance to have only two layers of membrane. Several measurements have been done on the same microcapsule, but these values have to be considered rather as an order of magnitude than an highly accurate measurement. However, the variations between the samples were sufficiently large to conclude. At this stage, it is not reasonable to conclude more precisely on the trend of the evolution of the membrane thickness with others formulation and process parameters. We can also deduce from Equation 4that the Young modulus of the membranes increased slightly from 1.6 to 3.2 MPa contrarily to G s that was multiplied by a factor 33. This order of magnitude of Young modulus is equivalent to the one measured by AFM on HSA microcapsules 38 . It means that the variation of G s can be mainly attributed to the variation of the membrane thickness and the membrane structure is almost homogeneous.

The surfacic shear elastic modulus G s increases also with the concentration of chitosan up to reach a plateau value for a critical concentration, about 0.5 % w/v (Figure 7). Above this critical concentration, there was so a sufficient quantity of chitosan molecules available to saturate the interface and the quantity of surfactant available for the complexation is the limiting factor, since the plateau value of G s depends on its concentration (Figure 7). Finally, we show that the elasticity of self-assembling biopolymer microcapsules can be easily controlled by the microfluidics process of fabrication in a large way (factor 10). Moreover, this process allowed us to get monodisperse suspensions of microcapsules with uniform elasticity that can be kept stable outside the system. Tuning the various properties (e.g. size, morphology, adhesion, mechanics) of microcapsules combined with costefficient route of farbrication is still a challenge and should open numerous innovations as in biotechnology 3 . For example, Gunes et al. 18 fabricated polyelectrolytes microcapsules with a flow-focusing microchip where the complexing molecules were dispersed separately in the oil and water phase according to their solubility, to avoid layer-by-layer construction. A similar system has been reproduced by Kaufman et al. 50 and in the present paper. More sophisticated microfluidic devices have been developed to estimate the elastic modulus simultaneously with the production of microcapsules [51][52][53] . These set-ups are useful to determine the kinetic process at the early stage of membrane formation 52 , but they are limited by the testing time-scale of a few seconds, and cannot give access to the final properties of the microcapsules whereas the typical time to form a stable membrane is in order of a few minutes. Moreover, the stability of the suspension outside the generation serum is not guaranteed 51 , i.e. the resulting osmotic shock may have too large consequences. In contrast to the previous workers, we have deliberately chosen to uncouple the formation of monodisperse droplets with surfactant/chitosan adsorption to prevent coalescence from the real formation of the membrane, to be able to control the complexation time (from 2 to 25 min), the growth of the membrane and therefore its elastic properties. The uncoupling fabrication of microcapsules with their mechanical characterization allowed us to guarantee the stability of the suspensions and also to measure a constitutive rheological parameter of the membrane, the surface elastic modulus, instead of a deformation under specific flow conditions making difficult to determine constitutive rheological parameters without numerical simulations 51 . Furthermore, high shear stress can be generated with the elongation flow chamber that have access to the non-linear interfacial rheological properties of chitosan/PFacidYN microcapsules.

Break-up, yield stress and plasticity of microcapsules

In the regime of large deformation (typically D ∞ >0.1), the nonlinear mechanical behavior of chitosan / PFacidYN microcapsules can also be tuned by the formulation (concentration, time) to obtain microcapsules that can present a yield stress associated to plastic deformations and/or that burst easily under an applied hydrodynamic stress. Few examples of these behaviors are presented in a quantitative way.

Chitosan / PFacidYN microcapsules with the thinnest membranes (< 100 nm) can be easily bursted under hydrodynamic stress with our elongation flow chamber (Figure 8-left). Beyond the small deformation regime, the evolution of the steady-state deformation tended to slow-down with the stress. The maximum steady deformation D ∞ was about 0.3 (L ∞ /R=1.3) for a stress of 36 Pa without membrane fracture (0.25% w/v chitosan, 0.066% w/w PFacidYN and 2 min complexation time). In this example, the deformation was elastic as microcapsules recovered their initial shapes at rest and stress larger than 40∼50 Pa were required to burst the microcapsule. The critical stress increased also with G s .

The process of break-up of chitosan / PFacidYN microcapsules occurred in the time scale of milliseconds (Figure 8-right). Whereas the long axis L increased with time, a pinch appeared at the center of the ellipsoidal microcapsule. The shape was unstationary and evolved to a dumbbell with a stronger central pinching due to the right and left lobs which stretched the capsule, reminiscent of recent numerical simulations 54 . Finally, a thin filament connected the two opposite lobs that moved away (see movie in supplementary materials). Contrary to Newtonian droplets, we have not been able in this study to obtain long and quite large threads. The elastic response of membrane protected efficiently the integrity of large lobs and prevented a leak of liquid from the lobs to the thread. This thin filament was unstable as expected from the Plateau-Rayleigh instability driven by the surface tension of fluid interfaces [55][56][57] which generated a pattern of pearls along the filament. These observations question the existence of a possible solid to liquid transition of the membrane at these rate of deformations that could be studied on plane interfaces by double-wall ring geometries [START_REF] Fuller | [END_REF] . Moreover, in elongation flow, droplets with surfactant also exhibited tip-streaming 58 which is not the case here at least in the range of studied parameters.

For the thickest Chitosan / PFacidYN membranes (> 200 nm), burst of microcapsules was not observed in the range of hydrodynamic stress available with our set-up. However, plasticity and relaxation times of few minutes were observed. The Taylor parameter at rest of microcapsules has been measured after applying different levels of hydrodynamic stress in interval of time between stopping the flow and measurement of 0.5∼1 min (Figure 9). In regime of small deformations there were no residual deformation that can be detected by our visualization method (D res < 0.008). However, the capsules that were exposed to larger stress did not return immediately to their spherical shape at rest. The residual deformation D res increased non-linearly with the value of applied stress on the microcapsules up to 0.1. The critical stress that separates elastic and plastic deformations can be quantified and varied with the formulation. From these observations, we cannot deduce the yield stress of the membrane. In fact, the interface being closed, the state of stress in the membrane is inhomogeneous when the capsule is deformed. Previous numerical simulations have shown that tension is maximal at the tips of the capsule 59 . We can so deduce that plastic deformations are probably localized near the tips of the deformed capsule. After the capsule has been highly deformed and the flow was stopped, we observed two phases of shape relaxation (Figure 10). In only few seconds, the deformation was much lower than the steady-state deformation D ∞ . On time scales of few minutes, the shape relaxed slowly. Microcapsules with large steady-state deformations had larger residual deformations, and seemed to require longer relaxation times (up to 40 min). Due to experimental limitations (e.g., sedimentation of capsules and working length of objective lens), longer relaxation time was not exerted.

In summary, extensional flow set-up allows us to exert large hydrodynamic stresses on the chitosan / PFacidYN microcapsules, for which deviation from the linear elasticity was observed, contrarily to most common methods like AFM or osmotic pressure test 12 . This method give insights for the first time on the nonlinear mechanical behavior of a kind polyelectrolyte microcapsules. Microcapsules presented remarkably non-linear behaviors, such as permanent deformation and membrane break-up. As the state of tension on the membrane is inhomogeneous, experiments on plane interfaces with double-wall ring geometries should be relevant to determine solid-liquid transitions of these kind of membranes 18,[START_REF] Fuller | [END_REF] . Chitosan / PFacidYN microcapsules stand out also of HSA capsules that can be stretched up to 180% of elongation without plasticity or break-up. 40 The difference between the two systems is the nature of bonds : electrostatic versus cross-linking. Chang and Olbricht 32,33 probably first conducted the study of the residual deformation and membrane break-up of nylon capsules with a size order of 2∼4 mm under an elongation flow and Couette flow. Other studies handled with the polysiloxane capsules 37,60 , with which permanent deformation and membrane fracture were also observed. For chitosan/PFacidYN microcapsules, probably it is also interesting to see that these different behaviors can be obtained the same complexing molecules by varying the membrane thickness. It has also been shown by Kaufman et al. 61 that complexing the polyelectrolytes with nanoparticles is a lever to control the plastic deformations for membranes of a few micrometers thick.

Conclusion

Self-assembling and self-healing 16,18 biopolymer chitosan/PFacidYN microcapsules exhibit a variety of linear and non-linear interfacial rheological properties. In the framework of interfacial properties of capsule membranes, this study generalizes the relevance of the extensional flow method 21,40,42 to probe the full scale of deformations of different kinds of microcapsules, whatever the cohesive strength of membrane materials from strong with crosslinking of HSA for example 38 to weak as shown in this study on polyelectrolyte multilayers.

If microfluidics is a way to produce monodisperse suspensions of droplets with a small standard deviation, a crucial step to make uniform suspensions, it does not prevent from having a larger deviation in interfacial rheological properties of capsule membranes. We have shown that combining the microfluidic tools with a cost-effective physico-chemical process involving a singlestep interfacial complexation of polyelectrolytes allows to generate various uniform suspensions of microcapsules. Indeed, it is remarkable to observe that with the same complexing molecules, the interfacial rheological properties can be tuned and various degrees of elasticity, yield stress, plasticity and break-up can be obtained. The surface elastic modulus ranged from 25 mN/m to 1 N/m and was correlated to the thickness of the membrane (40 -950 nm). This wealth of mechanical behavior can so be mainly attributed to the variation of the membrane thickness, as the order of magnitude of the Young modulus was almost constant whatever the concentration of complexing molecules was (1.6 -3 MPa, Eq. 4). As measurements have been realized in steadystate, the interfacial viscosity of these microcapsules was not determined. Experiments in shear flows are more suitable to determine this parameter by tracking the motion of the membrane with microparticles 41 , but the complexity is also to choose the accurate visco-elastic model 22,62 . Beyond the interests of microencapsulation in various domains of engineering, these objects can also be used as mastered model systems of red blood cells for micro-hydrodynamics experiments 63 . Finally, this paper is a step forward to 'engineer' the interfacial rheology of microcapsules for specific purposes.

Fig. 1

 1 Fig. 1 Fabrication of monodisperse chitosan / PFacidYN microcapsules. Left: Formation of aqueous droplets of chitosan dispersed in an oil phase containing the surfactant PFacidYN in a T-junction. The shell is formed around an aqueous phase droplet by electrostatic interactions between chitosan and PFacidYN molecules. Right: Optical micrograph of monodisperse chitosan / PFacidYN microcapsules.

Fig. 2

 2 Fig. 2 Measurement of membrane thickness by AFM. Left: Optical microscopy of a dried microcapsule on a silicium wafer. The microcapsule has collapsed on itself. Middle: Topographic image of a 30 x 30 µm 2 area of a dried microcapsule. The membrane thickness is estimated by measuring the height of the collapsed capsule flat regions. Right: Three dimensional reconstruction of the topographic image.

Fig. 3

 3 Fig. 3 Stretching of microcapsules. Left: Individual microcapsules are stretched in the center of the cross (1 x 1 mm 2 ) by an elongation flow. The deformation is characterized by the small S and large L axis of the microcapsule, and by the Taylor parameter D = (L -S)/(L + S). Right: Temporal evolution of the Taylor parameter D as a function of dimensionless time for different rates of elongation ε. The steady-state deformation D ∞ is determined as the plateau value of each curve. See Movie 1 in Supplementary Materials.

  pipes and any air bubbles to minimize the transient time when the flow was switched on. The visualization was achieved with an inverted microscope Olympus IX-71 with a 20× magnification. A high speed video camera (Photron Fastcam SA3) enables us to acquire up to 5000 frames per second. The image resolution was 0.851 µm by pixels. The temperature was controlled by a digital thermocouple with a precision of 0.1 • C. All the experiments were carried out at 22 ± 0.5 • C.

Fig. 4

 4 Fig.4 The steady-state Taylor parameter D ∞ increases linearly with the hydrodynamic stress σ for two different formulations of microcapsules. Each symbol corresponds to a microcapsule. The surface shear elastic modulus G s is determined from the slope of the curve for each microcapsule (dotted lines, Equation3). White symbols: 0.25 % w/v chitosan, 0.66 % PFacidYN, 2 min complexation time, 102 ± 1 µm radius. Black symbols : 0.11 % w/v chitosan, 0.36 % PFacidYN, 2 min complexation time, 102 ± 1 µm radius.

Fig. 5 GFig. 6

 56 Fig.5G s as a function of the radius R: 0.25%w/v chitosan, 0.36%w/w PFacidYN, 2 min complexation time. Comparison is made with HSA microcapsules (data from de Loubens et al.38 , 10% w/w HSA). G s is almost constant for chitosan / PFacidYN microcapsules, whereas it varies of several orders for HSA microcapsules. Each symbol is the mean of 4 microcapsules for chitosan /PFacidYN microcapsules and represents an individual microcapsule for HSA. Error bars are of the size of the symbol.

Fig. 7

 7 Fig. 7 Surface shear elastic modulus G s as a function of concentrion of complexing molecules. Top: G s increases linearly with the concentration of PFacidYN. 0.25%w/v chitosan, 2 min complexation time, 102 ± 1 µm radius. Bottom: G s saturates when the concentration of chitosan reaches a threshold of 0.5 %w/v. 0.36%w/w and 0.66%w/w PFacidYN, 2 min complexation time, 102 ± 1 µm radius. Error bars are the standard deviations.

D∞Fig. 8

 8 Fig. 8 Break-up of a chitosan / PFacidYN microcapsules in an elongation flow. Left: The steady-state Taylor parameter D ∞ increases non-linearly with the hydrodynamic stress σ above values of 0.1. Images show the evolution of the steady-state shape. Each symbol corresponds to a microcapsule. Right: When the break-up threshold is exceed (σ = 50.56 Pa), the microcapsule is stretched up to the formation of two lobs connected by a thin filament (t = 146 ms) that is then destabilized in droplets (t = 190 ms). The center of the flow chamber is indicated by +. 0.25%w/v chitosan, 0.066%w/w PFacidYN, 2 min complexation time, 102 1 radius. See Movie 2 in supplementary materials.

Fig. 9

 9 Fig. 9 Plastic deformations of chitosan / PFacidYN microcapsules. Top: The residual Taylor parameter D res increases with the hydrodynamic stress σ just after stopping the elongation flow (within 0.5 min). Bottom: The yield stress demarcates elastic and plastic deformation σ p is defined by D >0.01 and depends on the formulation of microcapsules. Lines are guides for the eyes. 0.36%w/w PFacidYN, 2 min complexation time, 102 ± 1 µm radius.

Fig. 10

 10 Fig. 10 Relaxation of the residual deformation. Top: A microcapsule that is spherical at rest is deformed up to D ∞ of 0.54. Just after the flow (1 min), a residual deformation D res remains at rest and decreases slowly with the time. Bottom: Taylor parameter D as a function of the time. The first and second parts of the curves correspond to the elongationnal flow and to the period at rest respectively. Microcapsules do not recover their initial shape (i.e. D res > 0.01 on the timescale of the observation (about 45 min). Circles: 0.11%w/v chitosan, 0.36%w/w PFacidYN, 2 min complexation time, G s = 0.065 N/m, σ = 153 Pa. Squares: 0.25%w/v chitosan, 0.36%w/w PFacidYN, 6 min complexation, G s = 0.230 N/m, σ = 744 Pa.

  25 % w/v chitosan, 0.66 % PFacidYN, 2 min complexation time, 102 ± 1 µm radius. Black symbols : 0.11 % w/v chitosan, 0.36 % PFacidYN, 2 min complexation time, 102 ± 1 µm radius.

tion D ∞ under elongation flow by

45

[1][2][3][4][5][6][7][8][9][10]

1-10 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

This work has benefited from financial support from the ANR Polytransflow (ANR-13-BS09-0015), from Labex MEC (ANR-11-LABX-0092) and from CNES. Also, K. Xie acknowledges the China Scholarship Council (CSC 201504490090) to fund his PhD scholarship. F. Dubreuil ackowledges the Nanobio-ICMG plateform (FR 2607) for the AFM facilities. The authors thank Cécile Soubeyrand from IMCD group for the gift of Wacker AP 1000 and Laurent Eymeoud from Plasgaard for the gift of lecithinYN (Palsgaard 4448, E442).

Supplementary Materials

The following files are available free of charge.

• Movie 1: Deformation of a chitosan / PfacidYN microcapsule in the regime of large deformations.

• Movie 2: Break-up process of a chitosan / PfacidYN microcapsule.

• Tables 1 &2: Formulations of microcapsules studied and surface shear elastic modulus G s .


Johannes Wolf

Konstantinos Nikolakopoulos

Patrice Longère

Jean-Philippe Crété

Jean-Marc Cadou

XFEM-CZM combination for the numerical treatment of ductile fracture

Keywords: Damage Mechanics, Fracture Mechanics, Computational Methods, XFEM, Cohesive Zone Model

teaching and research institutions in France or abroad, or from public or private research centers.

ABSTRACT

We are here interested in the numerical prediction using commercial finite element computation code of the response of large engineering structures made of ductile materials when submitted to accidental overloads. The challenge consists in reproducing within a unified methodology based on the finite element method (FEM) the different steps of failure as observed in the materials under consideration, namely diffuse microvoiding-induced damage, strain/damage localization then macro-cracking, if any.

The behavior of the material in presence of ductile damage is governed by the microporous plasticity Gurson model in the context of standard FEM and the kinematic consequences of the macro-crack propagation throughout the meshing are described in the context of the eXtended FEM (XFEM), see Crété et al. (2014). A key point for the prediction of ductile fracture is the numerical treatment of the continuous-to-discontinuous transition, intermediate stage, i.e. localization within a narrow band. To describe it, a model of 'cohesive strong discontinuity' and the XFEM are combined. Criteria for the (more or less) diffuse-damage-tolocalization and localization-to-macro-cracking transitions are as well proposed.

The feasibility of the methodology has been shown, see Wolf et al. (2017), via the implementation of user element subroutines into the commercial finite element computation code Abaqus, and boundary value problems are solved.


Konstantinos Nikolakopoulos

Jean-Philippe Crété

Patrice Longère

Patrice Longère Xfem

XFEM-CZM combination for the numerical treatment of ductile fracture

Keywords: Endommagement, Rupture ductile, Méthodes numériques, XFEM, Modèle de zone cohésive

Introduction

Nous nous intéressons ici au traitement numérique tri-dimensionnel de la ruine progressive de structures de grandes dimensions au moyen de la méthode des éléments finis. Cette ruine s'opère typiquement suivant 3 étapes : endommagement diffus, localisation de la déformation, fissuration. Dans la présente approche, la première étape d'endommagement diffus est régie par le modèle de GTN [START_REF] Tvergaard | Needleman Analysis of cupcone fracture in a round tensile bar[END_REF] et l'ouverture ultime de la fissure est décrite via l'utilisation de la méthode des éléments finis étendus (XFEM), comme proposée pour le cas des métaux ductiles par Crété et al. [START_REF] Crété | Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM[END_REF] . Pour la phase intermédiaire de localisation de la déformation, la méthode retenue est l'association d'un segment cohésif (CZM) et d'une discontinuité forte dans l'élément fini 3D dans le cadre de la formulation XFEM [START_REF] Wolf | Numerical modeling of strain localization in engineering ductile materials combining cohesive models and X-FEM[END_REF]. Enfin, la mise en oeuvre de la méthodologie est développée en tant que routine utilisateur (UEL) dans le code commercial de calculs par éléments finis ABAQUS.

Endommagement diffus et fissuration 2.1 Endommagement diffus GTN

Afin de décrire le comportement du matériau en présence de micro-vides, le modèle de Gurson [START_REF] Gurson | Continuum theory of ductile rupture by nucleation and growth : Part i -yield criteria and flow rules for porous ductile media[END_REF] est utilisé, tel que modifié par Tveergard et Needleman [START_REF] Tvergaard | Needleman Analysis of cupcone fracture in a round tensile bar[END_REF], dans le cadre de la méthode des éléments finis standard FEM.

Φ = σ eq σ y 2 + 2q 1 f cosh - 3 2 q 2 p m σ y -(q 3 f ) 2 -1 (1)
avec σ eq la contrainte équivalente, p m le pression, σ y la contrainte d'écoulement, f la fraction volumique des vides et q 1 , q 2 , q 3 des constantes.

Fissuration XFEM

Dans cette étude, la méthode utilisée pour décrire les conséquences cinématiques de la présence et de la propagation d'une fissure dans la structure est la méthode des éléments finis étendus (XFEM) [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF]. Les noeuds de chaque élément fini fissuré sont enrichis avec des degrés de liberté supplémentaires. Il a été montré que dans le cas d'un matériau ductile élasto-plastique fortement non linéaire la partie asymptotique peut être omise avec une perte de précision négligeable et entraînant une réduction importante du nombre de degrés de la liberté du problème [START_REF] Crété | Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM[END_REF]. En utilisant la variante de la base décalée [START_REF] Zi | New crack-tip elements for xfem and applications to cohesive cracks[END_REF] le champ de déplacement enrichi peut s'écrire comme suit :

u(x) = i∈I N i α i + j∈J (H(x) -H j )N j β j , où H j = -0.5 ou + 0.5 (2)
3 Phase intermédiaire de localisation

Combinaison XFEM-CZM

La combinaison d'une discontinuité forte avec une bande cohésive dans le cadre de la XFEM [START_REF] Wolf | Numerical modeling of strain localization in engineering ductile materials combining cohesive models and X-FEM[END_REF] peut se présenter sous la forme suivante, sans tenir compte des efforts extérieures :

Ω B T D t BdΩ Ω B T D t B * dΩ Ω B* T D t BdΩ Ω B* T D t B*dΩ + 4 Γ D N T C loc NdΓ da db = - Ω B T σdΩ Ω B* T σdΩ + Γ D N T T loc dΓ
(3) avec σ le vecteur contenant les contraintes, T loc le vecteur des efforts cohésifs, B la matrice des dérivées spatiales des fonctions d'interpolation, B* la matrice des dérivées spatiales des fonctions (H(x) -H j )N j et Γ D la bande de discontinuité. L'indice loc fait référence au repère local de la bande de localisation. Le schéma d'intégration choisi provisoirement dans cette étude est celui proposé dans [START_REF] Belytschko | A finite element with embedded localization zones[END_REF], dans lequel les 8 points de Gauss de l'élément sont tous utilisés deux fois, une pour la partie supérieure et une pour la partie inférieure de l'élément coupé. Il s'agit d'une approche simple et plus légère que celle proposée par Elguedj et al. [START_REF] Elguedj | Combescure Appropriate extended functions for x-fem simulation of plastic fracture mechanics[END_REF].

Modèle de CZM adopté

La loi de cohésion développée est une loi linéaire avec un plateau initial. Elle est caractérisé par une variable d'endommagement scalaire D, prenant des valeurs comprises entre 0 si zone saine et 1 si zone totalement endommagée. T loc (∆) est activée quand le critère de transition endommagement diffus-localisation est vérifiée et suit provisoirement une loi linéairement décroissante. ∆ est un déplacement équivalent [START_REF] Camanho | Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials[END_REF] qui prend en compte la composante normale δ 1 et les composantes de cisaillement δ 2,3 du saut de déplacement, exprimées dans le repère local de la bande de localisation.

∆ = << δ 1 >> 2 +δ 2 2 + δ 2 3 , où << δ 1 >>= max(0, δ 1 ) (4)
La traction initiale au moment de la localisation de l'élément est calculée via l'utilisation d'un tenseur de contrainte moyenné σ sur tous les points Gauss :

T 0 =< σ n > loc (5)
Le vecteur des efforts cohésifs s'écrit alors :

T loc = (1 -D) * T 0 , où D =      0 si ∆ ≤ κ 0 ∆-κ 0 κc si κ 0 < ∆ < κ 0 + κ c 1 sinon (6)
avec κ 0 le seuil de déplacement équivalent, à partir duquel l'adoucissement commence et κ c la valeur critique au-delà de laquelle nous considérons l'élément comme rompu. En dérivant, on obtient alors la matrice tangente de la loi cohésive :

C loc (i, j) = ∂T loc (i) ∂δ j = - ∂D ∂δ j * T 0 (i) = - 1 κ c δ j ∆ * T 0 (i), C loc (i, 1) = 0 si δ 1 < 0 (7)
Sur la Fig. 1, nous pouvons voir l'impact de la valeur κ c sur l'adoucissement. Cette méthode conduit à une loi cohésive qui nécessite de fixer une seule constante, à savoir κ c , car T 0 est déterminée automatiquement par le critère de transition endommagement diffus-localisation.

Transitions

-ENDOMMAGEMENT DIFFUS VERS LOCALISATION : Afin de détecter le plan de localisation dans le matériau, on réalise une analyse de bifurcation [START_REF] Rice | The localization of plastic deformation[END_REF]. Avec l'utilisation de l'opérateur tangent élastoplastique D t et l'hypothèse des petites déformations, le critère consiste en la recherche du plan n pour lequel :

det n T D t n = 0 (8)
On utilise une condition de porosité critique pour l'activation de la bande cohésive dans l'élément. Celle-ci n'est insérée que lorsque la porosité moyenne de tous les points de Gauss de l'élément dépasse une valeur critique f c .

-LOCALISATION VERS FISSURATION : Le critère pour la transition de la localisation à la fissuration de l'élément est atteinte lorsque la variable d'endommagement D dans (6) atteint la valeur 1.

Application

Le modèle numérique présenté dans les chapitres précédents a été implanté en tant qu'élément fini utilisateur UEL dans le code commercial par éléments finis Abaqus. Pour vérifier la méthodologie, une éprouvette en traction composée d'éléments hexaèdriques à intégration complète (C3D8) dans le cadre des petites déformations est utilisée (voir Fig. 2). Celle-ci a un centre plus étroit, ce qui oblige la localisation à apparaître au centre de l'éprouvette. Un déplacement est appliqué comme présenté sur la Fig. 2. Un maillage fin et un grossier ont été utilisés. Il en résulte que la partie adoucissante de la loi cohésive pour le maillage grossier a la même pente mais elle commence tardivement.

Conclusions

Il a été démontré que la méthodologie présentée peut décrire toutes les différentes étapes du phénomène de rupture ductile en 3D. Une légère dépendance au maillage a également été constatée, une réalité qui devra être traitée ultérieurement dans cette étude. Enfin, dans le maillage 1, nous pouvons également assister à un petit saut de la force de réaction qui est dû à un mauvais calcul du vecteur des forces internes au moment de la localisation.

Ces faiblesses sont en phase d'amélioration. En outre, notre objectif est de passer de l'utilisation de la théorie des petites déformations à la théorie des grandes déformations. Notre intérêt pour cette application provient du fait que la déformation plastique peut impliquer des changements importants dans la géométrie de la structure qui ne peuvent pas être négligés. Ce travail a été cofinancé par la DGA.

Références

FIGURE 1 -

 1 FIGURE 1 -Impact de κ c sur l'adoucissement

TABLE 1

 1 La raison pour laquelle une si grande valeur de κ 0 a été choisie tient au fait que cela contribue à la convergence. Les éléments de la zone de localisation entrent dans la phase de localisation à deux moments différents. D'une part, cela indique que notre modèle peut décrire les fissures qui se propagent en plusieurs incréments. D'autre part, cela conduit à ce que certains éléments localisés se trouvent dans la phase de plateau, tandis que d'autres sont déjà entrés dans la phase d' adoucissement. Cela peut entraîner des difficultés de convergence et l'utilisation d'une valeur de plateau plus grande permet à plus d'éléments d'entrer simultanément dans l'adoucissement. Dans les deux cas, l'analyse de bifurcation nous donne le même plan de localisation. Nous pouvons remarquer que le maillage grossier atteint la localisation à un moment légèrement différent de celui du maillage fin. Ceci est dû au fait que la porosité augmente avec un taux légèrement plus élevé dans le maillage fin, ce qui est également FIGURE 2 -(à gauche) Eprouvette, (à droite) affichage de l'effort de réaction pour les maillages 1 et 2 visible à partir de la courbe de force de réaction du maillage 1 légèrement plus basse que celle du maillage 2.

				-Données des deux calculs
	κ 0 κ c	f c	Maillage 1 (Fin)	Maillage 2 (Grossier)
	3	4 0.03 11704 éléments (4 dans l'épaisseur) 1336 éléments (2 dans l'épaisseur)