# On the benefits of the use of complementary numerical models for the structural analysis of a historic masonry monument: Application to a sexpartite vault of Notre-Dame de Paris cathedral

T. Parent, Pierre Morenon, P. Nougayrede, Paul Taforel, Frédéric Dubois, S. Morel

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

## To cite this version:

T. Parent, Pierre Morenon, P. Nougayrede, Paul Taforel, Frédéric Dubois, et al.. On the benefits of the use of complementary numerical models for the structural analysis of a historic masonry monument: Application to a sexpartite vault of Notre-Dame de Paris cathedral. International Journal of Solids and Structures, 2025, 310, pp.113224. 10.1016/j.ijsolstr.2025.113224. hal-04974857

HAL Id: hal-04974857

https://insa-toulouse.hal.science/hal-04974857v1

Submitted on 31 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

Distributed under a Creative Commons Attribution 4.0 International License
# On the benefits of the use of complementary numerical models for the structural analysis of a historic masonry monument : application to a sexpartite Vault of Notre-Dame de Paris Cathedral

Parent T. ¹*, Morenon P. ², Nougayrede P. ²,³, Taforel P. ³, Dubois F.⁴, Morel S. ¹

¹I2M, Université de Bordeaux, UMR5295, Bordeaux, France

²LMDC, Université de Toulouse, UPS, INSA, Toulouse, France

³GSA, Architecture School of Paris-Malaquais, PSL University, Paris, France

⁴MiMeTICS engineering, Montpellier, France

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

*Corresponding Author: thomas.parent@u-bordeaux.fr

## ABSTRACT:

The preservation of a historic masonry monument with structural pathologies requires the implementation of a multidisciplinary diagnostic approach including a mechanical assessment of the structure. To this end, a wide variety of models (such as Continuum Homogeneous Models (CHM), Block-Based Models (BBM) and geometry-based models (GBM)) are available in the literature to simulate the mechanical behavior of this type of structure. In this article, three different models (two BBM and one CHM) are used simultaneously to assess the structural behavior of a choir span of Notre-Dame de Paris Cathedral in Paris, with an aim to assist architects with diagnosis and repair operations within a time frame of approximately one year. The same set of assumptions (loading, geometry, mechanical parameters) was applied to all three models. Using these models simultaneously allowed the research team to build an optimized structural assessment based on constant interaction between the different modeling approaches. This enabled the team to leverage the strengths of each method, ensuring a high degree of confidence in the final results. The comparison of models also revealed potential improvements that could be made to each model, with the goal of reducing calculation times and/or enhancing the consistency of results, particularly in terms of failure mechanisms and behavior laws.

KEY WORDS: Masonry; DEM; FEM; NSCD; Failure analysis; Structural assessment; Damage model; FCZM model; Notre-Dame de Paris cathedral; cross vault
35

# 1. Introduction

36

37

38

39

40

Masonry, an ancient building material presents significant scientific challenges in understanding and modeling its nonlinear mechanical behavior, due to its multi-scale, anisotropic, and heterogeneous composition. Presently, extensive research, both in numerical [Roca, 2010], [d'Altri, 2019] and experimental fields [Sahrosis, 2016], is dedicated to this composite material for contemporary construction as well as the preservation of historical heritage

41

42

43

44

45

46

47

48

49

Beyond the scientific challenge of accurately modeling the mechanical behavior of masonry, the preservation of historical masonry constructions following the guidelines established by organizations like ICOMOS [Icomos, 1964], [Icomos, 2003], particularly those labeled as Patrimonial Monuments, is imperative. The fire of April 15, 2019 at Notre-Dame de Paris (NDP) Cathedral is there to remind us that the masonry heritage is fragile with regard to fire hazards [Praticò, 2020]. A monument with structural pathology requires a diagnostic assessment. It involves a multidisciplinary approach, including archival research, in-situ measurements, use of NDT (Non-Destructive Testing) methods, geometric surveys, mechanical characterization, modeling implementation, and numerical results analysis [Roca, 2010].

50

D'Altri recently proposed a classification of modelling into four strategies [d'Altri, 2019]:

51

52

(i) Block-Based Models (BBM) with a block-to-block description of the masonry considering the mechanical interactions between blocks.

53

54

(ii) Continuum Homogeneous Models (CHM) in which a homogenized masonry material is considered.

55

(iii) Geometry-Based Models (GBM) in which the structure is modelled as a rigid body.

56

57

(iv) Equivalent Frame Models (EFM) in which the structure is idealized into panel-scale structural components.

58

59

60

61

62

63

The CHMs use the finite element method (FEM) combined with non-linear damage models [Faria, 1998], [Berto, 2002], [Sellier, 2013], which have demonstrated effectiveness in assessing Gothic cross-ribbed vaults stability [Parent, 2017] and historical masonry structures under seismic actions [Roca, 2013] [Degli Abbati, 2019]. In contrast, the BBMs consider the actual stereotomy of the masonry structure and the mechanical properties of each constituent; considering the masonry's intrinsic heterogeneity and anisotropy due to joint beds [Malomo, 2023]. BBMs describe
failure and cracking patterns in detail at block-mortar interfaces, employing joint-elements in FEM codes as interface-element-based approach according to d'Altri classification [Lourenço, 1997], [Baraldi, 2017] or alternatively, contact laws in Discrete Element Method DEM codes as contact-based approach [Sahrosis, 2016], [Boukham, 2024], [Ferrante, 2024]. The distinct capabilities of DEM codes models lie in their ability to automatically detect new contact points and allow significant deformations, including complete separation. These characteristics render them especially well-suited for block-by-block simulations with large displacement. Furthermore, Geometry-Based Models (BBM) mainly rely on Limit Analysis-Based Solutions. Heyman was the first to apply limit theorems to masonry structures [Heyman, 1966]. The main hypotheses of this method are: (i) no tensile strength, (ii) infinite compressive strength, and (iii) no sliding between blocks." This approach has been adapted for the formulation of modern computational analysis using static theorem [Block, 2007], [Block 2014], and kinematic theorem [Milani, 2015].

A recent study compared advanced numerical modelling techniques for the mechanical behavior of masonry structures [Schiavoni, 2023]. The study compared two BBMs (using Non-Smooth DEM and Smooth DEM methods) and a finite element continuous model incorporating damage behavior, aiming to model seismic behavior in the Santa Maria Annunziata church. The comparison focused on computation times, crack profiles, and displacements at specific points within the structure. The primary objective of these study was to enhance understanding of the strengths and limitations of different models, and their potential applications in analysing historic masonry structures. The main conclusions underline that each model possesses its own strengths and weaknesses, and that the choice of a calculation method depends on the specific problem. Furthermore, the study highlighted that detailed BBMs may not be relevant if data regarding actions and mechanical properties of the structure are not well-identified.

In the present paper, we propose involving three models (two BBM strategies and one CHM approach) namely (i) a contact-based approach [Dubois, 2006] in a DEM code, (ii) an interface-element-based approach [Pegon, 2001], [Brocato, 2012], [Brocato, 2015] in a FEM code and (iii) a CHM approach with a coupled orthotropic plastic-damage model in a FEM code [Sellier, 2013], respectively. All of these approaches have previously been used by the authors to simulate old masonry structures and to design as well as predict behavior of new structures. In this work, the
models are used to assess the quasi-static mechanical behavior of a span of the Notre-Dame de Paris (NDP) cathedral which is comprised of a sexpartite cross-vault and their flying buttress abutments. This work was carried out in the context of assisting the architects in charge of diagnosis and repair operations over the course of approximately one year. The calculations were made to simulate the behaviour of the structure both before and after the fire. Only the pre-fire calculations are presented here; the post-fire calculations will be published later.

In this time-constrained study, the idea of using the three models mentioned above simultaneously was not to benchmark them, but to take advantage of their complementarity to meet the study deadlines, while ensuring the reliability of the results obtained. Therefore, without going into all the details of the calculation strategy used, the relevance of BBM methods in the modelling of masonry failure mechanisms (location, crack orientation, etc.) is used here in favor of the CHM in order to verify the failure mechanisms obtained by the latter. This is done so that the advantage of the CHM in terms of reduced time consumption during meshing generation and calculation. The different aspects of the structural evaluation, if necessary, can be checked through BBM for further accuracy.

The objective of this article is twofold: (i) to demonstrate the benefits of using various numerical strategies to rapidly and accurately evaluate the stability of a span of the NDP cathedral under quasi-static horizontal loading and (ii) to identify potential improvements to be made to each model to either shorten the calculation times and/or improve the accuracy of the results obtained.

After a brief presentation of the NDP cathedral, the studied span of the choir (architecture and materials) in particular in Part 2, Part 3 of this article addresses the common calculation hypotheses including models, geometries, an original computation methodology called "working point analysis", loading and boundaries conditions. Section 4 presents detailed results related to displacements, failures mechanisms, cracking patterns and behavior law of each structural element. Finally, Part 5 provides a comparison of the models, highlighting the strengths, weaknesses and potential areas of improvement for each model.

