Pierre-Jean Tisserand

Héloïse Rostagni

Cédric Giry

Thi Thanh

Huyen Nguyen

Rodrigue Desmorat

Frédéric Ragueneau

An orthotropic damage model with internal sliding and friction for masonry-like material

Keywords: quasi-brittle materials, orthotropic damage, internal sliding, masonry

Quasi-brittle materials like masonry exhibit mechanical properties and develop nonlinearities that are mainly driven by their joints. In the present work, a Continuum Damage Mechanics point of view is considered to describe the macroscopic behavior of masonry. Based on damage models developed for ceramic matrix composite materials, a fixed directional damage approach is proposed. From this formulation, unilateral effect, as well as internal sliding and friction coupled with damage, are introduced. Numerical examples of the response of the model for different loading cases involving cyclic and non-proportional loadings are carried out and compared to experimental results.

Conventions

The following equations are not written with Einstein's Convention. The stress and strain vector notations are defined using the Kelvin notation and Betcherew's basis:

σ σ σ =         σ 1 σ 2 σ 3 σ 4 σ 5 σ 6         =         σ 11 σ 22 σ 33 √ 2σ 23 √ 2σ 13 √ 2σ 12         ε ε ε =         ε 1 ε 2 ε 3 ε 4 ε 5 ε 6         =         ε 11 ε 22 ε 33 √ 2ε 23 √ 2ε 13 √ 2ε 12        
The elastic orthotropic compliance tensor is defined using the Kelvin notation and Betcherew's basis:

S =         S 11

Introduction

Running bond masonry is a widely used structural material in bridge and tunnel heritage, and a large amount of these are still in service. They are facing the regular increase in railway traffic, and natural hazards are the main unpredictable causes of their failure. In order to better prevent the impact of natural hazards, overall seismic ones, on masonry railway infrastructure, robust material models are needed. A running bond masonry is constituted of quasi-rectangular-shaped blocks linked by mortar joints. For a high aspect ratio of the blocks, it leads to an orthotropic elastic behavior of the masonry at the macroscale [START_REF] Pande | Equivalent elastic moduli for brick masonry[END_REF]. This quasi-brittle material develops cracks mainly in mortar joints ( [START_REF] Dhanasekar | The failure of brick masonry under biaxial stresses[END_REF], [START_REF] Page | The biaxial compressive strength of brick masonry[END_REF]). For cyclic uniaxial loading, one can observe unilateral effects due to the closure of the developed cracks. Under cyclic shear loading, hysteretic dissipation develops ( [START_REF] Silva | Experimental assessment of in-plane behaviour of three-leaf stone masonry walls[END_REF], [START_REF] Anthoine | Shear-compression testing and analysis of brick masonry walls[END_REF], [START_REF] Salmanpour | Displacement capacity of contemporary unreinforced masonry walls: An experimental study[END_REF], [START_REF] Naraine | Behavior of brick masonry under cyclic compressive loading[END_REF], [START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF]). This last phenomenon can be linked to internal sliding in the cracks generated at the mesoscale as for concrete [START_REF] Ragueneau | Damage model for concrete-like materials coupling cracking and friction, contribution towards structural damping: first uniaxial applications[END_REF].

As developed for instance in [START_REF] Ghiassi | Numerical Modeling of Masonry and Historical Structures: From Theory to Application[END_REF], [START_REF] Roca | Structural analysis of masonry historical constructions. classical and advanced approaches[END_REF], [START_REF] D'altri | Modeling strategies for the computational analysis of unreinforced masonry structures: review and classification[END_REF], different scales and different modeling strategies can be considered to describe the mechanical behavior of masonry. Among the first modeling strategy to evaluate the failure of masonry structure, one can find the work of Heyman [START_REF] Heyman | The stone skeleton[END_REF]. He applied the limit analysis method to masonry, considering simple hypotheses (i.e., no tensile strength and infinite compressive strength for masonry, and no possible sliding between blocks) in order to evaluate the strength of masonry arch. In this framework, several developments have been made later on to improve the representativeness of the masonry behavior like the derivation of the ultimate strength of masonry from a homogenization approach [START_REF] De Buhan | A homogenization approach to the ultimate strength of brick masonry[END_REF] or the definition of numerical tools to evaluate the limit states of a block assembly ( [START_REF] Livesley | Limit analysis of structures formed from rigid blocks[END_REF], [START_REF] Baggio | Limit analysis for no-tension and frictional three-dimensional discrete systems[END_REF]). In order to describe the failure of masonry structures and to describe explicitly the cracking in masonry, the discrete element method developed originally by Cundall [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] has shown its efficiency (e.g. [START_REF] Lemos | Discrete element modeling of masonry structures[END_REF]). In this general framework, some specific developments have been made, like, for instance, the description of interactions between bodies through the non-smooth contact dynamics method ( [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF], [START_REF] Jean | The non-smooth contact dynamics method[END_REF]).

For evaluating the response of masonry structures considering moderate loadings without provoking complete failure, homogeneous description of the masonry at the macroscale within a continuous description has shown its efficiency (see for instance [START_REF] Degli Abbati | Seismic assessment of interacting structural units in complex historic masonry constructions by nonlinear static analyses[END_REF], [START_REF] Limoge-Schraen | Toward a large-scale seismic assessment method for heritage building: vulnerability of masonry baroque churches[END_REF]). Different modeling strategies can be found to describe masonry as a homogeneous continuous medium. Due to its low resistance under tensile stress, a category of models proposes the perfectly no-tension material hypothesis [START_REF] Piero | Constitutive equation and compatibility of the external loads for linear elastic masonry-like materials[END_REF]. It allows to determine the maximum capacity of a structure; nevertheless, it needs dedicated numerical strategies [START_REF] Angelillo | A finite element approach to the study of no-tension structures[END_REF] and does not allow to investigate the softening response of a structure under seismic loadings. The nonlinear behavior of masonry can also be described by classical continuum theory as smeared cracks (e.g. [START_REF] Rots | Smeared and discrete representations of localized fracture[END_REF]), plasticity (e.g. [START_REF] Lourenço | Multisurface interface model for analysis of masonry structures[END_REF]), damage (e.g. [START_REF] Papa | A unilateral damage model for masonry based on a homogenisation procedure[END_REF], [START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF]) or the coupling of damage and plasticity (e.g. [START_REF] Addessi | A 2d cosserat finite element based on a damage-plastic model for brittle materials[END_REF] or [START_REF] Lee | Plastic-damage model for cyclic loading of concrete structures[END_REF] with application to masonry structures in [START_REF] Degli Abbati | Seismic assessment of interacting structural units in complex historic masonry constructions by nonlinear static analyses[END_REF]). As shown in the last example, numerous models in this continuum mechanics framework are inspired by models initially developed for concrete. Unless there is a robust and efficient description of nonlinear phenomena in quasi-brittle materials, these models generally miss the description of the anisotropic nature of the nonlinearities in masonry media. To assess this anisotropic nonlinear behavior of cracking phenomenon, an orthotropic damage continuum framework has been proposed by some authors ( [START_REF] Berto | An orthotropic damage model for masonry structures[END_REF], [START_REF] Pelà | A localized mapped damage model for orthotropic materials[END_REF]). In [START_REF] Berto | An orthotropic damage model for masonry structures[END_REF], the authors build a damage tensor for in-plane problems based on a combination of scalar damage variables associated with the normal and parallel directions of the bed joints. The influence of these scalar damage variables on shear response is derived from equilibrium at the macroscale using the effective stress tensor and friction angles. In [START_REF] Pelà | A localized mapped damage model for orthotropic materials[END_REF], using mapped tensors, the authors build an orthotropic damage model from an isotropic one. Both models are able to describe the progressive orthotropic degradation of masonry; nevertheless, no mechanism is introduced to describe the hysteretic loops observed for cyclic shear loadings.

