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

HAL
open science

# Bayesian interface based calibration of a novel rockfall protection structure modelled in the Non-smooth contact dynamics framework

Ritesh Gupta, Franck Bourrier, Vincent Acary, Stéphane Lambert

► To cite this version:

Ritesh Gupta, Franck Bourrier, Vincent Acary, Stéphane Lambert. Bayesian interface based calibration of a novel rockfall protection structure modelled in the Non-smooth contact dynamics framework. Engineering Structures, 2023, 297, pp.116936. 10.1016/j.engstruct.2023.116936. hal-04039263v2

HAL Id: hal-04039263

https://hal.science/hal-04039263v2

Submitted on 13 Oct 2023

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>

Distributed under a Creative Commons Attribution 4.0 International License
Bayesian interface based calibration of a novel rockfall protection
structure modelled in the non-smooth contact dynamics framework

Ritesh Gupta<sup>a,*</sup>, Franck Bourrier<sup>a,b</sup>, Vincent Acary<sup>b</sup>, Stéphane Lambert<sup>a</sup>

<sup>a</sup>Univ. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, 38000 Grenoble, France

<sup>b</sup>Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000 Grenoble, France

**Abstract**

This article presents the development and calibration of a numerical model simulating the response of a
novel rockfall protection structure subjected to localised dynamic loading. This structure is made of piled-
up concrete blocks interconnected via metallic components whose dynamics response under projectile
impact is examined via real-scale experiments. The corresponding numerical model is developed in
a python based open source software SICONOS which implements the Non-Smooth Contact Dynamics
(NSCD) method. The geometrical features and mechanical properties are incorporated in the model via
specific developments pertinent to the modelling requirements in SICONOS. Some parameters peculiar to
the numerical model are calibrated against the spatial-temporal measurements from two full-scale impact
experiments. The Bayesian interface statistical learning method aided by the polynomial chaos expansion
based meta-model of the NSCD model is deployed for the calibration. The additional understanding of
the model dynamics through the byproducts of the meta-model is highlighted. In the end, the NSCD
model is successfully calibrated against the spatial-temporal response of the experimental structure with
more than 90% accuracy for impact energies up to 1 MJ.

*Keywords:* rockfall, impact, meta-model, statistical learning, NSCD modelling, calibration, Bayesian
interface

# 1. Introduction

The Non-Smooth Contact Dynamics (NSCD) method developed by M. Jean and J.J. Moreau [1, 2]

is a modelling and numerical framework for the dynamics of multi-body systems in interaction through

frictional contact interfaces. The present work extends the use of NSCD modelling [3, 4] to study the

response to impacts of a complex structure made of individual concrete blocks interconnected via metallic

<footnote>*Corresponding author</footnote>

<footnote>Email address: ritesh.gupta@inrae.fr (Ritesh Gupta)</footnote>
components, thus forming an articulated structure. The structure is proposed by the company Géolithe and is intended to serve as passive protection against gravity-driven natural hazards. It is in particular designed to intercept rockfall, similarly to other massive structure types including embankments [5], concrete walls [6], concrete blocks cushioned with gabions [7], gabion structures [8] and sea containers [9]. This innovative technology offers the possibility to build massive vertical walls, with reduced footprint, high deformability, and versatility.

The development of this technology is based on small-scale experiments, to address the structure response considering different design options, and real-scale impact experiments to demonstrate the structure capacity in arresting rock blocks with kinetic energy exceeding 1 MJ [10]. These experiments revealed the intricacy of the impact response of structures made from interconnected concrete blocks. The impact induces sliding at the base and tilting of the structure, where the amplitude depends on the impact energy and distance to the impact point. In addition, the interconnection between concrete blocks improves the structure stability preventing excessive concrete block displacements in the impact vicinity.

The impact response of small- and real-scale structures has been previously investigated in the numerical modelling framework of the finite volume formulation code FLAC3D [10]. These models provided detailed insights into the energy dissipation capacities of the protection walls. However, the computation time for the real-scale structures ranging from 10 to 20 hours hinders the exhaustive investigations of the structure's mechanical and dynamic response. In particular, improving the design of such structures requires better quantifying their efficiency when exposed to rockfall. In line with some previous research on flexible barriers (e.g., [11, 12, 13]), this implies taking into account the variety of loading cases observed during real rockfalls and therefore performing a very large number of numerical simulations of the structure's response to the impact.

Given the above, we developed an NSCD-based model of the structure, using the SICONOS software [14]. SICONOS is an open-source scientific software primarily targeted at modelling and simulating non-smooth dynamical systems in C++ and in Python, including mechanical systems, switched electrical circuits, sliding mode control systems, and biology. The NSCD method guarantees that Signorini's condition at the velocity level is satisfied without introducing contact stiffnesses and that the coulomb friction, especially in the sticking mode is satisfied without resorting to viscous friction. This peculiarity of NSCD is enforced with the use of a specific implicit scheme (Moreau-Jean scheme) and yields a robust
numerical method, which consistently models the threshold phenomena (friction, contact) and the dissipation properties of the model in discrete time, in particular impact dissipation and energy properties [15].

Considering its intended use, the developed model is kept as simple as possible while allowing satisfactory mimicking of the whole structure's response to impact, as observed during the real-scale experiments. Such simplifications drastically reduce the computational time as each component of the structure is modelled as a rigid body instead of a deformable mesh resulting in a huge reduction in the degrees of freedom of the model. The energy dissipation due to plastic strains and breakages of the blocks during impact is modelled in a simple manner as energy dissipation at the contact points. Escallón et al. [16] reported a conceptually similar approach for the FEM model of flexible rockfall protection barriers and observed a significant reduction in the model computational time.

Further, the calibration of model parameters becomes necessary to validate the model performance. For this, an existing experimental database [10] is used as a reference and some identified model parameters are calibrated using the Bayesian interface statistical learning method.

A common requirement of any stochastic method is a large number of model computations which are impractical to implement due to a large computation cost. This gap is bridged by the use of a meta-model, representative of the actual model having a negligible computational cost and possessing high accuracy and reliability. Such meta-models have been used by many researchers, e.g., for dam engineering problems [17, 18] and for nuclear containment structures [19].

In the present work, the NSCD-based model is first described in detail, including specific developments to cope with the particular features of this structure type. A polynomial chaos expansion (PCE) based meta-model of the NSCD model is created for the stochastic analysis. The additional byproduct of PCE-based meta-models towards the relative influence of model parameters is discussed. The model parameters are calibrated based on the Bayesian interface approach considering the measurements obtained from two impact experiments. A discussion on model features in line with the perspectives for future research concludes the work.

## 2. Experimental structure and its impact response

Referring to [10] for further information, the experimental details use for the calibration of the NSCD
model are summarised hereinafter. The considered articulated structures consist of piled-up blocks, made
from concrete, and are reinforced with steel rebars, presented in Figure 1. The block extremities in the
horizontal plane are rounded and their upper and lower faces present empty spaces. The wall is made
up of 38 blocks and four half-blocks stacked in four rows. The blocks are arranged following a zig-zag
pattern (at 45° angle) to improve the structure stability against titling. The wall is 3.2 m in height and
about 14 m in length.

<figure><img src="image_3.png" /><figcaption>Figure 1: Experimental full-scale structure, recreated after [10]</figcaption></figure>

Each block is traversed by two cylindrical holes of diameter 154mm along the vertical axis to receive tubes and slings (see Figure 1(iv)). These latter metallic components, hereafter referred to as connectors, link blocks together. More precisely, a hollow steel tube of diameter 139.7mm connects each couple of adjacent layer blocks, preventing their relative lateral displacement. In addition, a cable runs from the structure base to its top, in the succession of tubes, and mainly prevents relative displacements along the vertical direction. From a conceptual point of view, these connectors provide the structure with mechanical continuity with the aim of increasing the number of blocks involved in the structure's response proportional to the impact loading. Another key feature in the design is the presence of various free spaces (hereby referred to as 'mechanical plays' or simply 'plays'): a 40 mm side-to-side distance between the same layer concrete blocks, about 14 mm difference between the external diameter of the metallic tubes and the inner diameter of the cylindrical holes and, a precisely non-quantifiable slack in the vertical cables. These plays give the structure a certain deformation capacity. All geometrical and mechanical characteristics of this structure are presented in Table 1 along with the ones of the impacting projectile.
<table><caption>Table 1: Geometrical and mechanical characteristics of the articulated concrete blocks structure</caption><thead><tr><th>Component</th><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="2">Ine</td><td>Length</td><td>1.56</td><td>m</td></tr><tr><td>Width</td><td>0.76</td><td>m</td></tr><tr><td rowspan="3">Concrete Block</td><td>Height</td><td>0.80</td><td>m</td></tr><tr><td>Mass (m<sub>block</sub>)</td><td>1800</td><td>kg</td></tr><tr><td>Hole diameter</td><td>154</td><td>mm</td></tr><tr><td>Ine</td><td>Diameter (ext.)</td><td>139.7</td><td>mm</td></tr><tr><td rowspan="2">Cylindrical bars</td><td>Diameter (int.)</td><td>123.7</td><td>mm</td></tr><tr><td>Length</td><td>0.8 (3 pcs.)</td><td>m</td></tr><tr><td></td><td></td><td>0.4 (2 pcs.)</td><td>m</td></tr><tr><td>Ine</td><td>Length (edge-to-edge)</td><td>≈14</td><td>m</td></tr><tr><td rowspan="2">Structure</td><td>Height</td><td>3.2</td><td>m</td></tr><tr><td>Pattern angle</td><td>45 (from y-axis)</td><td>°</td></tr><tr><td rowspan="4">Ine<br>Projectile</td><td>Side length (l<sub>proj</sub>)</td><td>1.1</td><td>m</td></tr><tr><td>Mass (m<sub>proj</sub>)</td><td>2600</td><td>kg</td></tr><tr><td>Impact angle</td><td>0</td><td>°</td></tr><tr><td>Impact point location</td><td>y ≈ 7.0, z = 1.7</td><td>m</td></tr></tbody></table>

Two identical real-scale walls were built and subjected to one impact each at different energy levels.

The impact experiments were performed with the pendulum testing facility of the Université Gustave Eiffel test site (Montagnole, France). It involved a 2600 kg reinforced concrete projectile whose shape is in accordance with the requirements for flexible barriers testing [20]. Its size is about one-third the structure's height (1.1 m). The velocity at impact is 20m/s and 28m/s during the first and second impacts, corresponding to kinetic energies of 520 kJ and 1020 kJ respectively. The angle of incidence and rotation velocity at impact were zero. The projectile impacted the structure at its mid-length and at about 1.7m from the ground. These impacts in particular induced wall sliding at its base, tilting, change
in conformation as well as concrete block damage.

The structure impact response with time is described based on measurements of displacement and acceleration along two vertical lines, as presented in Figure 2. They are termed ‘impact’ and ‘distant’ axis which are representative of the impact point and 3.5m to its left respectively. The experimental data used in this study concerned the displacement of concrete blocks at the *Top* and *Base* in the impact axis and points *C* and *D* in the distant axis. Besides, the data from two points *A* and *B* in the impact axis is also available which is used for the cross-reference purpose only.

<figure><img src="image_4.png" /><figcaption>Figure 2: Data acquisition locations in the experimented full-scale structure, recreated after [10]</figcaption></figure>

The displacement governs the structure sliding and tilting, both leading to structure failure if in excess. In this study, the database of displacement evolution with time at the impact vicinity and at distance is considered. The displacement data used in this study is plotted in Figure 3 for both impact and distant axes. A secondary y-axis is added for the distant axis albeit to the relatively low displacement amplitude compared to the impact axis.