## 2. Description of the studied choir span

### 2.1. Architectural elements
The NDP cathedral exemplifies the characteristics of early Gothic architecture in the Paris Basin. Its construction began in the year 1163 AD. The main vessel of the cathedral consists of the western nave and the eastern choir (Figure 1). The intersection of the transept connects the main vessel with both the northern and southern transepts, as illustrated in Figure 1(a). In this study, the X-axis corresponds to the direction transverse to the nave, the Y-axis corresponds to the axis of the nave, and the Z-axis represents the vertical axis.

<figure><img src="image_3.png" /><figcaption>Figure 1: (a) Plan view of NDP cathedral and localisation of the studied span of the choir and (b) sexpartite vault of the choir</figcaption></figure>

The central vessel of the cathedral features sexpartite ribbed vaults (Figure 1(b)), which are supported alternatingly by "strong" piers (receiving the diagonal and lateral transverse rib arches) and "weak" piers (receiving the central transverse rib arches), as depicted in Figure 2. An interior gallery on two floors (the first aisle) and an exterior single-level gallery (the second aisle) alongside the central vessel form the main architectural and structural elements of NDP cathedral. These aisles are covered by quadripartite ribbed vaults. The upper flying buttresses, which are the most slender elements, are supported by the longitudinal side walls of the central portion known as the 'Gouttereau' walls. The "Tas de charge" serves as the connection area between the vaults and the "Gouttereau" walls. The massive buttresses are abutments rising outside the cathedral. They are connected with the upper and lower flying buttresses. The numerical study presented in this paper will focus on the vault, the 'Gouttereau' wall, the upper flying buttress and the buttress.
<figure><img src="image_4.png" /><figcaption>Figure 2: 3D view of a typical span of Notre-Dame de Paris cathedral according to [Choisy, 1899]</figcaption></figure>

140

141

142

## 2.2. In situ materials

143

144

145

146

147

148

149

The ashlar masonry in place consists of blocks of limestone bonded with thin lime mortar of thickness generally <10 mm [Domède, 2024]. An experimental campaign provides approximate values for the mechanical properties of the limestones [Dabat, 2021]. Two types of limestones have been identified by the authors: (i) Soft limestones located principally in the vault webs, exhibiting a porosity of 26-43%, elasticity modulus of 5-17 GPa and a compressive strength of 5-19 MPa; and (ii) Hard limestones located in the ribs, and ‘Goutterau’ walls, exhibiting a porosity of 11-35%, an elasticity modulus of 12-42 GPa and a compressive strength of 9-50 MPa. A large in-situ sonic
auscultation carried out recently on the cathedral confirms the previously presented orders of magnitude for elasticity modulus and strength [Domede, 2024] and shows that the stones constituting the webs of all the cathedral stone vaults (sexpartite and quadripartite) are the softest and lightest, with a compressive strength of about 11 MPa and Young's modulus of about 13 GPa. The stones constituting the buttresses are the hardest and stiffest, with a compressive strength of about 42 MPa and Young's modulus of about 37 GPa.

Furthermore, a recent study had analyzed mortar extracted from the collapsed transverse arch of the nave and from the choir vault [Mechlin, 2022]. It concludes that the mortars were made from a hydrated lime with poor hydraulicity and a fine sand (0/2 to 0/4 mm in size). These elements will be used in selecting the mechanical parameters for the models presented in the following section, related to the calculation assumptions.

## 3. Calculation assumptions

As mentioned in Section 1, the present study was carried out in the context of assisting the architects in charge of diagnostic and repair operations of NDP cathedral. Given the limited time availability for the study during the diagnostic phase, it was decided to take advantage of the complementarity between existing calculation models (in particular between BBM and CHM approaches) to meet the study deadlines while ensuring the reliability of the results obtained. The study strategy therefore consisted of relying on the relevance of BBMs in the modelling of masonry failure mechanisms (location, crack orientation, etc.) to check the results obtained from the Continuum Homogeneous Model and thus take advantage of the CHM, which is less-time consuming in terms of meshing and calculation, to investigate different aspects of the structural evaluation before looping/checking (if necessary) with the BBM. Thus, the final results presented below result for the most part from constant interaction between BBM and CHM models and do not correspond at all to the conditions of a benchmark in which each model is used independently of the others.

### 3.1. *Models*

Three types of calculation models are used simultaneously in this work (Figure 3): (i) a "Discrete element Block-Based Model" noted DEM_BBM, (ii) a "Finite element Block-Based Model"
178 noted **FEM_BBM**, and (iii) a "Finite element Continuum Model" noted **FEM_CHM**. Throughout the
179 article, the three models are consistently represented by the following colours: green for
180 **DEM_BBM**, red for **FEM_BBM**, and blue for **FEM_CHM**. For a detailed description of the behavior
181 laws of the three models used in this study, the reader can refer to the previous work of the authors
182 [Parent, 2024].

183 **Discrete element Block-Based Model: DEM_BBM.** The discrete element approach is
184 implemented using the LMGC90 software [Dubois, 2006]. It follows the concept of hard contact,
185 which means there is no penetration between blocks, and it is solved using the non-smooth contact
186 dynamics (NSCD) algorithm (implicit method) [Dubois, 2018], forming the core of the LMGC90
187 software. The management of the contact between bodies is based on the so-called
188 candidate/antagonist approach [Taforel, 2015]. The 3D model is created by (i) discretizing the
189 masonry on a block-by-block basis; and (ii) applying an interaction law that characterizes the
190 contact between blocks by simulating the behavior of the mortar-block interface. Blocks are
191 considered elastically deformable and modelled by finite elements using extended blocks
192 [Lourenço, 1997]. A homogenized isotropic elastic behavior is used to take into account both the
193 elasticity of the joint and the stone. The frictional cohesive zone model FCZM contact laws
194 governing interactions between blocks describe the linear, damageable, and post-rupture frictional
195 characteristics of joint-blocks interfaces [Bisoffi, 2019], [Venzal, 2020], [Boukham, 2024].
<figure><img src="image_5.png" /><figcaption>Figure 3: The 3 models used: 2 BBM approaches and 1 CHM approach</figcaption></figure>

**Finite element Block-Based Model: FEM_BBM.** The second model is based on a block-to-block discretization of the masonry using the finite element method (Cast3M code). In this model, the stone blocks are represented by homogeneous, isotropic, elastic three-dimensional elements, while the mortar joints are represented by bi-dimensional interface elements exhibiting unilateral and nonlinear behavior [Brocato, 2012], [Brocato, 2015]. This behavior includes (i) preventing penetration (through penalization), (ii) allowing the decohesion (rupture in mode I), and permitting slipping (rupture in mode II) based on a Mohr-Coulomb criterion, with a non-associated flow rule. Young's modulus of the blocks is chosen according to the in-situ characterization, and the stiffness of the joint elements is adjusted to match the homogenized modulus of the masonry.

**Finite element Continuum Model: FEM_CHM.** The third model employs the finite element method (Code_Aster) using a nonlinear model that combines plasticity and orthotropic damage resulting from the induced rotating crack model [Sellier, 2013]. The initial homogenized properties in this model are isotropic. The constitutive law has already demonstrated its effectiveness in the analysis of the quasi-static stability of Gothic cross-ribbed vaults [Parent, 2017]. The law considers
the following aspects : (i) Asymmetry in tensile/compressive behavior of the material, as the tensile strength of masonry is typically significantly lower compared to its compressive strength; (ii) Anisotropic post-peak tensile damage due to softening behavior, leading to a decrease in apparent modulus and residual tensile strength; (iii) Calculation of crack openings resulting from plastic strains, along with the application of the Hillerborg method to regularize the dissipated energy in each principal tensile stress direction; (iv) Compressive-shear damage utilizing a Drucker-Prager criterion that accounts for the increase in shear strength with normal stress; and (v) Ability to reclose tensile cracks in case of unloading (reduction of applied tensile force) or compressive reloading of an element previously damaged in tension and shear-compression dilatancy (volume increase as the peak stress is approached).

The mechanical parameters used for each model are reported in Table 1 for the CHM model and in Table 2 for the BBM model as explained in [Parent, 2024]. They were chosen in order to obtain equivalent mechanical behaviors of the masonry in elastic domain until the peak in tension and compression. Note that in tension, a quasi no-tension approach (tensile strength of 0.1 MPa) was chosen. The homogenized modulus of masonry was chosen at $E = 1.3 \text{ GPa}$, a value derived from Italian codes [IMIT, 2009], [Borri, 2015], which was adopted by the engineers in-charge of the structural assessment of the Cathedral's masonry [DCE5, 2021].