Finally, the anisotropic behavior of the masonry can also be obtained through lower scale informed multiscale approaches. An efficient approach in this framework that includes various nonlinear local phenomena like friction without inducing large computational time is proposed by [START_REF] Marfia | Multiscale damage contact-friction model for periodic masonry walls[END_REF] using the Transformation Field Analysis (TFA) [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF].

The development of crack families associated with specific material direction can also be observed in ceramic matrix composites. To describe this nonlinear behavior and the effect of crack families on the response at the macroscale, the concept of fabric tensor is considered (e.g. [START_REF] Marcin | Development of a macroscopic damage model for woven ceramic matrix composites[END_REF]). The consistent decomposition of the degradation in a thermodynamical framework allows to develop a rigorous and numerically robust model. Furthermore, the coupling between damage and other mechanisms like plasticity for sliding in cracks (e.g. [START_REF] Desmorat | Continuum damage mechanics for hysteresis and fatigue of quasi-brittle materials and structures[END_REF]) can be easily performed in this formalism.

The present paper proposes to develop a novel approach with respect to existing continuous models for masonry by considering a description of the orthotropic damage through a decomposition per direction and associated crack family. Furthermore, coupling between damage and friction is introduced to describe the hysteric loop observed during cyclic shear loading. These aspects are of main importance when dealing with the dynamic response of structures subjected to medium earthquakes. This paper is divided into three parts. The first one concerns the theoretical formulation (section 2).

Different assumptions about the formulation and the constitutive equations are developed. First of all, the mechanical framework and the coupling between elasticity, damage and plasticity are explained. Then, the evolution of internal variables during loading/unloading conditions is presented. The formulation of an orthotropic damage model and the introduction of plasticity to describe friction under cyclic shear loading represent a new contribution to masonry modeling. The second part deals with the numerical implementation and algorithm (section 3). The last part highlights numerical validation and applications.

Theoretical formulation

The running bond masonry is described as a homogeneous material whose orthotropic behavior is defined by the orientation of the joints. Figure 1.a) gives a schematic description of the masonry with the orthotropic directions (1, 2, 3). The constitutive law associated with the elasticity is given by,

ε ε ε = S 0 : σ σ σ (1)
With σ σ σ, the Cauchy's stress tensor, ε ε ε, the strain tensor and S 0 the elastic compliance tensor. In the natural frame of the masonry, the tensor S 0 using the Kelvin notation is given in equation [START_REF] Angelillo | A finite element approach to the study of no-tension structures[END_REF].

S 0 =                  1 E 1 -ν 12 E 1 -ν 13 E 1 0 0 0 -ν 12 E 1 1 E 2 -ν 23 E 2 0 0 0 -ν 13 E 1 -ν 23 E 2 1 E 3 0 0 0 0 0 0 1 2.G 23 0 0 0 0 0 0 1 2.G 13 0 0 0 0 0 0 1 2.G 12                  (2)
E i , ν ij and G ij are elastic material parameters that can be identified by homogenization techniques ( [START_REF] Anthoine | Derivation of the in-plane elastic characteristics of masonry through homogenization theory[END_REF], [START_REF] Massart | Mesoscopic modeling of failure in brick masonry accounting for three-dimensional effects[END_REF]).

The measure of crack density proposed in equation (3) by [START_REF] Bristow | Microcracks, and the static and dynamic elastic constants of annealed and heavily coldworked metals[END_REF] is used to quantify the damage d i associated with each crack family linked to the orthotropic directions,

d i ∼ = 1 V n l 3 n,i (3)
where V is the representative volume and l n is the radius of the n th circular crack of the fracture plane of normal vector i.

Equation ( 3) is used in [START_REF] Kachanov | Elastic solids with many cracks and related problems[END_REF] to formulate cracks effect on elastic moduli. A simple formulation of the effective damage compliance tensor S ef f considering non-interacting orthogonal cracks is given in equation [START_REF] Anthoine | Shear-compression testing and analysis of brick masonry walls[END_REF].

S ef f = S 0 + i (d i A (i) : S 0 ) (4)
Where d i is a damage variable evolving from 0 for undamaged material to +∞ for totally damaged material and associated with the crack family i. A (i) is a fabric tensor (a fourth-order tensor), giving the influence of d i on each component of the elastic compliance tensor S 0 .

The damage development in the masonry is mainly governed by mortar joints, leading to orthogonal crack patterns. A physical representation of this damage mechanism is a system of three orthogonal fracture planes whose normals are related to orthotropic elastic directions (i.e., one family of cracks per axis of the natural frame of the masonry). In order to describe the damage mechanisms according to the formalism defined in equation ( 4), three damage variables (d 1 , d 2 , d 3 ) are introduced (Fig. 2). The effective compliance is defined by:

S ef f = S 0 + d 1 A (1) : S 0 + d 2 A (2) : S 0 + d 3 A (3) : S 0 (5)
with A (i) the fabric tensor associated with the direction i. The expressions of these fabric tensors are inspired by the models for composite materials (see for instance [START_REF] Marcin | Development of a macroscopic damage model for woven ceramic matrix composites[END_REF]). More particularly, these models introduce the concept of damage effect tensor [START_REF] Zheng | On damage effective stress and equivalence hypothesis[END_REF] to express the influence of damage on the effective compliance tensor. As a family of cracks influences directly the compliance in its associated direction, a coefficient equal to one is introduced in the fabric tensor for

A (i)
ii (i.e., the term associated with the normal direction of the considered crack family). Furthermore, a simple coupling is introduced by these tensors along the plane related to the normal direction of the family of cracks in order to account for the effect of damage on shear response. The fabric tensors are given in equations ( 6), [START_REF] Berto | An orthotropic damage model for masonry structures[END_REF] and [START_REF] Bristow | Microcracks, and the static and dynamic elastic constants of annealed and heavily coldworked metals[END_REF].

A (1) =         1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A (1) 55 0 0 0 0 0 0 A (1) 66         (6)
A (2) =         0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 A (2) 44 0 0 0 0 0 0 0 0 0 0 0 0 0 A (2) 66         (7)
A (3) =         0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 A (3) 44 0 0 0 0 0 0 A (3) 55 0 0 0 0 0 0 0         (8)
The expressions of the fabric tensors are chosen to describe an orthotropic damage behavior. A more complex anisotropic response like the one described for instance in [START_REF] Marfia | Multiscale damage contact-friction model for periodic masonry walls[END_REF] could be handled by introducing additional non-zero terms in the tensors A (i) .

The thermodynamic potential is Gibbs free enthalpy density,

ρΨ ⋆ = 1 2 σ σ σ : S ef f : σ σ σ = 1 2 σ σ σ : S 0 : σ σ σ + 1 2 3 i=1 d i [σ σ σ : (A (i) : S 0 ) : σ σ σ] (9)

Unilateral effects

In order to take into account the unilateral effect of damage between traction and compression by crack closure, the normal stress in the natural basis (i.e., σ i for the direction i) is written by splitting the stress as in [START_REF] Ladevèze | Sur une théorie de l'endommagement anisotrope[END_REF] (equation ( 10)).

i

∈ [1, 2, 3] σ i =< σ i > + + < σ i > - (10)
< x > + (= max(0, x)) takes the positive part of x and < x > -(= min(0, x)) the negative part of x. By using the equations ( 4), ( 6), ( 7), ( 8) and [START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF], the thermodynamic potential (equation ( 9)) can be developed as,

