3D hybrid modeling approach combining the finite and
discrete element methods: Validation based on masonry
shear wall tests

Ali Boukham, Vincent Venzal, Thomas Parent, Stephane Morel, Frédéric Dubois, Bernard Solbes

<figure><img src="image_1.png" /></figure>

To cite this version:

Ali Boukham, Vincent Venzal, Thomas Parent, Stephane Morel, Frédéric Dubois, et al. 3D hy-brid modeling approach combining the finite and discrete element methods: Validation based on masonry shear wall tests. International Journal of Solids and Structures, 2024, 289, pp.112638. 10.1016/j.ijsolstr.2023.112638 . hal-05012281

HAL Id: hal-05012281

https://cnrs.hal.science/hal-05012281v1

Submitted on 31 Mar 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Highlights

3D Hybrid Modelling Approach Combining the Finite and Discrete
Element Methods: Validation Based On Masonry Shear Wall Tests

Ali BOUKHAM, Vincent VENZAL, Thomas PARENT, Stephane MOREL, Frederic DUBOIS, Bernard SOLBES

* Limestone masonry walls are tested under shear-compression configuration.

* Experimental displacement and strain fields are measured using Digital Image Correlation.

* A new 3D hybrid modeling approach is developed, combining finite element method and discrete element method.

* An existing frictional cohesive zone model is improved by adding the linear stiffness of the mortar joint.

* Validation and calibration of the proposed modeling approach are conducted by comparing numerical results with experimental data from the shear wall.
# 3D Hybrid Modelling Approach Combining the Finite and Discrete Element Methods: Validation Based On Masonry Shear Wall Tests

Ali BOUKHAMª,*, Vincent VENZALᶜ, Thomas PARENTª, Stephane MORELª, Frederic DUBOISᵇ, Bernard SOLBESª

ªUniversité de Bordeaux, UMR5295, Institut de Mécanique et d'Ingénierie-Bordeaux (I2M), Dépt. Génie Civil et Environnemental (GCE), Bordeaux, 33000, France

ᵇLMGC, Université de Montpellier, CNRS, Montpellier, France

ᶜAIA Ingénierie, 10 rue Ariane, Bât C Mérignac, 33700, France

## Abstract

The masonry structures represent a remarkable portion of the world's built heritage. The preservation of these buildings requires a multidisciplinary diagnostic approach that includes a structural assessment to estimate their stability. However, this task remains difficult because of the high nonlinearity behavior of the masonry material and the complexity of the geometries of this type of building. For this reason, several robust modeling strategies have been developed. Among these strategies, block-based models (BBM) based on the discrete element method and more specifically the NSCD “Non Smooth Contact Dynamics” allows to take into account the actual geometrical texture of masonry. However, in this type of approach, blocks are generally rigid or elastic and nonlinearity appears only at the interfaces between blocks governed by cohesive frictional laws. In this study, a new hybrid approach is

<footnote>* Corresponding author</footnote>

<footnote>Email address: ali.boukham@u-bordeaux.fr (Ali BOUKHAM)</footnote>
proposed. It consists of integrating a model combining (i) a frictional cohesive zone model “FCZM” describing the block-mortar interface and (ii) a quasi-brittle damage model associated with the block while (iii) the elasticity of mortar joints is explicitly taken into account in the FCZM model by adding a linear spring in both normal and tangential directions. Thus, avoiding the use of extended blocks and corresponding homogeneous properties. The validation of this numerical approach is carried out on the basis of an experimental campaign on masonry shear wall test subjected to compression.

*Keywords:* Masonry, Damage, Non Smooth Contact Dynamics, Finite Element Method, Frictional Cohesive Zone Model

# 1. INTRODUCTION

Masonry is one of the oldest construction techniques still used today Aydenlou (2020). This technique includes any construction method based on the principle of stacking blocks. The variety of masonry types is due to the diversity of materials (blocks and mortar joints) and methods of arranging blocks with each other, that differ from one region to another Maldonado et al. (2019). In fact, construction techniques of masonry have often changed over time to meet specific needs. These historical and diversity aspects of masonry make that a great part of the building cultural heritage consists of monumental masonry structures such as towers, castles, churches, mosques, temples, etc.

In recent years, many historical masonry buildings were damaged either due to weathering, fire incidents or earthquakes (Colonna et al. (2017), Bertetto
et al. (2021), Çelebi et al. (2010), Saretta et al. (2021)). These incidents can have drastic consequences humanly, culturally or financially. In 2019, the famous Notre-Dame cathedral of Paris was hit by a devastating fire. The whole timber roof was burned and a part of the stone vaults collapsed (Praticò et al. (2020), Parent et al. (2023)). Thus, the preservation and restoration of such buildings requires a multidisciplinary diagnostic approach that includes a structural assessment to estimate their stability (Roca et al. (2010)). However, the heterogeneity, anisotropy and strong non-linearity of masonry make this task complex. It is even more complicated when historical monuments presents a complex geometry as well as a complex loading history. Hence it is necessary to develop robust numerical modeling strategies dedicated to masonry structures and able to take into account the non-linear phenomenon occurring in the history of the structure.