<table><caption>Table 1: Main mechanical properties of the homogenized material (CHM model)</caption><thead><tr><th>Mechanical properties</th><th>Notation</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Bulk Density</td><td>ρ</td><td>1950</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Homogenized Young modulus</td><td>E</td><td>1.3</td><td>GPa</td></tr><tr><td>Poisson ratio</td><td>ν</td><td>0.2</td><td>-</td></tr><tr><td>Masonry compressive strength</td><td>R<sub>c</sub></td><td>10.0</td><td>MPa</td></tr><tr><td>Masonry tensile strength</td><td>R<sub>t</sub></td><td>0.1 (interface)</td><td>MPa</td></tr><tr><td>Tensile Fracture energy</td><td>G<sub>t</sub></td><td>4.5</td><td>J.m<sup>-2</sup></td></tr></tbody></table>

<table><caption>Table 2: Main mechanical properties of the block/mortar interface (BBM models) [Parent, 2024]</caption><thead><tr><th colspan="2">Mechanical properties</th><th>Notation</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="2">Mode I<br>(Normal direction)</td><td>Maximum stress</td><td>σ<sub>I</sub><sup>e</sup></td><td>0.1</td><td>MPa</td></tr><tr><td>Cohesive energy</td><td>Gf<sub>I</sub></td><td>4.5</td><td>J.m<sup>-2</sup></td></tr><tr><td rowspan="2">Mode II<br>(Tangential direction)</td><td>Maximum stress</td><td>σ<sub>II</sub><sup>e</sup></td><td>0.5</td><td>MPa</td></tr><tr><td>Cohesive energy</td><td>Gf<sub>II</sub></td><td>300</td><td>J.m<sup>-2</sup></td></tr><tr><td>Frictional parameter</td><td>Frictional coefficient</td><td>μ</td><td>0.80</td><td>-</td></tr></tbody></table>
231

## 3.2. Working point analysis

232

233

234

235

236

237

Vault structures in churches and cathedrals are not subjected to vertical overload except for their own weight. Indeed, the vaults of the choir central vessel of NDP cathedral must be studied with respect to potential actions on the vessel structure and possibly leading to displacements of the vault piers (wind action, for example). Hence, the utilization of the Working Point Analysis method developed in [Parent, 2017] and explained in detail in [Parent, 2024] becomes essential for the structural evaluation of the sexpartite vaults studied.

238

239

240

241

242

243

244

245

246

247

248

249

The working point methodology involves the mechanical individual assessment for both the vaults and its supports (Figure 4). The supports are defined here as the structural elements designed to receive the horizontal thrust of the vaults. This methodology aims to determine the equilibrium point (or working point), for the combined system of the vault and its supports, under their dead load. In this method, the separation between the vault and its supports is carried out. Two distinct separations (or cut-off) are analysed (Figure 4(b)). In the cut-off n°1, the interface between the two substructures corresponds to (i) the vertical planes of the internal faces of the ‘Gouttereau’ walls, and (ii) the horizontal planes located at mid-height to the vertical interface between the upper flying-buttress and the ‘Gouttereau’ walls. In the cut-off n°2, the interface is defined by a vertical plan located at the external face of the ‘Gouttereau’ wall. This second separation was chosen to group the ‘Gouttereau’ wall and the ‘Tas-de-charge’ with the sexpartite vault in the same substructure.

250

251

252

253

254

255

256

257

258

259

At this interface between the two substructures (vault and support), the resulting horizontal force acting on the vault-support interface is denoted as $\vec{F}$ and its associated displacement is denoted as $\vec{\delta}$. For each substructure, a behavior law is computed using numerical simulation. A mechanical behavior law is defined as the relationship between the horizontal force $F$ and the horizontal displacement $\delta$. Consequently, two curves are obtained independently, each corresponding to one of the substructures (Figure 4(c)). Hence, the working point $(\delta^w, F^w)$ is defined as the intersection point between the two curves and corresponds to the equilibrium point (paired values of horizontal force – associated displacement) at the vault-support interface. Note that this point is expected to be different for “strong” and “weak” piers, due to the specific typology of sexpartite vaults.
<figure><img src="image_6.png" /><figcaption>Figure 4: (a) General view of the cut-off needed for the working point analysis (b) Zoom on the zone of interest for two different cut-off (c) Schematic of the Force-Displacement law expected for the support and the vault.</figcaption></figure>

### 3.3. Geometry

Due to the symmetries of the choir span, as shown in Figure 5(a), and to limit computation times, it has been decided to model only one quarter of the span comprising of one-quarter of the sexpartite vault and half of weak and strong supports. This implies that the applied displacements are also symmetrical. The geometries shown in Figure 5 (b) and (c) exhibit the two sub-structures for the cut-off n°2. In green, the supports are composed of the buttress and the flying buttresses while in blue, the substructure comprises of the quarter sexpartite vault, the “Tas-de-charge” and the ‘Gouttereau’ walls. With regard to cut-off n°1, the only difference concerns the ‘Gouttereau’ wall and the “Tas-de-charge”, which are attributed to the green substructure of the supports (as depicted in Figure 4(a)). To optimize computational time, the first and second aisle levels with their quadripartite vault are not considered in this study. The displacements at the base of ‘Gouttereau’ wall are blocked.

The geometric envelope of the structures has been defined by close fitting to the point-cloud data collected before the fire in a laser-scanner survey carried out by Andrew Tallon
[Sandron, 2013]. This envelope defines the geometry of the model for both CHM and BBM approaches. In the case of CHM, the volume inside the envelope is occupied by a homogeneous, isotropic material. In the case of BBMs, information on the geometry is enriched by a segmentation of the volume into blocks (representing several real blocks) as illustrated in Figure 5(b) and (c). Regarding the ribs of the vaults and the flying buttress, the discretization of the blocks in the model corresponds to the actual stereotomy of the structure. The other parts of the structure (Buttress, webs and gouttereau wall) have been discretized more coarsely, making a compromise that allows for a sufficiently detailed description to capture the elastic and fracture anisotropic behavior of the masonry, but not too detailed to limit computational time. For those structural elements, the width-to-height ratio of the blocks has been preserved, and size ratios between the real blocks and the blocks in the model were chosen approximately as 2 for the abutment, 3 for the vault ribs, and 2 for the 'Gouttereau' walls. Nevertheless, this block-to-block cutting of the geometry matches the in-situ stereotomy logic of the cathedral.

<figure><img src="image_7.png" /><figcaption>Figure 5: Geometry of the vault and its supports (a) Symmetries considered (b) and (c) Geometry of the Span quarter including sexpartite vault with its support for the block-based.</figcaption></figure>

## 3.4. Boundaries conditions and loads

### 3.4.1. Loading sequence
296 **Supports.** The loading sequence applied to the supports in the case of the cut-off n°1 is as

297 follows:

298 • Phase A: Application of the self-weight.

299 • Phase A<sup>bis</sup>: Application of the self-weight of the vault on vault/support interface through a vertical force F<sub>vault</sub> (75 kN and 150 kN for weak and strong piers respectively. These values have been obtained in a preliminary FEM calculation on the vault).

303 • Phase B: Application of a horizontal displacement on the horizontal "Tas-de-Charge" surface shown with red segment in Figure 4(b) and Figure 6. The displacement is applied outward and inward to obtain complete force-displacement laws F<sub>x</sub> vs. δ.

307 The loading sequence applied to the supports in the case of the cut-off n°2 is as follows.

308 • Phase A: Application of the self-weight.

309 • Phase B: Application of a horizontal displacement on the vertical surface located at the head of the buttress shown with red segment in Figure 4(b) and Figure 6. The displacement is applied outward and inward to obtain the complete force-displacement laws F<sub>x</sub> vs. δ.

313 **Vaults.** The loading sequence applied to the ¼ Vault is as follows.

314 • Phase I: Application of the self-weight until the working point of the weak and strong piers corresponding to the mechanical equilibrium point between the vault and its supports is reached. In this phase, the previous force-displacement law of the weak and strong supports is gradually applied on the red segment at the vault-support interface shown in Figure 4(b), simulating a synchronous formwork removal of the vault and the upper flying buttress.

320 • Phase II: From the working point, an identical horizontal outward and inward increment of displacement to both weak and strong interfaces is applied until reaching a vault failure mechanism.

323 **3..4.2. Boundary conditions**

324 **Supports.** The boundary conditions applied to the supports for cut-off n°1 are as follows
(Figure 6):

326 • Buttress base (yellow surface) is clamped (all the displacement are blocked)

327 • 'Gouttereau' wall base (blue surface):

328 ○ Phase A: displacement u<sub>x</sub> and u<sub>y</sub> are blocked, vertical displacement u<sub>z</sub> is progressively imposed downward on this surface to accompany the settlement observed at the top of the buttress when the self-weight of the
structure is applied. This displacement $u_z$ corresponds to the settlement of the top of the buttress at the end of self-weight application, performed in a preliminary calculation.

* Phase A<sup>bis</sup> and B: blocking of displacements $u_x, u_y$ and $u_z$

* Vault/support interface:

* Phase A: blocking of displacements $u_x, u_y$ and $u_z$

* Phase A<sup>bis</sup>: blocking of displacements $u_x$ and $u_y$ and application of the vault self-weight through a vertical force $\overrightarrow{F_{z,vault}}$