2ρΨ ⋆ = < σ 1 > 2 + E1 (1 + d 1 ) + < σ 1 > 2 - E1 + < σ 2 > 2 + E2 (1 + d 2 ) + < σ 2 > 2 - E2 + (11)
< σ 3 > 2 + E3 (1 + d 3 ) + < σ 3 > 2 - E3 -2 ν 21 σ 1 σ 2 E1 -2 ν 31 σ 1 σ 3 E1 -2 ν 32 σ 3 σ 2 E2 + σ 2 4 2G 23 (1 + d 2 A (2) 44 + d 3 A (3) 44 ) + σ 2 5 2G 13 (1 + d 1 A (1) 55 + d 3 A (3) 55 ) + σ 2 6 2G 12 (1 + d 1 A (1) 66 + d 2 A
(2) 66 )

Internal sliding and friction

Quasi-brittle materials exhibit hysteretic dissipation and permanent strains due to frictional sliding at the surface of the cracks. The hysteretic loops can be characterized thanks to cyclic loadings (see for instance [START_REF] Mazars | The unilateral behaviour of damaged concrete[END_REF] for concrete). For masonry material, this phenomenon is mainly observed along with the shear stress component. In the present work, the internal sliding and friction are introduced only for the shear stress component and are coupled to damage according to the proposal of [START_REF] Desmorat | Continuum damage mechanics for hysteresis and fatigue of quasi-brittle materials and structures[END_REF] for the isotropic case.

The frictional sliding phenomenon is introduced with sliding strain components denoted ε π 4 , ε π 5 and ε π 6

and their associated stresses σ π 4 , σ π 5 and σ π 6 . The Gibbs free enthalpy can be decomposed into two parts: ρΨ ⋆ = ρΨ ⋆ e -w s with ρΨ ⋆ e the elastic energy density and w s the energy density stored by internal sliding and friction.

For the sake of clarity, the elastic energy density is decomposed also into two parts : ρΨ ) ⋆ the elastic energy density associated with the normal stress components. The full Gibbs free enthalpy density is written by,

ρΨ ⋆ = ρ(Ψ normal e ) ⋆ + ρ(Ψ shear e ) ⋆ -w f riction s (12) 2ρ(Ψ normal e ) ⋆ = (P N : σ σ σ) : S 0 : (P N : σ σ σ) + 3 i=1 d i < σ i > 2 + E i ( 13
)
where P N is an operator of projection that selects the diagonal part of σ σ σ in the natural basis (i.e., the normal components in the orthotropic basis).

In order to introduce the internal sliding, ρ(Ψ shear e ) ⋆ is written in accordance with the proposal of [START_REF] Desmorat | Continuum damage mechanics for hysteresis and fatigue of quasi-brittle materials and structures[END_REF] for each shear component:

ρ(Ψ shear e ) ⋆ = 1 2 (σ 6 -σ π 6 ) 2 (1 + A (1) 66 d 1 + A (2) 66 d 2 ) 2G 12 + (σ π 6 ) 2 2G 12 g 12 (d 1 , d 2 ) + 1 2 (σ 5 -σ π 5 ) 2 (1 + A (1) 55 d 1 + A (3) 55 d 3 ) 2G 13 + (σ π 5 ) 2 2G 13 g 13 (d 1 , d 3 ) + 1 2 (σ 4 -σ π 4 ) 2 (1 + A (2) 44 d 2 + A (3) 44 d 3 ) 2G 23 + (σ π 4 ) 2 2G 23 g 23 (d 2 , d 3 ) (14)
The function g ij (d i , d j ) is assumed to give the same properties to σ π ij as in [START_REF] Desmorat | Continuum damage mechanics for hysteresis and fatigue of quasi-brittle materials and structures[END_REF], when the material is undamaged,

d i = d j = 0, σ π ij = 0 =⇒ g(d i , d j ) = 0 and when it is totally damaged, g(d i , d j ) = 1.
It leads to the formulation in equations ( 15) to [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF].

g 12 (d 1 , d 2 ) = A (1) 66 d 1 + A (2) 66 d 2 1 + A (1) 66 d 1 + A (2) 66 d 2 (15) g 13 (d 1 , d 3 ) = A (1) 55 d 1 + A (3) 55 d 3 1 + A (1) 55 d 1 + A (3) 55 d 3 (16) g 23 (d 2 , d 3 ) = A (2) 44 d 2 + A (3) 44 d 3 1 + A (2) 44 d 2 + A (3) 44 d 3 (17)
The energy density stored by internal sliding and friction is also modified in order to accommodate the kinematic hardening plasticity like friction stress (equation ( 18)).

w f riction s = 1 2 b 4 α 2 4 + b 5 α 2 5 + b 6 α 2 6 ( 18
)
b i are material parameters, their identification by experiments is needed; α i are kinematic hardening variables.

From the state potential ρΨ ⋆ , the state laws can be derived for total strains, shear sliding strains, kinematic back stresses and damage release rate (see Appendix A for the full derivation of the equations),

ε i = ρ ∂Ψ ⋆ ∂σ i , ε π k = -ρ ∂Ψ ⋆ ∂σ π k , X k = -ρ ∂Ψ ⋆ ∂α k , Y j = ρ ∂Ψ ⋆ ∂d j (19)
2.3. Evolution and loading/unloading conditions

Damage

The main phenomenon to initiate damage is linked to the tension. The extension of damage criteria for each damage parameter has to take also into consideration the shear. Different formulations for the extension indicator can be formulated. The one chosen for this model is a normal extension and shear strain combination. As three orthogonal damaging planes have been exhibited, the variable governing the damage d i is denoted εi and can be formulated from a projection on these planes.

The variable εi can be calculated for each damage variable using the equations ( 20) to [START_REF] Jean | The non-smooth contact dynamics method[END_REF].

ε1 = < ε 1 > 2 + + 1 2 β 12 ε 2 6 + 1 2 β 13 ε 2 5 (20) ε2 = < ε 2 > 2 + + 1 2 β 21 ε 2 6 + 1 2 β 23 ε 2 4 (21) ε3 = < ε 3 > 2 + + 1 2 β 31 ε 2 5 + 1 2 β 32 ε 2 4 ( 22
)
Where β ij and β ik are introduced to express the influence of shear strain on the damage. These material parameters have to be identified with experiments. The equivalent strain variable εi can be used as a quantity of interest for a non-local method of regularization as presented in section 3.4.

The threshold function is given in equation ( 23),

f i = χ i (ε i -K i ) -ln (1 + d i ) K i εi ≤ 0 ( 23
)
where K i is the initial extension threshold for undamaged material and χ i is the damage sensibility, a parameter controlling post-peak damage brittleness. These parameters can be identified on uniaxial tensile tests for each direction.

The loading conditions f i = 0 and ḟi = 0 drive the explicit evolution law for d i (equation ( 24)).

d i = εmax i K i e χi(ε max i -Ki) -1, εmax i = max τ ≤t εmax i (τ ) (24)
Finally, the damage evolution of the three damage variables is described with three independent explicit scalar expressions.

Internal sliding and friction

Criterion functions f π 4 , f π 5 and f π 6 are defined as follows to govern the loading/unloading conditions in shear (equation [START_REF] Kachanov | Elastic solids with many cracks and related problems[END_REF]). The influence of the stress components associated with the normal of the fracture planes is introduced through a Mohr-Coulomb-like criterion.

f π 4 =|σ π 4 -X 4 | + µ 4 (< σ 2 > -+ < σ 3 > -) ≤ 0 (25) f π 5 =|σ π 5 -X 5 | + µ 5 (< σ 1 > -+ < σ 3 > -) ≤ 0 (26) f π 6 =|σ π 6 -X 6 | + µ 6 (< σ 1 > -+ < σ 2 > -) ≤ 0 (27)
Each µ i is a material parameter corresponding to a friction coefficient and can be identified by biaxial shear tests.