Measures along the impact axis are derived from video records while measures from cable extensometers are considered along the distant axis. Due to the absence of this later data for the 520 kJ impact test at point *C*, the evolution obtained from simulations presented in [10] are considered. This exception is thought to have a minor influence due to the very small amplitude in displacement observed in this specific case. The dynamic response of the structure is clearly visible through the relative initiation time lag between observations corresponding to impact and distant axes. Besides, for the 520kJ test, nearly equal displacement magnitudes are observed at the end of the test for the points in the same axis, indirectly presenting nearly no-permanent post-impact tilt in the structure and small relative displacement between adjacent layer blocks. However, for the 1020kJ test, a significant difference is observed for both axes, inferring a combination of permanent tilt and relative displacement between the same axis blocks.
<figure><img src="image_5.png" /><figcaption>Figure 3: Post-processed experimental database for the time evolution of displacement for (a) 520 kJ and (b) 1020kJ impact test</figcaption></figure>

The observed response at the point *Top* for the 1020kJ test is an exception due to the severe damage in concrete blocks in the impacted vicinity, up to fracture. This resulted in the loss of mechanical continuity within the structure, which was thought to cause a relatively larger block displacement at the wall top than only due to the combination of a permanent tilt and relative block movements in the structure.

## 3. Numerical model

The rockfall protection structure is modelled in a python-based software package named SICONOS, which implements the NSCD method with complex geometries [14]. In this section, the NSCD framework is first introduced, followed by the numerical model development and finally the model response evaluation aiming to reproduce the whole structural displacement response upon impact by a projectile.

### 3.1. NSCD main principles and relevance

The NSCD method was developed to solve multi-body multi-contact problems with rigid and/or deformable bodies [1, 2, 3] with rigid contact laws such as the Signorini's model of unilateral contact and the Coulomb's dry friction without any kind of regularization (viscous friction) or compliance. It finds many modelling applications in the field of civil engineering and geomechanics such as; masonry and stone structures exposed to static, cyclic, and dynamic loading [21]; cohesive and non-cohesive granular materials [22]; standard Finite Element Method (FEM) from the early stage of its development [23]; rockfall propagation on slopes [24] in particular and mechanisms and linkages [25].
In this study, the system is modelled as a collection of rigid blocks connected by unilateral constraints with Coulomb friction. The finite-freedom dynamics of rigid bodies with unilateral constraints is known to be non-smooth, in the sense that the velocities of the system possess jumps when a contact is closing with a positive relative velocity. These velocity jumps are described through the introduction of an impact law. Thereby, the NSCD method is able to perform the numerical time integration of the multi-body system under impacts (of a projectile, in our case) in a dynamic condition. Some details and notation on the NSCD method are given in Appendix A to support the following discussion.

Thus, with the NSCD method, robust and efficient simulations of the dynamics and the quasi-statics of multi-body systems with contact, impact, and Coulomb friction are possible. For each body, the Newton-Euler equations are discretized in time using an implicit scheme:

$$M(v_{k+1} - v_k) - hF_{k+\theta} = G^\top(q_{k+1})P_{k+1} \quad (1)$$

here, $M$ is the inertia matrix, $v$ is the body twist (translational and rotational velocity) vector, $q$ is representative of body position and orientation and $G$ is the Jacobian of the gap function (see Appendix A.2 for details), $h$ is the timestep. $F_{k+\theta}$ refers to the external forces applied to the bodies (limited to gravity forces, in this context) except contact interactions. $P_{k+1}$ refers to all contact interactions that are modelled as percussions.

The percussions $P_{k+1}$ allow modelling impact involving Coulomb friction in the context of the study. The percussion calculations are governed by the user-defined friction ($\mu$) and restitution ($e$) coefficients. The tangential component of percussion calculations reflects the energy dissipation due to friction and the normal component accounts for energy dissipation due to inelastic collision which, in other words, represents the damage in the model body.

Notably, in this work, the restitution coefficient of magnitude less than one (i.e., inelastic collision) globally accounts for the energy dissipation due to plasticity and fracture in the structure. Nonetheless, a complex representation of plasticity, damage, and fracture processes is possible in the NSCD framework. For example, modelling blocks as deformable bodies using Finite Element based approaches involving block plasticity models (e.g., [26]) and cohesive zone models in-between the finite elements (e.g., [27]) where fracture could occur. However, such a complex block model shall substantially increase the computational duration and hence it is not adapted in the present work.
The non-impulse terms of the dynamics are approximated with a $\theta$-method. For $\theta = 1$, the scheme is fully implicit and maximizes the numerical dissipation. It is therefore well adapted to quasi-static evolutions, or for efficiently reaching a static equilibrium with large time-steps. For $\theta = 1/2$, the mechanical energy (the sum of the kinetic and potential energy) is conserved. Furthermore, the dissipation is always positive and is given by a second order approximation of the actual dissipation. These energy and dissipation properties detailed in [15] render the scheme robust and stable with quite large time-steps, and are of particular interest for the analysis of the dissipative processes in protection structures. These are the main advantages with respect to classical smooth DEM approaches [28, 29] where the use of explicit scheme renders difficult the energy analysis and the contact compliance generates spurious oscillations that prevent to obtain quickly static equilibrium without artificial damping. With a non-smooth dynamics, direct higher order approximations are possible, but some care has to be taken. The use of standard higher scheme — Newmark, HHT, generalised-$\alpha$, implicit Runge-Kutta to name a few — is not possible and we need to rely on specific schemes such as non-smooth Newmark or generalised-$\alpha$ schemes [30, 31, 32], or time discontinuous Galerkin methods [33, 34].

## 3.2. NSCD structure model development in SICONOS

By contrast with common applications of NSCD, new challenges arise when modelling such articulated-concrete-blocks structures. These relate to the shape of the concrete blocks, the presence of structural components passing through these blocks, the interaction between these later components, and the modelling of the connectors.

In SICONOS, each object (i.e., block and connector) of the wall is modelled as an assembly of many rigid body components with simple geometrical shapes (referred to as 'primitive shapes') for contact detection. These shapes present computational time efficiency and accuracy in contact detection in comparison to a typical numerical mesh. The interaction between different objects is user-controlled through so-called 'contactors', inferring that in its absence, two physical objects placed adjacently, would be invisible to each other. The model configuration and its components are shown in Figure 4. The modelling strategy that is developed to cope with these specific features is presented hereafter.

### 3.2.1. Blocks modelling

The concrete blocks geometry feature from the experimental structure is presented in Figure 5, inclusive of the effective contact area between two adjacent layer blocks in the wall. The block's surface
<figure><img src="image_6.png" /><figcaption>Figure 4: The NSCD model of the articulated concrete block structure and its components</figcaption></figure>

comprising of the upper, lower and lateral side faces, is modelled as a combination of ‘primitive’ shapes.

<figure><img src="image_7.png" /><figcaption>Figure 5: Experimental concrete block design and effective contact surfaces between adjacent layer blocks</figcaption></figure>

Different combinations of primitive shapes are considered for geometry and contact area relevance with the experimental block, and for numerical efficiency and computational accuracy in the NSCD framework. This resulted in the best block model option, presented in Figure 6.

The user-controlled contactor definition feature enables to place cylinder and box superimposed to collectively constitute the block model. Subsequently, the interaction between layered concrete blocks in the wall is enabled through cylinder primitives and interaction between the projectile and the wall through both cylinder and cuboidal box primitives, detailed later in the section 3.2.4.

Several contact points are set along the contact surface to model the interaction between adjacent layer blocks. The location of the contact point is defined in order to maximize the contact surface between the interacting bodies. Notably, as the friction and restitution coefficients control the interaction between bodies, the difference in the physical contact area between adjacent layer numerical blocks being different
<figure><img src="image_8.png" /><figcaption>Figure 6: Block design in the NSCD model</figcaption></figure>

from the experimental counterpart becomes irrelevant.

The block model is assigned with the same mass as of the real block and the inertia of the model
blocks is obtained through a 2D surface mesh using the Convex Hull [35] feature in SICONOS. Notably,
modelling the block cylindrical holes is ignored and a specific strategy is developed for accounting for the
interaction between the block and the connectors, detailed in the next section.

3.2.2. Connectors modelling

The connector between blocks (see Figure 1.iv), are modelled as a single component accounting for the
influence on the structure response of both the sling and tubes. This required defining a mechanical model
governing the connector-block interaction (which is insured by the tubes) and the interaction between
connectors (which is in particular controlled by the slings). The collective connector model illustrative
profile is presented in Figure 7 for before and after impact situations.

It consists of four cylindrical primitive shape components, named as ‘bars’, at locations corresponding
to the real-structure tube positions and possessing the same length (3 bars with 0.8m and 1 bar with
0.4m), diameter, inertia and mass as the tube. Notably, a fifth bar is added at the top block having
similar geometrical characteristics as the short-length bar, however with a negligible mass. This is done
to complete the configuration and the bounding limits of the sling in the real structure.

The mechanical plays existing in the structure, in both the horizontal and vertical directions, are
accounted for in the model. The horizontal play in the real structure results from the difference between
the external diameter of the tube (139.7mm) and the diameter of the hollow cylinder in the block (154mm).
The vertical play in the real structure results from the relative looseness in the sling. The detailed
connector modelling strategy to accommodate horizontal and vertical plays is presented hereafter.

The user control for assigning the contact detection enables, the connector bars passing through the
<figure><img src="image_9.png" /><figcaption>Figure 7: Schematic representation of the deformation of the structure upon projectile impact showing plays in action</figcaption></figure>

concrete blocks without a direct interaction between them. However, in reality, the steel bars inserted in cylindrical block holes, interact with the concrete blocks. Hence, this interaction is assigned through additional components named as hollow circular disks added to the block model object, as illustrated in Figure 8. These disks are modelled as 2D mesh and have the same internal diameter as that of the hole in the real block, hence ensuring the presence of the requisite horizontal play of about 7mm. Each block is equipped with four hollow disks meaning that each of the two bars passing through a block interacts with this block via two disks. The position of these disks with respect to the lower/upper block faces and to the extremities of the bar is characterized by parameter $d_z$, which is kept alike for simplicity. The positioning of these disks can not be explicitly defined and hence it is characterized as a model parameter.

The vertical play is considered equally distributed along the structure's height as six local vertical plays ($v_{p,l} = v_p/6$) whose development is described as follows. Each bar object is equipped with three contacting components to represent the influence of the sling, presented in Figure 9.

These components are comprised of a cuboid box of length equal to 2.5 times the bar diameter (placed above the bar), a sphere of diameter equal to the bar diameter (placed below the bar) and another cuboid box of length equal to 1.5 times the bar diameter (placed within the bar) named as 'bar-bar penetration stopper'. The vertical position of these components is assigned such that the sphere from the adjacent
<figure><img src="image_10.png" /><figcaption>Figure 8: Concept of the connector modelling incorporating horizontal play</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Concept of the connector modelling incorporating vertical play</figcaption></figure>

bar detects contact with the boxes of the current bar. The distance between the box placed above the bar and the bar surface, as well as the distance between the sphere and the bar surface, is half the local vertical play. This way, the box-sphere combination from two adjacent bars collectively accounts for a local vertical play. The adjacent bars are assigned invisible to each other and their interpenetration is avoided via contact detection between this bar-bar penetration stopper and the sphere from the adjacent bar.

In the NSCD model, all components of each model object has no relative movement between them meaning that the circular disks in the block model object move according to the corresponding block. Similarly, all three components of the bar object move according to the bar movement. This feature enables to indirectly implement the block-bar interaction and sling functionality representative of the real structure. An illustrative representation of the connector model (inclusive of the plays) in action
247

248

249

250

251

252

is added in Figure 7 for understanding. Here, the interaction of the vertical play constituents (box and sphere) can be clearly seen. The distribution of vertical play (as six local vertical plays) enables mimicking the looseness in the sling as an indirect chain made from multiple linkages. Once the complete vertical play is consumed, the setup of bars mimics the tension in the sling. Such a state is reached depending (mostly) on the assigned looseness in the vertical play and the energy of the impacting projectile. This makes the vertical play ($v_p$) a model parameter.