* Phase B: Displacement $u_x = \delta$ is imposed outward ($\delta \le 0$) and inward ($\delta \ge 0$) to obtain behavior law supports. $u_z$ is free and $\overrightarrow{F_{vault}}$ is maintained.

In the case of cut-off n°2, the differences lie in the fact that: (i) there is no need to apply vertical displacement at the base of the ‘Gouttereau’ wall nor the self-weight of the vault $\overrightarrow{F_{vault}}$, and (ii) the displacement $\delta$ is directly applied at the head of the buttress arch.

<figure><img src="image_8.png" /><figcaption>Figure 6: Boundary conditions applied on the support in the cases of cut-off n°1 and cut-off n°2</figcaption></figure>

**Vaults.** Figure 7 displays the boundary conditions applied to the vault sub-structure. The vault-support interfaces are shown in red (strong pier) and green (weak pier). Each pier presents two interfaces: (i) the ‘Formeret’ interface representing the vertical interface separation between the lateral web (referred as ‘Formeret’ in French) and the ‘Gouttereau’ wall and (ii) the “Tas-de-charge” interface symbolizing the horizontal separation between the vault and the "Tas-de-charge".
351

In case of the cut-off n°1, the boundary condition applied on the quarter sexpartite vault substructure is as follows:

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

* Lateral and central transverse ribs, as well as longitudinal web ridge, exhibit symmetries by restricting the displacement normal to their surfaces.

* 'Formeret' interface : free

* « Tas-de-charge » interface :

* Phase I : The displacement $u_y$ and $u_z$ are blocked. In the direction x, progressive application of the behavior law support is applied.

* Phase II : The displacement $u_y$ and $u_z$ are blocked. Displacement $u_x = \delta$ is imposed outward ($\delta \le 0$) or inward ($\delta \ge 0$) to obtain behavior laws of the vault on each pier.

The objective of cut-off n° 2 was to shift the boundary conditions of the vault beyond its direct interface shared with the "Tas-de-charge" and the 'Gouttereau' wall. So, in this case, the only differences with the cut-off n°1 lies in the fact that: (i) the transverse blue symmetries are extended to the "Tas-de-charge" and the 'Gouttereau' walls, (ii) the 'Gouttereau' wall base is clamped and (iii) the $u_x = \delta$ displacement is imposed at the lower part interface between the flying buttress and the "Gouttereau" wall displayed in purple on Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7 : Boundary conditions applied on the quarter sexpartite vault substructure</figcaption></figure>

368

369

370

## 4. Numerical modelling Results

371

### 4.1. Supports
The **Figure 8** presents the mechanical behavior laws *F* vs. *δ* obtained by the three simulations for cut-off n°2. In this case, there is no distinction between strong and weak piers in the support laws since the buttress arches and their buttresses are identical for both piers. At the end of phase A, corresponding to the application of the self-weight with *δ* blocked, thrusts of 76 kN, 72 kN, and 77 kN are obtained for the **FEM_CHM**, **DEM_BBM**, and **FEM_BBM** models, respectively (Table 3). A limited damage or cracking at the top of the flying buttress is observed at the end of this loading phase attesting that the structure is well-designed.

<table><caption>Table 3: Thrust at the end of phase A corresponding to the self-weight application</caption><thead><tr><th colspan="2">Thrust</th><th>FEM_CHM</th><th>DEM_BBM</th><th>FEM_BBM</th><th>Mean (COV%)</th></tr></thead><tbody><tr><td rowspan="2">Cut-off n°1</td><td>Strong pier</td><td>113 kN</td><td>103 kN</td><td>109 kN</td><td>108 kN (4,6 %)</td></tr><tr><td>Weak pier</td><td>74 kN</td><td>67 kN</td><td>71 kN</td><td>71 kN (5,0 %)</td></tr><tr><td>Cut-off n°2</td><td></td><td>76 kN</td><td>72 kN</td><td>77 kN</td><td>75 kN (3,5 %)</td></tr></tbody></table>

Phase B of loading reveals an asymmetric behavior of the support. With an outward displacement *δ* applied, the horizontal force *F* strongly increases with respect to the displacement *δ*. The active thrust state of the support is obtained. Failure occurs at the buttress with a diagonal crack (A) at mid-height shown in Figure 9, generated by a combination of normal force and bending moment in the abutment. In this failure mode, the buttress arch is minimally damaged and effectively transmits the thrust forces from the vaults to the buttresses.

When an inward displacement *δ* is imposed, the passive state is achieved. An asymptotic force of about 60 kN is established regardless of the level of horizontal displacement. A triple hinges failure mechanism (associated with cracks (B) in Figure 9) is observed on the upper flying buttress, while the buttress remains undamaged (Figure 9). Note that the three models show similar responses in active and passive states.

A 3D kinematic limit analysis was conducted to assess the magnitude of passive and active thrusts, taking into account the predicted collapse mechanisms of the three approaches under Heyman's hypotheses (no tensile strength, infinite compressive and sliding strength). This analysis was performed using the numerical tool described in Nougayrede, 2024. For both passive and active thrust conditions, the collapse mechanisms considered in the kinematic approach involve three hinges, as illustrated in Figure 10(a). The same block geometry and block density were
considered in this calculation. The magnitudes of the forces associated with these hinges are represented by the green segments in Figure 10(b). The horizontal component of thrust in the active and passive states at hinge n°1 are compared to the asymptotic values of the force-displacement relationship for the support, shown in Figure 8. These values are indicated by black dotted lines. The passive thrust of 49.2 kN and the active thrust of 464.1 kN are of the same order of magnitude as the asymptotic values obtained from the three different models. It is important to note that this type of limit analysis may overestimate these limit values due to the assumption of infinite sliding strength. Indeed, in passive state, the results between the limit analysis and the three models are consistent, as the failure mechanism corresponds to the formation of the three hinges. In contrast, the observed difference between the limit analysis and the three models in active state is related to the fact that the failure peak is more associated with a shear failure of the abutment (diagonal crack (A) in Figure 9) involving sliding between blocks.

<figure><img src="image_10.png" /><figcaption>Figure 8: Force-Displacement law of the support in the case of cut-off n°2</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Failure mechanism of the support in active and passive states: normal crack opening fields</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Limit analysis of the support: (a) considered failure mechanisms with 3 hinges, of the support in active and passive states, (b) associated forces on hinges</figcaption></figure>
Figure 11 illustrates the comparison of the support behavior laws between cut-off n°1 and n°2. Each of the three graphs corresponds to one model. The solid lines depict the behavior law obtained previously in the case of cut-off n°2. The dashed lines represent the distinct responses of the strong and weak piers for cut-off n°1.

<figure><img src="image_13.png" /><figcaption>Figure 11 : Force-Displacement law of the support: comparison of the two cut-off strategies</figcaption></figure>

At the end of phase A corresponding to the self-weight application (including the weight of the vaults), the thrusts reported in Table 3 highlight the asymmetry between the two piers. An average horizontal force of 108 kN and 71 kN are obtained in the strong and weak pier, respectively. This asymmetry is due to the eccentricity with respect to the axis of the ‘Gouttereau’ wall of the volume of the “Tas-de-charge”, which is larger in the strong pier than in the weak pier (Figure 7); as well as the vertical load $\vec{F}_{vault}$ (150 kN on the strong pier and 75 kN on the weak pier). The eccentricities of these two vertical forces increase the bending moment in the ‘Gouttereau’ wall and consequently increase the inward horizontal thrust at the support. It should be noted that this difference in the thrusts between the two piers remains relatively constant
regardless of the applied displacement $\delta$, since the forces and eccentricities involved do not vary significantly with the imposed displacement $\delta$.

It was found that the support behavior laws obtained in the case of cut-off n°1 and n°2 in both passive and active state exhibit a similar overall shape. However, the stiffness of the linear portion and the peak rupture values in active state are higher in cut-off n°1 due to the role of the ‘Gouttereau’ wall in carrying the horizontal forces.

## 4.2. Vault

Figure 12 presents the behavior laws of the strong and weak piers for the vault in cases of cut-off n°1 and cut-off n°2 for the three models. In the cut-off n°1, (Figure 12(a)), the horizontal force corresponds to the support reaction at the base of the vault “Tas de charge” interface surface. In the cut-off n°2, the horizontal force plotted in Figure 12(b) represents the sum of the forces passing through the ‘Formeret’ and the “Tas-de-charge” interfaces.

At the end of phase I corresponding to self-weight application, the working points ($\delta^w, F^w$) for each pier reported in Table 4 are reached. In the case of cut-off n°2, the displacement $\delta$ considered is the one at point A indicated in Figure 7. Similar trends are observed in both cut-offs; the thrust in the strong pier is higher (approximately 110 kN) than in the weak pier (approximately 65 kN). This clearly demonstrates the asymmetric behavior of the sexpartite vault between strong and weak pier due to its geometry.