The evolution law derives from a dissipation potential F π through the normality rule in a non-associated framework. The internal sliding model is non associated regarding the expression based on [START_REF] Armstrong | A Mathematical Representation of the Multiaxial Bauschinger Effect[END_REF] and by considering the flow developed only in the shear direction (equation ( 28)).

F π i = |σ π i -X i | + 1 2 a i X 2 i ( 28
)
a i are material parameters that influence the shape of the hysteretic loops for shear cyclic loadings. The equation [START_REF] Lemaitre | Mécanique des matériaux solides-3eme édition[END_REF] shows the normality rule by introducing the internal sliding multiplier λπ ,

επ i = λπ i ∂F π i ∂σ π i = λπ i σ π i -X i |σ π i -X i | (29) αi = -λπ i ∂F π i ∂X i = λπ i σ π i -X i |σ π i -X i | -a i X i State variables V k Associated variables A k ε i σ i ε π i -σ π i α i X i d i -Y i

Variables and parameters summary:

The variables and material parameters of the model are summarized in Tables 1 and2 respectively. The parameters can be obtained from the literature, by experimental test or by virtual test on masonry samples. The response of a masonry sample for a tensile test along the direction i is mainly influenced by the parameters E i , ν ij , ν ik , K i and χ i . The elasticity parameters can be identified from the initial linear response. The parameters K i define the threshold for the damage initiation. Experimentally, they correspond to the strain level observed along the direction i for the initiation of the first crack. Finally, the parameters χ i is linked to the energy dissipated by cracking phenomena from damage initiation up to failure. As it is in practice complex to perform a direct tensile test on a masonry wall, this response is efficiently obtained thanks to virtual testing (e.g. [START_REF] Oliver-Leblond | A beam-particle model to identify constitutive laws for quasi-brittle materials under complex loading: From concrete to masonry[END_REF]). The response of a masonry sample for a shear test along l tangential component is mainly influenced by the parameters

G ij , β ij , A (i) ll , A (j)
ll , µ l , a l and b l . The elasticity parameters can be obtained from the initial linear response. The parameters a l , b l and µ l are linked to the size of a hysteretic loop at a fixed damage state. This last parameter introduces the effect of a confinement stress on the shear response. The parameters β ij weight the contribution of the shear strain components on the development of the damage variables. Finally, the parameters

A (i) ll and A (j)
ll define the influence of the damage variables on the shear response. As damage and friction are coupled, a pragmatic approach to identify these parameters can be done by comparing the response of the model with experimental data for a cyclic shear loading. Cycle at fixed damage state needs to be performed to identify a l and b l . Different levels of confinement need to be considered to evaluate the parameters µ l . After this first set of identification, the coefficients of coupling (β ij , A (i) ll , A (j) ll ) can be identified according to monotonic shear loading along different directions.

Condition for thermodynamic admissibility

The thermodynamic potential previously presented (equations ( 13), ( 14) and ( 18)) is convex as a sum of quadratic terms. The use of a non-standard framework for damage involves that the thermodynamic consistency is proven by the positivity of the dissipation D. According to [START_REF] Lemaitre | Mécanique des matériaux solides-3eme édition[END_REF], the dissipation D is expressed with state variables (V k in table 1) and associated variables (A k in table 1). Following [START_REF] Lubliner | On the thermodynamic foundations of non-linear solid mechanics[END_REF], the dissipation can also be derived from the Gibbs free enthalpy density (equation ( 30)).

D = ρ Ψ⋆ -ε ε ε • σ σ σ = k V k Ȧk (30)
D can be split according to the different dissipation mechanisms (equation ( 31)).

D = D damage + D sliding (31)
The dissipation due to sliding D sliding is build on the same framework as standard plasticity. As a consequence, this term is always positive by construction.

The damage dissipation D damage (equation ( 32)) is deduced from the table (Table 1) and the definition (equation ( 30)).

D damage = Y 1 ḋ1 + Y 2 ḋ2 + Y 3 ḋ3 (32)
The study of D damage exhibits some load paths which could involve negative dissipation (due to sliding) if no attention is paid. Assuming no normal stress, Y 1 is negative when

ε 2 6 ≤ (ε 6 -ε π 6 ) 2 or ε 2 5 ≤ (ε 5 -ε π 5 ) 2
, and ḋ1 ≥ 0. For Y 2 and Y 3 , the same statement can be written. This thermodynamic violation can be theoretically and numerically avoided including explicitly the condition (equation ( 33)) in the modeling.

Y i ḋi ≥ 0 ( 33
)
When Y i < 0, ḋi = 0 With these conditions, which are respected in this formulation, the intrinsic dissipation is positive (D ≥ 0)

for any loading, uniaxial or not, proportional or not. So the constitutive equations presented in this model respect the thermodynamic principles.

The initial no-completion of the dissipation positivity condition is due to the coupling between the damage and the internal sliding. When g ij (d i , d j ) is an increasing function, the physical behavior (i.e., the increase of shear internal sliding with the damage increasing) is well-represented, but the second principle of thermodynamic is violated if d i is allowed to grow at negative Y i . Nevertheless, when g ij (d i , d j ) is decreasing, the second principle is respected, but the sliding is no more relied to damage increasing. In addition, negativity of the dissipation has never been registered during all validation tests.

Numerical implementation

In order to implement the model in a finite element code (i.e., strain-driven computation), additional tools need to be introduced. At time τ + ∆τ , an input strain ε ε ε τ +∆τ = ε ε ε τ + ∆ε ε ε is applied to the model with ε ε ε τ the strain at time τ and ∆ε ε ε the strain increment. The strain is expressed in the natural coordinate axis of the masonry by using rotation operator R (equation 34), which links the strain and stress in the general basis (• ⋆ ) to the strain and the stress in the natural basis of masonry (•).

ε ε ε = Rε ε ε ⋆ R T ( 34
)
σ σ σ = Rσ σ σ ⋆ R T (35) R = R X R Y R Z (36)

Damage

For damage, an implicit numerical scheme is used with a direct solution in a closed-form with respect to the total strain level. The equivalent strain is determined for each damage direction according to equa-tions [START_REF] Marfia | Multiscale damage contact-friction model for periodic masonry walls[END_REF] to [START_REF] Mazars | The unilateral behaviour of damaged concrete[END_REF].

ετ+∆τ 1 = (< ε 1 > τ +∆τ + ) 2 + 1 2 β 12 (ε τ +∆τ 6 ) 2 + 1 2 β 13 (ε τ +∆τ 5 ) 2 (37)
ετ+∆τ 2 = (< ε 2 > τ +∆τ + ) 2 + 1 2 β 21 (ε τ +∆τ 6 ) 2 + 1 2 β 23 (ε τ +∆τ 4 ) 2 (38)
ετ+∆τ 3 = (< ε 3 > τ +∆τ + ) 2 + 1 2 β 31 (ε τ +∆τ 5 ) 2 + 1 2 β 32 (ε τ +∆τ 4 ) 2 (39)
The damage evolution law is implemented according to equation [START_REF] Milani | Homogenised limit analysis of masonry walls, part i: Failure surfaces[END_REF].

d τ +∆τ i = max τ ετ+∆τ i K i e χi(ε τ +∆τ i -ki) -1, d τ i (40)

Internal sliding and friction

For the internal sliding and friction, a return mapping algorithm [START_REF] Ortiz | An analysis of a new class of integration algorithms for elasto-plastic constitutive relations[END_REF] is used to determine the stress state and the evolution of the internal variables. A damage-elastic step will be proceeded to calculate (σ π 4 ) trial , (σ π 5 ) trial and (σ π 6 ) trial (equations ( 41) to ( 43)).