253

254

255

256

257

258

259

This intricate setup allows us to model the large displacements in the model structure courtesy to the accompanying horizontal play ($h_p$). Its presence controls the maximum relative lateral displacement between the box and sphere components of the vertical play. Thereby, a possible contact between the vertical play components from adjacent bars is always ensured, which makes the connector model's functionality robust. An additional box and a sphere are added at the bottom and topmost layer of block objects respectively. This is to respectively mimic the anchorage of the sling with the structure's bottom and to make sure the bars do not fly out of the structure's top surface upon impact.

260

### 3.2.3. *Projectile modelling*

261

262

263

264

265

The experimental equivalent projectile is modelled as a convex 2D mesh of manually assigned vertices. The projectile shape and size ($1/3^{rd}$ of the wall height) features are in accordance with the requirements for flexible barriers testing [20]). The mass of the projectile is assigned 2600 kg, and the inertia is assigned through the convex hull (as used for the block model). The input kinetic energy corresponding to the two impact tests (i.e., 520kJ and 1020kJ) is derived to assign the input velocity to the projectile model.

266

### 3.2.4. *Interaction laws*

267

268

269

270

The model objects (i.e., blocks, connectors, projectile, and ground) are defined as a combination of different contactors, collectively making the model equivalent to the real structure. These contactors are defined with a uniquely assigned collision group to identify them and differentiate one from another. The collision groups of the contactors in the model bodies are listed in Table 2.

271

272

273

274

275

The interaction between model bodies is assigned through Newton impact friction non-smooth law [3].

This law governs the interaction between a pair of contactors (identified through their respective collision groups) via a user-defined coefficient of friction ($\mu$) and coefficient of restitution (e). The assigned interaction laws in our NSCD model are listed in Table 3. Here, the various interactions in the numerical model bodies are chronologically listed, analogous to the construction of the real structure.
<table><caption>Table 2: Model bodies constitution for interaction</caption><thead><tr><th>ine</th><th>Collision group</th><th>Contactor</th><th>Collision group</th><th>Contactor</th></tr></thead><tbody><tr><td>ine 0</td><td></td><td>Block (cuboid portion)</td><td>5</td><td>Hollow disk</td></tr><tr><td>1</td><td></td><td>Block (curvature portion)</td><td>6</td><td>Box (vertical play)</td></tr><tr><td>2</td><td></td><td>Bar (odd)</td><td>7</td><td>Sphere (vertical play)</td></tr><tr><td>3</td><td></td><td>Bar (even)</td><td>8</td><td>Projectile</td></tr><tr><td>4</td><td></td><td>Bar-Bar penetration stopper</td><td>9</td><td>Ground</td></tr><tr><td>ine</td><td></td><td></td><td></td><td></td></tr></tbody></table>

The frictional interaction between the projectile and the blocks and that between blocks are governed by the concrete-concrete friction coefficient ($\mu_{cc}$). Similarly, the interaction of each block of the wall with the ground is governed by the concrete-soil friction coefficient ($\mu_{es}$). The restitution coefficient ($e$) is kept the same for both of these interaction types.

Besides, all interactions between steel material and blocks and ground are assigned a constant friction coefficient of a relatively low magnitude of 0.2 for a minimum energy dissipation through these frictional interactions in the overall model system. Also, the restitution coefficient with zero magnitude is assigned presuming all such contacts being perfectly inelastic (i.e., no rebound after contact). Moreover, as the sling is a continuous element and bars do not interact with each other, the friction coefficient between the boxes and sphere of connector elements (vertical play and penetration stopper) is assigned zero. This is to indirectly represent consistent stretching in the sling as local vertical plays start to be fully consumed and subsequently tension in the sling in the event of complete consumption of the available vertical play. Also, $e = 0$ in the bar-to-disk contact assures no relative rebound once the horizontal play is consumed and the bar comes in contact with the disk (analogues to the tube and concrete block hole interaction in the real structure.)

### 3.2.5. Simulation database

The structure model is impacted at the wall center height and in the middle of its footprint length. The model wall response is computed for a one-second duration with a 0.25 ms time step for which the simulation duration of about 20 minutes is reported. A demonstrative post-impact deformed structural
<table><caption>Table 3: Assigning the interaction between model bodies contactors</caption><thead><tr><th>line Construction</th><th>Interacting pair (Colliding groups)</th><th>Friction</th><th>Restitution</th></tr></thead><tbody><tr><td rowspan="2">line<br>Foundation</td><td>Block (1) - Ground (9)</td><td>μ<sub>cs</sub></td><td>e</td></tr><tr><td>Odd bar (2) - Ground (9)</td><td>0.2</td><td>0</td></tr><tr><td rowspan="2">Stacking blocks</td><td>Block (1) - Block (1)</td><td>μ<sub>cc</sub></td><td>e</td></tr><tr><td>Odd bar (2) - Disk (5)</td><td>0.2</td><td>0</td></tr><tr><td rowspan="2">Inserting tubes</td><td>Even bar (3) - Disk (5)</td><td>0.2</td><td>0</td></tr><tr><td>Box (6) - Sphere (7)</td><td>0</td><td>0</td></tr><tr><td rowspan="2">Inserting sling</td><td>Penetration stopper (4) - Sphere (7)</td><td>0</td><td>0</td></tr><tr><td>Block (0,1) - Projectile (8)</td><td>μ<sub>cc</sub></td><td>e</td></tr><tr><td>Impacting the wall</td><td></td><td></td><td></td></tr></tbody></table>

conformation is presented in Figure 10. Here, the initial (pre-impact) configuration is also presented (similar to Figure 4) where transparent blocks reveal the underlying connector model characteristics.

<figure><img src="image_12.png" /><figcaption>Figure 10: Qualitative representation of initial and post-impact model conformations and highlighted connector model characteristics</figcaption></figure>

The output files of the model impact tests provided detailed information on the displacement, rotation, and velocities of each block and thereby enabled direct comparison with the experimentally obtained results at the corresponding locations. Besides, a detailed output database recorded for the contact forces also enabled the additional post-processing of the model structure behavior extending the limits
of the experimentally acquired database to a computationally advanced level.

### 3.3. Evaluation of the NSCD model response

The presented NSCD model comprises of many input parameters, most of which are directly taken from the real structure description (Table 1). By contrast, five model parameters couldn't be given a precise value from this description. The first two take into account the imposed constraints for the block-connector ($d_z$) and connector-connector ($v_p$) interaction while the other three concern the interaction laws between contacting bodies ($\mu_{cc}$, $\mu_{cs}$ and $e$).

The relative disk position ($d_z$) could not be precisely derived from the block and tube geometry. The vertical play ($v_p$) magnitude is variable from place to place in the structure and it is difficult to measure precisely. The coefficients governing friction between blocks ($\mu_{cc}$) and that between the blocks and the soil ($\mu_{cs}$) are difficult to measure precisely as dealing with the dynamic loading of an articulated structure where relative movement between bodies include translation and rotation. Last, the restitution coefficient ($e$) is a model parameter whose precise value couldn't be derived from some mechanical or geometrical features of the real structure.

<table><caption>Table 4: Model parameters with their considered range and mean value</caption><thead><tr><th>ine Parameter</th><th>Possible range</th><th>Mean value</th><th>Unit</th></tr></thead><tbody><tr><td>ine Bar-block contactor disk position ($d_z$)</td><td>5 - 10</td><td>7.5</td><td>cm</td></tr><tr><td>Vertical play ($v_p$)</td><td>1 - 10</td><td>5.5</td><td>cm</td></tr><tr><td>Friction coefficient concrete-concrete ($\mu_{cc}$)</td><td>0.25 - 0.55</td><td>0.4</td><td>(-)</td></tr><tr><td>Friction coefficient concrete-soil ($\mu_{cs}$)</td><td>0.3 - 0.6</td><td>0.45</td><td>(-)</td></tr><tr><td>Restitution coefficient ($e$)</td><td>0 - 0.3</td><td>0.15</td><td>(-)</td></tr><tr><td>ine</td><td></td><td></td><td></td></tr></tbody></table>

These parameters are listed in Table 4 with their range of possible values. The range for the vertical play is derived from observations of the real structure. The range for the disk position is determined as a fraction of the block height. The concrete-concrete friction angle is considered with the applicable range information retrieved from the literature [36]. The range for the friction angle between concrete and soil is established assuming slightly higher values than between concrete and concrete. The restitution
coefficient magnitude is considered ranging from perfectly inelastic collision (i.e., $e = 0$) upto a magnitude

towards perfectly elastic collision. This upper limit is selected based on multiple simulations such that

the rebound effect does not cause excessive inter-blocks displacement.