Furthermore, the horizontal displacements obtained at the end of self-weight application (end of phase I) at the « Tas-de-charge » interface are generally positive in both the strong and weak piers, indicating that the vault's equilibrium point is located inside the central nave. This numerical result is validated by the measurement of the tilts of the internal colonnettes attached to the side walls [BESTREMA, 2021], obtained from the point cloud data recorded before the fire by [Tallon,2013]. This displacement is more pronounced in the weak pier. Further, the inward equilibrium of the nave can be explained by the extreme thinness of the choir's vaults (between 12 and 15 cm), which contributes to limit the thrust of the vaults on their supports. A similar calculation conducted at the level of a nave bay with significantly thicker vaults (between 22 and 25 cm) showed that these equilibrium points were located outside the central vessel. However, it can be noticed that these extremely limited displacements of about 1 mm remain and the quasi-non-
existent damage obtained at the end of this loading phase testify the design knowledge of the gothic builders.

462

<figure><img src="image_14.png" /><figcaption>Figure 12: Force-Displacement law of vault: comparison of the two cut-of strategies.</figcaption></figure>

463

464

<caption>Table 4: Working point of the vaults: Horizontal force F and their corresponding displacement at the end of phase I corresponding to the self-weight application of the vault.</caption>

465

466

<table><thead><tr><th></th><th></th><th></th><th>FEM_CHM</th><th>DEM_BBM</th><th>FEM_BBM</th><th>Mean (COV%)</th></tr></thead><tbody><tr><td rowspan="4">Cut-off<br>n°1</td><td rowspan="2">Strong Pier</td><td>F<sup>w</sup>(kN)</td><td>112</td><td>101</td><td>107</td><td><strong>106 kN (5%)</strong></td></tr><tr><td>δ<sup>w</sup> (mm)</td><td>0,47</td><td>0,18</td><td>0,26</td><td><strong>0,3 mm (45%)</strong></td></tr><tr><td rowspan="2">Weak Pier</td><td>F<sup>w</sup>(kN)</td><td>68</td><td>61</td><td>64</td><td><strong>64 kN (6%)</strong></td></tr><tr><td>δ<sup>w</sup> (mm)</td><td>0,98</td><td>0,63</td><td>0,97</td><td><strong>0,9 mm (27%)</strong></td></tr><tr><td rowspan="4">Cut-off<br>n°2</td><td rowspan="2">Strong Pier</td><td>F<sup>w</sup>(kN)</td><td>108</td><td>128</td><td>105</td><td><strong>114 kN (11%)</strong></td></tr><tr><td>δ<sup>w</sup> (mm)</td><td>-0,45</td><td>0,58</td><td>0,04</td><td><strong>0,05mm</strong></td></tr><tr><td rowspan="2">Weak Pier</td><td>F<sup>w</sup>(kN)</td><td>73</td><td>64</td><td>64</td><td><strong>67 kN (7%)</strong></td></tr><tr><td>δ<sup>w</sup> (mm)</td><td>0,88</td><td>0,61</td><td>0,65</td><td><strong>0,71 mm (20%)</strong></td></tr></tbody></table>

467

468

469

Starting from these working points, inward and outward displacements are applied at the vault-support interface to obtain the behavior laws of the vaults in the strong and weak piers (Figure
470
471

[12]). In parallel, Figure 13 presents, in bottom view, the failure mechanisms in passive and active state where crack openings are plotted on a scale ranging from 0 to 1 mm.

472
473
474
475
476
477
478
479

In the case of cut-off n°1 (Figure 12 (a)), it is first observed that all three models provide very similar responses. Furthermore, the asymmetry in the behavior laws between the passive state and the active state is analogous to what was observed in section 4.1 for the supports: in active thrust, a significant increase in force with displacement is observed. In passive thrust, an asymptotic behavior is reached for relatively small displacement. It should be noted that in the weak pier, the estimated force peak of approximately 80 kN is relatively low compared to the thrust obtained at the working point of 65 kN. This calculation highlights a relatively limited margin of safety for this loading case.

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

480
481
482

<caption>Figure 13: Failure mechanism of the vaults viewed from below in active and passive states: normal crack opening fields</caption>

483
484
485

In the case of cut-off n°2 (Figure 12 (b)), trends in behavior laws similar to those of cut-off n°1 are observed. The failure mechanisms associated with passive and active states can be analyzed in Figure 13. In passive state, two main cracks can be distinguished. Crack (A)
propagates at the ridge line of the central vaults. This crack is typical for a vault experiencing an outward supports displacement. Crack (B), known as the Sabouret crack, concerns the separation of the lateral vault from the 'Gouttereau' wall (represented by the 'Formeret' interface in Figure 7). This type of crack, typical for sexpartite or quadripartite rib-vaults, forms naturally and does not necessarily constitute a pathology compromising the stability of the vault [Sabouret, 1934] [Abraham, 1934]. In active thrust, crack profiles are more diffuse. However, simulations using all three models indicated marked cracking patterns at the base of the vault (area (C) in Figure 13). It should be noted that the analysis limit performed previously on the support to compare the results is difficult to perform on vaults, as the failure mechanisms of these complex 3D structures are hard to identify

When comparing behavior laws of cut-off n°1 and cut-off n°2 in Figure 12, it can be observed that the laws in passive state in the weak and strong piers are relatively similar. This can be explained by the creation of the Sabouret crack in the case of cut-off n°2, which leads to a the boundary condition at the 'Formeret' interface equivalent to the free boundary condition in this interface applied in cut-off n° 1. In active thrust, the laws between the two cut-offs are very different in the weak pier. Cut-off n°2 shows a greater peak force. This is due to the contribution of the 'Formeret' interface, which transmits the compression/shear forces from the lateral web to the 'Gouttereau' wall. On the contrary, this is not the case in cut-off n°1 where this interface is free to move and thus not constrained by the 'Gouttereau' wall.
<figure><img src="image_16.png" /><figcaption>Figure 14 : Force-Displacement law of vault: contribution of ‘Formeret’ and “Tas-de-charge” in the horizontal force bearing, according to the definitions of interfaces in Figure 7</figcaption></figure>

To illustrate this phenomenon, Figure 14 shows, in the case of cut-off n°2, the repartition of the horizontal forces of the vault between the ‘Formeret’ and "Tas-de-charge" interfaces (interfaces shown in Figure 7) in carrying the horizontal forces of the vault in the case of weak pier (Figure 14(a)) and weak piers (Figure 14(b)). For clarity, only the DEM_BBM model is plotted, but the other two models show similar trends. It is confirmed that in both piers, in passive state, the thrusts are quickly and entirely carried by the "Tas-de-charge" interface, which is explained by the formation of the Sabouret crack. Inversely, in active state, the contribution of the ‘Formeret’ interface in carrying the thrusts is significant and even becomes predominant in the weak pier. In conclusion, considering realistic boundary conditions at the ‘Formeret’ interface of sexpartite rib-vaults is essential in cases where the equilibrium points or loadings are directed inward the central nave. In this case, the cut-off n°2 is essential to accurately model the transfer of forces from the lateral web to the ‘Gouttereau’ wall.

## 5. Comparison between the three models
Although the three models were in constant interaction during the study in order to take
advantage of the complementarities between them, the study carried out here allows us to compare
the models in terms of (i) geometry and mesh, (ii) computation time, and (iii) obtained results
(behavior laws and failures mechanisms) and to deduce the potential improvements to be made to
each model in order to reduce calculation times and/or increase the relevance of the results
obtained, as described in detail in this section.

5.1. Geometry and mesh

In the case of complex 3D structures involving double gaussian curvature geometry such as the
sexpartite vaults of NDP cathedral, which are not easily accessible for direct surveys, the
conventional method involves the usage of a point cloud obtained from a 3D laser scan of the
monument to generate the geometry. Based on this survey, two solutions can be used.

The first method involves designing an ideal mathematical geometry and manually fitting it as
closely as possible to the point cloud. This method requires the modeler to have a good
understanding and knowledge of the construction and architectural principles of ancient masonry
structures. Secondly, an interesting method involves the use of algorithms that semi-automatically
generate meshes from point clouds, such as the Cloud2FEM software [Castellazzi, 2015]
[Castellazzi, 2022]. This strategy allows for the generation of a homogeneous finite element mesh
and is not currently suitable for BBMs. Therefore, it has been chosen to implement the first method,
which is better suited to all models proposed in this study.

Fitting the geometric model envelope with the point cloud is common to all three models
and requires a high level of expertise; especially when dealing with geometrically complex
structures. From the envelope, meshing the geometry for the **FEM_CHM** model is relatively direct.
In contrast, block-to-block methods require block segmentation, which can be challenging when
dealing with complex geometry. Segmentation of structures with double curvature requires
discretization methods such as the Chebyshev network [Butzer, 1999]. Additionally, in the
**FEM_BBM** approach, it is necessary to ensure that each block shares common nodes with its
neighbors at the interface. This procedure is essential to automatically generate the bi-dimensional
interface element describing the behavior of the joint between blocks. As for the **DEM_BBM** model,
the candidate points of each block must coincide with their corresponding surface of neighboring
550