(σ π 4 ) trial τ +∆τ = (ε τ +∆τ 4 -(ε π 4 ) τ )2G 23 g 23 (d τ +∆τ 2 , d τ +∆τ 3 ) ( 41
) (σ π 5 ) trial τ +∆τ = (ε τ +∆τ 5 -(ε π 5 ) τ )2G 13 g 13 (d τ +∆τ 1 , d τ +∆τ 3 ) (42)
(σ π 6 ) trial τ +∆τ = (ε τ +∆τ 6 -(ε π 6 ) τ )2G 12 g 12 (d τ +∆τ 1 , d τ +∆τ 2 ) (43)
The iteration process in the context of return mapping uses a Newton-Raphson scheme to determine the solution of the nonlinear equations, following the Taylor's series expansion of sliding thresholds: f π 4 = 0, f π 5 = 0 and f π 6 = 0 (equation ( 44)).

f π i | k+1 = 0 ≈ f π i | k + ∂f π i ∂σ π i k ∆σ π i | k+1 + ∂f π i ∂X i k ∆X i | k+1 (44)
During the iteration process, the total strain no more evolves. According to the state laws and the normality rules, the variation of the stress friction is defined by equations ( 45) to [START_REF] Page | The biaxial compressive strength of brick masonry[END_REF].

∆σ π 4 | k+1 = -∆λ π 4 | k+1 (2G 23 g 23 (d 2 , d 3 )) ∂F π 4 ∂σ π 4 k (45) ∆σ π 5 | k+1 = -∆λ π 5 | k+1 (2G 13 g 13 (d 1 , d 3 )) ∂F π 5 ∂σ π 5 k (46) ∆σ π 6 | k+1 = -∆λ π 6 | k+1 (2G 12 g 12 (d 1 , d 2 )) ∂F π 6 ∂σ π 6 k (47)
According to the state laws and the normality rules, the variation of the kinematic back stress X is defined by equation 48.

∆X i | k+1 = -∆λ π i | k+1 b π i ∂F π i ∂X i k (48)
The expression of the internal sliding multiplier is defined by equation 49.

∆λ π 4 | k+1 = f π 4 | k ∂f π 4 ∂σ π 4 k (2G 23 g 23 (d 2 , d 3 )) ∂F π 4 ∂σ π 4 k + ∂f π 4 ∂X 4 k b π 4 ∂F π 4 ∂X 4 k (49) ∆λ π 5 | k+1 = f π 5 | k ∂f π 5 ∂σ π 5 k (2G 13 g 13 (d 1 , d 3 )) ∂F π 5 ∂σ π 5 k + ∂f π 5 ∂X 5 k b π 5 ∂F π 5 ∂X 5 k (50) ∆λ π 6 | k+1 = f π 6 | k ∂f π 6 ∂σ π 6 k (2G 12 g 12 (d 1 , d 2 )) ∂F π 6 ∂σ π 6 k + ∂f π 6 ∂X 6 k b π 6 ∂F π 6 ∂X 6 k (51) (52)
The iterative scheme illustrated in the algorithm 1 stops when

f π i | k+1 f π i |0
≤ 10 -8 . |f π i | 0 | corresponds to the absolute value of the threshold function at the first iteration.

Total stress calculus and unilateral effect

The computation of the global stress σ σ σ should be obtained by inverting the state's laws. But the chosen formulation is impossible to be analytically inverted because of the positive/negative part splitting of the stress.

A solution is to proceed with an iterative scheme. A prediction of the stress state is obtained with,

σ 1 = C 11 ε 1 + C 12 ε 2 + C 13 ε 3 (53)
σ 2 = C 12 ε 1 + C 22 ε 2 + C 23 ε 3 (54)
σ 3 = C 13 ε 1 + C 23 ε 2 + C 33 ε 3 (55)
Where C ij are components of the elastic stiffness tensor, which can be obtained with the components of the compliance tensor (i.e., C ij = S -1 ij ):

C 11 =
S 11 1+d1 E1 1 E1 1+d1 E1 1+d1 E1 1 E1 1 E1 1+d1 E1 1 E1 S 22 1+d2 E2 1+d2 E2 1 E2 1+d2 E2 1 E2 1+d2 E2 1 E2 1 E2 S 33 1+d3 E3 1+d3 E3 1+d3 E3 1 E3 1+d3 E3 1 E3 1 E3 1 E3
Table 3: Description of the different cases for S ij

Mesh dependency, damage localization and energy regularization

To avoid or limit mesh dependency in finite element method due to strain-softening, regularization techniques have to be considered. For the case-studies considered in this work a simple energetic regularization [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] is considered. For non-local methods, the non-local equivalent strains (built from local εi ) are good candidates for the definition of non-local quantities that drive the damage variable d i ( [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]). The surface energy dissipation associated with the whole development of the fracture process zone is defined by G f (equation ( 62)). For energetic regularization, the parameters of the damage evolution law are calibrated in order to dissipate the same amount of energy in the volume of a finite element as the surface energy G f .

G f = h i +∞ 0 σdε (62)
with h i the characteristic size of the element along the normal direction of the crack. In order to derive an analytical expression for the parameters associated with the damage evolution law, this implies that σ σ σ is integrable up to total failure. By considering uniaxial formulation for each damage variable d i , one can get a relationship between the brittleness coefficient χ i , G f and h i (equation ( 63)).

χ i = E i K i h i G fi -h i 1 2 E i (K i ) 2 (63)
This energetic regularization formulation is available for the mode I fracture. G fi is a material parameter which has to be identified on uniaxial tensile test. It should be pointed out that energetic regularization formulation suffers mesh bias (such as dependency to mesh orientation) as shown for instance in [START_REF] Jirásek | Evaluation of directional mesh bias in concrete fracture simulations using continuum damage models[END_REF] or [START_REF] Mosler | Embedded crack vs. smeared crack models: a comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias[END_REF]. The choice of this method here is mainly driven by its easiness of use and by the fact the computational time is not increased.

Response of the model at the local scale

This section presents local results allowing to bring out the robustness of the formulation according to numerical applications, with the physical representation of nonlinear mechanisms. One case study at the structural scale shows that the modelling is able to describe strong nonlinear behaviors for diverse loading states.