<figure><img src="image_13.png" /><figcaption>Figure 11: Input sample of size 300 (based on Sobol' sampling method) comprising of five parameters following uniform distribution in their respective ranges of variability</figcaption></figure>

As a first step, the overall response of the NSCD model is assessed through a sample set of different combinations of the five input parameters. This sample set is statistically generated using sampling methods based on Monte-Carlo (MC) and quasi MC (such as Sobol', Halton) [37] implemented in the UQlab input module [38]. The Sobol' sampling method is used to define a set of 300 combinations of these five input parameters, as presented in Figure 11. This sampling method provides a random selection of the input parameters within their predefined ranges such that the input space is well-filled. Besides, without any prior bias to the magnitude of the five model parameters, uniform distributions are assigned for all.

The global minimum and maximum range of outcomes from these 300 different simulations are presented in Figures 12a and 12b together with the NSCD model outcome from the mean value set of the input parameters. Globally, the range of numerical outcomes envelops the experimental observations at all measurement points. However, the range of model outcomes for the point Top remains below the experimental observations for the 1020 kJ impact. The loss of mechanical continuity (as reported in section 2) introduces a difference between the real structure and its numerical counterpart where no
<figure><img src="image_14.png" /><figcaption>Figure 12: Comparison of the experimental observations against the range of numerical computation response for displacement for (a) 520kJ and (b) 1020 kJ impact energy tests. The response obtained from the deterministic set of mean values (*) of model parameters is also highlighted.</figcaption></figure>

concrete block fracture is possible.

In spite of the difference at the point *Top* for 1020kJ test, the overall model-to-experiment correspondence is good considering all the data describing the structure response. This suggests that the NSCD model has the capacity to model the structure global response, assuming no block fracture. Also, the envelope infers that there could exist at least one set of five parameters that allow the NSCD model to possibly replicate the experimental response. This demands an extensive model calibration procedure presented hereinafter.

4. Calibration of the NSCD model

The aim of the calibration is to find the best set of the five model parameters (Table 4) such that the simulated spatio-temporal structure impact response mimics that of the real structure considering
four observation points ( i.e., points *Top*, *Base*, *C* and *D* as shown in Figure 2). The time evolution of the displacement along the x-axis is accounted for considering three representative time instants: the initial stage, the maximum, and the rest, referred to as *init*, *max* and *rest* respectively. The displacement measured at these three-time instants is considered to present a global narrative of the structure's displacement, and thus performance from the point of view of its practical application. Notably, the *init* time instant, where the displacement evolves rather linearly, is aimed at reflecting the dynamics of the structure, is arbitrarily fixed to 0.1 s after the impact beginning in the impact axis (i.e., for points *Top* and *Base*) and 0.25 s in the distant plane (i.e., for points *C* and *D*).

### 4.1. Calibration strategy

The five model parameters are calibrated considering the displacement at the four measuring points, at the three different times and for the two impact experiments. The calibration is performed using the Bayesian interface statistical learning procedure on the meta-model of the articulated model structure, following the conceptual pathway presented in Figure 13.

<figure><img src="image_15.png" /><figcaption>Figure 13: Conceptual workflow of the NSCD model calibration using a combination of stochastic methods</figcaption></figure>

This complex pathway is motivated by the fact that reaching this ambitious calibration objective demands a very large number of model computations. This is practically impossible with each model computation taking about 20 minutes. This limitation is averted using surrogate of NSCD model, representative of the structure response for different points at various time instances, through a collection of meta-models. Here, we use polynomial chaos expansion (PCE) based meta-modelling technique for surrogate model development, presented hereinafter.
## 4.2. Meta-model creation and accuracy assessment

A meta-model is a mathematical modelling tool that bridges the gap between the current version of the numerical model and the requisite large number of computations for stochastic analysis. Here, the small statistically created input sample of the model parameters (presented in Figure 11) is used to define the mathematical expression incorporating the model behaviour with the desired accuracy.

A polynomial chaos expansion (PCE) based meta-modelling technique implemented in the framework of UQ[py]Lab (a python library of UQLab [39]) is used in this work (see Appendix B for details). A total of 24 meta-models are created from the displacement response at three time instances - for four different points - for two different energies (i.e., $3 \times 4 \times 2 = 24$). For each meta-model creation, a mathematical relation (as presented in Equation B.2), is established between the 300 distinct sets of input parameters and the corresponding displacement output set. The mathematical relation here can be analogically referred to as obtaining a regression for a 2D database.

The accuracy of a meta-model is accessed through the leave-one-out (LOO) error (Equation B.5) in the order of $10^{-2}$ to $10^{-1}$ for all the computational cases. Thereby, as the NSCD model is highly nonlinear and non-smooth, we accept the meta-model accuracy and the underlying uncertainty. Such a level of accuracy with the added advantage of a single model computation within micro-seconds is deemed sufficient for the upcoming Bayesian interface based calibration method.

Further, the accuracy of the meta-models is assessed by comparing the meta-models predictions with the NSCD model computations for the 520 kJ and 1020 kJ impact tests (Figures 14a and 14b). Notably, the set of 300 NSCD model simulations that was used to create the meta-model is reused as a validation. The choice of not using an independent set is made based on the order of magnitude of LOO and the relative difference between NSCD and meta-model outputs for the same set of input parameters.

The cloud of points is rather well aligned with the diagonals in these plots, indicating that the predicted response by the meta-models fits with the corresponding NSCD simulations results. The reliability of the predictions along the impact axis is relatively better than in the distant one for both impacts. The distant axis points move relatively less and present more divergence from the diagonal line which highlights the significance of the variability in model dynamics farther from the impact axis as the model interaction parameters vary.
<figure><img src="image_16.png" /><figcaption>Figure 14: Comparison of displacements evaluated at four points (Top, Base, C and D) each at three time instances (initial, maximum and rest) from 300 NSCD model runs, and the corresponding evaluation by the PCE based meta-models with LOO error of order $10^{-2}$ to $10^{-1}$ for (a) 520kJ and (b) 1020 kJ impact tests</figcaption></figure>

### 4.3. Relative influence of the model parameters

Before proceeding with the calibration process, the influence of each parameter on the spatio-temporal displacement response of the structure is investigated through the Sobol sensitivity method also known as analysis of variance [40]. This method decomposes the variance of the output parameters as the sum of the contributions of the different input parameters including the possible interaction between input parameters. Each contribution is characterized by the ratios of the partial variance to the total variance, called Sobol sensitivity indices.

The accurate Sobol indices computation demands a large number of model computations (of order $10^6$). This is highly impractical in the present study if the NSCD model is to be used directly. Nonetheless, the meta-model of the NSCD model can be directly used to compute Sobol indices at zero cost [41, 42]. The UQlab sensitivity analysis module [43] is used for such computations. The total Sobol indices for the three time instances, at all considered points and for both impact tests are presented in Figure 15.

The Sobol indices are observed to significantly differ from one parameter to the other. The higher values are observed for the vertical play ($v_p$) and concrete-to-soil friction coefficient ($\mu_{cs}$) e.g., at points
<figure><img src="image_17.png" /><figcaption>Figure 15: Total Sobol indices of the spatio-temporal displacement response for the sensitivity analysis of the NSCD model parameters for (a) 520kJ and (b) 1020kJ impact test</figcaption></figure>

D and C respectively. Even though other parameters globally exhibit lower Sobol indices, their overall influence cannot be neglected. Besides, we observe that the Sobol index for a given parameter varies over time and is different from one impact energy to the other, e.g., '$v_p$' ranges from 0.01 to 0.7, and 'e' ranges from 0.04 to 0.22.

Overall, these results aid in understanding the structure response highlighting the mechanisms at work with time, complementing the insights to the experimental data discussed in Section 2. The influence of the vertical play ($v_p$) on the displacement of the upper blocks during both impact tests decreased with time. This highlighted the significant upward displacement (i.e., positive z-axis) of the blocks initiated in the early stage after the impact happened. By contrast, the influence of the concrete-to-soil friction coefficient ($\mu_{cs}$) showed a global increasing trend with time, which is attributed to the progressive base sliding of the wall. Furthermore, the Sobol index of the concrete-to-concrete friction coefficient, $\mu_{cc}$ suggested that the relative displacement between blocks is much higher in the upper part of the structure (points Top and D).

Besides, the sum of indices for a given situation exceeds one, expressing the interaction between input parameters. In other words, changing the value of a given parameter modifies the structure dynamic response which consequently alters the influence of other parameters. Overall, non-zero magnitude of all parameters affirms to consider their variability for calibrating the structure response. It also validates the identifiability of all parameters [44] for the Bayesian interface based studies presented hereinafter.
### 4.4. Bayesian interface for model calibration

The Bayesian interface is inspired by the Bayes' theorem [45] - a representation of the changing beliefs - simply demonstrating that probability of a 'hypothesis' being correct becomes more reliable with supporting 'evidence'. In our work, the 'hypothesis' states that the 300 distinct sets of input parameters give a correct range of displacement output. The idea of reliability increase in on our hypothesis, given the experimental database as 'evidence', means that there exist at least one set of input parameters such that the recorded evidence is reproduced, in other words, the model is calibrated. Hence, our goal is to find the optimal values of the input parameters that allow one to best fit the model predictions to the observations. A brief description of the Bayesian interface is presented in Appendix C.

#### 4.4.1. Representation

The available displacement experimental data at the points of interest and at three time instances are taken as a benchmark to calibrate the NSCD model parameters to obtain the requisite set of parameters. Here, we use the PCE-based meta-models of NSCD model as the forward model or 'prior' (see Equation C.2) to accelerate Bayesian computations, an approach extensively used by many researchers [46, 47, 48]. The meta-model based calibrated set of input parameters i.e., 'posterior' (see Equation C.3) is the representative of the NSCD model calibration.

In the UQlab Bayesian interface framework, the uncertainty in the model prediction is assigned via added Gaussian discrepancy (see Equation C.5) when correlated with the recorded experimental data. In the present work, as the model calibration relies on a single corresponding experimental measurement, deemed precise, the discrepancy with a known residual variance of order $10^{-20}$ is manually assigned.

#### 4.4.2. Implementation

The Bayesian interface module in UQLab [49] is implemented to obtain a set of input parameters considering three approaches referred to as 'point', 'energy', and 'all'. These three approaches are inspired by the scope of the user's interest to reproduce the experimental outcome which ranges from a particular time instance calibration to overall model behavior calibration in both space and time.

The first approach amounts to giving priority to a given point, meaning that the model is calibrated to provide precise predictions over time for this specific point courtesy of the multiple model output feature (demonstrated in Equation C.8). The second approach, 'energy' is in line with strategies consisting
in calibrating a model at a given impact energy before using it for another impact energy. The ‘all’ approach aims at ensuring that the model provides reliable predictions at any point and whatever the impact energy. The second and third approaches are implemented courtesy of the multiple forward model feature (demonstrated in Equation C.10).

<table><caption>Table 5: Calibrated sets of model parameters following three different approaches</caption><thead><tr><th rowspan="3">Parameter</th><th colspan="12">Approaches</th></tr><tr><th colspan="4">Point</th><th colspan="4">1020kJ</th><th colspan="2">Energy</th><th rowspan="2">All</th><th rowspan="3">Units</th></tr><tr><th colspan="4">520kJ</th><th colspan="4">1020kJ</th><th>520kJ</th><th>1020kJ</th></tr><tr><th>Top</th><th>Base</th><th>C</th><th>D</th><th>Top</th><th>Base</th><th>C</th><th>D</th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>ine d<sub>z</sub></td><td>5.7</td><td>9.6</td><td>5.3</td><td>5.4</td><td>9.5</td><td>6.2</td><td>6.8</td><td>5.0</td><td>7.2</td><td>6.8</td><td>6.8</td><td>cm</td></tr><tr><td>v<sub>p</sub></td><td>5.2</td><td>3.8</td><td>6.7</td><td>6.5</td><td>9.0</td><td>8.0</td><td>7.6</td><td>7.0</td><td>6.5</td><td>7.1</td><td>7.1</td><td>cm</td></tr><tr><td>μ<sub>cc</sub></td><td>0.251</td><td>0.419</td><td>0.474</td><td>0.265</td><td>0.378</td><td>0.512</td><td>0.510</td><td>0.342</td><td>0.351</td><td>0.323</td><td>0.316</td><td>(-)</td></tr><tr><td>μ<sub>cs</sub></td><td>0.311</td><td>0.304</td><td>0.340</td><td>0.346</td><td>0.306</td><td>0.351</td><td>0.301</td><td>0.301</td><td>0.320</td><td>0.308</td><td>0.307</td><td>(-)</td></tr><tr><td>e</td><td>0.190</td><td>0.194</td><td>0.138</td><td>0.033</td><td>0.152</td><td>0.290</td><td>0.283</td><td>0.200</td><td>0.198</td><td>0.220</td><td>0.222</td><td>(-)</td></tr></tbody><tfoot><tr><td colspan="13">ine</td></tr></tfoot></table>

The calibrated sets of input parameters are reported in Table 5 obtained as a ‘maximum a posteriori’ and ‘mean’ point estimates (see Equation C.15) for ‘point’ approach and for both other approaches respectively. Depending on the calibration approach, significant differences in parameter values are observed. For example, $v_p$ varied from 3.8 to 9 cm depending on the calibration strategy. A ratio of about two between the minimum and the maximum values attributed to a parameter is also observed for $d_z$ and $e$ in particular. Interestingly, the set of values calibrated following the ‘all’ approach is very close to that obtained following the ‘energy’ approach while considering the 1020kJ impact. Notably, the results obtained from the basic single-time instance approach are not presented in this paper due to limited interest.

### 4.4.3. Interpretation

The displacement response comparison between experimental results and NSCD model results obtained considering the magnitude of these input parameters is presented in Figure 16. The interpretation of these observations is presented hereinafter.
<figure><img src="image_18.png" /><figcaption>Figure 16: Bayesian inversion based calibration of (a) 520kJ and (b) 1020kJ impact test. Curves in *black* show the experimental data. Curves in *red*, *green* and *blue* show the NSCD simulation results corresponding to approaches ‘point’, ‘energy’ and ‘all’, respectively.</figcaption></figure>

The first approach (*i.e.*, 'point' approach) results in a total of four calibrated sets of input parameters for each impact test. This approach allows us to represent the influence of model dynamics at different points as the best set of input parameters is reported distinct for each point and for both energies. The experimental and numerical response, respectively in *black* and *red* in Figure 16, are in a very good quantitative and qualitative agreement with each other justifying the usage of stochastic methods for model calibration. This very good agreement emphasizes the success of choosing only three representative time instances over the whole impact duration. The significant difference observed at the point *Top* for the 1020kJ test is attributed to the loss of mechanical continuity as discussed in section 2. The obtained eight distinct sets of parameters (Table 5) highlight the significance of local dependency of the model constitution to accurately replicate the experimental response.

The second approach (*i.e.*, 'energy') collectively takes into account the model outputs at three time instances for all four measurement points for a given impact test. This way, two calibrated sets of input parameters are obtained from the Bayesian interface one each for the 520kJ and 1020kJ impact test. This approach presents the influence of the impact energy on the global model response. The corresponding response obtained from the NSCD model run (*in green*) is compared against the experimental observation.
Similar to the ‘point’ approach, a good qualitative agreement is observed between the experimental and numerical responses. Interestingly, the recorded ‘energy’ response using a single set of input parameters presented a similar response to the four distinct sets of these parameters in the ‘point’ approach for each energy level, highlighting the successful implementation of the multiple forward model feature of UQLab.

The third approach (i.e., ‘all’) collectively takes into account the model outputs at three time instances for all four measurement points for both impact tests. Hence, a single set of calibrated input parameters is obtained representing the global response of the NSCD model, with time and space. Here again, the comparison between the corresponding response obtained from the NSCD model (in blue) is in good qualitative and quantitative agreement with the experimental observation.

### 4.4.4. Application

Finally, the counterpart of the input sample distribution before the calibration (see Figure 11), the posterior input sample distribution obtained from the ‘all’ approach is presented in Figure 17. The data point corresponding to the calibrated set is also highlighted. Notably, from the histogram of the outputs, all input parameters present a locally emphasised magnitudes corresponding to the good fits with the desired output. However, the final set of parameters deviates slightly from the most probable magnitude of a particular input parameter. This difference is negligible for all parameters except $\mu_{cc}$ where the final selected parameter magnitude is slightly larger than the most probable one. Besides, a very narrow range of $\mu_{cs}$ parameter indirectly highlights its relatively high sensitivity to the model response (as also highlighted in Figure 15).

Further, as the counterpart of the displacement response before the calibration process (see Figure 12), the post-calibration response (from the ‘all’ approach) is presented in Figure 18. Overall, a good qualitative agreement between experimental and numerical observations is reported favouring the calibrated set of input parameters. From the quantitative comparison of the collective experimental and numerical (posterior mean) displacement measurements at three time representatives, a root mean square error (RMSE) of 9.9cm and 2.3cm is reported for impact and distant axis respectively which account for about 8.6% and 5.6% of representative displacements from the respective axis. A relatively better accuracy to reproduce the distant axis response favours the model’s practical utility. Continuing with the loss of mechanical continuity argument for the point Top during the 1020-kJ test (refer to section 2), the reported relatively large quantitative and qualitative difference is accepted and also is not considered in
<figure><img src="image_19.png" /><figcaption>Figure 17: Output sample obtained after Bayesian calibration corresponding to ‘all’ approach. Here, diagonal plots present histograms of the five model parameters and other plots present their distribution, where colours reflect the density of points. The red ⋆ and line present the calibrated value of the respective model parameter.</figcaption></figure>

RMSE estimation.

Moreover, the range of possible model outputs is also presented using all sets of input parameters from the posterior distribution. The relatively narrow width of the band in comparison to the prior sample validates the computational efficiency of the Bayesian interface. Besides, in co-relation to the input sample distribution (see Figure 18), most of the output trends in the band are likely to be closer to the posterior-mean response. Nonetheless, such co-relation is not made as the eventual goal is to obtain a deterministic set of input parameters for the further utility of the NSCD model.

Globally, the NSCD model is well-calibrated as a positive consequence of the Bayesian interface approach where a large number of model computations are favoured by the usage of the PCE-based meta-models. As an extension, for better and more reliable model calibration, the diversity in the experimental data both in time and space is favourable.

# 5. Discussion

## 5.1. *Benefits of the calibration approach*

This research proposed a complex approach making use of Bayesian interface accelerated with PCE-based meta-models for calibrating the developed NSCD model.
<figure><img src="image_20.png" /><figcaption>Figure 18: Comparison of the experimental displacement response against the NSCD model response obtained from the best calibrated set of input parameters along with the band of possible posterior sample based outcomes for (a) 520kJ and (b) 1020kJ impact tests respectively</figcaption></figure>

The model calibration concerned five input parameters whose values could not be precisely measured or determined a priori. The proposed stochastic analysis required an investment of 300 NSCD model runs to retrieve the essential database. By comparison, a simple calibration method, for example based on Monte Carlo simulations, would require a much larger number of simulations for obtaining a reliable set of model parameters. The combination of the meta-model and Bayesian interface thus presented a practical route to surpass the induced computational time constraint. On that subject, it is worth highlighting that the average computational time for the simulation of one structure impact response is of about 20 minutes with the developed NSCD model. By comparison, a typical duration of 10 hours was necessary with the finite difference model proposed by Furet et al. [10]. This highly significant difference further
539

highlights the practical implementation potential of the NSCD model given the output of interest.

540

The five model parameters were calibrated considering 24 outputs describing the structure displacement at four points and three-time instances and for two impact energies. The best set of model parameters was determined minimizing the deviation with the experimental data for this set of 24 outputs. A simple calibration method would place an emphasis on a limited number of outputs (thus considering a smaller number of points and time instances) or would rely on a qualitative evaluation of the simulated response with time, without any quantification of the deviation. For these reasons, a simple calibration method would fail in satisfactorily describing the complex spatio-temporal response of the structure.

541

542

543

544

545

546

547

The proposed calibration approach thus appears to result in a much more reliable set of model parameters. It could advantageously be used for calibrating any model where many parameters cannot be precisely determined a priori.

548

549

550

## 5.2. *On the optimum calibration strategy*

551

552

553

554

555

556

557

558

559

560

561

562

The results presented in section 4.4 revealed that the calibrated value of a particular model parameter significantly depends on the chosen calibration strategy. Ratios of about two are observed between extreme values of a given model parameter from one strategy to the other. This in particular indicates that a calibration based on the recorded structural response at a single point only may result in a set of parameters that may not be reliable for estimating the structure response at another location. This is demonstrated in Figure 19, where the focus is placed on the model calibrated against the displacement of the point *Base* during the 520kJ energy impact test. It is clearly seen that, despite a very good qualitative and quantitative agreement of the predictions by the model with the experimental response at this specific point, the model fails to correctly reproduce the responses at nearly all remaining locations for the 520kJ-energy test and at all four locations for the 1020 kJ-energy test. This is in line with the basic idea that the intended use of the model should be taken into account when defining the model parameters calibration strategy.

563

564

565

566

567

By contrast, the ‘all’ calibration approach appeared to result in displacement predictions at a specific point as good as the prediction by the model calibrated for this specific point according to the ‘point’ approach. This demonstrates that calibrating the model considering a large number of data is not at the expense of a reduction in prediction capacity as compared to models calibrated for a more specific purpose. In the end, the significance of using all spatio-temporal database for calibrating the model is
<figure><img src="image_21.png" /><figcaption>Figure 19: Comparison of experimental displacement response with the NSCD model response obtained from the point Base calibration set from 520kJ for (a) 520 kJ and (b) 1020kJ impact test</figcaption></figure>

fully justified.

### 5.3. Energy dissipation

Estimating the energy dissipative capacities of the protective structures exposed to a localized dynamic loading is of paramount importance for understanding their mechanical response and improving their design. Numerous research related to rockfall protection has considered the energy issue based on simulation results, with application to flexible barriers [50, 51, 52], embankments and walls [5, 10] and dissipative materials or components [53, 54]. The provided data mainly aimed at explaining how the incident block kinetic energy is transferred and dissipated in the impacted structure and its foundation.

As for articulated concrete block structures, finite difference numerical simulations presented in [10] suggested that 50% of the projectile kinetic energy is dissipated by concrete plasticization due to shear and tension. The NSCD model proposed in the present study accounts for energy dissipation via $e, \mu_{cc}$ and $\mu_{cs}$. The former indirectly accounted for damage to concrete blocks, while the two latter accounted for dissipation by friction between structure components.

Nevertheless, the significant difference in model parameters from one calibration strategy to the other, including the parameters governing energy dissipation ($e, \mu_{cc}$ and $\mu_{cs}$), suggests that energy dissipation estimated from the corresponding simulations may lead to different conclusions, in particular in terms of the respective contribution of each dissipative mechanism.
In the end, this highlights that simulation results interpreted in terms of energy dissipation should be considered with caution. The contribution by the various dissipative mechanisms (plasticization, friction...) strongly depends on the modelling strategy, on the constitutive laws and on the value attributed to the model parameters. This comes in addition to the classical energy conservation issue associated with some numerical schemes. These general but critical comments undoubtedly hold for any numerical model of structures exposed to dynamic loading, such as rockfall protective structures (e.g. [5, 10, 50, 52]). Further investigation is necessary for evaluating the influence of the model parameters on the contribution of each mechanism to energy dissipation.

# 6. Conclusions and perspectives

In this work, a numerical model of a structure consisting of concrete blocks connected one to the others via metallic components and designed to intercept rockfall is presented and successfully calibrated. The model is kept simple and developed in SICONOS software package based on NSCD approach mainly for the computation time-saving purpose. NSCD revealed efficiency for modelling such a complex structure where many components interact one with each other and some are invisible to each other. Each component in the real structure is directly and indirectly reproduced in the NSCD model considering the interaction mechanisms and the design framework of objects with contactors in SICONOS. An average computation time of 20 mins is reported for the NSCD model in comparison to highly expensive FEM computations lasting for about 10 hours for the same computation.

Five model parameters required calibration. Three parameters governed the mechanical interaction between the various model components and two concerned the structure description. The calibration of these five model parameters is conducted based on the displacement with time at four different locations in the structure measured during real-scale experiments considering two projectile kinetic energies. The calibration is conducted considering the Bayesian interface statistical learning method, accompanied by the meta-modelling techniques. This is done in view of accounting for the spatial and temporal displacement response of the structure upon a projectile impact. The meta-modelling techniques presented a surrogate of the NSCD model which represented a similar response albeit negligible computation time (microseconds) in comparison to the approximated 20 minutes time for one NSCD model computation. Subsequently, a large number of model computations are made possible by meta-models for the efficient
613

614

615

usage of the Bayesian interface for model calibration. Besides, the Sobol sensitivity analysis is made possible through the surrogate model which presented the relative influence of one parameter to the other both in space and time.

616

617

618

619

620

621

In the calibration process different calibration approaches, ranging from the local point to the global structure's response are presented. It is evidenced that the model parameters value significantly depended on the calibration strategy, meaning on the number and variety of data used in this purpose. The final set of parameters value is reported to nearly mimic the spatial-temporal response of the real structure for both impact tests. An overall quantitative deviation of numerical results with the experimental evolution is reported to be 8.6% for the impact axis and 5.6% for the distant axis.

622

623

624

As a perspective, the developed NSCD model will be used to investigate the response of articulated concrete blocks structures when exposed to impacts under different conditions, in view of quantifying their real efficiency, in a similar approach as in Mentani et al. [11], Toe et al. [12] and Lambert et al. [13].

625

## Acknowledgment

626

627

628

629

This work is conducted in the frame of the SMART PROTECT research project, led by Géolithe with INRAe, INRIA and Myotis as partners and is funded by the Auvergne-Rhône-Alpes region. The authors would like to thank Dr Emilie Rouzies for her valuable assistance and helpful comments on statistical learning methods.

630

## References

631

632

633

[1] M. Jean, J. J. Moreau, Dynamics in the presence of unilateral contacts and dry friction: a numerical approach, in: G. Del Pietro, F. Maceri (Eds.), Unilateral problems in structural analysis. II, CISM 304, Springer Verlag, 1987, pp. 151–196.

634

635

636

[2] M. Jean, The non smooth contact dynamics method, Computer Methods in Applied Mechanics and Engineering 177 (1999) 235–257. Special issue on computational modeling of contact and friction, J.A.C. Martins and A. Klarbring, editors.

637

638

[3] F. Dubois, V. Acary, M. Jean, The Contact Dynamics method: A nonsmooth story, Comptes Rendus Mécanique 346 (2018) 247–262.
639

640

641

[4] V. Acary, B. Brogliato, Numerical methods for nonsmooth dynamical systems. Applications in mechanics and electronics., Lecture Notes in Applied and Computational Mechanics 35. Berlin: Springer. xxi, 525 p. , 2008.

642

643

[5] C. Ronco, C. Oggeri, D. Peila, Design of reinforced ground embankments used for rockfall protection, Natural Hazards and Earth System Science (2009).

644

645

646

[6] A. C. Y. Yong, C. Lam, N. T. K. Lam, J. S. Perera, J. S. H. Kwan, Analytical solution for estimating sliding displacement of rigid barriers subjected to boulder impact, Journal of Engineering Mechanics (2019).

647

648

[7] R. Green, J. S. Finlan, Rapid Design of a Modular Rockfall Protection Wall in Response to the 2016 KaiKoura Earthquake, 2021, pp. 151–165. doi:10.1061/9780784483688.015.

649

650

[8] S. Lambert, F. Bourrier, P. Gotteland, F. Nicot, An experimental investigation of the response of slender protective structures to rockfall impacts, Canadian geotechnical journal (2019).

651

652

[9] C. Williams, J. Morkeh, K. Dorfschmidt, C. Poon, P. Matlashewski, J. Carvalho, Innovative rockfall solutions based on calibration and field testing, Mining, Metallurgy & Exploration (2021).

653

654

655

[10] A. Furet, S. Lambert, P. Villard, J. P. Jarrin, Experimental and numerical impact responses of an innovative rockfall protection structure made of articulated concrete blocks, Rock Mechanics and Rock Engineering (2022).

656

657

658

659

[11] A. Mentani, L. Govoni, G. Gottardi, S. Lambert, F. Bourrier, D. Toe, A new approach to evaluate the effectiveness of rockfall barriers, Procedia Engineering 158 (2016) 398–403. VI Italian Conference of Researchers in Geotechnical Engineering, CNRIG2016 - Geotechnical Engineering in Multidisciplinary Research: from Microscale to Regional Scale, 22-23 September 2016, Bologna (Italy).

660

661

662

[12] D. Toe, A. Mentani, L. Govoni, F. Bourrier, G. Gottardi, S. Lambert, Introducing meta-models for a more efficient hazard mitigation strategy with rockfall protection barriers, Rock Mechanics and Rock Engineering (2018).

663

664

[13] S. Lambert, D. Toe, A. Mentani, F. Bourrier, A Meta-Model-Based Procedure for Quantifying the On-Site Efficiency of Rockfall Barriers, Rock Mechanics and Rock Engineering 54 (2021) 487–500.
665

666

667

[14] V. Acary, F. Perignon, Siconos: A Software Platform for Modeling, Simulation, Analysis and Control of Nonsmooth Dynamical Systems, SIMULATION NEWS EUROPE, ArgeSIM/ASIM 17 (2007) 19–26.

668

669

670

[15] V. Acary, Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact, Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 96 (2016) 585–603.

671

672

[16] J. Escallón, C. Wendeler, E. Chatzi, P. Bartelt, Parameter identification of rockfall protection barrier components through an inverse formulation, Engineering Structures 77 (2014) 1–16.

673

674

[17] M. A. Hariri-Ardebili, B. Sudret, Polynomial chaos expansion for uncertainty quantification of dam engineering problems, Engineering Structures 203 (2020) 109631.

675

676

[18] X. Guo, D. Dias, C. Carvajal, L. Peyras, P. Breul, Reliability analysis of embankment dam sliding stability using the sparse polynomial chaos expansion, Engineering Structures 174 (2018) 295–307.

677

678

679

[19] R. Gupta, D. Rossat, X. Dérobert, J. Baroth, M. Briffaut, G. Villain, F. Dufour, Blind comparison of saturation ratio profiles on large rc structures by means of ndt and sfe—application to the vercors mock-up, Engineering Structures 258 (2022) 114057.

680

[20] EOTA, Ead 340059-00-0106: Falling rock protections kits, 2018.

681

682

[21] A. Rafiee, M. Vinches, C. Bohatier, Application of the nscd method to analyse the dynamic behaviour of stone arched structures, International Journal of Solids and Structures 45 (2008) 6269–6283.

683

684

685

[22] M. Jean, V. Acary, Y. Monerie, Non-smooth contact dynamics approach of cohesive materials, Philosophical Transactions of the Royal Society of London Series A Mathematical and Physical Sciences (1934-1990) 359 (2001) 2497–2518.

686

687

[23] M. Wronski, M. Jean, Some computational aspects of structural dynamics problems with frictional contact, Contact Mechanics (1995) 137–144.

688

689

690

[24] F. Bourrier, V. Acary, Predictive Capabilities of 2D and 3D Block Propagation Models Integrating Block Shape Assessed from Field Experiments, Rock Mechanics and Rock Engineering 55 (2022) 591–609.
691

692

693

[25] N. Akhadkar, V. Acary, B. Brogliato, Multibody systems with 3D revolute joints with clearances: an industrial case study with an experimental validation, Multibody System Dynamics 42 (2018) 249–282.

694

695

[26] V. Acary, F. Bourrier, B. Viano, Variational approach for nonsmooth elasto-plastic dynamics with contact and impacts, Computer Methods in Applied Mechanics and Engineering 414 (2023) 116156.

696

697

698

[27] N. Collins-Craft, F. Bourrier, V. Acary, On the formulation and implementation of extrinsic cohesive zone models with contact, Computer Methods in Applied Mechanics and Engineering 400 (2022) 115545.

699

700

[28] P. A. Cundall, O. D. Strack, A discrete numerical model for granular assemblies, geotechnique 29 (1979) 47–65.

701

[29] N. Bićanić, Discrete element methods, Encyclopedia of computational mechanics (2004).

702

703

704

[30] Q. Z. Chen, V. Acary, G. Virlez, O. Brüls, A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints, International Journal for Numerical Methods in Engineering 96 (2013) 487–511.

705

706

707

708

[31] O. Brüls, V. Acary, A. Cardona, On the Constraints Formulation in the Nonsmooth Generalized-α Method, in: S. I. Publishing (Ed.), Advanced Topics in Nonsmooth Dynamics. Transactions of the European Network for Nonsmooth Dynamics, 2018, pp. 335–374. URL: https://hal.inria.fr/ hal-01878550. doi:10.1007/978-3-319-75972-2\_9.

709

710

711

[32] O. Brüls, V. Acary, A. Cardona, Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme, Computer Methods in Applied Mechanics and Engineering 281 (2014) 131–161.

712

713

714

[33] T. Schindler, S. Rezaei, J. Kursawe, V. Acary, Half-explicit timestepping schemes on velocity level based on time-discontinuous Galerkin methods, Computer Methods in Applied Mechanics and Engineering 290 (2015) 250–276.

715

716

717

[34] T. Schindler, V. Acary, Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook, Mathematics and Computers in Simulation 95 (2014) 180–199.
[35] C. B. Barber, D. P. Dobkin, H. Huhdanpaa, The quickhull algorithm for convex hulls, ACM Trans. Math. Softw. 22 (1996) 469–483.

[36] M. Taklas, M. Leblouba, S. Barakat, A. Fageeri, F. Mohamad, Concrete-to-concrete shear friction behavior under cyclic loading: experimental investigation, Scientific Reports 12 (2022) 9451.

[37] H. Niederreiter, 1. Monte Carlo Methods and Quasi-Monte Carlo Methods, 1992, pp. 1–12. URL: https://epubs.siam.org/doi/abs/10.1137/1.9781611970081.ch1. doi:10.1137/1.9781611970081.ch1. arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9781611970081.ch1.

[38] C. Lataniotis, S. Marelli, B. Sudret, UQLab user manual – The Input module, Technical Report, Chair of Risk, Safety & Uncertainty Quantification, ETH Zurich, 2015. Report UQLab-V0.9-102.

[39] S. Marelli, B. Sudret, UQLab: A Framework for Uncertainty Quantification in Matlab, 2014, pp. 2554–2563. doi:10.1061/9780784413609.257.

[40] I. M. Sobol, A screening design for factorial experiments with interactions, Mathematical and Computer Modelling 1 (1993) 407–414.

[41] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering & System Safety 93 (2008) 964–979. Bayesian Networks in Dependability.

[42] E. Rouzies, C. Lauvernet, B. Sudret, A. Vidard, How is a global sensitivity analysis of a catchment-scale, distributed pesticide transfer model performed? Application to the PESHMELBA model, Geoscientific Model Development 16 (2023) 3137–3163.

[43] S. Marelli, C. Lamas, K. Konakli, C. Mylonas, P. Wiederkehr, B. Sudret, UQLab user manual – Sensitivity analysis, Technical Report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. Report UQLab-V2.0-106.

[44] S. Dobre, T. Bastogne, C. Profeta, M. Barberi-Heyob, A. Richard, Limits of variance-based sensitivity analysis for non-identifiability testing in high dimensional dynamic models, Automatica 48 (2012) 2740–2749.

[45] T. Bayes, N. Price, LII. An essay towards solving a problem in the doctrine of chances. by the late Rev. Mr. Bayes, F. R. S. communicated by Mr. Price, in a letter to John Canton, A. M. F. R. S, Philosophical Transactions of the Royal Society of London 53 (1763) 370–418.
745

746

[46] Y. M. Marzouk, D. Xiu, A stochastic collocation approach to Bayesian inference in inverse problems, Communications in Computational Physics (2009).

747

748

[47] Y. M. Marzouk, H. N. Najm, Dimensionality reduction and polynomial chaos acceleration of bayesian inference in inverse problems, Journal of Computational Physics 228 (2009) 1862–1902.

749

750

[48] L. Yan, L. Guo, Stochastic collocation algorithms using 1-1-minimization for bayesian solution of inverse problems, SIAM Journal on Scientific Computing 37 (2015) A1410–A1435.

751

752

753

[49] P.-R. Wagner, J. Nagel, S. Marelli, B. Sudret, UQLab user manual – Bayesian inversion for model calibration and validation, Technical Report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. Report UQLab-V2.0-113.

754

755

[50] H. Xu, C. Gentilini, Z. Yu, X. Qi, S. Zhao, An energy allocation based design approach for flexible rockfall protection barriers, Engineering Structures 173 (2018) 831–852.

756

757

[51] L. Castanon-Jano, E. Blanco-Fernandez, D. Castro-Fresno, Design of a new energy dissipating device and verification for use in rockfall protection barriers, Engineering Structures 199 (2019) 109633.

758

759

[52] D. Di Giacinto, L. Grassia, G. Capriello, E. Ruocco, A novel steel damping system for rockfall protection galleries, Journal of Constructional Steel Research 175 (2020) 106360.

760

761

762

[53] M. Previtali, M. O. Ciantia, S. Spadea, R. P. Castellanza, G. B. Crosta, Multiscale modelling of dynamic impact on highly deformable compound rockfall fence nets, Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 174 (2021) 498–511.

763

764

[54] S. Yan, Y. Wang, D. Wang, S. He, Application of EPS geofoam in rockfall galleries: Insights from large-scale experiments and fdem simulations, Geotextiles and Geomembranes 50 (2022) 677–693.

765

766

[55] G. D. Saxce, Z. Q. Feng, New inequality and functional for contact with friction: The implicit standard material approach, Mechanics of Structures and Machines 19 (1991) 301–325.

767

768

[56] V. Acary, M. Brémond, O. Huber, Advanced Topics in Nonsmooth Dynamics., Acary, V. and Brüls. O. and Leine, R. (eds). Springer Verlag, 2018, p. 375–457.

769

770

[57] S. Marelli, B. Sudret, Uqlab user manual - polynomial chaos expansions, 2015. doi:10.13140/RG.2.1.3778.7366.
## Appendix A. NSCD method

The constitutive framework of the NSCD method is described hereinafter.

### Appendix A.1. Newton-Euler equations

The configuration of each rigid body in three dimensions is described by the position of its center of mass $x_g \in \mathbb{R}^3$ and the orientation of the body-fixed frame with respect to a given inertial frame $R \in SO^+(3)$. The velocity of the body is given by the velocity of the center of mass $v_g = \dot{x}_g \in \mathbb{R}^3$ and the angular velocity of the body expressed in the body-fixed frame $\Omega \in \mathbb{R}^3$. A possible formulation of the Newton-Euler equations of motion for each body is

$$ 
\left\{ 
\begin{aligned}
m \dot{v}_g &= f(t, x_g, v_g, R, \Omega) \\
I \dot{\Omega} + \Omega \times I \Omega &= M(t, x_g, v_g, R, \Omega) \\
\dot{x}_g &= v_g \\
\dot{R} &= R \tilde{\Omega}
\end{aligned}
\right. 
\qquad (A.1) 
$$

where $m > 0$ is the mass, $I \in \mathbb{R}^{3 \times 3}$ is the matrix of moments of inertia around the center of mass and the axis of the body-fixed frame, and $f(\cdot) \in \mathbb{R}^3$, respectively $M(\cdot) \in \mathbb{R}^3$, are the total forces and respectively torques with respect to the center of gravity applied to the body. The matrix $\tilde{\Omega} \in \mathbb{R}^{3 \times 3}$ is given by $\tilde{\Omega}x = \Omega \times x$ for all $x \in \mathbb{R}^3$. In the implementation of SICONOS, the orientation matrix $R$ is parameterized by a unit quaternion $p$ such that $R = \Phi(p)$. In the equations of motion (Equation A.1), the occurrences of $R$ are substituted by their corresponding expressions in terms of $p$. The differential equation of Lie type $\dot{R} = R\tilde{\Omega}$ is replaced by $\dot{p} = \Psi(p)\Omega$ in the unit quaternion space. The closed formulae for $\Phi$ and $\Psi$ can be found in any textbook on rigid body dynamics. We denote by $q$ the vector of coordinates of the position and the orientation of the body, and by $v$ the body twist. In matrix notation, the relation between the body twist $v$ and the time derivative of $q$ is

$$ q := \begin{bmatrix} x_g \\ p \end{bmatrix}, \quad v := \begin{bmatrix} v_g \\ \Omega \end{bmatrix}, \quad \dot{q} := \begin{bmatrix} \dot{x}_g \\ \Psi(p)\Omega \end{bmatrix} = \begin{bmatrix} I & 0 \\ 0 & \Psi(p) \end{bmatrix} \quad v := T(q)v \qquad (A.2) $$

with $T(q) \in \mathbb{R}^{7 \times 6}$. A compact matrix form of the Newton-Euler equations is

$$ 
\left\{ 
\begin{aligned}
\dot{q} &= T(q)v, \\
M\dot{v} &= F(t, q, v),
\end{aligned}
\right. 
\qquad (A.3) 
$$
where $M \in \mathbb{R}^{6 \times 6}$ is the total inertia matrix and $F(t, q, v) \in \mathbb{R}^6$ collects all the forces and torques applied to the body given by

$$ M := \begin{pmatrix} mI_{3 \times 3} & 0 \\ 0 & I \end{pmatrix}, \quad F(t, q, v) := \begin{pmatrix} f(t, x_g, v_g, R, \Omega) \\ I\Omega \times \Omega + M(t, x_g, v_g, R, \Omega) \end{pmatrix}. \quad (A.4) $$

In the sequel, we assume that the vectors $q \in \mathbb{R}^{7n}$ and $v \in \mathbb{R}^{6n}$ collect the configuration parameters and the velocities of the $n$ bodies in the model. The equation of motion (Equation A.3) is rearranged accordingly.

## Appendix A.2. Contact modelling as unilateral constraints and Coulomb friction

For two contacting bodies A and B, we assume that we are able to define one or more contact pair composed of two candidate contact points $C_A$ and $C_B$ and a local frame at contact ($C_A, N, T_1, T_2$) (see Figure A.20 for details). In this local frame, the normal gap function $g_N$ is defined as the signed distance between the point $C_A$ and $C_B$, that is $g_N = (C_B - C_A) \cdot N$. The unilateral contact is given by Signorini's condition as

$$ 0 \le g_N \perp r_N \ge 0, \quad (A.5) $$

where $r_N$ is the local contact normal force. The symbol $g_N \perp r_N$ means that $g_N r_N = 0$. If the contact is open $g_N > 0$, the contact force $r_N$ has to be equal to zero. If the contact is closed $g_N = 0$, then the contact force is nonnegative $r_N \ge 0$. As we said before, the dynamics of finite dimensional system requires the definition of an impact law. Let's consider the relative velocity $u_N$ defined by $u_N = \dot{g}_N$. One of the simplest impact is the Newton impact law given by

$$ u_N^+ = -eu_N^- \text{, if } g_N = 0 \text{ and } u_N^- \le 0, \quad (A.6) $$

where $e \in [0, 1]$ is the coefficient of restitution, $u_N^+$ is the post-impact velocity and $u_N^-$ the pre-impact velocity.

For the tangential behaviour of the interface, we consider the Coulomb friction that can be formulated with the Signorini's condition at the velocity level and in terms of the complete reaction forces $r = [r_N, r_T]^\top$
<figure><img src="image_22.png" /><figcaption>Figure A.20: A contact pair and a local contact frame</figcaption></figure>

810 and the relative velocity $u = [u_N, u_T]^T$ as

$$ egin{cases} r = 0 & \text{if } g_N > 0 \quad (\text{no contact}) \\ r = 0, u_N \ge 0 & \text{if } g_N = 0 \quad (\text{take-off}) \\ r \in K, u = 0 & \text{if } g_N = 0 \quad (\text{sticking}) \\ r \in \partial K, u_N = 0, \|u_T\|_{r_T} = \|r_T\|_{u_T} & \text{if } g_N = 0 \quad (\text{sliding}) \end{cases} \quad (A.7) $$

811 where $K = \{r \in \mathbb{R}^3, \|r_T\| \le \mu r_N\}$ is the usual Coulomb friction cone. Using the modified local relative velocity introduced by De Saxcé and Feng [55, 56], $\hat{u} := u + \mu \|u_T\|_{\mathbf{N}}$ and the dual cone of K, i.e., 812 $K^* = \{z \in \mathbb{R}^3 | z^T x \ge 0 \text{ for all } x \in K\}$, the contact model is equivalent to 813

$$ K^* \ni \hat{u} \perp r \in K. \quad (A.8) $$

## Appendix A.3. Equations of motion with contact and friction

The gap function $g_N$ for a contact is generally a function of the configuration of the bodies $q$. Let us consider that we have a set of $m$ contact pairs given by the unilateral constraints:

$$ g_N^\alpha(q) \ge 0, \alpha \in \mathcal{I} \subset \mathcal{N}, |\mathcal{I}| = m. \quad (A.9) $$

817 For a contact $\alpha \in \mathcal{I}$, the relative normal velocity is related to $v$ through the relation

$$ u_N^\alpha = \dot{g}_N^\alpha(q) = J_{g_N}^\alpha(q) \dot{q} = J_{g_N}^\alpha(q) T(q) v := G_N^\alpha(q) v, \quad (A.10) $$

818 where $J_{g_N}^\alpha(q)$ is the Jacobian of $g_N^\alpha$ with respect to $q$. The same type of relation can be written for the 819 tangential relative velocity $u_T$ as $u_N = G_T^\alpha(q)v$ and we get for $u^\alpha$

$$ u^\alpha := G^\alpha(q)v. \quad (A.11) $$
By duality, the total force generated by the contact $\alpha$ is

$$G^{\alpha, \top}(q)r^{\alpha} := G_{N}^{\alpha, \top}(q)r_{N}^{\alpha} + G_{T}^{\alpha, \top}(q)r_{T}^{\alpha}. \qquad (A.12)$$

Altogether the equation of motion with contact and Coulomb friction is given by

$$\left\{
\begin{array}{ll}
\dot{q} = T(q)v, \\
\\
M\dot{v} = F(t, q, v) + G^{\top}(q)r, \\
\\
r^{\alpha} = 0, & \text{if } g_{N}^{\alpha}(q) > 0, \\
\\
K^{\alpha,*} \ni \hat{u}^{\alpha} \perp r^{\alpha} \in K^{\alpha}, & \text{if } g_{N}^{\alpha}(q) = 0, \\
\\
u_{N}^{\alpha,+} = -e^{\alpha}u_{N}^{\alpha,-}, & \text{if } g_{N}^{\alpha}(q) = 0 \text{ and } u_{N}^{\alpha,-} \le 0
\end{array}
\right.
\quad \alpha \in \mathcal{I},$$

## Appendix A.4. Time-discretization

Let us consider a time discretization $t_0 < \dots < t_{k-1} < t_k < \dots < T$ with a constant time step

$h = t_{k+1} - t_k$. The Moreau-Jean scheme for the system is

$$\left\{
\begin{array}{ll}
q_{k+1} = q_k + hT(q_{k+\theta})v_{k+\theta} & \\
\\
M(v_{k+1} - v_k) - hF_{k+\theta} = G^{\top}(q_{k+1})P_{k+1}, & \qquad (A.13) \\
\\
P_{k+1}^{\alpha} = 0, & \alpha \notin \mathcal{I}_k \\
\\
K^{\alpha,*} \ni \hat{u}_{k+1}^{\alpha} + e^{\alpha}u_{N,k}^{\alpha} N \perp P_{k+1}^{\alpha} \in K^{\alpha} & \alpha \in \mathcal{I}_k.
\end{array}
\right.$$

where the notation $x_{k+\theta} = (1-\theta)x_k + \theta x_{k+1}$ is used for $\theta \in [0, 1]$ and the set $\mathcal{I}_k$ is the set of contact activated at the velocity level

$$\mathcal{I}_k = \{\alpha \in I \mid g_{N,k}^{\alpha} + \gamma u_{N,k}^{\alpha} \le 0\} \text{ with } \gamma \in [0, \frac{1}{2}]. \qquad (A.14)$$

In the time-stepping method, the unknown $P_{k+1}$ is an approximation of the impulses of the contact reaction measure $d_i$ over the time interval, that is

$$\int_{(t_k, t_{k+1}]} d_i \approx P_{k+1}. \qquad (A.15)$$

The contact reaction measure is related to the contact force when the motion is smooth enough by $d_i = r(t)dt$. The system (Equation A.13) is a second-order cone complementarity problem that is solved by SICONOS using a Gauss-Seidel method with projection [56].
## Appendix B. PCE based meta-model

Consider the articulated structure model represented by a $M(\mathbf{X})$ as an equivalent mathematical model. Here, $\mathbf{X} \in \mathbb{R}^M$ is a random vector with independent components described by the joint probability density function (PDF) $f_{\mathbf{X}}$. Consider also a finite variance computational model as a map $Y = M(\mathbf{X})$, with $Y \in \mathbb{R}$ such that:

$$ \mathbb{E}[Y^2] = \int_{D_{\mathbf{X}}} M(\mathbf{x})^2 f_{\mathbf{X}}(\mathbf{x}) d\mathbf{x} < \infty \quad (B.1) $$

Then, under the assumption of Equation B.1, the PCE of $M(\mathbf{X})$ is defined as:

$$ Y = M(\mathbf{X}) = \sum_{\alpha \in \mathbb{N}^M} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) \quad (B.2) $$

where, the $\Psi_{\alpha}(\mathbf{X})$ are multivariate polynomials orthonormal with respect to $f_{\mathbf{X}}$, $\alpha \in \mathbb{N}^M$ is a multi-index that identifies the components of the multivariate polynomials $\Psi_{\alpha}$ and the $y_{\alpha} \in \mathbb{R}$ are the corresponding coefficients. In practical applications, the sum in Equation B.2 needs to be truncated to a finite sum by introducing the truncated polynomial chaos expansion:

$$ \mathcal{M}(\mathbf{X}) \approx \mathcal{M}^{PC}(\mathbf{X}) = \sum_{\alpha \in \mathcal{A}} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) \quad (B.3) $$

where, $\mathcal{A} \subset \mathbb{N}^M$ is the set of selected multi-indices of multivariate polynomials.

In this work, the least-angle regression (LARS) method is used to create the PCE meta-model truncated to the maximum polynomial degree ($p$) ranging from 1 to 20, and using hyperbolic truncation scheme ($q$) ranging from 0.5 to 1.

$$ \mathcal{A}^{M,p,q} = \{ \alpha \in \mathcal{A}^{M,p} : \|\alpha\|_q \le p \}, \quad \text{where, } \|\alpha\|_q = \left( \sum_{i=1}^{M} \alpha_i^q \right)^{1/q} \quad (B.4) $$

The accuracy of the constructed PCE is estimated by computing the leave-one-out (LOO) cross-validation error ($\epsilon_{LOO}$). It consists in building $N$ meta-models $M^{PC\setminus i}$, each one created on a reduced experimental design $\mathcal{X} \setminus \mathbf{x}^{(i)} = \{\mathbf{x}^{(j)}, j = 1, \dots, N, j \neq i\}$ and comparing its prediction on the excluded point $\mathbf{x}^{(i)}$ with the real value $y^{(i)} = M(\mathbf{x}^{(i)})$ [57]. The leave-one-cross-validation error can be written as:

$$ \epsilon_{LOO} = \frac{\sum_{i=1}^{N} \left( M(\mathbf{x}^{(i)}) - M^{PC\setminus i}(\mathbf{x}^{(i)}) \right)^2}{\sum_{i=1}^{N} \left( M(\mathbf{x}^{(i)}) - \hat{\mu}_Y \right)^2} \quad (B.5) $$
where, $\hat{\mu}_Y$ is the mean of the experimental design sample.

## Appendix C. Bayesian interface for model calibration

Consider the computational model $\mathcal{M}$ that allows the analyst to predict certain quantities of interest gathered in a vector $\mathbf{y} \in \mathbb{R}^{N_{out}}$ as a function of input parameters $\mathbf{x}$:

$$ \mathcal{M}: \mathbf{x} \in \mathcal{D}_{\mathbf{X}} \subset \mathbb{R}^M \mapsto \mathbf{y} = \mathcal{M}(\mathbf{x}) \in \mathbb{R}^{N_{out}} \quad (C.1) $$

The Bayesian interface for model calibration focuses on identifying the input parameters of a computational model to recover the observations in the collected output data set. It comprises of a computational forward model $\mathcal{M}$, a set of input parameters $\mathbf{x} \in \mathcal{D}_{\mathbf{X}}$ that need to be inferred, and a set of experimental data $\mathcal{Y}$. Here, $\mathcal{Y} \stackrel{\text{def}}{=} \{\mathbf{y}_1, \dots, \mathbf{y}_N\}$ is a global data set of $N$ independent measured quantities of interest ($\mathbf{y}_i$).

The forward model $\mathbf{x} \mapsto \mathcal{M}(\mathbf{x})$ is a mathematical representation of the system under consideration. The lack of knowledge on the input parameters is modelled by considering them as a random vector, denoted by $\mathbf{X}$ which is assumed to follow a so-called prior distribution (with support $\mathcal{D}_{\mathbf{X}}$), as presented in Figure 11 in the present work.

$$ \mathbf{X} \sim \pi(\mathbf{x}) \quad (C.2) $$

The Bayesian statistics combine this prior knowledge of the parameters with the few observed data points to obtain a statistical model called posterior distribution ($\pi(\mathbf{x} | \mathbf{y})$) of the input parameters, using Bayes' theorem [45], expressed as:

$$ \pi(\mathbf{x} | \mathbf{y}) = \frac{\pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x})}{\pi(\mathbf{y})} \quad (C.3) $$

Now, considering the available data set ($\mathcal{Y}$) as independent realizations of $\mathbf{Y} | \mathbf{x} \sim \pi(\mathbf{y} | \mathbf{x})$, the collected measurements result in the definition of the likelihood function $\mathcal{L}(\mathbf{x}; \mathcal{Y})$, which is a function of input parameters $\mathbf{x}$:

$$ \mathcal{L}: \mathbf{x} \mapsto \mathcal{L}(\mathbf{x}; \mathcal{Y}) \stackrel{\text{def}}{=} \prod_{i=1}^{N} \pi(\mathbf{y}_i | \mathbf{x}) \quad (C.4) $$
This implicitly assumes independence between individual measurements in $\mathcal{Y}$. Intuitively the likelihood function for a given $\boldsymbol{x}$ returns the relative likelihood of observing the data at hand, under the assumption that it follows the prescribed parametric distribution $\pi(\boldsymbol{y} | \boldsymbol{x})$.

As all models are simplifications of the real world, a discrepancy term ($\epsilon$) is introduced to connect real-world observations ($\mathcal{Y}$) to the predictions by the model. In practice, the discrepancy term represents the effects of the measurement error and model inaccuracy. The discrepancy term introduced here reads:

$$ \boldsymbol{y} = \mathcal{M}(\boldsymbol{x}) + \epsilon \qquad (C.5) $$

Here, the $\epsilon$ is assumed as an additive Gaussian discrepancy [49] with a zero mean and given covariance matrix ($\boldsymbol{\Sigma}$):

$$ \epsilon \sim \mathcal{N}(\epsilon | \mathbf{0}, \boldsymbol{\Sigma}) \qquad (C.6) $$

Taking insights from the discrepancy term definition, a particular measurement point ($\boldsymbol{y}_i \in \mathcal{Y}$, is a realization of the Gaussian distribution with mean value $\mathcal{M}(\boldsymbol{x})$ and covariance matrix $\boldsymbol{\Sigma}$. This distribution is named as discrepancy model and is expressed as:

$$ \pi(\boldsymbol{y} | \boldsymbol{x}) = \mathcal{N}(\boldsymbol{y} | \mathcal{M}(\boldsymbol{x}), \boldsymbol{\Sigma}) \qquad (C.7) $$

In application, the discrepancy model defines the connection between the supplied data ($\mathcal{Y}$) and the forward model. In the present work, as the model calibration relies on a single experimental measurement, the discrepancy model with known residual variance is assigned.

Afterwards, the $N$ independent available measurements gathered in the data-set (i.e., $\mathcal{Y} = \{\boldsymbol{y}_1, ..., \boldsymbol{y}_N\}$) are used to define the likelihood function as:

$$ \begin{aligned} \mathcal{L}(\boldsymbol{x}; \mathcal{Y}) &= \prod_{i=1}^{N} \mathcal{N}(\boldsymbol{y}_i | \mathcal{M}(\boldsymbol{x}), \boldsymbol{\Sigma}) \\ &= \prod_{i=1}^{N} \frac{1}{\sqrt{(2\pi)^{N_{\text{out}}} \det(\boldsymbol{\Sigma})}} \exp \left( -\frac{1}{2} (\boldsymbol{y}_i - \mathcal{M}(\boldsymbol{x}))^{\top} \boldsymbol{\Sigma}^{-1} (\boldsymbol{y}_i - \mathcal{M}(\boldsymbol{x})) \right) \end{aligned} \qquad (C.8) $$

Moreover, in the present work, the experimental data is retrieved through various points and for different impact energies. The Bayesian interface allows incorporating all these data points together for the model calibration by arranging the elements of $\mathcal{Y}$ in disjoint data groups and defining different likelihood functions for each data group [49].
Denoting the g-th data group by $\mathcal{G}^g = \{\mathbf{y}_i\}_{i \in \mathbf{u}}$, where $\mathbf{u} \subseteq \{1, ..., N\}$, the full data-set can be combined by

$$ \mathcal{Y} = \bigcup_{g=1}^{N_{gr}} \mathcal{G}^{(g)} \qquad (C.9) $$

Here, each of the $N_{gr}$ data groups contains measurements collected at the same measurement point. In the context of the present work, these measurements are experimental data recorded at 'init', 'max' and 'rest' time instances at a particular point. This makes it evident to have a different likelihood function $\mathcal{L}^{(g)}$ describing the experimental conditions that led to measuring $\mathcal{G}^{(g)}$. Assuming the Independence between $N_{gr}$ measurement conditions, the full likelihood function can then be written as

$$ \mathcal{L}(\mathbf{x}_{\mathcal{M}}, \mathbf{x}_{\epsilon}; \mathcal{Y}) = \prod_{g=1}^{N_{gr}} \mathcal{L}^{(g)}(\mathbf{x}_{\mathcal{M}}, \mathbf{x}_{\epsilon}^{(g)}; \mathcal{G}^{(g)}) \qquad (C.10) $$

where, $\mathbf{x}_{\epsilon}^{(g)}$ are the parameters of the g-th discrepancy group. The different model output groups are assigned through a model output map (MOMap) vector [49].

Thereby, following Bayes' theorem, the posterior distribution $\pi(\mathbf{x} | \mathcal{Y})$ of the parameters ($\mathbf{x}$) given the observations in $\mathcal{Y}$ can be written as:

$$ \pi(\mathbf{x} | \mathcal{Y}) = \frac{\mathcal{L}(\mathbf{x}; \mathcal{Y})\pi(\mathbf{x})}{Z} \qquad (C.11) $$

Here, $Z$ is a normalizing factor, known as the marginal likelihood or evidence, is added to ensure that this distribution integrates to 1:

$$ Z \stackrel{\text{def}}{=} \int_{\mathcal{D}_{\mathbf{x}}} \mathcal{L}(\mathbf{x}; \mathcal{Y})\pi(\mathbf{x})d\mathbf{x} \qquad (C.12) $$

The closed-form solutions do not exist in practice the posterior distribution is obtained through Markov chain Monte Carlo (MCMC) simulations. In the present work, among many proposed algorithms in [49], the Adaptive Metropolis (AM) algorithm is used and 100 parallel chains with 1000 steps are assigned to the MCMC solver. The start of the covariance adaptation in AM algorithm is assigned at the 100th step (see Wagner et al. [49] for more details).

The output ($\mathbf{y}$) predictive capabilities of the Bayesian interface is assessed through the comparison of prior ($\pi(\mathbf{y})$) and posterior ($\pi(\mathbf{y} | \mathcal{Y})$) output distributions as:
$$ \pi(\mathbf{y}) = \int_{D_{\mathbf{x}}} \pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x})d\mathbf{x} \qquad (C.13) $$

$$ \pi(\mathbf{y} | \mathcal{Y}) = \int_{D_{\mathbf{x}}} \pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x} | \mathcal{Y})d\mathbf{x} \qquad (C.14) $$

Lastly, in the present studies, the purpose of the Bayesian interface is to obtain the ‘best set of input parameters’. Given the posterior distribution ($\pi(\mathbf{x} | \mathcal{Y})$), we are interested in finding a suitable set among the posterior computed set i.e., $\mathbf{X} | \mathcal{Y}$. This is done through a point estimator ($\hat{\mathbf{x}}$) computed from:

$$ \pi(\mathbf{y} | \mathcal{Y}) \stackrel{\text{def}}{=} \pi(\mathbf{y} | \hat{\mathbf{x}}) \qquad (C.15) $$

This point estimator can be a mean or mode (maximum a posteriori ‘MAP’) [49] of the posterior distribution as per the user’s choice.