blocks. This formalism implies rigor and precision for constructing block-to-block meshes, which
takes time and a considerable level of expertise.

551

552

## 5.2. Computation Time

553

Table 5 outlines the main characteristics necessary to evaluate the computation times involved in numerical simulations of the support under horizontal displacement (Phase B). It should be noted that these calculation times are only indicative. The simulations have been carried out within the time constraints related to the cathedral reconstruction project (1 year) and have not been optimized for each method. Consequently, the time comparison may not be considered as trully objective although trends may be visible. In addition, it is noteworthy that the **FEM_CHM** and **DEM_BBM** methods have an equivalent order of magnitude of number of finite elements, while the **FEM_BBM** method relies on much finer meshes. Additionally, the number of blocks is identical for the two block-to-block methods (266 blocks) as they are based on a common geometry. However, it is worth mentioning that the number of interface elements in the **FEM_BBM** method is approximately ten times higher than the number of candidate points in the **DEM_BBM** method.

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

For each of the three models employed, the computing stations exhibit distinct characteristics, which may impact the objective comparison of computing times. Nevertheless, it is observed that the computing time per number of cores effectively used in the calculation, denoted as $t_{cores}$, is lowest for the **FEM_CHM** method $t_{CHM} = 236$ min. Block-to-block methods consume more computational resources, with computing times of $124 \times t_{CHM}$ for the **DEM_BBM** method and $266 \times t_{CHM}$ for the **FEM_BBM** method. When considering the restitution calculation time, i.e., the real time required to finalize the numerical simulations: $t_{eff} = t_{cores}/$number of cores, it is observed that the computing time of the **FEM_CHM** model remains unchanged as it employs single-core computation. The trend reverses between the two block-to-block models, as the **FEM_BBM** model can parallelize computations more effectively. A total of 12 cores in parallel were used for the **FEM_BBM** method, compared to a maximum of 3 cores for the **DEM_BBM** method. As expected, the **FEM_CHM** method is much faster and less resource intensive than block-to-block models. Those calculation times in terms of CPU time should be interpreted with caution even if the trends thus outlined are consistent.
<table><caption>Table 5: Main characteristics of the time calculations for the support under horizontal displacement (phase B): number of the finite elements and the blocks, characteristics of the scientific computing station and time calculations</caption><thead><tr><th></th><th></th><th>FEM_CHM</th><th>DEM_BBM</th><th>FEM_BBM</th></tr></thead><tbody><tr><td rowspan="2">Volumic Finite Element</td><td>Type</td><td>Tetrahedron/ Hexahedron</td><td>Tetrahedron</td><td>Tetrahedron</td></tr><tr><td>Number</td><td>15 668/ 23 296</td><td>14 056</td><td>581 098</td></tr><tr><td rowspan="2">Block interactions</td><td>Type</td><td>-</td><td>Candidat contact point</td><td>Triangular Joint element</td></tr><tr><td>Number</td><td>-</td><td>3 460</td><td>36 182</td></tr><tr><td colspan="2">Number of blocks</td><td>-</td><td>266</td><td>266</td></tr><tr><td rowspan="2">Characteristics of the scientific computing station</td><td>RAM</td><td>32 Go</td><td>512 Go</td><td>128 Go</td></tr><tr><td>Number of cores (Frequency)</td><td>8 (3,7 GHz)</td><td>72</td><td>28 (2,4 GHz)</td></tr><tr><td colspan="2">Number of cores used</td><td>1</td><td>3</td><td>12</td></tr><tr><td rowspan="2">t<sub>cores</sub><br>Computational time × Number of cores</td><td></td><td>236 min<br>(1 × t<sub>CHM</sub>)<br>(03 h 56 min)</td><td>29 362 min<br>(124 × t<sub>CHM</sub>)<br>(20 d 09 h 22 min)</td><td>28 032 min<br>(277 × t<sub>CHM</sub>)<br>(45 d 10 h 08 min)</td></tr><tr><td>t<sub>eff</sub><br>Restitution calculation time</td><td>236 min<br>(1 × t<sub>CHM,eff</sub>)<br>(03 h 56 min)</td><td>9787 min<br>(41 × t<sub>CHM,eff</sub>)<br>(6 d 19 h 7 min)</td><td>2 336 min<br>(10 × t<sub>CHM,eff</sub>)<br>(1 d 14 h 56 min)</td></tr></tbody></table>

### 5.3. Behavior laws and failures mechanisms

The comparison of results in terms of behavior laws and failure mechanisms reveals the following. At the end of the self-weight application (end of loading Phase A), damages are limited for all models, attesting the knowledge of Gothic builders in designing masonry historical buildings. The horizontal forces obtained are very close between the three models for both the vault and the supports, as evidenced by the coefficients of variation below 10% indicated in Table 3 and Table 4. It can be concluded that for moderate loading rates equivalent to the Serviceability Limit States of the Eurocode (SLS), the three approaches yield very similar results. This is consistent since, in the case where non-linearity or damage are limited, the behavior of the structure is far from the breakup and, in this case, the thrusts depend mainly on (i) the geometric envelope of the structure, (ii) the material's elasticity modulus, and (ii) the applied self-weight. These three assumptions are common to the three proposed approaches. A primary conclusion of this study is that for loadings which maintain the structure far from the breakup, one can prefer the CHM methods, which are easier to implement and with faster computational time.
During the application of horizontal displacements at the vault/support interface (phase B), the three models exhibit cracking patterns at approximately the same locations. More noticeable differences between models arise from the nonlinearities of mechanical behavior laws which are specific to each model. Indeed, it is noticed firstly that in the active thrust of the support (flying buttress), the **FEM_CHM** model exhibits a lower and less ductile active limit thrust than the BBMs. This is possibly explained by the observed difference in cracking profile at the abutment level (Figure 9). The double crack obtained in the **FEM_CHM** model does not correspond to the real stereotomy of the masonries, unlike the BBMs. Moreover, pre-peak damage in compression at the external compressed fiber of the buttress at mid-height also contributes to this behavior difference. Secondly, in the active thrust of the weak pier vault in the case of cut-off n° 2 (Figure 12 (b)), the behavior law of the continuous model shows a lower active thrust than the one obtained by discrete models. The origin of this difference can be explained by (i) the use of a crack reclosure law with a finite energy and a softening behavior in compression-shear in the **FEM_CHM** model that differ of the frictional behavior of cracked joints (infinite energy dissipate by friction) on the BBM and/or (ii) the consideration of the damageable nature of the “Tas de charge” (in the **FEM_CHM**), which tends to crack in the continuous method while it is modeled by a single elastic block in the BBMs

### 5.4. Benefits of using multiple calculation methods in parallel

This work allowed the research team to challenge the calculation results obtained by each method, leading to continuous re-assessment. The specificity and complementarity of each model highlighted a wide variety of phenomena affecting the mechanical behavior of the structure. This finding would not have been possible with the use of a single method. Finally, the relatively similarity of simulated responses obtained from each model enabled the team to validate their results and conclusions, providing a sufficient level of confidence to the architects in charge of diagnostic and repair operations of the cathedral Notre-Dame. Additionally, the comparison of calculation results from this work highlighted the advantages and limitations of each of the three models, thereby identifying issues of improvements and numerical development for each model.

For the **FEM_CHM** model, there is scope for future development in the consideration of the initial masonry anisotropy (in both elasticity and failure criterium). Currently, only the damage is anisotropic. It should be noted that this anisotropy of elasticity and failure is naturally accounted for
in BBMs with the stereotomy description of masonry. Secondly, efforts will be made to describe
friction for large displacements in pre-damaged and cracked elements.

For the **DEM_BBM** model, the calculation times related to automatic contact detection at
each time step, which can be relatively long and may prove prohibitive for large-scale structures.
One improvement option involves using joint finite elements as used in the **FEM_BBM** approach,
as long as the joints remain in the linear domain. This would limit the number of contacts that need
to be automatically detected. Then, interfaces undergoing relative damage could be replaced by
cohesive elements with interfaces to finely simulate failures with automatic detection of new
contacts and large relative displacements between blocks.

For the **FEM_BBM** model, the improvement path would involve implementing an FCZM-
type law [Venzal, 2020] in the joint elements, which would allow for a finer modeling of the post-
failure behavior of mortar joints (taking into account the mix of modes and the coupling between
cohesion and friction in Mode II).

Finally, in BBM models, it is envisaged to integrate the damage model of [Sellier, 2013] to
model the nonlinear behaviour of the block in tension and compression. A hybrid model is currently
under development [Morenon, 2022], [Boukham, 2024] on the LMGC90 calculation code.

## 6. Conclusion