These tests are performed on one linear element (CUB8) in the finite element solver Cast3M (http://www-cast3m.cea.fr, [START_REF] Verpeaux | A modern approach of large computer codes for structural analysis[END_REF]). Table 4 summarizes the material parameters used for these three tests. This set of parameters has been deliberately chosen in order to strengthen the orthotropic effects. Furthermore, they are chosen to be representative of masonry behavior. For these local tests, the parameters associated with the influence of shear on damage (i.e., shear component of fabric tensor A (i) and coefficient β ij ) are taken equal to 0.5. This choice is motivated by the fact a smaller coupling with the shear component than the normal one can be expected for the fabric tensor. This has been observed for composites (both ceramic and organic composites [START_REF] Marcin | Development of a macroscopic damage model for woven ceramic matrix composites[END_REF]). For the coefficient β ij , it allows to introduce a larger value for shear strain with respect to normal strain that is needed to initiate damage. Three tests are carried out in order to illustrate the main nonlinear mechanisms introduced in the model. The first test is a unidirectional tensile/compressive test to provoke damage and unilateral effect. The second test is a cyclic shear test with no pre-compression to illustrate the coupling between damage and friction mechanisms. The third test is a cyclic shear test with 0.1MPa of pre-compression to show the model response for a non-proportional loading and highlight the influence of confinement on the shear response, particularly the friction mechanism. For this purpose, the damage variables D i are plotted. These variables progress from 0 for a virgin material to 1 for a fully damaged material. These variable are obtained according to equation (64).

D i = d i 1 + d i ( 64
)
Uniaxial cyclic tensile/compressive test. This test (Fig. 3) illustrates the damage evolution and the unilateral of the tensile strain ε 11 can be observed. Furthermore, due to the stress decomposition (i.e., equation ( 10)), while stress becomes negative, the original stiffness is recovered.

Cyclic shear test. This test (Fig. 4) illustrates the damage development with shear and the coupling with friction. One can observe a progressive degradation of the shear modulus due to damage. Furthermore, the stress confinement is observed with the width of the loop. The higher the stress confinement, the larger the hysteretic loop. Similar evolutions are observed with or without confinement as the total strain drives the damage variables.

Structural case study

In order to evaluate the model capacity to describe the masonry behavior at the structural scale, the response of a confined hollow wall test in shear [START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] and a confined solid wall in cyclic shear [START_REF] Anthoine | Shear-compression testing and analysis of brick masonry walls[END_REF] are studied. This case study is considered to provide a qualitative comparison with respect to experimental data regarding the global response and the main failure mechanisms.

. Initial calibration of the model

To get the order of magnitude of the material parameters influencing the shear response, a calibration is made on tests performed on walls [START_REF] Vermeltfoort | Deformation controlled tests in masonry shear walls, part 2[END_REF] with the same masonry constituents (i.e. Joosten solid clay bricks and mortar) as the one used for the hollow wall [START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] modeled in the next section. The walls have the following geometrical characteristics: 1 x 1 x 0.098 (m x m x m). This test is labelled J6D in [START_REF] Vermeltfoort | Deformation controlled tests in masonry shear walls, part 2[END_REF]. The same loading system as the one shown in figure 8 is used. A confinement of 0.120 N/mm 2 is applied along the vertical direction to the solid panel. It should be stressed out that this part is not dedicated to the definition of a general method to identify the parameters of the model. To identify in a proper manner the whole set of parameters, further specific developments need to be done, using for instance methods proposed for peculiar problems, like the identification of failure criteria (e.g. [START_REF] Mistler | In-plane and out-of-plane homogenisation of masonry[END_REF] or [START_REF] Milani | Homogenised limit analysis of masonry walls, part i: Failure surfaces[END_REF]) or for damage model (e.g. [START_REF] Luciano | Homogenization technique and damage model for old masonry material[END_REF]). As strong coupling between damage and friction as well as anisotropic interactions are introduced in the present model, numerical experiments with a mesoscale model need to be carried out to optimize the identification of the whole set of material parameters (see as an example [START_REF] Chisari | Identification of mesoscale model parameters for brick-masonry[END_REF] for a procedure to identify a mesoscale model from experimental data).

In the present feasibility study, the set of parameters of the model associated to shear has been identified thanks to a simple minimization of the difference between the experimental response and the numerical one. Table 5 gives the values of the parameters obtained from the identification. Figure 6 between the global experimental and numerical responses for monotonic shear loading.

Mesh sensitivity analysis

In order to analyse the sensitivity to the mesh size, a dedicated study is performed on the same structural test (i.e., confined wall under monotonic shear loading). Only the global response is compared as it has been shown and observed local regularization technique suffered mesh bias (see for instance [START_REF] Jirásek | Evaluation of directional mesh bias in concrete fracture simulations using continuum damage models[END_REF] or [START_REF] Mosler | Embedded crack vs. smeared crack models: a comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias[END_REF]). A structured mesh with N x N parallelepipeds is considered, where N is the number of elements along the width and the height of the wall. For the analysis, the following values of N are considered: 5, 10, 20, 40, 50 and 100 elements. It corresponds to an average mesh size of: 0.2m, 0.1m, 0.05m, 0.025m, 0.02m and 0.01m. Figure 7 illustrates the results obtained regarding the mesh sensitivity for the shear wall tested by [START_REF] Vermeltfoort | Deformation controlled tests in masonry shear walls, part 2[END_REF]. For mesh size smaller than 0.05 m, the response tends to show the same maximum force and, the first softening point. Nevertheless, for large values of damage, a part of the dissipation is produced by the friction in the model. As this part is not directly affected by the energetic regularization, small differences are observed regarding the residual force.

To tackle this problem and to get a full mesh objectivity of the results (i.e. damage path and global response), the authors suggest to use advanced non-local regularization methods introducing an internal length in the constitutive euqations ( [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], ...) with example of application for the masonry in [START_REF] Peerlings | A thermodynamically motivated implicit gradient damage framework and its application to brick masonry cracking[END_REF] and [START_REF] Toti | Nonlocal damage propagation in the dynamics of masonry elements[END_REF].

Confined hollow wall under monotonic shear loading

The same material parameters as for the test on the solid panel are used here, excepted for the initial threshold K 1 and K 2 and the damage brittleness G f1 and G f2 . Indeed these modifications were motivated by the fact a too brittle behavior was observed when using the same set of parameters as the one of the previous subsection. The new values are: K 1 = 2.5•10 -4 , K 1 = 0.9•10 -4 , G f1 = 450N/m and G f2 = 350N/m. The characteristic lengths h i of the element for each damage variable d i are evaluated prior to the computation considering the size of the element along the normal direction of the crack family d i . The main objective of this case study is to appreciate the ability of the constitutive model to describe the nonlinear response of masonry at the structural scale with complex loadings and geometry.

In figure 8, one can see the experimental set-up as well as the mechanical boundary conditions and the loading. Finite elements with linear interpolation functions are considered. The average mesh size for the computation is equal to 2cm. The results are analyzed up to a damaged state of the wall that corresponds to the full development of diagonal shear cracks, as illustrated in figure 9. After this stage, the description in a continuous framework tends to be no more relevant. Furthermore, at the local scale, the model response reaches the residual shear stress. The sum of local contributions leads to a residual shear force at the global scale and so to no more loss of strength. A way to describe the failure after this high damage level could be to introduce a transition to strong discontinuity like in [START_REF] Kakarla | Coupled continuous-discrete formulation based on microplane and strong discontinuity models for representing non-orthogonal intersecting cracks[END_REF]. 8: Boundary conditions and loading for the numerical test of the hollow wall in shear with a normal stress [START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF] Comparisons between experimental and numerical load/displacement curves are presented in figure 10. Regarding the intrinsic variability in the results illustrated by the difference between the two experimental results in the figure 10, the order of magnitude of the peak force is captured by every computation. The lelepiped mesh. More significant differences are observed for the other quantities like the displacement at the peak force or the global dissipated energy. A non-local formulation could improve these results, as discussed in the section 5.1.2). Indeed the energetic formulation chosen works well for fractures in mode I, and this experiment shows a combination of mode I and other fracture modes. Furthermore, this regularization applies only for dissipation due to damage and may suffer from a mesh dependency in the damage map (i.e., the damage path is determined by the orientation of the mesh). To tackle this problem, the authors suggest to use more advanced methods introducing internal length ( [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], ...). are caught by the model for structures sustaining shear load. The main diagonal cracks are described in maps d 1 and d 2 . Furthermore, the horizontal cracks due to global bending observed experimentally on the lateral boundaries are also seen in map d 2 . As the continuous model is intended to model the crack pattern in a diffuse manner/diffusely, the local cracking and some peculiar cracks like the ones observed above the window are not directly described. In order to get this fine information (i.e., crack path, crack opening), one can, for instance, consider post-processing analysis using Discrete Element Method like in [START_REF] Oliver-Leblond | Non-intrusive global/local analysis for the study of fine cracking[END_REF]. Finally, the differences observed experimentally on the global response and the crack patterns may be caught numerically by introducing a spatial variability in the model parameters.

At last, the cyclic behavior, which is of major importance for the calculation of structure under seismic solicitations, will be investigated in the next section.

Confined solid wall under cyclic shear loading