Several modelling strategies exist in the literature and can be separated into four major families (Smoljanović et al. (2018), Rafiee et al. (2008), D'Altri et al. (2019), Roca et al. (2010)) (i) block based models with or without joints (Abdulla et al. (2017), Petracca et al. (2017)) (ii) continuous models where masonry is considered as a homogeneous continuous material based on different homogenisation techniques (Anthoine (1995), Cavalagli et al. (2011), Di Nino and Luongo (2019), Zucchini and Lourenço (2009), Parent et al. (2017)) (iii) macro-models that model the structure by a system of panels with a non-linear mechanical response (Di Nino et al. (2017)) (iv) geometry-based models that consider the structure as a rigid solid and assume that the masonry has an infinite compressive strength, a zero tensile strength and prohibition
39

40

of block’s sliding (this approach uses only the geometric and loading data of the structure).

41

42

43

44

45

46

47

48

Among the block based models, the NSCD (Non Smooth Contact Dynamics) is one of the discrete element methods developed to model the dynamics of a set of distinct bodies (rigid or deformable) in interaction (Dubois and Jean (2006)). The particularities of this method compared to other discrete element methods are an implicit time integration, an implicit resolution of the contact problem and the possibility of taking into account non smooth contact laws. Thus, several NSCD-based modeling strategies can be adapted for masonry structures.

49

50

51

52

53

54

55

56

57

The different modelling approaches are often validated by comparing the developed model response with the experimental response of in plane confined shear tests of masonry walls under quasi-static monotonic and cyclic loading (Zucchini and Lourenço (2009), D'Altri and de Miranda (2022), Orduña (2017)). This type of test is particularly adapted to mimic the masonry behaviour subjected to wind loads or a wall supporting the thrust of an arch. To a lesser extent, it can also be used to simulate, in a quasi-static manner, the structural responses caused by seismic forces acting within the wall plane.

58

59

60

61

62

63

Moreover, the shear/compression test of masonry walls is considered as an adequate structural test case since it includes all typical masonry failure mechanisms (block-mortar bond tensile failure, block-mortar bond shear sliding, diagonal masonry cracking, masonry crushing and block and mortar tensile cracking as described by Lourenco (1996)). Its response highly
64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

depends on the typology of masonry and on the mechanical behaviour of each masonry constitutive element (blocks, mortar, block/mortar interface) (Tomaževič (2009)). Thus, it allows to assess the capacity of the model to reproduce the complex non linear response of masonry and to capture the observed failure mechanisms. Several authors have conducted compression/shear tests on different types of masonry walls and proposed different test set ups to evaluate the shear mechanical behaviour of masonry walls (Godio et al. (2019), Petry and Beyer (2014), Rezaie et al. (2020)). These studies were intended to assess the impact of the loading type (monotonic or cyclic), the boundary conditions (fixed-fixed boundary condition, cantilever boundary condition), the shear span or the vertical load ratio on the failure mechanisms, stiffness, strength and drift capacity of the wall. Several analytical models were then proposed to describe experimental responses of masonry walls to enlight the community on the limitations of the normative codes (Petry and Beyer (2015), Vanin et al. (2017)). In this work, the scope of the shear wall test is to validate the hybrid modelling approach.

The work presented here is based on the work initially carried out by Venzal et al. (2020) to model the mechanical behaviour of masonry walls under shear loading. In this previous study, a simplified approach based on the NSCD method was implemented. The masonry was modelled in 2D as a stack of homogenized isotropic extended elastic blocks. The stiffness of the mortar joints was integrated in the blocks using linear homogenization. The interaction between these extended blocks was described by a frictional cohesive contact law called FCZM. This frictional cohesive law is adapted to
model the mechanical response of the block mortar interface. It takes into account the progressive damage of the interface in tensile (Mode I), shear (Mode II) and tensile-shear mixed mode. In addition, in the case of a load that combines compression and shear, FCZM allows to gradually take into account friction with respect to the damage evolution. The identification of the mechanical parameters of this FCZM law was carried out by means of an experimental characterization campaign using a direct tensile test on masonry duos and a compression/shear test on masonry triplets (Bisoffi-Sauve et al. (2019), Venzal et al. (2020)).

The limitation of this approach is its inability to reproduce the damage of the blocks observed in walls subjected to shear due to the high compressive stresses and in masonry buildings subjected to a thermo-mechanical loading as in the case of a fire. Neglecting this aspect can highly impact the mechanical response of the masonry wall. Thus, in such a case it is necessary to model the blocks with a non linear damage model rather than a simple linear elastic one. Also, homogenization of such non linear models is a difficult challenge. To overcome both issues, a novel modelling approach is proposed in this work where masonry is modelled by a stacking of non extended, non-homogenized and non linear blocks (a model combining damage and plasticity is used for the blocks (Sellier (2018), Morenon et al. (2022))), that is, keeping their own dimensions and their own material properties. The contribution of the mortar joints on the elasticity of the masonry is then considered by adding linear springs (in both normal and tangential directions) to the FCZM "Frictional Cohesive Zone Model" law. The particularity of the proposed modelling
strategy compared to other hybrid models that combines finite and discrete element methods (Sarhosis and Lemos (2018), Smoljanović et al. (2013)) is the incorporation of a damage-plasticity model in the blocks while solving the contact problem in a NSCD formulation according to FCZM.

The present paper will be divided into two main sections. The first section will describe the experimental program as well as the experimental results obtained. While the second section will focus on the calibrating procedure of the proposed modelling approach, and its validation by comparing the simulated responses to the experimental ones.

## 2. Experimental results

The purpose of this section is to summarize the results of the experimental program carried out by Venzal (2020) on masonry walls subjected to compressive/shear loading.

### 2.1. *Experimental program*

In this experimental campaign, 9 walls were tested under compressive/shear loading. The walls were built outside the test stand on a metallic profile (reinforced HEA300) (Figure 1). They were left to cure for 90 (±7) days under laboratory conditions and were watered for 15 days in order to minimise shrinkage effects. Three different normal compressive stresses (0.3MPa, 0.4MPa and 0.5MPa) were applied by means of a vertical distributed load (three walls for each compressive stress) to quantify its impact on the wall
behaviour.

138

The walls were 1, 85m wide, 1.2m tall and 0.15m thick and were composed of 11 courses and they were constructed by an experienced mason. They were made of soft limestone (E = 11.1 GPa; Rc = 9.5 MPa) collected from Fron-tenac quarry. A natural hydraulic lime mortar often used in ancient masonry building restoration was used. Lime mortar is constituted by NHL3.5 and 0/2mm sand according to mass ratios of water:NHL:sand=1.5:1:6.5 chosen to optimize the workability of the mortar.

146

Concerning the mechanical boundary conditions of the walls, a cantilever configuration is adapted even though it does not represent realistically the boundary conditions of a wall in a building (Petry and Beyer (2014)). This choice is made to facilitate the modeling of this mechanical test, since a cantilever boundary condition is easier to reproduce numerically. In this condition, the height of zero moment $H_0$ remained at a constant value of $H_0 = H$ throughout the entire duration of the shear test. The wall was glued to two metallic supports (reinforced HEM300) using an epoxy resin (Figure 3). The lower metallic support was clamped to the test bench while the upper metallic support authorises the in-plane rotation and vertical displacements. The vertical loading was first applied and kept constant during the test, then an horizontal displacement was imposed to the upper metallic support. Figure 2 illustrates the test stand in which the shear-compression tests were carried out.

147

148

149

150

151

152

153

154

155

156

157

158

159

160
<figure><img src="image_3.png" /><figcaption>Figure 1: Protocol of walls fabrication.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Picture of the test stand.</figcaption></figure>
## 2.2. Measurements

During the tests, the walls were instrumented in order to measure the following quantities (figure 3) : (a) the vertical and horizontal forces applied, which were measured respectively by two load cells with an accuracy of 0.5%; (b) the rotation of the upper metallic support measured by means of two lvdt's placed vertically on the bottom surface of the metallic support; (c) the imposed horizontal displacement measured by an lvdt placed at the level of application of the horizontal force (location of the different load cells and lvdts are shown in the figure 3).

In addition to the above mentioned hard-wired measurements, the 3D displacement and strain fields of one face of the wall were measured using

<figure><img src="image_5.png" /><figcaption>Figure 3: Experimental setup of shear under constant vertical load tests of masonry walls.</figcaption></figure>
“Digital Image Correlation” (DIC) technique. This technique allows to obtain the displacement and strain fields of all the area of interest and thus to accurately analyze the failure mechanisms using the different measured fields. It allows also to compare the measured fields to the ones obtained from numerical simulations. The optical measurement consisted of a stereo-camera system placed on the front of the panel.

The system acquired high-resolution grey-scale images. The typical procedure used was: (1) preparing by lightening the region of interest. Speckles pattern was not realized on the specimen thanks to the natural texture and contrast of the masonry panel (this approach was validated by measuring the noise level on a set of static images before each test); (2) setting up, pointing and focusing the PIKE F-241B (2048 ×2048 pixels) cameras with their SCHNEIDER XENOPLAN objectives (focal length = 29,3 mm) on the specimen; (3) calibrating the camera system by using a standard calibration target plate; (4) running the test and acquiring images; (5) computing the 3D displacement field by using the software VIC-3D. In DIC, it is known that the post-processed fields are highly influenced by the subset size, step size, and filter size parameters (Belda et al. (2020), Yaofeng and Pang (2007)). In this study, these parameters were chosen for each tested wall to obtain an optimal trade-off between spatial precision and smoothness. To resume, the main characteristics of the stereo correlation system used in this study are summarized in Table 1.
<table><caption>Table 1: Materials and DIC parameters used for the stereo correlation.</caption><thead><tr><th colspan="2">Stereo camera vision system</th></tr></thead><tbody><tr><td>Cameras</td><td>2 PIKE F-421B (image size = 2048x2048)</td></tr><tr><td>Lens</td><td>2 Schneider Xenoplan (focal length = 29,3 mm)</td></tr><tr><td>Lighting</td><td>5 LED projectors</td></tr><tr><td>Calibration target</td><td>planar calibration target 1,20x0,84 m2</td></tr><tr><td>Scale factor</td><td>1 Pixel/mm</td></tr></tbody></table>

## 2.3. Experimental results

### 2.3.1. Shear Wall responses

Figure 4 shows the “$F_x - U_x$” experimental responses of all the walls (where $U_x$ is the applied horizontal displacement while $F_x$ is the corresponding horizontal reaction force). These curves were obtained by taking the envelope of the cyclic responses of each tested wall (the cyclic responses of the walls are not the aim of this article). Each of the nine curves exhibit the same shape with (i) an elastic domain, (ii) a pre-peak nonlinear branch and (iii) a post-peak non-linear branch. Based on the “$F_x - U_x$” curves, each wall’s strength and initial stiffness were estimated and given in table 2. They correspond, respectively, to the maximum load of the ‘$F_x - U_x$’ curve and the initial linear portion’s slope. The table 2 highlights that the higher the applied compressive stress the higher the wall strength. This is explained by the fact that a high number of the wall bed (or horizontal) joints are under compressive/shear loading configuration. Also, since the shear strength of the block/mortar interfaces increases with respect to the compressive stresses due to friction, the wall resists to higher forces when the
normal loading is higher. In a similar manner, the stiffness increases with
increasing compressive stress. In fact, one can make the assumption that
even at laboratory conditions, the interfaces are predamaged (low effective
adherent surface compared to the apparent surface) during the fabrication
process and the cure due to the mortar shrinkage. Thus, at low imposed dis-
placement, some interfaces where the shear total stress is below the friction
sliding threshold are blocked. This phenomena can explain the increase of
stiffness at higher compressive stresses Venzal et al. (2020). However, since
the crack surfaces susceptible to slide are not predominant at low imposed
displacements and differs from wall to wall (the predamage level of each wall
is not controlled), the latter conclusion may not always hold. This explains
the high dispersion of the results.

Note that in this linear phase, the wall is developing no damage and no
residual displacement is yet observed on the cyclic unloading branch. The
initial stiffness of the wall exhibited by this linear phase is an important
property to fit the mechanical parameters of the constitutive materials of the
wall as discussed further in the paper.

2.3.2. Shear wall failure mechanisms

The DIC results of the nine walls exhibits two distinct failure mechanisms
depending on the applied normal stress and correspond to those described in
Figure 5. For more detail, the typical displacement and strain fields obtained
for each tested wall are shown in Appendix B (the DIC fields of wall “3 -
0.5Mpa” are displayed in figure 6). The first failure mechanism takes place
for low applied compressive stresses (see figure 5a). The horizontal crack (a
<table><caption>Table 2: Stiffness and strength of the tested walls.</caption><thead><tr><th>σ<sub>N</sub></th><th>Parameter</th><th>wall 1</th><th>wall 2</th><th>wall 3</th><th>Mean value</th><th>CV</th></tr></thead><tbody><tr><td rowspan="2">0.3 MPa</td><td>Stiffness (×10<sup>7</sup>) [N/m]</td><td>5.7</td><td>7.5</td><td>6.1</td><td>6.4</td><td>12%</td></tr><tr><td>Strength [kN]</td><td>69.1</td><td>63.5</td><td>64.4</td><td>65.7</td><td>4%</td></tr><tr><td rowspan="2">0.4 MPa</td><td>Stiffness (×10<sup>7</sup>) [N/m]</td><td>6</td><td>5.2</td><td>9.8</td><td>7</td><td>29%</td></tr><tr><td>Strength [kN]</td><td>70.9</td><td>78.6</td><td>80.6</td><td>76.7</td><td>5%</td></tr><tr><td rowspan="2">0.5 MPa</td><td>Stiffness (×10<sup>7</sup>) [N/m]</td><td>6.3</td><td>7.4</td><td>8.5</td><td>7.4</td><td>12%</td></tr><tr><td>Strength [kN]</td><td>98.3</td><td>99.5</td><td>92.9</td><td>96.9</td><td>3%</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: The mechanical response $F_x - U_x$ of the tested walls for the three applied compressive stresses.</figcaption></figure>
is initiated at the lower left corner on the first bed mortar interface (where the tensile stresses are the highest due to the maximal bending moment at the wall base) and then propagates horizontally along the interfaces of the first bed joint over a greater or lesser distance with respect to the level of the applied compressive stress. Thus, for low applied compressive stress, the crack (a) propagates along the entire length of the first bed joint.

In the second failure mechanism for higher compressive stress (see figure 5b), the crack (a) propagates for a short distance and stops due to high compressive stress ahead of the tip. Then, a second diagonal cracks (b) initiates at the middle of the wall where shear and tensile stresses are the highest and then propagates to form a stair-step crack along the diagonal. Moreover, the $\varepsilon_{xx}$, $\varepsilon_{yy}$ and $\varepsilon_{xy}$ strain fields (where $\varepsilon_{xx}$ and $\varepsilon_{yy}$ are the axial strains along the X and Y axes and $\varepsilon_{xy}$ is the shear strain in the XY plane) show that the diagonal crack initiates first due to the mode I opening of the head joints followed by the mode II sliding of the bed joints at the middle of the wall. Finally, when the diagonal crack reaches the area of the lower right corner of the wall, the lower right block (block (c) in figure 5b) crushes in compression.

Based on the “$F_x - U_x$” curves and on the observed failure mechanisms, four different phases can be distinguished on the shear wall responses. Figure 6 displays the four phases of the response of one tested wall (0.5MPa) and the corresponding displacement and strain fields at different imposed horizontal displacement of each phase (only necessary fields are plotted since it is sufficient to have a good insight into the wall’s crack pattern). The
<figure><img src="image_7.png" /><figcaption>Figure 5: Failure mechanisms of the tested walls : (a) for low applied compressive stresses (b) for high applied compressive stress.</figcaption></figure>
first phase, corresponds to the linear behavior of the wall. This phase is characterized by the initial stiffness. This first phase is followed by a non-linear phase where $F_x$ keeps increasing, the nonlinearity is caused mostly by the horizontal crack (a) initiation and propagation at the block-mortar interfaces. Additionally, another contributing factor in this nonlinear phase is the blocks' pre-peak compressive damage, which contributes to the overall nonlinearity of the wall, as discussed later in this article. The third phase exhibits the initiation and propagation of diagonal crack (b). The last observed phase is characterized by a sharp decrease in $F_x$ due to the crushing of the blocks highly loaded in compression (in the right lower corner of the block) caused by the horizontal and diagonal cracks propagation and so the decrease of the effective surface submitted to the normal compressive load.

To resume, figure 7 displays the horizontal displacement of the nine walls at the end of the mechanical test. As mentioned, the cracking pattern depends highly on the applied compressive stress. It can be also observed that, for each normal stress level, the crack patterns vary between each wall. This is due to the dispersion of the cohesive mechanical properties of the block-mortar interface.
<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: The mechanical response $F_x - U_x$ of the wall 3 (0.5MPa) with the corresponding DIC displacement and strain fields.</figcaption></figure>
## 3. Numerical modelling

This section aims to present the novel modelling approach and its validation by comparing the numerical results to the experimental ones.

### 3.1. Modelling approach

The proposed approach is a 3D simplified micro modelling approach based on the so-called “non smooth contact dynamics method” or NSCD (Dubois and Jean (2006)). In view of the discontinuous character of the masonry formed by stacking of blocks, this block-by-block numerical approach has advantages over other modeling methods given its ability to capture physical phenomena occurring locally at the block/mortar interfaces. The masonry is modeled here by stacking non-extended and non-homogenized blocks (keeping their own dimensions and mechanical properties). Furthermore, mortar joint and block/mortar interface is modeled by an interface with an initial gap (Figure 8). On this basis, (i) the interactions between blocks (estimated by the Discrete Element Method “DEM”) are described by a new frictional cohesive zone model in which is integrated the normal and tangential stiffness of mortar joints and (ii) the damage and plasticity of the finite elements blocks are described by a quasi-brittle damage model. The added linear springs in the contact law allow to explicitly take into account the elastic mechanical behavior of the mortar joint in the model. This modelling procedure have three main benefits: (i) avoiding complex nonlinear homogenization of the homogenized material incorporating the mechanical behavior of the block and the mortar joint. (ii) avoiding modelling the mortar as 3d continuum material with finite elements since a micro detailed approach is more time
<figure><img src="image_10.png" /><figcaption>Figure 7: DIC horizontal displacement fields ($U_x$) at the failure of each tested wall.</figcaption></figure>
consuming than a micro simplified approach (as the one developed in this article) (iii) the ability to adjust the normal and tangential stiffness of the mortar joint, to accurately model vertical and shear stiffness at the scale of the wall. In the following, the new contact law between blocks accounting for stiffnesses of joints is described followed by a quasi-brittle damage model describing the nonlinear behavior of the blocks.

### 3.1.1. Contact law of the interface

The contact law of the interface includes the mechanical behavior of both mortar and the block/mortar interface. This is achieved by adding an elastic spring law in both normal and tangential directions (describing the normal and tangential stiffnesses of mortar joint) to a cohesive frictional law called FCZM initially proposed by Venzal et al. (2020) (Figure 9).

The stress-separation law ($\sigma_I - \delta_I$) in tension and compression (mode I) are defined respectively by equations 1 and 2:

$$ \begin{array}{l} \text{Mode I} \\ \text{tension} \\ (\delta_I \ge 0) \end{array} \quad \left\{ \begin{array}{l} \delta_I = \delta_I^{\text{coh}} + \delta_I^{\text{spring}} \\ \sigma_I = \sigma_I^{\text{coh}} = \sigma_I^{\text{spring}} \\ \sigma_I^{\text{coh}} = (1 - d_{\text{coh}}) K_I^{\text{coh}} \delta_I^{\text{coh}} \\ \sigma_I^{\text{spring}} = (1 - d_{\text{spring}}) K_I^{\text{spring}} \delta_I^{\text{spring}} \\ d_{\text{coh}} = \max_{\text{history}} \left( \max \left( 1 - \frac{\sigma_I^e}{K_I^{\text{coh}} \delta_I^{\text{coh}}} e^{\phi_I \left( \frac{\sigma_I^e}{K_I^{\text{coh}}} - \delta_I^{\text{coh}} \right)} \right); 0 \right) \end{array} \right. \quad (1) $$
$$
\begin{equation}
\begin{aligned}
& \text{Mode I} \\
& \text{compression} \\
& (\delta_I < 0)
\end{aligned}
\quad
\left\{
\begin{aligned}
\delta_I &= \delta_I^{\text{spring}} \\
\delta_I^{\text{coh}} &= 0 \\
\sigma_I &= \sigma_I^{\text{coh}} = \sigma_I^{\text{spring}} = (1 - d_{\text{spring}}) K_I^{\text{spring}} \delta_I^{\text{spring}}
\end{aligned}
\right.
\tag{2}
\end{equation}
$$

Where $\delta_I$ the total normal opening of the interface, $\delta_I^{\text{spring}}$ the elongation (opening) of the normal elastic spring, $\delta_I^{\text{coh}}$ the normal opening of the FCZM law, $\sigma_I = \sigma_I^{\text{coh}} = \sigma_I^{\text{spring}}$ the actual normal stress of the total law which is the same for the normal spring and FCZM law in so far as the elastic spring and the FCZM law are connected in series. Concerning the cohesive part of the contact law, $d_{\text{coh}}$ is the damage variable of the cohesive part that can only increase (the damage law is defined based on the assumption that the post-peak branch of the cohesive law follows an exponential softening Venzal et al. (2020)), $\phi_I$ a parameter computed to ensure that the dissipated energy is equal to the mode I fracture energy $G_{fI}$ : $\int_0^\infty \sigma_I^e e^{\phi_I(\delta_I^e - \delta_I^{\text{coh}})} d\delta_I^{\text{coh}} = G_{fI}$ and $(\sigma_I^e, K_I^{\text{coh}}, G_{fI})$ are the input parameters of the law in mode I defining respectively the elastic limit stress [Pa], the interface stiffness [Pa/m] and the cohesive energy [J/m²]. Concerning the normal elastic spring which describes the normal stiffness of the mortar joints in the wall, this one is calibrated from the theoretical one $K_I^{\text{spring}} = \frac{E_{\text{limemortar}}}{e_j}$ (335)

For mixed mode shear+compression, the stress-sliding law ($\sigma_{II} - \delta_{II}$) is
defined by the following set of equations 3 :
$$
\begin{equation}
\begin{aligned}
\text{Mixed} \quad \left\{
\begin{aligned}
\delta_{II} &= \delta_{II}^{\text{coh}} + \delta_{II}^{\text{spring}} \\
\sigma_{II} &= \sigma_{II}^{\text{FCZM}} = \sigma_{II}^{\text{coh}} + \sigma_{\mu} = \sigma_{II}^{\text{spring}} \\
\sigma_{II}^{\text{spring}} &= (1 - d_{\text{spring}}) K_{II}^{\text{spring}} \delta_{II}^{\text{spring}} \\
\sigma_{II}^{\text{coh}} &= (1 - d_{\text{coh}}) K_{II}^{\text{coh}} \delta_{II}^{\text{coh}} \\
|\sigma_{\mu}| &\le \sigma_{\mu_c} = \mu d_{\text{coh}}^p |\sigma_N| \\
d_{\text{coh}} &= \max_{\text{history}} \left( \max\left(1 - \frac{\sigma_{II}^e}{K_{II}^{\text{coh}} \delta_{II}^{\text{coh}}} e^{\phi_{II} \left( \frac{\sigma_{II}^e}{K_{II}^{\text{coh}} \delta_{II}^{\text{coh}}} - \delta_{II}^{\text{coh}} \right)} ; 0 \right) \right)
\end{aligned}
\right\} \\
\text{Mode}
\end{aligned}
\tag{3}
\end{equation}
$$

In a similar way to mode I, $\delta_{II}$ corresponds to the total tangential displacement of the interface, $\delta_{II}^{\text{spring}}$ the elongation of the tangential elastic spring, $\delta_{II}^{\text{coh}}$ the tangential sliding of the FCZM law, $\sigma_{II} = \sigma_{II}^{\text{FCZM}} = \sigma_{II}^{\text{spring}}$ the actual tangential stress of the total law which is the same for the tangential spring and FCZM law. Concerning the cohesive part, $d_{\text{coh}}$ is the damage variable already introduced in the description of the Mode I law, $\phi_{II}$ a parameter computed to ensure that the dissipated energy is equal to $G_{fII} : \int_0^\infty \sigma_{II}^e e^{\phi_{II}(\delta_{II}^e - \delta_{II}^{\text{coh}})} d\delta_{II}^{\text{coh}} = G_{fII}$, $\sigma_N$ is the compressive stress applied and $(\sigma_{II}^e, K_{II}^{\text{coh}}, G_{fII}, p, \mu)$ are the input parameters of the law in mode II defining respectively the elastic limit stress [Pa], the FCZM stiffness [Pa/m], the cohesive energy [$J/m^2$], the exponent of the frictional stress evolution and the coefficient of friction. Concerning the tangential elastic spring which describes the tangential stiffness of the mortar joints in the wall, this one is calibrated from the theoretical one $K_{II}^{\text{spring}} = \frac{G_{\text{limemortar}}}{e_j}$ (where $G_{\text{limemortar}} = \frac{E_{\text{limemortar}}}{2(1+\nu)}$ is the shear modulus of lime mortar) using the same damage parameter $d_{\text{spring}}$ as in mode I.
① Non linear finite elements

② Interface with an initial gap

<figure><img src="image_11.png" /><figcaption>Figure 8: Modelling approach adapted for masonry : (a) masonry constituents and (b) modelling strategy of the assembly</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Rheological model of the proposed contact law : (a) in mode I (b) in mode II
(c) masonry constituents.</figcaption></figure>
Noting that in tension, an exponential cohesive law is used for the FCZM law (eq 1). While in compression, a signorini law is used since a physical interface should exhibit an infinite stiffness in compression (eq 2). For this reason, the stiffness of the total law in compression is the same as of the elastic spring. While for shear+compression mode, the FCZM law considers that the shear stress is the sum of an exponential cohesive stress (mode II) and frictional stress, where this latter increases with the evolution of $d_{coh}$ (eq 3).

Figure 10a illustrates the mechanical behavior of this contact law in both pure mode I. The response of the spring (represented by the red curve) is consistently linear elastic in both compression and tension, with the same stiffness considered for both regimes. The tensile response of the cohesive part (represented by the blue curve) exhibits an initial linear elastic phase as long as the stress $\sigma_I$ does not exceed the elastic limit $\sigma_I^e$, i.e., $\sigma_I < \sigma_I^e$. During this phase, the interface is undamaged ($d = 0$). Once the stress reaches the value $\sigma_I^e$, the interface starts to undergo damage, which is reflected by a softening phase of the curve. Indeed, this softening part is accompanied by a decrease in the stiffness. The stress continues to decrease until it tends toward 0, which means that the interface is fractured and no longer transmits stress ($d = 1$). It should be noted that the linear part of the curve depends on the initial stiffness and the elastic stress, whereas the softening part is governed by the dissipated cohesive energy, which corresponds to the area under the curve. For the compressive response of the cohesive part, an infinite stiffness is considered (Signorini law) and there is no damage evolution
380 in compression.

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

Based on the response of both the linear spring and the cohesive parts, the total response of the contact law (represented by the black curve) exhibits a similar behavior in tension as the cohesive part. The linear spring affects only the total initial stiffness and the ratio of the stored elastic energy to the total dissipated energy. Indeed, the total dissipated energy is equal to the cohesive energy since the linear spring does not contribute to the dissipation of energy). Furthermore, in compression, the total response is similar to the spring's response (since the interface stiffness is infinite in compression, the total stiffness is equal to the spring's stiffness).

The same observations holds for mode II (Figure 10b). However, as mentioned above under a combined load of compression and shear, the damaged surface experiences friction. Hence, it is necessary to consider the contribution of the frictional stress. According to the FCZM model, the shear stress is expressed as the sum of the cohesive stress (mode II) $\sigma_{II}^{coh}$ and a frictional stress $\sigma_{\mu}$. As it can be seen from figure 10b, the frictional cohesive response (depicted by the blue solid curve) results from the summation of the cohesive and frictional curves. Noting that the frictional stress increases with damage. For $d = 0$ the frictional stress is null and thus the response is purely cohesive, while for $d = 1$ the frictional stress is maximal and the cohesive stress is null, and thereby causing the response to become primarily frictional. Similar to the mode I, the linear spring affects only the total initial stiffness $K_{II}^{eq}$ and the ratio of the stored elastic energy to the total dissipated
405 energy.

406
407 The material’s parameters of the cohesive frictional part of the interface
408 were identified from triplet shear test and duo direct tensile test Venzal et al.
409 (2020). Moreover, the stiffness of the normal and tangential springs (K<sub>I</sub><sup>spring</sup>
410 ; K<sub>II</sub><sup>spring</sup>) were determined from the elastic modulus, shear modulus and the
411 thickness of the lime mortar joints (E<sub>limemortar</sub> = 4.3 GPa ; G<sub>limemortar</sub> =
412 1.8 GPa ; e<sub>j</sub> = 7 mm) according to the equations : (i) K<sub>I</sub><sup>spring</sup> = E<sub>limemortar</sub>/e<sub>j</sub>
413 (ii) K<sub>II</sub><sup>spring</sup> = G<sub>limemortar</sub>/e<sub>j</sub>. Where E<sub>limemortar</sub>, G<sub>limemortar</sub> were determined by
414 Bisoffi-Sauve (2016). To resume, the law contact parameters are shown in
415 table 3.
<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>Figure 10: Mechanical response of the proposed contact law : (a) in mode I (b) in mode II.</figcaption></figure>
<table><caption>Table 3: Law contact (FCZM+spring) parameters used. (Venzal et al. (2020))</caption><tr><td rowspan="3">FCZM (mode I)</td><td>σ<sup>e</sup><sub>I</sub> (MPa)</td><td>0.05</td></tr><tr><td>K<sup>coh</sup><sub>I</sub> (Pa/m)</td><td>1.21 × 10<sup>11</sup></td></tr><tr><td>G<sub>fI</sub> (J/m<sup>2</sup>)</td><td>3</td></tr><tr><td rowspan="3">FCZM (mode II)</td><td>σ<sup>e</sup><sub>II</sub> (MPa)</td><td>0.27</td></tr><tr><td>K<sup>coh</sup><sub>II</sub> (Pa/m)</td><td>2.8 × 10<sup>11</sup></td></tr><tr><td>G<sub>fII</sub> (J/m<sup>2</sup>)</td><td>206</td></tr><tr><td rowspan="2">Linear spring</td><td>K<sup>spring</sup><sub>I</sub> (Pa/m)</td><td>4.3 × 10<sup>11</sup></td></tr><tr><td>K<sup>spring</sup><sub>II</sub> (Pa/m)</td><td>1.8 × 10<sup>11</sup></td></tr></table>
3.1.2. Quasi-brittle damage model for blocks

For the blocks, a model called “ENDO3D” coupling damage and plastic-
ity was used (Morenon et al. (2022), Sellier (2018)). Specifically, this model is
suitable to simulate the mechanical behaviour of quasi-brittle materials such
concrete or stones since it is capable of capturing the asymmetric response
in tension/compression, dilatancy, energy regularization in tension (with the
Hillerborg method) and cracks reclosure in compression. This is achieved by
incorporating a multi-surface plasticity approach and different damage ten-
sors for tension and compression.

425

Figure 11 illustrates the mechanical response of the ENDO3D model in
the case of a isostatic cube (modeled by one finite element) subjected to
a uniaxial tensile/compression load. At the beginning, the stress increases
linearly up to the tensile strength (point A), insofar as pre-peak damage is
disabled. Then a Rankine criterion is activated, which leads to the devel-
opment of plastic strain and post-peak damage which causes the softening
of the law. When unloading is applied (point B), the initial stiffness is re-
stored at a certain negative stress equal to the crack reclosing stress (point
C). Indeed, at point C the reclosing criterion of cracks is activated (similar-
ily to the Rankine criterion). From point C to D, the phase corresponds to
a pre-peak compressive damage. At the compressive strength (point D), a
Drucker-Prager criterion activates and causes post-peak damage. This phase
is accompanied by an inversion the sign of variation in volumetric strain
(green curve) attesting to the dilating behavior of the ENDO3D model after
the compression peak. After unloading once again (point E), the stress in-

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440
creases until it reaches the residual tensile stress previously obtained at the end of the first tensile cycle (Point F). In parallel one reactivates the damage post-peak in tension until the complete rupture in tension (point G). Note that the post-peak damage is always accompanied by plasticity since the post-peak damage depends on the plastic strain, which is computed from the effective stress tensor.

The main model parameters are summarized in Table 4, as follows: E (Young's modulus), ν (Poisson's ratio), $R_t$ (tensile strength), $\varepsilon_{pic_t}$ (strain at tensile peak stress), $G_{ft}$ (tensile surface energy dissipation), $R_c$ (compressive strength), $\varepsilon_{pic_c}$ (strain at compressive peak stress), $\delta$ (coefficient representing the effect of confinement pressure on compressive strength), $\beta$ (coefficient controlling dilation), and $\varepsilon_{kdc}$ (characteristic strain controlling the compressive post-peak behavior). The block's mechanical properties showed in Tables 4 were taken from (Bisoffi-Sauve (2016), Parent (2015)) due to the similarities between the materials tested in these studies and the ones used for the masonry walls. Note that Figure 11 illustrates the damage-plasticity model response under uniaxial tensile-compression loading with the parameters given in Table 4.

Finally, the proposed contact law as well as the damage-plasticity model were implemented in the discrete element software LMGC90. Thus providing a hybrid model merging continuum and discrete approaches dedicated to the modeling of the mechanical behaviour of stone masonry.
<figure><img src="image_15.png" /><figcaption>Figure 11: Model response of a cube under a uniaxial tensile/compression test along z axis with parameters calibrated on the soft Frontenac limestone.</figcaption></figure>
<table><caption>Table 4: ENDO3D parameters used</caption><tbody><tr><td rowspan="2">Elastic</td><td><i>E</i> (GPa) Bisoffi-Sauve (2016)</td><td>11.1</td></tr><tr><td><i>ν</i> Bisoffi-Sauve (2016)</td><td>0.2</td></tr><tr><td rowspan="3">Tensile</td><td><i>R</i><sub>t</sub> (MPa)</td><td>1.2</td></tr><tr><td><i>ε</i><sub>pic</sub></td><td>0.11 × 10<sup>-3</sup></td></tr><tr><td><i>G</i><sub>ft</sub> (J/m<sup>2</sup>)</td><td>18.8</td></tr><tr><td rowspan="5">Compression</td><td><i>R</i><sub>c</sub> (MPa) Bisoffi-Sauve (2016)</td><td>9.5</td></tr><tr><td><i>ε</i><sub>pic</sub></td><td>0.13 × 10<sup>-2</sup></td></tr><tr><td><i>δ</i> Parent (2015)</td><td>0.4</td></tr><tr><td><i>β</i> Parent (2015)</td><td>0.2</td></tr><tr><td><i>ε</i><sub>kdc</sub></td><td>0.15 × 10<sup>-2</sup></td></tr></tbody></table>
3.2. Boundary conditions and loading

In LMGC90 software, the 3D geometry of the masonry wall was created using a micro modeling technique. Masonry units were meshed using cubic elements (2 elements along x,y axes and one element along z axis), while the mortar and block/mortar interface were modeled by a 2D contact elements with an initial gap (4 contact points by face element) (figure 12). The number of mesh elements and of contact points were determined in order to optimize time calculation while preserving accurate results. Rigid elements were used to model the upper and lower metallic supports due to their high stiffness compared to the wall's stiffness. Figure 12 displays the discrete element model of the masonry wall with the corresponding block's finite element meshes and contact points. Boundary conditions similar to the experimental ones were applied to the numerical model. The first and last courses of blocks are respectively clamped to the lower and upper rigid bodies (metallic supports) through a contact law not allowing relative displacements between elements. All the degrees of freedom of the lower rigid body are constrained, while for the upper rigid body the in-plane rotation is allowed while a vertical load and an horizontal displacement are imposed (figure 12). Regarding the load history: (i) gravity, (ii) normal loading and (iii) imposed horizontal displacement were applied successively and gradually using tangent hyperbolic functions to avoid sudden variations that can lead to dynamic effects.

The simulations presented in section 3.3 were run in parallel on two In-
tel(R) Xeon(R) Gold 6248R 3.00GHz CPUs using OpenMP for parallelization.
A time step of dt = 1.25 × 10⁻⁴ seconds was used and the run time of
<figure><img src="image_16.png" /><figcaption>Figure 12: Numerical model of the wall.</figcaption></figure>

simulations with nonlinear blocks was 1200 minutes compared to 850 minutes
with linear blocks.

3.3. *Numerical results*

In this section, the different numerical results obtained from the numeri-
cal wall model with linear elastic blocks as well as with nonlinear blocks are
presented. This comparison aims to highlight the model's capability with lin-
ear blocks in reproducing the experimental results while also shedding light
on its limitations. Additionally, the enhancements brought about by the use
of the hybrid model are shown. Within the same section, the calibration
procedure for the different models is presented.

The "$F_x - U_x$" responses obtained from the **model with linear blocks**,
using the parameters derived from the characterization tests (shown in table 3) and without taking into account any pre-damage in the spring and the cohesive law i.e. with a $d_{spring_0} = d_{coh_0} = 0$ are plotted in figure 14a. It can be seen that the numerical model is not able to reproduce the experimental responses (for the three normal stress level), the obtained wall's stiffness and strength are highly overestimated. This can be attributed to the fact that the mechanical parameters employed for the mortar and the interface may not be applicable at the structural scale and might be overestimated. In the case of the specimen, curing conditions and fabrication protocols are carefully optimized and well-controlled, which is not necessarily the case at the scale of a wall.

To obtain the correct horizontal wall stiffness (figure 14b), an isotropic damage (isotropic means here that same a damage value is shared for both normal and tangential directions) $d_{spring_0} = 0.968$ is applied on the elastic spring part of the contact law (representing the mortar joint). The applied damage can be seen as a simple decrease of the mortar joint stiffness compared to the stiffness expected of elastic properties of lime mortar estimated experimentally. Note that $d_{spring_0} = 0.968$ is considered at the onset of simulation and this value is kept constant during computation. The impact of this initial isotropic damage $d_{spring_0}$ on the mode I response of intact contact law is depicted in the figure 13. As it can be seen, only the stiffness of the linear elastic spring (red curve) is decreased while the cohesive part (blue curve) stay unchanged. Therefore, the stiffness of the total contact law (black curve) is decreased as well, resulting in an increase in the ratio of
elastic energy to dissipated energy. These observations hold true for mode II as well. Furthermore, it is possible to fit the wall's shear stiffness by calibrating the elastic properties of the block. But this choice is unrealistic, since the elastic properties of a limestone test sample and a limestone block are the same and known a priori (due to the dimensions of the representative elementary volume of limestone). Also, the block/mortar interface is not a bulk and thus having a theoretical infinite stiffness and it was shown in a previous work Venzal (2020) that it is impossible to fit the wall's shear stiffness by calibrating the interface stiffnesses. This discrepancy imply that the difference in the observed and the numerical stiffness is conditioned by the mortar joint stiffness. Thus, fitting the shear wall stiffness by calibrating the elastic spring's stiffnesses appear to be more physical. For these reasons, the linear springs were added to the cohesive contact law.

541

Although the stiffness is correctly captured by the numerical model, the nonlinear phase is still not correctly reproduced and the wall's strength is still overestimated (figure 14b). This is a consequence of the fact that the elastic cohesive stress and the cohesive energy of mode I and II are overestimated and were not decreased by the parameter $d_{spring_0}$. To overcome that, another initial damage value $d_{coh_0}$ were considered for the frictional cohesive part of the contact law (representing the block/mortar interface behavior). Figure 13 shows the impact of $d_{coh_0}$ on the intact cohesive component (blue curve) and on the total law (black curve) in mode I. The elastic cohesive stress, interface stiffness and the cohesive energy are decreased. The same observation hold true for pure mode II. In a compressive-shear mode the impact of $d_{coh_0}$

542

543

544

545

546

547

548

549

550

551

552
is explained in Venzal et al. (2020). In Appendix A the impact of an initial damage $d_{coh_0}$ on the shear wall response curve is depicted. In these figures, the set of $d_{coh_0}$ tested were chosen so that they reduce the initial cohesive strength of mode II by 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80% and 90%. The best fit of the experimental curve is obtained with an initial damage $d_{coh_0} = 0.99944$ (figure 14c), corresponding to a 60% decrease of the shear strength $\sigma_{II}^e$ and 85% of the tensile one $\sigma_I^e$. The corresponding numerical results gives satisfying results and the three first phases of the wall's response mentioned previously were well captured. However, since the model with linear blocks is not capable of reproducing the compressive and tensile damage of the blocks, it was not able to capture the sharp decrease of the horizontal force observed experimentally due to the crushing of the blocks highly loaded in compression. Note that, the initial damage $d_{coh_0}$ applied on the cohesive part of the contact law (part describing the behavior of the block/mortar interface) is necessary for the accurate modeling of the failure mechanisms, this initial damage does not affect the initial stiffness of the wall (unlike the initial damage of the joint $d_{spring_0}$).
<figure><img src="image_17.png" /><figcaption>Figure 13: Impact of $d_{spring_0}$ and $d_{coh_0}$ on the total, FCZM and spring mode I laws</figcaption></figure>

Figure 14d illustrates the response of the model with nonlinear blocks. Finally, the model gives more accurate results using the same $d_{coh_0}$ and $d_{spring_0}$ value as in the model with linear blocks. It successfully captures all the phases of the wall's response including the sharp decrease of the horizontal force. Furthermore, it effectively reproduces all the discussed failure mechanisms including cracks at both interfaces and blocks as shown in figure 15. It can be seen that the same crack patterns observed experimentally are reproduced, as well as the displacement fields $U_x$ and $U_y$ obtained by DIC as shown in figures 15, 16.

Furthermore, the applied initial damages $d_{coh_0}$ and $d_{spring_0}$ are high and may seem nonphysical. However, applying such high initial damage values is due to numerous physical factors (assumptions not yet validated) related to
cure conditions and scale effects between a masonry sample and a masonry
structure:

* The process of building a masonry wall is more complex than making a masonry sample. In fact, with all taken precautions the mortar joints made by the mason would not be as perfect as the joints of a masonry sample. With this in mind, we can assume that the mortar's defect density at the scale of the wall is higher compared to the masonry sample. These defects are voids in the apparent mortar joint's volume which leads to lower apparent mechanical properties. Also, since the joint's thickness is only 10 mm while a 0.2 mm sand is used the representative elementary volume of a mortar joint can be considered smaller than the elementary volume representative of a 4x4x16cm3 test specimen. Thus, the lime mortar material of a thin joint is different from that of a 4x4x16cm3 specimen.

* The limestone blocks used in this study exhibits high porosity (20-30%) and high permeability. In this case, the hydraulic curing of lime mortar in a joint between two limestone blocks is disrupted near the joint-block interfaces by the absorption of water contained in the mortar by the blocks. This water absorption strongly degrades the mechanical quality of the mortar near the interfaces. Note that in the case of low-porosity stone, this phenomenon of mechanical degradation of interfaces (due to mortar hydration failure) is probably less significant than in the case of high-porosity stone.

* In a structure (in this case, a wall), hydraulic shrinkage of the joints
during curing has an impact on the mechanical properties of the joint
mortar and block/mortar interface due to the arrangement of the blocks
in the masonry structure. For example, horizontal shrinkage of verti-
cal joints occurs freely without a rigid movement of the adjacent blocks
(assembled by the vertical joints). This leads to a relative displacement
at the level of the interface and as a consequence to lower cohesive prop-
erties of the interface. This is not the case for the vertical shrinkage of
the horizontal joints since the weight of the adjacent blocks allows for
a rigid movement of these blocks and thus a null relative displacement
at the level of the interface. In the same way, the horizontal shrinkage
of horizontal joints and vertical shrinkage of vertical joints are dis-
turbed (especially near block corners). This phenomenon of shrinkage
of mortar joints, which takes place in any masonry structure, leads to a
deterioration in the mechanical properties of the lime mortar joint and
of the block/mortar interface compared to the mechanical properties
expected from the test specimens (in which these phenomena do not
occur).

In addition, it should be noted that the limestone blocks tend to exhibit
premature damage in compression at low imposed horizontal displacements,
as depicted in Figure 15d. The effect of this premature damage in compres-
sion is observed in figure 11 : the branch CB exhibits non-linearity from 50%
of Rc to the peak strength. This premature damage has a slight unnoticeable
effect on the stiffness and strength of the wall model with linear blocks (if
the same d$_{spring}$$_{0}$ and d$_{coh}$$_{0}$ calibrated for the model with non linear blocks
are used for the model with linear blocks). In cases where the extent of pre-
633 peak compressive damage in limestone is more significant, its impact on the

634 stiffness and strength of the wall become more pronounced. Consequently,

635 the assumption of linear elastic blocks is no longer appropriate for masonry

636 walls constructed with soft limestone exhibiting strong pre-peak compressive

637 damage, even under moderate loading conditions.
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d)</figcaption></figure>

Figure 14: Comparison of the numerical and experimental response curves Fx-Ux (blocks with a damage-plasticity model) : (a) $d_{coh_0} = d_{spring_0} = 0$ ; linear blocks (b) $d_{spring_0} = 0.968$ ; $d_{coh_0} = 0$ linear blocks (c) $d_{spring_0} = 0.968$ ; $d_{coh_0} = 0.99944$ ; linear blocks (d) $d_{spring_0} = 0.968$ ; $d_{coh_0} = 0.99944$ ; non linear blocks.
<figure><img src="image_22.png" /><figcaption>Figure 15: Different damage fields of the numerical wall with nonlinear blocks at ($\sigma_N$ = 0.5MPa) : (a) compressive damage in the blocks ($U_x$ = 8mm) (b) tensile damage in the blocks ($U_x$ = 8mm) (c) interface damage ($U_x$ = 8mm) (d) compressive damage in the blocks ($U_x$ = 0.5mm).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 16: Comparison of the numerical and experimental displacement field for test wall 2 - 0.5MPa at an imposed displacement of 6 mm : (a) horizontal displacement $U_x$ (b) vertical displacement $U_y$</figcaption></figure>

4. Conclusions and perspectives

In this work, a new modelling strategy was developed in which damage can **initiate** in both interfaces and blocks. This new approach, al-
lows to capture all the typical failures mechanisms experienced by masonry structures even at large displacements. Also, since the appropriate mechanical parameters found experimentally are put in the blocks and in the joint+interface/block assembly, no homogenization step is needed. This is achieved by introducing a new contact law in which a normal and a tangential elastic spring describing the normal and tangential stiffnesses of mortar joints are added to a frictional cohesive zone model describing the block-mortar interface behaviour.

The importance of considering the nonlinearity in blocks becomes evident when comparing the numerical results of two wall models: one using linear elastic blocks and the other incorporating non-linear blocks. The *model with linear blocks* produces satisfactory results, successfully capturing the failure mechanisms where cracks initiate and propagate along the joint/block interface. However, by introducing damage and plasticity in the blocks, a more accurate estimation of the wall's stiffness and strength is achieved. The inclusion of damage and plasticity allows for a better representation of the wall's behavior. The early pre-peak compressive damage in the blocks slightly affects the wall's stiffness and strength, this impact becomes no longer neglected with higher prepeak compressive damage. This implies that the *model with nonlinear blocks* reproduces more accurately the shear wall's response. Furthermore, the *model with nonlinear blocks* is capable of capturing the force decrease observed at high displacements, which is attributed to the crushing of the lower right block. This failure mechanism could not be captured by the *model with linear blocks*, and it may be significant in accurately
predicting the behavior of the wall under extreme loading conditions. Thus, the model with nonlinear blocks can provide a deeper understanding of the structural behavior and accurately assess the condition of masonry buildings.

In the scope of this work, the thermal damage of masonry constituents will be incorporated in the contact law and in the law behavior of the blocks in order to take into account the impact of the temperature on masonry. The proposed modeling approach will be applied specifically to assess the stability of the cathedral of Notre Dame of Paris, following the fire incident. An experimental campaign is also planned to properly characterize the thermomechanical behavior of the materials used in the structure under investigation.

## 5. Acknowledgements

The authors express their gratitude to AIA Ingénierie company, particularly Alain Popinet, Agency Director, for their financial support in conducting this study. The authors thank also Philippe TARIS for his expertise in mechanical assemblies and invaluable assistance in developing the test setup.

The French National Research Agency (ANR) is also acknowledged for funding the DEMMEFI project (Grant ANR-20-CE22-0004). Special thanks are extended to the CESA, TMH, and DOYEUX Sablières companies for providing the necessary materials, including sand and lime, for the research.

The experimental phase of the study was done using the Xyloplate platform equipment from the Equipex XYLOFOREST (ANR-10-EQPX-16), and the authors appreciate the support received from the French National Research
<footnote n="690">690 Agency (ANR) in this regard.</footnote>
<figure><img src="image_24.png" /><figcaption>(a) 0.3MPa</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b) 0.4MPa</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure A.17: Comparison of the numerical and experimental response curves Fx-Ux with different values of initial damage $d_{coh_0}$ (model with linear blocks).</figcaption></figure>

692

Appendix B.
<figure><img src="image_27.png" /><figcaption>Figure B.18: In plane displacement and strain fields of the wall 1 (0.3MPa) at three different imposed displacements..</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure B.19: In plane displacement and strain fields of the wall 2 (0.3MPa) at three different imposed displacements.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure B.20: In plane displacement and strain fields of the wall 3 (0.3MPa) at three different imposed displacements.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure B.21: In plane displacement and strain fields of the wall 1 (0.4MPa) at three different imposed displacements.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure B.22: In plane displacement and strain fields of the wall 2 (0.4MPa) at three different imposed displacements.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure B.23: In plane displacement and strain fields of the wall 3 (0.4MPa) at three different imposed displacements.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure B.24: In plane displacement and strain fields of the wall 1 (0.5MPa) at three different imposed displacements.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure B.25: In plane displacement and strain fields of the wall 2 (0.5MPa) at three different imposed displacements.</figcaption></figure>

## References

Abdulla, K.F., Cunningham, L.S., Gillie, M., 2017. Simulating masonry wall behaviour using a simplified micro-model approach. Engineering Structures 151, 349–365. URL: https://www.sciencedirect.com
com/science/article/pii/S0141029616314833, doi:https://doi.org/
10.1016/j.engstruct.2017.08.021.

699
700
701
Anthoine, A., 1995. Derivation of the in-plane elastic characteristics of masonry through homogenization theory. International Journal of Solids and Structures 32, 137–163.

702
703
704
705
706
707
Aydenlou, R.M., 2020. Chapter three - types of existing buildings: detailed introduction and seismic rehabilitation, in: Aydenlou, R.M. (Ed.), Seismic Rehabilitation Methods for Existing Buildings. Butterworth-Heinemann, pp. 193–553. URL: https://www.sciencedirect.com/science/article/pii/B9780128199596000038, doi:https://doi.org/10.1016/B978-0-12-819959-6.00003-8.

708
709
710
711
712
Belda, R., Megías, R., Feito, N., Vercher-Martínez, A., Giner, E., 2020. Some practical considerations for compression failure characterization of open-cell polyurethane foams using digital image correlation. Sensors 20. URL: https://www.mdpi.com/1424-8220/20/15/4141, doi:10.3390/s20154141.

713
714
715
Bertetto, A.M., D'Angella, P., Fronterre, M., 2021. Residual strength evaluation of notre dame surviving masonry after the fire. Engineering Failure Analysis 122, 105183.

716
717
718
719
Bisoffi-Sauve, M., 2016. Etude des ouvrages maçonnés en pierre par la méthode des éléments discrets : caractérisation et modélisation du comportement cohésif des joints. Theses. Université de Bordeaux. URL: https://tel.archives-ouvertes.fr/tel-01373072.
Bisoffi-Sauve, M., Morel, S., Dubois, F., 2019. Modelling mixed mode fracture of mortar joints in masonry buildings. *Engineering Structures* 182, 316–330. doi:https://doi.org/10.1016/j.engstruct.2018.11.064.

Cavalagli, N., Cluni, F., Gusella, V., 2011. Strength domain of non-periodic masonry by homogenization in generalized plane state. *European Journal of Mechanics - A/Solids* 30, 113–126. doi:10.1016/j.euromechsol.2010.10.009.

Çelebi, M., Bazzurro, P., Chiaraluce, L., Clemente, P., Decanini, L., DeSorinis, A., Ellsworth, W., Gorini, A., Kalkan, E., Marcucci, S., et al., 2010. Recorded motions of the 6 april 2009 mw 6.3 l’aquila, italy, earthquake and implications for building structural damage: overview. *Earthquake Spectra* 26, 651–684.

Colonna, S., Imperatore, S., Zucconi, M., Ferracuti, B., 2017. Post-seismic damage assessment of a historical masonry building: the case study of a school in teramo, in: *Key Engineering Materials, Trans Tech Publ.* pp. 620–627.

D’Altri, A., Sarhosis, V., Milani, G., Rots, J., Cattari, S., Lagomarsino, S., Sacco, E., Tralli, A., Castellazzi, G., Miranda, S., 2019. Modeling strategies for the computational analysis of unreinforced masonry structures: Review and classification. *Archives of Computational Methods in Engineering* 27. URL: https://link.springer.com/article/10.1007/s11831-019-09351-x, doi:10.1007/s11831-019-09351-x.

D’Altri, A.M., de Miranda, S., 2022. Prediction of flexural drift capacity
in masonry walls through a nonlinear truss-based model. International Journal of Solids and Structures 243, 111593. URL: https://www.sciencedirect.com/science/article/pii/S0020768322001330, doi:10.1016/j.ijsolstr.2022.111593.

Di Nino, S., D'Annibale, F., Luongo, A., 2017. A simple model for damage analysis of a frame-masonry shear-wall system. International Journal of Solids and Structures 129, 119–134.

Di Nino, S., Luongo, A., 2019. A simple homogenized orthotropic model for in-plane analysis of regular masonry walls. International Journal of Solids and Structures 167, 156–169.

Dubois, F., Jean, M., 2006. The non smooth contact dynamic method: recent lmgc90 software developments and application. Analysis and simulation of contact problems, 375–378.

Godio, M., Vanin, F., Zhang, S., Beyer, K., 2019. Quasi-static shear-compression tests on stone masonry walls with plaster: Influence of load history and axial load ratio. Engineering Structures 192, 264–278.

Lourenco, P., 1996. Computational strategies for masonry structures.

Maldonado, N.G., Martín, P., del Solar, G.G., Domizio, M., 2019. Historic masonry, in: Heritage. IntechOpen.

Morenon, P., Sellier, A., Domède, N., 2022. Computational performances optimization of a non-linear mechanical behaviour model for geomaterials. Academic Journal of Civil Engineering 40, 207–210. URL: https://
journal.augc.asso.fr/index.php/ajce/article/view/3115, doi:10.
26168/ajce.40.1.51.

Orduña, A., 2017. Non-linear static analysis of rigid block models for struc-
tural assessment of ancient masonry constructions. International Journal
of Solids and Structures 128, 23–35.

Parent, T., 2015. Méthodologie de Diagnostic de Structures Maçonnées
Anciennes. Theses. Université Toulouse 3. URL: https://tel.archives-ouvertes.fr/tel-01403780.

Parent, T., Brocato, M., Colas, A.S., Domede, N., Dubois, F., Garnier, D., Gros, A., Mindeguia, J.C., Morel, S., Morenon, P., Nougayrede, P., Taforel, P., 2023. A multi-model structural analysis of the vaults of notre-dame de paris cathedral after the 2019 fire and a proposal for a hybrid model merging continuum and discrete approaches. Journal of Cultural Heritage URL: https://www.sciencedirect.com/science/article/pii/S1296207423000626, doi:https://doi.org/10.1016/j.culher.2023.05.009.

Parent, T., Domede, N., Sellier, A., 2017. Multi-scale mechanical behavior of a gothic monument composed of ashlar masonry. application to the design of a reinforcement technique. International Journal of Architectural Heritage 11, 399–414. URL: https://doi.org/10.1080/15583058.2016.1238970, doi:10.1080/15583058.2016.1238970, arXiv:https://doi.org/10.1080/15583058.2016.1238970.

Petracca, M., Pelà, L., Rossi, R., Zaghi, S., Camata, G., Spacone, E., 2017.
788 Micro-scale continuous and discrete numerical models for nonlinear analysis of masonry shear walls. Construction and Building Materials 149, 789 296–314. doi:10.1016/j.conbuildmat.2017.05.130.
790

791 Petry, S., Beyer, K., 2014. Influence of boundary conditions and size effect 792 on the drift capacity of urm walls. Engineering structures 65, 76–88.

793 Petry, S., Beyer, K., 2015. Force-displacement response of in-plane-loaded 794 urm walls with a dominating flexural mode. Earthquake Engineering & 795 Structural Dynamics 44, 2551–2573.

796 Praticò, Y., Ochsendorf, J., Holzer, S., Flatt, R.J., 2020. Post-fire 797 restoration of historic buildings and implications for notre-dame de 798 paris. Nature Materials 19, 817–820. URL: https://doi.org/10.1038/ 799 s41563-020-0748-y, doi:10.1038/s41563-020-0748-y.

800 Rafiee, A., Vinches, M., Bohatier, C., 2008. Application of the NSCD method 801 to analyse the dynamic behaviour of stone arched structures. Interna- 802 tional Journal of Solids and Structures 45, 6269–6283. URL: https:// 803 linkinghub.elsevier.com/retrieve/pii/S0020768308003065, doi:10. 804 1016/j.ijsolstr.2008.07.034.

805 Rezaie, A., Godio, M., Beyer, K., 2020. Experimental investigation of 806 strength, stiffness and drift capacity of rubble stone masonry walls. Con- 807 struction and Building Materials 251, 118972.

808 Roca, P., Cervera, M., Gariup, G., Pelà, L., 2010. Structural analysis of ma- 809 sonry historical constructions. classical and advanced approaches. Archives
810 of Computational Methods in Engineering 17, 299-325. doi:10.1007/s11831-010-9046-1.
811

812 Saretta, Y., Sbrogio, L., Valluzzi, M.R., 2021. Seismic response of masonry buildings in historical centres struck by the 2016 central italy earthquake. calibration of a vulnerability model for strengthened conditions. Construction and Building Materials 299, 123911.
813
814
815

816 Sarhosis, V., Lemos, J., 2018. A detailed micro-modelling approach for the
817 structural analysis of masonry assemblages. Computers & Structures 206,
818 66-81. doi:10.1016/j.compstruc.2018.06.003.

819 Sellier, A., 2018. Anisotropic Damage and Visco-Elasto-Plasticity Applied to
820 Multiphasic Materials. Research Report. LMDC - Laboratoire Matériaux
821 et Durabilité des Constructions de Toulouse ; Université de Toulouse III -
822 Paul Sabatier ; INSA de Toulouse. URL: https://hal.insa-toulouse.
823 fr/hal-01710289.

824 Smoljanović, H., Živaljić, N., Nikolić, , 2013. A combined finite-discrete
element analysis of dry stone masonry structures. Engineering Structures
52, 89–100. doi:10.1016/j.engstruct.2013.02.010.
825
826

827 Smoljanović, H., Živaljić, N., Nikolić, , Munjiza, A., 2018. Numerical
828 analysis of 3D dry-stone masonry structures by combined finite-discrete
element method. International Journal of Solids and Structures 136-
137, 150–167. URL: https://linkinghub.elsevier.com/retrieve/
830 pii/S0020768317305486, doi:10.1016/j.ijsolstr.2017.12.012.
831
832 Tomaževič, M., 2009. Shear resistance of masonry walls and eurocode 6: shear versus tensile strength of masonry. Materials and structures 42, 889-907.
833
834

835 Vanin, F., Zaganelli, D., Penna, A., Beyer, K., 2017. Estimates for the stiffness, strength and drift capacity of stone masonry walls based on 123 quasi-static cyclic tests reported in the literature. Bulletin of Earthquake Engineering 15, 5435-5479.
836
837
838

839 Venzal, V., 2020. Modélisation discrète du comportement mécanique desouvrages maçonnés en pierre. Aspects expérimentaux - Analyse énergétique. Theses. Université de Bordeaux. URL: https://tel.archives-ouvertes.fr/tel-03003557.
840
841
842

843 Venzal, V., Morel, S., Parent, T., Dubois, F., 2020. Frictional cohesive zone model for quasi-brittle fracture: Mixed-mode and coupling between cohesive and frictional behaviors. International Journal of Solids and Structures 198, 17-30. URL: https://www.sciencedirect.com/science/article/pii/S0020768320301402, doi:https://doi.org/10.1016/j.ijsolstr.2020.04.023.
844
845
846
847
848

849 Yaofeng, S., Pang, J.H., 2007. Study of optimal subset size in digital image correlation of speckle pattern images. Optics and Lasers in Engineering 45, 967-974. URL: https://www.sciencedirect.com/science/article/pii/S0143816607000425, doi:https://doi.org/10.1016/j.optlaseng.2007.01.012.
850
851
852
853

854 Zucchini, A., Lourenço, P.B., 2009. A micro-mechanical homogenisation
<footnote n="855">855 model for masonry: Application to shear walls. International Journal of</footnote>

<footnote n="856">856 Solids and Structures 46, 871–886.</footnote>