Contrary to a typical benchmarking study wherein the methods being evaluated are used
completely independently of each other, the calculation strategy proposed in this study is to take
advantage of the complementarities between the models. The capabilities of Block-Based Models
(BBM) to describe the failure mechanisms of masonry structures were therefore used first to verify
the results obtained by the Continuous Homogeneous Model (CHM) model. This is then used to
take advantage of the CHM, which is less-time consuming in terms of meshing and calculation, to
investigate different aspects of the structural evaluation before looping/checking with the BBM
based on the necessity. This strategy not only made it possible to carry out the study within the
imposed deadline, but also enabled comparison of the results obtained by the different models and
thus to confer a certain level of confidence to the structural analysis of the sexpartite vault of the
choir of Notre-Dame to assist architects with diagnosis and repair operations.
In addition, comparing the three models on the same case study (i.e., wherever possible, from identical geometry, loading phases, boundary conditions, and material parameters) has helped to highlight several significant points both regarding the mechanical behavior of a typical span of the choir of the Notre-Dame and the comparison of the models. Key findings on the cathedral's mechanical behavior are as follows:

* The flying buttress and its abutment exhibit asymmetric behavior in active and passive thrust.

* The equilibrium thrusts in weak and strong supports are well-distinguished in a sexpartite vault, with approximately a thrust ratio of 3-5 between the two types of supports.

* The horizontal displacement of the vault supports under self-weight in weak and strong supports are located inside the central nave, which is related to the extreme thickness of the choir vaults and correlated to in situ tilts measurements.

* In active thrust of the vault, the 'Formeret' interface between the lateral rib of the vault and the 'Gouttereau' wall bear the majority of the vault's efforts. It is therefore necessary to model it accurately by prioritizing a cut-off at the head of the flying buttress. (cut-off n°2 in Figure 4)

* In passive thrust of the vault, the creation of Sabouret cracks leads to a transmission of thrust between the vault and the supports, mostly through the "Tas de charge".

* The limitation of damage during the application of self-weight to the vault and its supports with minimal displacements of the working points highlights the incredible skill of Gothic builders in designing such complex structures.

The following conclusions can be drawn based on the comparison between the two models

* CHM models are simpler to implement in terms of geometry as they do not require the consuming step of block-by-block partitioning of the structure. Meshing is therefore more direct.

* The restitution calculation time is lower for the continuous method than for the block-based methods, with significant ratios of 1 to 40 between the FEM_CHM and DEM_BBM models. Furthermore, the FEM_BBM model has a lower calculation time than the DEM_BBM model (ratio of 1 to 4), primarily due to its better ability to parallelize computations.

* If the structure remains in a quasi-linear domain, the different approaches yield equivalent results. This is because, in the limitation of damage and non-linearities, the thrusts obtained
primarily depend on three assumptions common to all three models: (i) geometry, (ii) material density, and (iii) material elasticity modulus.

* The localization of damage and cracks is relatively equivalent when CHM method and the BBM methods were used, although the CHM method cannot represent the realistic path of cracks due to its intrinsic continuous nature.

* Beyond the linear domain, during crack propagation, more pronounced differences in behavior laws were observed, which are likely related to distinct dissipative mechanisms between the different modeling approaches. The friction behavior of a damaged crack is simulated through a clack reclosure finite energy in the continuous approach which is not well adapted in modeling the friction in case of important shear displacement experienced by an interface.