To analyse the model ability to describe the hysteretic response at the structural scale, cyclic tests are performed on a solid panel and compared with the experimental campaign conducted by [START_REF] Anthoine | Shear-compression testing and analysis of brick masonry walls[END_REF]. Experimental conditions are reproduced, and the panel geometry, loading and boundary conditions are given in figure 13. A uniformly distributed vertical pre-compression load p of 0.6MPa is first applied and kept constant during the test, followed by a cyclic horizontal displacement u applied at the panel top. Two or three imposed cycles are performed at each amplitude [START_REF] Anthoine | Shear-compression testing and analysis of brick masonry walls[END_REF]. Regarding the boundary conditions, the bottom side of the wall is entirely blocked, while the top side is prevented from rotating.

The numerical study is conducted using cubic finite elements with linear interpolation functions (CUB8) with a size of 0.125m. The identification of material parameters is carried out on a monotonic shear test, allowing the envelope curve of the cyclic test to be matched.

The numerical cyclic response of the low panel is plotted in figure 13 and compared with the experimental one. The numerical analyses show good agreement in reproducing the force-displacement curve: the general behavior of the hysteresis mechanism is well represented. Due to friction between the crack surfaces, the hysteresis phenomenon develops with the progressive degradation of the masonry. The loops size enlarges with the increase of damage, which is in accordance with the rise of the surface area that can slide with the development of new crack surfaces. As the loading progresses, the shear modulus degrades, causing the horizontal tilt of the loops. At the end of the test, the stiffness drops and the last loops are very large, indicating that the specimen is close to failure. Figure 14 shows the damage field D 2 obtained at the end of the loading for the numerical test and the

Concluding remarks

In this paper, an original constitutive model for a masonry-like material has been presented. The model has been derived to deal with analysis involving multiaxial cyclic loadings. The theoretical framework of the thermodynamics of irreversible processes is considered to define the nonlinear processes, their evolution and their interactions with other processes. To describe the orthotropic nature of damage, a decomposition by crack families associated with each natural direction of the masonry is proposed. Crack families are defined by a scalar damage variable and a fabric tensor that define the effect of a crack family on the effective compliance tensor. The unilateral effect observed with cyclic loading is modeled thanks to a decomposition of the stress tensor in positive and negative parts. Hysteretic effects and permanent strains for the shear components are described with a direction-by-direction damage-sliding coupling. The numerical implementation of the model has been extensively described. The choice of a decomposition direction-by-direction allows to simplify the implementation of the model and provides a simplicity and robustness to the model. To reduce the classical mesh dependency observed for softening media, a simple energetic regularization approach has been considered. More complex approach like non-local methods can be derived easily by considering as non-local quantities the equivalent strains that drive each damage variables. This development will be carried out in future works. It should be mentioned the set of parameters that drive the nonlinear process has been defined and linked to specific phenomena. As a consequence, they can be obtained thanks to experimental or virtual tests. The capacity of the model to describe the response of masonry-like material for complex loadings has been illustrated with local tests. More particularly, the model achieves to reproduce hysteretic loops observed for cyclic shear loadings. Furthermore, the effect of confinement on the friction developed with shear loading is well-reproduced by the model. Finally, tests at the structural element level have been performed illustrating the capacity of the model to describe the response of masonry structures. Global quantities, as well as local ones like damage fields, are qualitatively well reproduced compared with experimental data. To conclude, this work provides a first step to later evaluate the seismic vulnerability of masonry structures like bridges or buildings for low to moderate earthquakes.

Appendices

510

A. State laws

From the state potential ρΨ ⋆ the state laws can be derived: for total strains,

ε 1 = ρ ∂Ψ ⋆ ∂σ 1 = < σ 1 > + E 1 (1 + d 1 ) + < σ 1 > - E 1 - ν 12 σ 2 E 1 - ν 13 σ 3 E 1 ( 65
)
ε 2 = ρ ∂Ψ ⋆ ∂σ 2 = < σ 2 > + E 2 (1 + d 2 ) + < σ 2 > - E 2 - ν 12 σ 1 E 1 - ν 23 σ 3 E 2 (66)
ε 3 = ρ ∂Ψ ⋆ ∂σ 3 = < σ 3 > + E 3 (1 + d 3 ) + < σ 3 > - E 3 - ν 13 σ 1 E 1 - ν 23 σ 2 E 2 (

Figure 1 .

 1 b) shows the associated homogenized material.

Figure 1 :

 1 Figure 1: Definition of the orthotropic directions (i.e., natural frame) of the masonry: a) explicit description, b) homogenized description of the masonry

Figure 2 :

 2 Figure 2: Effect of crack families on the compliance tensor

)

  ⋆ with ρ(Ψ shear e ) ⋆ the elastic energy density associated with the shear stress components and ρ(Ψ normal e

215effect.

  The following loading is applied: ε 11 starts from 0 to 2.5•10 -4 then decreases to -2•10 -4 , and increases up to the 6 • 10 -4 . With the evolution of damage, a progressive decrease of the stiffness with the increase (a) Stress-strain response (b) Evolution of damage variable D 1

Figure 3 :

 3 Figure 3: Unidirectional tensile/compressive test along direction 1

Figure 4 :Figure 5 :

 45 Figure 4: Cyclic shear test

5. 1 .

 1 Confined wall under monotonic shear loading 5.1.1

Figure 6 :

 6 Figure 6: Calibration of the model for a monotonic shear test on a solid panel (experimental result from [61])

Figure 7 :

 7 Figure 7: Mesh sensitivity analysis for the global response of a shear wall

Figure

  Figure8: Boundary conditions and loading for the numerical test of the hollow wall in shear with a normal stress[START_REF] Vermeltfoort | Shear tests on masonry walls[END_REF]

Figure 9 :

 9 Figure 9: Maximum damage value d 1 reached along the diagonal between the bottom left corner of the wall and the window.

Figure 10 :

 10 Figure 10: Comparison between force/displacement curves of numerical and experimental [60] tests.

  Figures 11 and 12 highlight a comparison between numerical damage maps and experimental crack maps (bottom) for d 1 (right) and the vertical displacement field (left), and d 2 (right) and the horizontal displace-ment field (left) respectively. The labels J2G and J3G in figures 11 and 12 correspond to labels of two experiments made in [60]. Figures 10, 11 and 12 show that global damage mechanisms at structural scale

Figure 11 :

 11 Figure 11: Comparison between the experimental vertical failure pattern (down)([60] and [33]) and the horizontal component of the displacement field (top left) and the damage field d 1 (top right) at 20 mm horizontal displacement results.

Figure 12 :

 12 Figure 12: Comparison between the experimental horizontal failure pattern (down)([60] and [33]) and the vertical component of the displacement field (top left) and the damage field d 2 (top right) at 20 mm horizontal displacement results.

Figure 13 :

 13 Figure 13: (left) Geometry of low panel [4] (right) Comparison between experimental (dot) and numerical (solid) forcedisplacement response curve under cyclic loading for the low panel

Figure 14 :

 14 Figure 14: Comparison between (left) experimental crack pattern [4] and (right) numerical damage map for the low panel

Table 1 :

 1 Variables summary table

	Mechanism	Parameter name	Parameter symbol
		Young modulus in i-direction	E i
	Elasticity	Poisson ratio in ij-plane	ν ij
		Shear modulus in ij-plane	G ij
		Shear participation of ε ij in i-extension's-direction	β ij
	Damage threshold and evolution	Extension threshold in i-direction	K i
		Damage brittleness in i-direction	χ i
	Shear damage	Participation of d i in σ l Participation of d j in σ l	(i) ll (j) A A ll
		Sliding coefficient in i shear component	µ i
	Sliding threshold and evolution	Hardening parameter for σ π i Sliding intensity parameter for σ π i	a i b i

Table 2 :

 2 Parameter summary table

  S 11 S 22 S 33 -S 11 S 2 23 -S 22 S 2 13 -S 33 S 2 12 + 2S 12 S 23 S 13 . The normal components of σ σ σ are positive or negative. Each sign combination leads to different cases for the compliance S ii values as shown in table 3. There are 8 cases, and for each one, S 11 , S 22 , S 33 have different values. The normal components of σ σ σ are computed with equations (53) and the new values for S

	sign(σ 1 ) sign(σ 2 ) with ∆S = Case sign(σ 3 )	1 + + +	C 12 = C 13 = C 22 = C 23 = C 33 = 2 -+ +	S 22 S 33 -S 2 23 ∆S S 13 S 23 -S 12 S 33 ∆S S 12 S 23 -S 13 S 22 ∆S S 11 S 33 -S 2 13 ∆S S 12 S 13 -S 23 S 11 ∆S S 22 S 11 -S 2 12 ∆S 3 4 5 + + --+ -+ -+	6 -+ -	7 + --	8 ---	(56) (57) (58) (59) (60) (61)

ij . This iteration starts if the sign of one of the normal components of σ σ σ is different from the previous converged time step. This iterative scheme is summarized in the algorithm 1.

Table 4 :

 4 Material parameters summary table for one element tests

  shows the comparison

	Elastic	Values	Elastic	Values	Elastic	Values
	parameters		parameters		parameters
	E1	5 GPa	G12	1.8 GPa	ν12	0.15
	E2	3 GPa	G13	1.6 GPa	ν13	0.15
	E3	5 GPa	G23	1.6 GPa	ν23	0.15
	Damage	Values	Internal sliding	Values	Internal sliding Values
	parameters		parameters		parameters
	G f 1	100 N/m	µ1	0.1	b6	10 5
	G f 2	80 N/m	µ2	0.1	b5	10 5
	G f 3	100 N/m	µ3	0.1	b4	10 5
	K1	4 × 10 -4	a6	10 -6
	K2	1.3 × 10 -5	a5	10 -6
	K3	4 × 10 -4	a4	10 -6
	A (1) 12	0.8
	A (2) 12	0.55

Table 5 :

 5 Material parameters summary table for one element tests

fracture energy regularization provides similar peak load and damage field results with a structured paral-

Acknowledgments

The authors wish to express their most grateful thanks to SNCF RESEAU for its financial and technical support.

This work was performed using HPC resources from the "Mésocentre" computing center of Centrale-Supélec and École normale supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http: //mesocentre.centralesupelec.fr).

B. Algorithm

Algorithm 1 General algorithm of the model 1. Inputs: ε ε ε τ , ∆ε ε ε, σ σ σ τ , d τ 1 , d τ 2 , d τ 3 , (ε π 4 ) τ , (ε π 5 ) τ , (ε π 6 ) τ , X τ 4 , X τ 5 , X τ 6 .

2. Actualization of the strain tensor: ε ε ε τ +1 = ε ε ε τ + ∆ε ε ε 25)) (c) Actualization of the internal sliding and friction variables:

trial and no evolution of internal variables • if f π i > 0 → Return mapping i. Computation of the internal sliding multiplier: ∆λ π i k+1 (Eq. ( 49)) ii. Actualization of the variables: σ π i k+1 (Eq. ( 45)), X i k+1 (Eq. ( 48)) iii. Computation of the internal sliding threshold function: f π i k+1 (Eq. ( 25)) iv. Convergence test: if f π i k+1 ≤ 0 → exit, else → go back to i.


Paul Mazoyer

Salma Yahiaoui

Julien Tissot

Evaluation of wind measurement uncertainty of Windcube in complex terrain

come


Kévin Filippini

Julien Tissot

Understanding Remote Sensor Calibration Uncertainties, The case of a Nacelle-mounted Pulsed LiDAR

come   L'archive ouverte pluridisciplinaire

To simplify, we can identify and regroup the different sources of uncertainties by answering the following questions:

 How accurate is the reference anemometer?

Nacelle-mounted Lidars are increasingly used for power performance verifications.

The IEC 61400-50-3 standard, to be published in 2020, will support the development of these practices with a clear and shared methodology, from calibration and installation of the sensor to uncertainty calculation and final reporting.

The advised method for calculating calibration uncertainties can seem hard to put into practice for most actors of the wind energy. The type of reconstruction algorithm will decide the method for calculating the propagation of uncertainties [3]. In the case of the WI PO 4Beams which reconstructs the wind with a model accounting for shear, the non-linearity of the equations used will require the use of the Monte Carlo Method.

Monte Carlo is a statistical simulation using a great number of random samplings as inputs. When applied to the propagation of uncertainties of a lidar, it consists in: This method is extremely powerful and flexible but is well known to be costly in term of processing time if implemented badly, e.g. using loops in the code.

However, with ingenious use of linear algebra with modern programming language, i.e. python, we were able to reach a total processing time below 1 second for 10 6 samples per RWS. In other words, the propagation of the uncertainties can be done in less than a second for one set of WFC.

This very high number of samples was specifically chosen in order to neglect additional uncertainty due to the MCM itself. Results and comparison to current practices

The table here shows the comparison of the expanded uncertainties (k=2) on the HWS (uHWS) per bin obtained with the MCM and the one as given by the third party certifier method.

Currently most of the certifiers use the following formula :

𝑢 𝐻𝑊𝑆 = σ 𝑖 4 𝑢 𝑅𝑊𝑆,𝑖
4cos(𝛼)cos(𝛽) Where α and β are the half opening horizontal and vertical angles of the WI PO scanning pattern. The final lidar uHWS calculated with the MCM are smaller or equal to the one given by the current method. 95% of all measurements are inside the interval defined as the "expanded uncertainties", i.e. the standard uncertainty multiplied by k = 2 that we give here. In addition to provide equivalent results than an analytical propagation of the uncertainty, the Monte Carlo Method has the unique advantage to be applicable to other WFC, such as Wind Direction and Shear. This method can give the dependence of the uncertainty to any variable on the wind field. For instance, in the graph below, we show the variation of the expanded uHWS in m.s -1 as a function of the HWS and the WD.

Uncertainties associated to lidar come from the calibration process itself.

MCM propagation for WI 4B can be implemented efficiently using uRWS from calibration as an input.

 Resulting uHWS are smaller when using this method.

 It can easily give uncertainties on multiple WFC, giving a better understanding of the lidar behavior.

  a "White box" calibration, the Radial Wind Speed (RWS) measured by each of the lidar's beams are individually compared to a reference wind speed measured by a certified cup anemometer [1].



  Are both sensors measuring at the same location?  What is the statistical deviation of the RWS measured by WI PO? All of these sources are the added in quadrature [1] in compliance with GUM methodology [

  1. A. Borraccino, M. Courtney et R. Wagner, Generic Methodology for Field Calibration of Nacelle-Based Wind Lidars, 2016. 2.Guide to the Expression of Uncertainty in Measurement 3.M. G. Cox et P. M. Harris, Software specifications for uncertainty evaluation, 2006.

Sensor Calibration Uncertainties The case of a Nacelle-mounted Pulsed LiDAR Kévin Filippini, Julien Tissot

  It uses a numeric method to calculate the uncertainties on Wind Field Characteristics (WFC) reconstructed such as Horizontal Wind Speed (HWS) and Wind Direction (WD): the Monte Carlo Method (MCM). Our objective here is to present how to propagate them on the WFC with the MCM and what are its advantages. To do so we will use the case of Leosphere's Wind Iris PO (WI PO) nacelle-mounted Lidar calibration.

	PO.130
	Leosphere, France

Abstract uRWS -Calibration Uncertainty on Radial Wind Speeds Conclusions uHWS -Uncertainty propagation using Monte Carlo Understanding Remote MEET US AT C3-A11