The comparison of the calculations highlighted the differences in behavior between models. So, in perspective of this research, it is envisaged to more finely analyze the mechanisms of the models and their associated mechanical parameters that are the source of these differences. In this context, it would be appropriate to conduct a comparison of the methods at smaller scales (scale of the interface in Mode I and Mode II, scale of a Representative Elementary Volume (REV), and then the scale of a small structural element. Furthermore, it is planned as a future work to conduct a sensitivity analysis of (i) the various mechanical parameters of the models (particularly the parameters characterizing the masonry stiffness) and (ii) the level of discretization of the block-by-block models on the structural response.

Moreover, the comparison highlighted the strengths and weaknesses of each of the three models, leading to the identification of improvement points. They are currently being addressed.

* For the **FEM_CHM** model: Development of a model that accounts for elastic and failure anisotropy and better consideration of friction in pre-damaged interfaces.

* For the **DEM_BBM** model: Usage of 2D finite elements to model the interface between blocks for undamaged joints to optimize calculation times.

* For the **FEM_BBM** model: Implementation of an enriched contact law in the joint elements that accounts for damage, mixed-mode behavior, and cohesive-frictional coupling observed in mortar joints.

Finally, a hybrid model is currently under development [Morenon, 2022], [Boukham, 2024].

It involves integrating the damage model of [Sellier, 2013] into the LMGC90 calculation code. The
objective of this model is to simulate both the non-linear behavior observed in the blocks and to optimize calculation at the structural scale; considering some areas with entire substructures in homogenized material (CHM) and others requiring a finer description in BBMs. Finally, a study on the thermo-mechanical behavior of the cathedral vaults under fire loading (April 2019 fire) is underway to evaluate the residual behavior of the fire-damaged sexpartite vaults.

717

718 Acknowledgements

The authors thank the Public Establishment in charge of the conservation and Restoration of Notre-Dame de Paris cathedral (EPRNDP) for its financial support of this study (convention R&D « Evaluation post-incendie des structures de la cathédrale Notre-Dame de Paris »). The authors would also like to thank BESTREMA and STUDIO COMES engineering and the architects in charge of diagnosis and repair operations of Notre-Dame de Paris cathedral for constructive technical discussions relating to the definition of the different hypotheses of the study and to the results obtained.

726 Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this work, the authors used Open AI ChatGPT in order to improve the English language. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication

730 7. Bibliography

[Abraham, 1934] Abraham, P., 1934. Viollet-Le-Duc et le rationalisme médiéval.

[Baraldi, 2017] Baraldi, D., Cecchi, A., 2017. A full 3D rigid block model for the collapse behaviour of masonry walls. European Journal of Mechanics - A/Solids 64, 11–28.
https://doi.org/10.1016/j.euromechsol.2017.01.012

[Berto, 2022] Berto, L., Saetta, A., Scotta, R., Vitaliani, R., 2002. An orthotropic damage model for masonry structures. International Journal for Numerical Methods in Engineering 55, 127–157.
https://doi.org/10.1002/nme.495

[BESTREMA, 2021] Annexe particulière 4 – Relevé des déformations du vaisseau central. Dossier de consultation des entreprise, Maçonneries incendiées, Grand comble, Achèvement hors d'eau, (chœur et nef). Décembre 2021.

[Bisoffi, 2019] Bisoffi-Sauve, M., Morel, S., Dubois, F., 2019. Modelling mixed mode fracture of mortar joints in masonry buildings. puctures 182, 316–330.
https://doi.org/10.1016/j.engstruct.2018.11.064
744 [Block,2007] Block P, Ochsendorf J (2007) Thrust network analysis: a newmethodology for
745 three-dimensional equilibrium. J Int AssocShell Spat Struct 48(3):167-173

746 [Block,2014] Block P, Lachauer L (2014) Three-dimensional (3d) equilibriumanalysis of gothic
747 masonry vaults. Int J Archit Herit8(3):312-335

748 [Borri, 2015] Borri, A., Corradi, M., Castori, G., De Maria, A., 2015. A method for the analysis
749 and classification of historic masonry. Bull Earthquake Eng 13, 2647-2665.
750 https://doi.org/10.1007/s10518-015-9731-4

751 [Boukham,2024] Boukham, A., Venzal, V., Parent, T., Morel, S., Dubois, F., Solbes, B., 2024. 3D
752 hybrid modelling approach combining the finite and discrete element methods: Validation based
753 on masonry shear wall tests. International Journal of Solids and Structures 289, 112638,
754 https://doi-org/10.1016/j.ijsolstr.2023.112638

755 [Brocato, 2012] Brocato, M., Mondardini, L., 2012. A new type of stone dome based on Abeille's
756 bond. International Journal of Solids and Structures 13, 1786-1801.
757 https://doi.org/10.1016/j.ijsolstr.2012.03.036

758 [Brocato, 2015] Brocato, M., Mondardini, L., 2015. Parametric analysis of structures from flat
759 vaults to reciprocal grids. International Journal of Solids and Structures 54, 50-65.
760 https://doi.org/10.1016/j.ijsolstr.2014.11.007

761 [Butzer, 1999] Butzer, P., Jongmans, F., 1999. P. L. Chebyshev (1821-1894): A Guide to his
762 Life and Work. Journal of Approximation Theory 96, 111-138.
763 https://doi.org/10.1006/jath.1998.3289

764 [Castellazzi, 2015] Castellazzi, G., D'Altri, A.M., Bitelli, G., Selvaggi, I., Lambertini, A., 2015. From Laser
765 Scanning to Finite Element Analysis of Complex Buildings by Using a Semi-Automatic Procedure.
766 Sensors 15, 18360-18380. https://doi.org/10.3390/s150818360

767 [Castellazzi, 2022] [Castellazzi, G., Lo Presti, N., D'Altri, A.M., de Miranda, S., 2022. Cloud2FEM: A
768 finite element mesh generator based on point clouds of existing/historical structures. SoftwareX
769 18, 101099. https://doi.org/10.1016/j.softx.2022.101099

770 [Choisy, 1899] Choisy, A., 1899. Histoire de l'architecture. G. Baranger.

771 [Dabat, 2021] T. Dabat, J.-D. Mertz, L. Leroux, Y. Melinge, 2021. Rapport n°1R Voûte
effondrée, murs bahuts et pignons. Caractérisation des carottes de pierre prélevées après
l'incendie du 15 Avril 2019.

774 [D'altri, 2019] D'Altri, A.M., Sarhosis, V., Milani, G., Rots, J., Cattari, S., Lagomarsino, S., Sacco, E.,
775 Tralli, A., Castellazzi, G., de Miranda, S., 2019. Modeling Strategies for the Computational
776 Analysis of Unreinforced Masonry Structures: Review and Classification. Arch Computat Methods
777 Eng. https://doi.org/10.1007/s11831-019-09351-x

778 [DCE, 2021] Dossier de consultation des entreprises, 2021. Maçonneries incendiées, grand comble,
779 achèvement hors d'eau (Chœur et Nef). Annexe particulière 6 Studio Comes – Projet de
780 reconstruction et consolidation des structures

781 [Degli Abbati, 2019] Degli Abbati, S., D'Altri, A.M., Ottonelli, D., Castellazzi, G., Cattari, S., de
782 Miranda, S., Lagomarsino, S., 2019. Seismic assessment of interacting structural units in complex
783 historic masonry constructions by nonlinear static analyses. Computers & Structures 213, 51-71.
784 https://doi.org/10.1016/j.compstruc.2018.12.001

785 [Domede, 2024] Domede, N., Parent, T., Guenser, C., Boukham, A., Morenon, P., Issa-Ibrahim,
786 A.S., 2024. Mechanical characterisation of the stones of Notre-Dame de Paris by in situ acoustic
787 velocity measurement. International Journal of Rock Mechanics and Mining Sciences 175,
788 105671. https://doi.org/10.1016/j.ijrmms.2024.105671
789 [Dubois, 2006] Dubois, F., Jean, M., 2006. The non smooth contact dynamic method: Recent LMGC90 software developments and application. Lecture Notes in Applied and Computational Mechanics 2006, 375-378. https://doi.org/10.1007/3-540-31761-9_44

792 [Dubois, 2018] Dubois, F., Acary, V., Jean, M., 2018. The Contact Dynamics method: A nonsmooth story. Comptes Rendus Mécanique, The legacy of Jean-Jacques Moreau in mechanics / L'héritage de Jean-Jacques Moreau en mécanique 346, 247-262. https://doi.org/10.1016/j.crme.2017.12.009

796 [Faria, 1998] Faria, R., Cervera, M., Oliver, J., 1998. A strain-based plastic viscous-damage model for massive concrete structures. Int. J. Solids and Structures. Erreur ! Référence de lien hypertexte non valide. https://doi.org/10.1016/S0020-7683(97)00119-4

799 [Ferrante, 2024] Ferrante, A., Dubois, F., Morenon, P., 2024. Comparison of Continuous and Discrete Modeling Strategies for the Structural Assessment of a Masonry Vault Under Dynamic Seismic Loading. International Journal of Architectural Heritage 0, 1-13. https://doi.org/10.1080/15583058.2024.2377297

803 [Heyman, 1966] Heyman J (1966) The stone skeleton. Int J Solids Struct 2(2):249-279

804 [Icomos, 1964] ICOMOS, 1964. International Charter for the conservation and restoration of monuments and sites.

806 [Icomos, 2003] Recommendations for the analysis, conservation and structural restoration of architectural heritage.

808 [IMIT, 2009] Circ. 02.02.2009, n. 617: Istruzioni per l'applicazione delle Nuove Norme Tecniche per leCostruzioni di cui al decreto ministeriale 14 gennaio 2008. Italian Ministry of Infrastructures and Transportation, Rome, Italy (in Italian)

811 [Lourenço, 1997] Lourenço P.B., Rots J.G., 1997. Multisurface Interface Model for Analysis of Masonry Structures. Journal of Engineering Mechanics 123, 660-668. https://doi.org/10.1061/(ASCE)0733-9399(1997)123:7(660)Schiavoni, M., Giordano, E., Roscini, F., Clementi, F., 2023. Numerical modeling of a majestic masonry structure: A comparison of advanced techniques. Engineering Failure Analysis 149, 107293. https://doi.org/10.1016/j.engfailanal.2023.107293

817 [Malomo, 2023] Malomo, D., Pulatsu, B., 2024. Discontinuum models for the structural and seismic assessment of unreinforced masonry structures: a critical appraisal. Structures 62, 106108. https://doi.org/10.1016/j.istruc.2024.106108

820 [Mechling, 2022] Mechling, J.-M., Baquet, M., 2022. Physicochemical characterization of the mortars from Notre-Dame de Paris cathedral - Preliminary findings and reconstitution of the fresh mix. Journal of Cultural Heritage. https://doi.org/10.1016/j.culher.2022.10.006

823 [Milani, 2015] Milani G (2015) Upper bound sequential linear programmingmesh adaptation scheme for collapse analysis of masonry vaults. Adv Eng Softw 79:91-110

825 [Morenon, 2022] 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.

828 [Nougayrede, 2024] Nougayrede, P., Ciblac, T., 2024, July. Stabilithos: a numeric tool for general limit analysis of stone structures. In 18th International Brick and Block Masonry Conference.

830 [Parent, 2017] 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.

833 [Parent, 2024] 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., 2024. A multi-model
835 structural analysis of the vaults of Notre-Dame de Paris Cathedral after the 2019 fire and a
836 proposal for a hybrid model merging continuum and discrete approaches. Journal of Cultural
837 Heritage. https://doi.org/10.1016/j.culher.2023.05.009

838 [Pegon, 2001] Pegon, P., Pinto, A.V., Géradin, M., 2001. Numerical modelling of stone-block
839 monumental structures. Computers and Structures 22–25, 2165–2181.

840 [Praticò, 2020] Praticò, Y., Ochsendorf, J., Holzer, S., Flatt, R.J., 2020. Post-fire restoration of
841 historic buildings and implications for Notre-Dame de Paris. Nat. Mater. 19, 817–820.
842 https://doi.org/10.1038/s41563-020-0748-y

843 [Roca, 2010] Roca, P., Cervera, M., Gariup, G., Pela, L., 2010. Structural Analysis of Masonry
844 Historical Constructions. Classical and Advanced Approaches. Archives of Computational
845 Methods in Engineering 17, 299–325.

846 [Roca, 2013] Roca, P., Cervera, M., Pelà, L., Clemente, R., Chiumenti, M., 2013. Continuum
847 FE models for the analysis of Mallorca Cathedral. Engineering Structures 46, 653–670.
848 https://doi.org/10.1016/j.engstruct.2012.08.005

849 [Sahrosis, 2016] Sarhosis, V., De Santis, S., Felice, G., 2016. A review of experimental
850 investigations and assessment methods for masonry arch bridges. Structure and Infrastructure
851 Engineering 12, 1439–1464. https://doi.org/10.1080/15732479.2015.1136655

852 [Sellier, 2013] A. Sellier, G. Casaux-Ginestet, L. Buffo-Lacarrière, and X. Bourbon, "Orthotropic
853 damage coupled with localized crack reclosure processing. Part I: Constitutive laws," Engineering
854 Fracture Mechanics, vol. 97, no. 0, pp. 148–167, Jan. 2013.

855 [Schiavoni, 2023] Schiavoni, M., Giordano, E., Roscini, F., Clementi, F., 2023. Numerical modeling
856 of a majestic masonry structure: A comparison of advanced techniques. Engineering Failure
857 Analysis 149, 107293. https://doi.org/10.1016/j.engfailanal.2023.107293

858 [Taforel, 2015] Taforel, M. Renouf, F. Dubois, C. Voivret, "Finite Element-Discrete Element
859 Coupling Strategies for the Modelling of Ballast-Soil Interaction", International Journal of Railway
860 Technology, 4(2), 73-95, 2015. doi:10.4203/ijrt.4.2.4

861 [Sabouret, 1934] Sabouret, V., 1934. Les voûtes d'arêtes nervurées, Rôle simplement décoratif
862 des nervures 205–209.

863 [Sandron, 2013] Sandron, D., Tallon, A., 2013. Notre-Dame de Paris. Parigramme, Paris.

864 [Sarhosis, 2016] Sarhosis, V., Bagi, K., Lemos, J., Milani, G., 2016. Computational Modeling of
865 Masonry Structures Using the Discrete Element Method. https://doi.org/10.4018/978-1-5225-
866 0231-9

867 [Venzal, 2020] Venzal, V., Morel, S., Parent, T., Dubois, F., 2020. Frictional cohesive zone
868 model for quasi-brittle fracture: Mixed-mode and coupling between cohesive and frictional
869 behaviors. International Journal of Solids and Structures 198, 17–30.
870 https://doi.org/10.1016/j.ijsolstr.2020.04.023