R Boulaud

C Douthe

K Sab

Modelling of curtain eect in rockfall barrier with the dynamic relaxation

Keywords: Discrete element method, numerical stability, sliding cable, dynamic relaxation, ring element, anti-submarine net

Flexible rockfall barriers are protection systems against risks of falling rocks. Their complex behaviour is still not well understood and the development of relevant models allowing for quick calculations is a need toward the optimisation of such structures. The key issue is the modelling of the so-called "curtain eect", where the cable net slides along the supporting cables. The present paper proposes hence a model of sliding cable submitted to concentrated forces. It addresses the stability conditions of the centred nite dierence scheme, especially in the framework of the Dynamic Relaxation Method where estimates are found for the best possible ctitious masses at nodes. Expressions are derived in the general case, but also for monotonically curved cables and looped cable elements or rings. Elementary problems are shown for validation, before a case study based on experiments conducted during the French National project C2ROP demonstrates the accuracy and the reliability of the proposed methodology. The present study is however limited to quasi-static loadings and therefore needs further developments to be extended to realistic dynamic cases which are discussed in the conclusion.

1. Introduction

General context

The hazard of landslide, mainly in mountain areas, compromises the safety of inhabitants.

The need to protect them, their properties and infrastructures against this risk requires the setting of protective structures. Rockfall barrier is an often used alternative because of its low weight and its high capacity to absorb energy. It can be installed, by specialised workers, on hard-to-reach areas. These structures are complex and involve non-standard connections in civil engineering. Schematically, their behaviour can be described as follows: a wire net (shown in Figure 1, but not necessarily) intercepts the rock trajectory. It then deforms, sliding along the supporting cables which are attached to the posts or the cli by dissipating devices that undergo irreversible transformations by yielding, rubbing or tearing. The dissipation of the rock kinetic energy is thus insured by the friction and plastic deformation of the net itself and by the brakes distributed along the cables. The sliding of the rings along each other and along the supporting cables is essential to the behaviour of the barrier: it allows the members to be rearranged in such a way that the structure can withstand the impact without stress concentration and adapt to any location of the impact.

The proper modelling of the sliding should give more insight into a very non-linear problem, and help to improve the results of companies trying to get an agreement for their products.

It was indeed shown in [START_REF] Chanut | Dynamic behavior of rock fall protection net fences: a parametric study[END_REF] that one out of two tests at the French experimental testing facility in Montagnole was negative.

Preprint submitted to International Journal of Solids and Structures

November 29, 2019 Indeed, 2008, a European Technical Agreement (ETA) was set for exible rockfall barriers in order to align the design standards between the EU countries. This agreement is in line with the recommendations of the technical guide ETAG 027, which is the reference publication for the design of barriers [START_REF] Eota | Etag 027[END_REF]. To receive it, the structure must withstand the fall of a block with a given energy with controlled deformation. The full-scale test procedure, proposed by the ETAG 027, is largely driven by the previous experiments on complete structures carried out by manufacturers and academics. We can notably mention the experiments of D. [START_REF] Peila | Evaluation of Behaviour of Rockfall Restraining Nets by Full Scale Tests[END_REF], where the block is guided by a trolley before the impact on the net [START_REF] Peila | Evaluation of Behaviour of Rockfall Restraining Nets by Full Scale Tests[END_REF] and the experiments of H. [START_REF] Grassl | Steel-net rockfall protection: Experimental and numerical simulation[END_REF] in partnership with the Swiss company Geobrugg [START_REF] Grassl | Steel-net rockfall protection: Experimental and numerical simulation[END_REF]. A product is hence approved for one critical energy level. Despite the development of powerful computational tools in the last decades [START_REF] Hambleton | Perforation of exible rockfall barriers by normal block impact[END_REF][START_REF] Volkwein | Numerische Simulation von exiblen Steinschlagschutzsystemen[END_REF][START_REF] Bertrand | Full-scale dynamic analysis of an innovative rockfall fence under impact using the discrete element method: From the local scale to the structure scale[END_REF][START_REF] Grassl | Steel-net rockfall protection: Experimental and numerical simulation[END_REF], this experiment is today the only recognised way to validate the performance of a exible barrier.

The Dynamic Relaxation Method

The behaviour of the exible barriers reveals strong geometric and material non-linearities, so that their modelling is complex and their calculations are time consuming, especially when the full dynamic simulation is conducted. The development of a modelling strategy that would allow quick calculations and parametric studies would thus be of great help in the comprehension of their complex behaviour and toward the development of new barriers.

The key issue is the modelling of sliding cables which is addressed here, with a focus on a stability condition with the Dynamic Relaxation Method (DRM). The dynamic relaxation method is a discrete numerical method used to determine the equilibrium state of a system submitted to mechanical loads which is very suitable for large displacement [START_REF] Day | An Introduction to Dynamic Relaxation[END_REF]. Its principle consists in linearising the dierential equation of the Newton's rst law by an explicit scheme of numerical integration. The static position of a structure is hence considered as the result of a smoothed dynamical process.

If the forces are conservative, the dynamical behaviour of the system has no inuence on its equilibrium state. The dynamic parameters (damping and mass) can thus be chosen to optimise the convergence of the algorithm. M. Papadrakakis proposed a general method to identify the critical parameters of a structure [START_REF] Papadrakakis | A method for automatic evaluation of the dynamic relaxation parameters[END_REF], but it is relatively complex and does not signicantly outperform other methods [START_REF] Veenendaal | An overview and comparison of structural form nding methods for general networks[END_REF]. That's why an alternative articial damping, the "kinetic damping" proposed by P. Cundall, is used in the algorithm presented in this paper [START_REF] Cundall | A computer model for simulating progressive large scale movements in block rock systems[END_REF]. It is based on the principle of mechanical energy conservation during an undamped dynamic process. Considering that in this case, each minimum of potential energy coincides with a maximum of kinetic energy, setting articially the kinetic energy to zero at every peak should lead progressively to the equilibrium state. This method is very popular among others [START_REF] Bel Hadj | Analysis of clustered tensegrity structures using a modied dynamic relaxation algorithm[END_REF][START_REF] Bel Hadj | A dynamic-relaxation formulation for analysis of cable structures with sliding-induced friction[END_REF][START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF], so that in the rest of this paper, numerical simulations will be carried out with a DR algorithm using kinetic damping.

This quasi-static modelling rely on the idea that the block will concentrate most of the inertial forces and of the kinetic energy, as the mass of the block is generally ve times greater than the mass of the net and that initially only the block is moving. It also relies on the fact that the brakes serve as fuses for the anchorages and surrounding cables or supporting members, cutting o peaks that may occur in usual impact problems. Therefore, a quasistatic calculation of the equilibrium state should give a quick and reasonable estimate of the equilibrium conguration after impact and hence provides to the engineer some useful global information of the barrier, such as estimates of the residual height, the brakes stroke or the depth of the net. The loading equivalence between static and dynamic loadings would then be based on the dissipated energy, not on forces. These assumptions build on preliminary results of the C2ROP research project [START_REF] Olmedo | Extended experimental studies on rockfall exible fences[END_REF] where similar quasi-static and dynamic experiments were conducted and evidenced similarity in the barrier behaviour (see the load-displacement curves of gure 3, note that for the dynamic test the load is calculated from the product mass time measured acceleration of the block.) However the net itself is not protected by the brakes and their fuses role. Therefore it needs to be underligned here that the dynamic aspects of the impact can not be predicted by the proposed quasi-static analysis which would fail to estimate maximum forces in the net due to back lash as well as local stress concentration around the block, which is known as "bullet eect" and tightly linked with the form and speed of the block (see for instance [START_REF] Hambleton | Perforation of exible rockfall barriers by normal block impact[END_REF]).

"Curtain eect" and sliding cable

As just said in the previous section, the modelling of the "curtain eect" is key to the understanding of the structure behaviour: when a block impacts the interception structure, the net slides along the support cables and tends to concentrate around the impacted zone.

The geometric reorganisation, created by the sliding of the net, allows to increase the motion (in the vertical direction) of the complete structure, hence minimising the internal forces.

This phenomenon, illustrated in Figure 4, is called "curtain eect" and is commonly modelled with a "sliding cable" [START_REF] Grassl | Steel-net rockfall protection: Experimental and numerical simulation[END_REF], whose tension depends on the variation of its whole length. The tension is hence constant in the cable and does not depend on the relative positions of successive nodes. We can nd a few models of such elements in the literature. The list of works presented here is not exhaustive. B. Zhou developed in 2004 a 3-nodes sliding cable model for an application to parachute systems [START_REF] Zhou | Finite element formulation for modeling sliding cable elements[END_REF]. A generalisation is proposed by C. Chen for a digital implementation in a FEM software in the purpose of studying the resistance of a dome [START_REF] Chen | Formulation and application of multi-node sliding cable element for the analysis of Suspen-Dome structures[END_REF].

L. Ghoussoub [START_REF] Ghoussoub | Analyse de quelques éléments du comportement des écrans de lets pare-blocs[END_REF] also proposed a generalisation of the sliding model of B. Zhou to n-nodes cable [START_REF] Zhou | Finite element formulation for modeling sliding cable elements[END_REF].

An alternative approach of the sliding of the cable through the post head was also proposed using dedicated Master and Slave nodes by Gentilini et al [START_REF] Gentilini | Three-dimensional numerical modelling of falling rock protection barriers[END_REF] and Tran et al [START_REF] Tran | Prototype of a wire-rope rockfall protective fence developed with three-dimensional numerical modeling[END_REF],

while Albaba et al proposed a rigid closed loop connections, formed by four nodes, and through which the cable passes [START_REF] Albaba | DEM Modeling og a Flexible Barrier Impacted by a Dry Granular Flow[END_REF].

Another part of its eld of use is also the case of cables passing through pulleys. This system is widely used in the industry to transfer loads on long and complex path. The particularity of pulleys is that they can apply a friction force on the cable so that tension is dierent on each side of it. M. Aufaure proposed a rst nite element model of a cable passing through one pulley to describe the deformation of a power line [START_REF] Aufaure | A nite element of cable passing through a pulley[END_REF][START_REF] Aufaure | Three-node cable element ensuring the continuity of the horizontal tension; A clamp-cable element[END_REF]. This work was then extended by F. Ju to identify the tangent stiness matrix of a cable passing through n pulleys. While staying in the framework of pulleys system, K. Hincz conducted a nonlinear analysis with the DR method to determine the geometry of an arch-supported cable net roof [START_REF] Hincz | Nonlinear Analysis of Cable Net Structures Suspended From Arches with Block and Tackle Suspension System, Taking into Account the Friction of the Pulleys[END_REF] and A. Bel Hadj et al proposed a model taking into account a friction force [START_REF] Bel Hadj | Analysis of clustered tensegrity structures using a modied dynamic relaxation algorithm[END_REF][START_REF] Bel Hadj | A dynamic-relaxation formulation for analysis of cable structures with sliding-induced friction[END_REF]. J. Coulibaly used also this approach to propose a model adapted to rockfall barriers and included friction forces [START_REF] Coulibaly | Sliding cable modeling: An attempt at a unied formulation[END_REF]. In Coulibaly et al [START_REF] Coulibaly | Sliding cable modeling: An attempt at a unied formulation[END_REF] as well as in Bel Hadj et al [START_REF] Bel Hadj | A dynamic-relaxation formulation for analysis of cable structures with sliding-induced friction[END_REF], the parametric studies show that friction might have a signicant inuence on the static equilibrium and dynamic of lightweight structures.

However, in the specic case studied here, the experiments on the rockfall barriers conducted in [START_REF] Olmedo | Extended experimental studies on rockfall exible fences[END_REF] have shown that the sliding occurs at the very beginning of the impact, so that it occurs on a conguration with low angles between cable successive elements (θ) and with low forces in cables (T ). Therefore one might estimate that the friction evaluated from the friction coecient µ as T (exp(µθ) -1) is low and that the associated dissipated energy is negligible compared to the energy dissipated by the brakes. Also the fact that the sliding cables are always linked with brakes at each end ensures that the magnitude of the tension in the cable end is xed by the brakes threshold which is quickly reached during the impact and almost independently of the friction coecient. A frictionless sliding cable should thus be able to capture the behaviour of the barrier at leading order.

Summary of contributions

The next two sections are hence dedicated to the discrete modelling of frictionless sliding cables and its numerical stability within the framework of a dynamic explicit scheme, with application to n-nodes monotonic supporting cables and ring nets. Some computations are carried in the fourth section on simple examples to assess the accuracy and stability of the numerical models. Finally, a case-study on a quasi-static experiment on a whole barrier prototype is proposed and demonstrates the interest of the frictionless proposed model. Some perspectives on the limitation of the present study and possible extensions will conclude the paper.

Modelling of a sliding cable element

As presented in introduction the "curtain eect" is currently modelled by means of a sliding cable element. Cables are connected to the rings along each edge of the net. In view of the diameters of cables and rings, a single contact point can be considered between these two elements. Such cables can either be modelled as continuous elements submitted to concentrated forces, or discretely as an ordered set of nodes linked together by spring interactions. The main purpose of this section is to identify the tangent stiness matrix of a frictionless sliding cable, which is necessary to determine a stability condition of the centred nite dierence scheme [START_REF] Papadrakakis | A method for automatic evaluation of the dynamic relaxation parameters[END_REF]. In a rst subsection, the expression of the internal forces vector will rst be identied as a function of displacements from the equilibrium of a continuous cable submitted to n concentrated forces. This straightforward approach is wellknown and will allow us to correctly introduce the notations required to dene the discrete internal forces vector. Moreover, it will be demonstrated that, with this case of concentrated forces, results do not depend on the description of the problem (continuous or discrete). The second subsection will be dedicated to the identication of the tangent stiness matrix of a sliding cable. From its expression, a stability condition of the nite dierences scheme will be proposed in a third section, easily calculable, in order to reduce computation costs.

2.1. Equilibrium of a sliding cable submitted to concentrated forces

Denition

We propose to write the equilibrium of a cable loaded by n concentrated forces (which may represent n rings connected to the cable in a rockfall barrier). A cable of length l is represented as a uni-dimensional element and parametrised by its arc-length, so that ∀ s ∈ [0; l], the position of a point s along the cable is given by: s → x(s). We also consider the n forces f i applied on the points of curvilinear abscissa: s 1 < ... < s n with s 1 = 0 and

s n = l. x i x i+1 x i-1 x 1 x n t 1 -t i-2 t i -t i -t n-1 f 1 f i-1 f i f i+1 f n -t i-1 t i-1 t i+1 -t 1 x 1 t 2 x n-1 f n-1 t n-1 -t n-2
t i = -t(s - i+1 ) = t(s + i ))
.

In the rest of this section, the following assumptions will be considered:

• The function s → x(s) is piece-wise dierentiable. The tangent vector is hence dened on each interval where f is dierentiable by: t : s → dx(s) ds .

• Since the bending and compression stiness are considered as null, the cable can only undergo tensile forces.

• The tensile force is uniform along the cable and given by: T (s) = T t(s) with s ∈ ]0, l[ when t is dened.

• Considering the large displacements of the cable, all the equilibrium equations are established on the deformed geometry.

Let

ds i ∈ 0, min 1 2 ||x(s i ) -x(s i-1 )|| , 1 2 ||x(s i+1 ) -x(s i )||
. Firstly, we write the equilibrium of a cable segment between two consecutive forces. The segment

[s i + ds i , s i+1 -ds i ]
is not loaded by any external force and its equilibrium is simply given by:

T t(s i + ds i ) = T t(s i+1 -ds i ) (1)
The cable is thus necessarily straight between two concentrated forces and

t(s + i ) = t(s - i+1 ) = x(s i+1 ) -x(s i ) ||x(s i+1 ) -x(s i )|| (2)
Then, the equilibrium of any segment [s ids i , s i + ds i ] on which the force f i is applied gives the change of tangent direction at the node:

T t(s + i ) -t(s - i ) = f i (3)
Finally, to describe the equilibrium of the whole cable, it still remains to write the equilibrium of its two extremities:

-T t(s + 1 ) = f 1 and T t(s - n ) = f n (4)
The deformed geometry of the cable is hence piece-wise linear and it is completely given by the position of the points on which the forces are applied (see Figure 5). The equilibrium of the whole cable F int = F ext writes hence:

T                 t(s + 1 )
. . .

t(s + i ) -t(s - i ) . . . -t(s - n )                 +                 f 1 . . . f i . . . f n                 = 0 (5)
The transition from the equilibrium of a continuous cable under multiple concentrated forces to the discrete model is thus obvious. The discrete model is completely determined by the ordered set of the n nodes positions [x(s 1 ) . . . x(s i ) . . . x(s n )] of the continuous cable which now will be simply noted [x 1 . . . x i . . . x n ] for a discrete cable. The tangent vectors can always be dened in the same way but for the sake of clarity, only the right side tangent will be considered:

t i = ∆x i l i with ∆x i = x i+1 -x i and l i = ||x i+1 -x i || (6)
The components of the vector t i in the direct orthogonal system (e 1 , e 2 , e 3 ) are respectively noted: t 1,i , t 2,i , t 3,i .

We assume then that the behaviour of the element is linear elastic and characterised by Hooke's law. Noting respectively E , S and , the Young's modulus, the cross section and the uniform strain of the cable, the tensile force is simply given by:

T = ES (7)
where can be evaluated from the deformed length of each segment and the global rest

length l 0 : = 1 2   n-1 i=1 l i l 0 2 -1   with l 0 = n-1 i=1 l rest i (8)
Considering the previous notations, the expression of the internal forces vector (5) in frictionless cables can be rewritten on the nodal displacements basis:

F int = ES                         -t 1,1 . . . t 1,i-1 -t 1,i t 2,i-1 -t 2,i t 3,i-1 -t 3,i . . . t 3,n-1                         (9)

Remarks

Before proceeding with the determination of the tangent stiness matrix, it is interesting to discuss the form of the internal forces vector [START_REF] Chen | Formulation and application of multi-node sliding cable element for the analysis of Suspen-Dome structures[END_REF].

Eulerian vs Lagrangian. As explained in the introductory paragraph of this section, similar studies to identify the tangent stiness matrix of a sliding cable have been already carried out. B. Zhou [START_REF] Zhou | Finite element formulation for modeling sliding cable elements[END_REF], C. Chen [START_REF] Chen | Formulation and application of multi-node sliding cable element for the analysis of Suspen-Dome structures[END_REF] and L. Ghoussoub [START_REF] Ghoussoub | Analyse de quelques éléments du comportement des écrans de lets pare-blocs[END_REF] have established an expression of the tangent stiness matrix for a 3 or n nodes sliding cable respectively; the developments of C. Chen and L. Ghoussoub being a generalisation of the problem dealt by B. Zhou.

Through their demonstration, they naturally identify an internal forces vector F int ante which can be compare to the present internal force vector F int present , such that:

F int present = l f inal /l initial 2 F int ante ( 10
)
The expressions of the internal forces vector ( 9) developed in the present paper thus dier from the one dened by C. Chen and L. Ghoussoub by the square of the lengths ratio. The study of the three demonstrations indicates that, while writing the virtual work of the element, a confusion between the Eulerian and the Lagrangian strains had been made, as well as in the choice of the reference conguration for the Cauchy stress tensor, with the consequence for the whole scheme not to be consistent. Indeed the Cauchy stress tensor σ lies in the deformed (Eulerian) conguration while the deformation considered by previous authors corresponds to the Lagrangian deformation tensor which is built on the initial conguration. The strain energy considered is thus an hybrid product of one Eulerian object and one Lagrangian object which is not consistent. To be correct, it requires to take into account the metric changes between Eulerian and Lagrangian deformation ds ds 2 = l initial /l f inal 2 which would cancel the additionnal term in [START_REF] Coulibaly | Sliding cable modeling: An attempt at a unied formulation[END_REF]. It is thus claimed that the actual formulation is the correct one. Moreover the consistency of the present model is highlighted by the fact that, in the case of concentrated forces, the discrete model becomes fully equivalent to the continuous description.

Discrete normal vectors. The main assumption of a sliding cable model involves that, the external forces are orthogonal to the cable tangent in the deformed geometry. The demonstration is relatively simple in the case of distributed forces and was done above for concentrated forces. The cable being not continuously dierentiable, it is impossible to dene the normal vector to the curve wherever a concentrated force is applied, but a discrete normal vector can be dened. Actually one of the usual approaches is to set it as the direction given by the bisector of two successive tangents (see Figure 5):

n i = t i+1 -t i ||t i+1 -t i || (11)
The norm of t i+1 -t i depends on the discrete turning angle θ, which is the angle between t i+1 and t i :

||t i+1 -t i || = 2| sin θ 2 | (12)
Comparing ( 11) with ( 9), we notice that the intensity of the internal forces depends on the local discrete normals of the cable. Therefore, in order to simplify notations in the rest of this section, even if the internal forces at the points x 1 and x n are not directed by the discrete normal (which can not be dened at these points), we will note N the vector such that: [START_REF] Chen | Formulation and application of multi-node sliding cable element for the analysis of Suspen-Dome structures[END_REF], so that N is dened by:

F int = ES N in
N T = [-t 1,1 , . . . , t 1,i -t 1,i-1 , t 2,i -t 2,i-1 , t 3,i -t 3,i-1 , . . . , t 3,n-1 ] (13)

Identication of the tangent stiness matrix

The expression of internal forces in a sliding cable as a function of its nodal displacements have been identied in the previous section. It should be noted that the description of the cable equilibrium was only a way of achieving this result and the expression can be extended to any kinematically admissible displacement. This section is dedicated to the identication of the sliding cable tangent stiness matrix. It is currently calculated from the expression of the internal forces vector (9) by:

K = grad U (F int ) ( 14
)
where U is the vector of the actual or eulerian nodal displacements. F int and U are vectors of size 3n and, to avoid any confusion, their components will be indexed with the index j ∈ [1, 3n], whereas the nodes of the discrete cable are indexed with the index i ∈ [1, n].

For example, considering that j = 3i and using the expression (9) it will be noted:

F int j = ES N j = ES (t 1,i-1 -t 1,i ) , F int j-1 = ES (t 3,i-2 -t 3,i-1 ) , F int j+1 = ES (t 2,i-1 -t 2,i ) (15)
The same notations are considered for the nodal displacements, and the components of the matrix K can now be expressed:

K jk = ∂F int j ∂U k = ES N j ∂ ∂U k + ∂N j ∂U k (16)
In order to simplify the rest of the development, K is considered as the sum of two square matrices: the elastic stiness matrix and the geometric stiness matrix, respectively noted

K E and K G : K E,jk = ESN j ∂ ∂U k , and K G,jk = ES ∂N j ∂U k (17)
It remains now to calculate the two partial derivatives. Firstly, we determine the expression of K E . The uniform strain of the cable is given by [START_REF] Chanut | Dynamic behavior of rock fall protection net fences: a parametric study[END_REF], which depends on the deformed length of each segment. According to the notations introduced in ( 6), it can be expressed as:

l i = ∆x 2 1,i + ∆x 2 2,i + ∆x 2 3,i (18)
We can hence write:

∂ ∂U k = l l 2 0 n-1 i=1 1 2l i ∂∆x 2 1,i ∂U k + ∂∆x 2 2,i ∂U k + ∂∆x 2 3,i ∂U k (19)
Now we use the nodal displacement indexing and Kronecker symbol δ jk to develop this expression:

∂ ∂U k = l l 2 0 3n-3 j=1 x j+3 -x j l j ∂x j+3 ∂U k - ∂x j ∂U k = l l 2 0 3n-3 j=1 x j+3 -x j l j δ (j+3)k -δ jk (20)
Finally, the coecients of the K E matrix are dened as follows:

K E,jk = ES l 0 N j                -l l 0 x k+3 -x k l k if k ∈ [1, 3] l l 0 x k -x k-3 l k-3 -x k+3 -x k l k if k ∈ [4, 3n -3] l l 0 x k -x k-3 l k-3 if k ∈ [3n -2, 3n] (21)
The expression of the discrete normal N is recognised in the expression after the bracket, so that:

K E = ES l 0 l l 0 N ⊗ N = K 0 l l 0 N ⊗ N (22)
The second part of the stiness matrix K G depends on the coecients ∂N j ∂U k . The global indexing is still used and, to clarify the rest of the developments, we consider that j is a multiple of 3, even if there are no diculties of computation in the other cases. We consider the general case where j ∈ [4, 3n -3] and once again dierent sub-cases must be considered.

∂N j ∂U k = ∂ ∂U k x j -x j-3 l j-3 - x j+3 -x j l j (23)
After developments,

∂N j ∂U k = δ jk l j-3 - ∂l j-3 ∂U k x j l 2 j-3 - δ (j-3)k l j-3 - ∂l j-3 ∂U k x j-3 l 2 j-3 - δ (j+3)k l j - ∂l j ∂U k x j+3 l 2 j + δ jk l j - ∂l j ∂U k x j l 2 j (24)
Using [START_REF] Grassl | Steel-net rockfall protection: Experimental and numerical simulation[END_REF], the coecients are given by:

• k = j: ∂N j ∂U j = l 2 j-3 -(x j -x j-3 ) 2 l 3 j-3 + l 2 j -(x j+3 -x j ) 2 l 3 j (25)
• k = j -3:

∂N j ∂U j-3 = (x j -x j-3 ) 2 -l 2 j-3 l 3 j-3 (26)
• k = j + 3:

∂N j ∂U j+3 = (x j+3 -x j ) 2 -l 2 j l 3 j (27)
• k = j -2:

∂N j ∂U j-2 = (x j -x j-3 )(x j+1 -x j-2 ) l 3 j-3 (28)
• k = j -1:

∂N j ∂U j-1 = (x j -x j-3 )(x j+2 -x j-1 ) l 3 j-3 (29)
• k = j + 1:

∂N j ∂U j+1 = - (x j -x j-3 )(x j+1 -x j-2 ) l 3 j-3 - (x j+3 -x j )(x j+4 -x j+1 ) l 3 j ( 30
)
• k = j + 2:

∂N j ∂U j+2 = - (x j -x j-3 )(x j+2 -x j-1 ) l 3 j-3 - (x j+3 -x j )(x j+5 -x j+2 ) l 3 j (31)
• k = j + 4:

∂N j ∂U j+4 = (x j+3 -x j )(x j+4 -x j+1 ) l 3 j ( 32
)
• k = j + 5:

∂N j ∂U j+5 = (x j+3 -x j )(x j+5 -x j+2 ) l 3 j (33) • k / ∈ [j -3, j + 5] ∂N j ∂U k = 0 (34)
The extension of previous calculations to the next two lines allows to construct the following 3 × 3 matrix, whose coecients are written with the index i running from 1 to n.

A i = - 1 l 3 i         ∆x 2 1,i -l 2 i ∆x 1,i ∆x 2,i ∆x 1,i ∆x 3,i ∆x 2 2,i -l 2 i ∆x 2,i ∆x 3,i sym ∆x 2 3,i -l 2 i         (35)
The matrix K G of size 3n × 3n is thus sparse and writes in accordance with the sub-matrix

A i as follows:

K G = ES                         A 1 -A 1 0 • • • 0 A 1 + A 2 -A 2 0 • • • 0 . . . . . . . . . • • • . . . A i-1 + A i -A i 0 • • • 0 sym . . . . . . • • • . . . A n-2 + A n-1 -A n-1 A n-1                         (36)
This nishes the determination of the tangent stiness matrix of a n-nodes sliding cable.

2.3. Stability condition of the Dynamic relaxation scheme 2.3.1. Denition of a general condition M. Papadrakakis proposed a general method to identify the critical parameters of a structure made of springs in the framework of the dynamic relaxation algorithm [START_REF] Papadrakakis | A method for automatic evaluation of the dynamic relaxation parameters[END_REF]. This method has been taken over by many authors, including M. Barnes who formulated a simple expression insuring the stability of the algorithm with kinetic damping [START_REF] Barnes | Form Finding and Analysis of Tension Structures by Dynamic Relaxation[END_REF]:

m i ≥ ∆t 2 2 K i max ( 37
)
where m i is the ctitious mass of node i , ∆t the time step and K i max corresponds to the highest stiness value of the spring elements connected with the node i.

About segments of vanishing length

As previously described, in the deformed geometry the external forces are directed toward the bisectors of the angles formed by two successive tangents. It is then easy to prove from this equilibrium condition that if two successive concentrated forces have the same direction, in the deformed geometry, their two points of application coincide (see Figure 6). Yet, as it will be proved further (see expression (51)), the critical ctitious mass of nodes depends, for a sliding cable element, on the inverse of the minimal length between two consecutive points. To avoid numerical instability due to the vanishing of a segment, a series of springs, which depend on the length of each segment, are introduced. The stiness of those springs is chosen in the form of a barrier function equal to zero when the length is above 5% of the initial segment length, and then equal to a logarithmic function of the length for lower values of the actual length. This numerical artefact helps convergence, avoids the collapse of any segment and visually does not aect the equilibrium shape (which can not distinguish between a 100% sliding and a 95% sliding).

Development of the stability condition

Hence, the stiness matrix K is dened as the summation of two real symmetrical matrices. K is thus diagonalisable and all its eigenvalues are real. In order to establish a simple analytic expression of K i max , the spectral radius of K E and K G , ρ K E and ρ K G respectively, are separately determined. Indeed, in the particular case of a real and symmetrical matrix, the spectral radius is equal to the matrix 2-norm, so that the triangle inequality leads to:

f f 2f (a) (b)
ρ K = ρ K E + K G ≤ ρ K E + ρ K G (38)
To evaluate these two spectral radii we rst establish that K E is a positive denite matrix.

Let X be a non-zero vector of size (3N). We may write:

X T • N ⊗ N T • X = X T • N N T • X = < N , X > 2 (39)
N being a strictly positive vector, K E is real and positive denite. Because all its eigenvalues are positive we have the obvious following upper bound:

ρ K E < tr K E (40)
From the analytic expression of N T (13), we obtain:

tr(N ⊗N T ) = 2+ n-1 i=2 ∆x 1,i l i - ∆x 1,i-1 l i-1 2 + ∆x 2,i l i - ∆x 2,i-1 l i-1 2 + ∆x 3,i l i - ∆x 3,i-1 l i-1 2 (41)
By developing the sum:

tr(N ⊗ N T ) = 2 + 2 n-1 i=2 1 - ∆x 1,i-1 l i-1 ∆x 1,i l i + ∆x 2,i-1 l i-1 ∆x 2,i l i + ∆x 3,i-1 l i-1 ∆x 3,i l i (42)
Considering the series (θ i ) of angles between two consecutive vectors along the cable (see Figure 7), we may write the previous expression as follows:

tr(N ⊗ N T ) = 2 n -1 - n-1 i=2 cos(θ i ) = 2 n -1 - n-1 i=2 1 -2 sin 2 θ i 2 (43)
Because the sinus values range from -1 to 1:

tr(N ⊗ N T ) ≤ 2 1 + 2 n-1 i=2 sin θ i 2 (44)
The expression (44) corresponds to the general case, where the nodes of the cable have three degrees of freedom in translation. And as it is, it can be implemented in the algorithm. In order to reduce the computational time, it is also possible to make simplifying assumptions and therefore to establish a simpler expression. In practice, the sliding cable is almost planar and its rotation angle θ is almost continuously increasing.. The value of its total curvature is thus lower than π (see Figure 7). Considering these assumptions, we know that all values of θ i are in the range of 0 to π and hence all values of sin (θ i ) are positive. Moreover, a current upper bound for the sinus is its argument. Therefore the expression (44) becomes in this particular case:

tr(N ⊗ N T ) ≤ 2 1 + n-1 i=2 θ i (45)
Then by using the regularity of the cable's curvature, we establish a simple upper bound of the spectral radius of the K E matrix: Now it remains to determine an upper bound for the spectral radius of K G dened in (36). To do this, we may use the Greshgorin circle theorem, which allows to calculate an upper bound of each eigenvalue of a diagonalisable matrix from its coecients:

ρ K E ≤ 2K 0 l l 0 (1 + π) (46) x 1 x n x i x i+1 θ i
ρ K G ≤ K 0 l 0 max i∈[1;N -1] n j=1 A j (47)
which can be developed using (36):

ρ K G ≤ K 0 l 0 max i 3 k=1 2 A i,(jk) + 2 A i-1,(jk) (48)
Then, using the expression (35) of the A i matrix, we get:

ρ K G ≤ K 0 l 0 max i 2 2 l i + 2 l i-1 (49)
Figure 8: Experimental tensile test [START_REF] Nicot | Etude du comportement mécanique des ouvrages souples fr protection contre les éboulements rocheux[END_REF] and discrete ring model. bending occurs. In the second stage, the length of the ring changes and axial stiness of the rope is activated. The bending of the ring can be reproduced thanks to diagonal struts as in [START_REF] Volkwein | Numerische Simulation von exiblen Steinschlagschutzsystemen[END_REF] but the strut stiness is dicult to evaluate. Therefore we prefer here to integrate directly the bending stiness in the tensile stiness of the ring which will hence become nonlinear as in [START_REF] Coulibaly | Nonlinear Discrete Mechanical Model of Steel Rings[END_REF]. This kind of element will be called "ring element" in the rest of this paper and its characteristics will be identied in section 5.3.2.

The calculations allowing the determination of the ring cable tangent stiness matrix are practically the same as those carried out in the case of the sliding cable. The only dierence with [START_REF] Day | An Introduction to Dynamic Relaxation[END_REF] is the expression of the normal vector N for the rst and last node which becomes:

N T = [t 1,n -t 1,1 , . . . , t 1,i-1 -t 1,i , t 2,i-1 -t 2,i , t 3,i-1 -t 3,i , . . . , t 3,n-1 -t 3,n ] (54)
with:

t n = x 1 -x n ||x 1 -x n || (55)
Taking into account this modication, the new expression of the ∂Nt ∂U matrix becomes:

∂N t ∂U =                         A n + A 1 -A 1 0 • • • 0 A 1 + A 2 -A 2 0 • • • 0 . . . . . . . . . • • • . . . A i-1 + A i -A i 0 • • • 0 sym . . . . . . • • • . . . A n-2 + A n-1 -A n-1 A n-1 + A n                        
where the A i matrices are dened as previously in (35).

The upper bound of K G spectral radius is the same as previously. However, due to changes in the expression of N , the calculation of the upper bound of ρ(K E ) is a little dierent. By keeping the same notations, we may rewrite the trace of N ⊗ N T :

tr(N ⊗ N T ) = 2 n-1 i=2 (1 -cos θ i ) (57)
There is no obvious simplifying assumption in the case of a n-nodes ring (a looped sliding cable with n nodes). In particular, the ring curvature may not be monotonic and the ring's convexity isn't ensured. Hence, the stability condition of a ring element depends on the number of its nodes:

tr(N ⊗ N T ) ≤ 8n (58)
Therefore, as previously the spectral radius of the tangent stiness matrix may be estimated of length h holding a beam of length d. The initial length of the sliding cable (which can be prestressed) is l 0 . We assume that the members stiness of the module are much higher than the one of the sliding cable, and thus the module's member can be considered non-deformable.

The equilibrium state of this structure can thus be analytically determined. The equilibrium of one of the sliding points is given by:

P = 2T sin(θ) (61)
Moreover, the tensile load in the cable can be expressed according to its deformation and prestress by:

T = K∆l + T 0 ( 62
)
where K is the stiness of the cable, ∆l its elongation and T 0 the initial prestress. By writing the cable deformation according to the angle θ (cf. Figure 9): and (P 2bis , A bis , P 1bis ) for the left and the right sides respectively. In addition to these four cables there are two symmetrical bars (A, B) and (B bis , A bis ), drawn in black.

l = f (θ) (63)
The two cables pairs, which constitute the structure, have the same stiness. The stiness of the two bars is variable in the study: it actually controls the relative sliding between the cables. The equilibrium state of this system is solved on one hand with the dynamic relaxation tool and on the other hand with a non-linear solver implemented in the soft-ware matlab R and using a trust-region algorithm (see J. Moré 1983 [21]). The problem is controlled in displacement by the vertical motion of the point C. The position of points P 1 , P 2 , P 2,bis and P 2,bis are xed. The resolution of the problem by the trust-region algorithm consists simply in solving a four-dimensional non-linear system, given by the equilibrium equations of points A and B and in which the unknown parameters are the positions of these two points. The relaxed lengths are given by the initial geometry which is described below, by taking into account the symmetry of the problem:

P 1 A P 2 B B bis P 2,bis A bis P 1,bis A B C B bis A bis 2T 1 T 1 T 1 + x T 1 + x T 2 T 1 + T 2
P 1 C = 7.5 m , P 2 C = 2.5 m , AC = 4.25 m , BC = 1.75 m ( 66
)
The four cables have the same bi-linear elastic behaviour with two stinesses respectively before and after the threshold: k c1 = 106 kN.m -1 and k c2 = 3, 6 kN.m -1 . They have also the same threshold: lim = 2.3%, so that the lateral cables reach it for a twice lower tensile stress. These bi-linear behaviours mimic the side cables in rockfall barrier which are generally linked with dissipating energy devices conferring to the set "cable + brakes" this material non-linearity. The bar's behaviour is elastic with a stiness of: k b = 0, 8 kN.m -1 for the rst two computation series whose results are presented in the gures 12 and 13. It illustrates the role of the cable net in the control of the sliding cable geometry.

Comparison of results

This simplied model of rockfall barriers presents some key aspects of their behaviour, especially the connection of two sliding cables via the net (here modeled by the bar AB).

Graphs in Figure 12 present the comparison of computation results obtained with the trustregion and the DR algorithms. We can observe that the kinematics of the structure given by the four unknowns of this static problem (Graph (b) and (c) of Figure 12) are identical with an error smaller than 0.1%. For very small vertical displacement of the structure (Y c < 20cm), we observe that motions of A and B are larger than during the rest of the loading path. It can also be noted that A moves away from the loaded zone during this rst displacement, while it moves toward the loaded zone afterwards.

Another signicant result concerns the whole behaviour of the structure (Graph (a) of Figure 12): the second phase of the behaviour, after passing the threshold of the central cable, is remarkably linear. Figure 13 highlights this impact of the yielding of the two central cables. The stiness of the whole structure increases until the central cables reach their threshold. After the threshold and during this linear phase, the apparent stiness remains constant. Thus we observe that, even if the tensile force in lateral cables is not negligible compared to the one in central cables (almost 50%), its inuence on the whole behaviour remains limited as it changes little the structure's geometry.

Convergence of the DR algorithm

The purpose of this last paragraph is to present some issues about the convergence and the stability of the DR algorithm. The calculations of the ctitious masses are given by the expression (52) and the current safety coecient g = 1.05 ([4, 14]). The parameters of the rst simulations series are used. We consider the equilibrium of the structure for a given vertical displacement of the point C: Y C = 5 m. The calculations are carried out by applying to the structure a given displacement or a given force. In each loading case, the load is reached in two dierent ways: in a series of 25 "incremental steps" (the step length is 0.2m for the study in displacement and 5.0kN for the study in force) and in a single step.

The results of these four series of simulations are presented in Table 1 and compared with the result given by the trust region method. We see that whatever the loading mode, the calculation converges toward the same solution with a high accuracy which highlights the relevance of dynamic relaxation when only one equilibrium positions is looked for. The results are presented in Figure 14. We see that the increase of the bar stiness has signicant inuence on the equilibrium geometry but relatively little on the reaction force.

F resultant (kN) Y C (m) X A (m) Y A (m) X B (m) Y B (
Points A and B slide towards each other, so that the lengths of the two central sliding cables increase as the lengths of the bars decrease.

Note that, with the parameters (stiness and geometry) chosen here, the range under 10 kN is not accessible with the two algorithms. Indeed, the positions of points P 1 , P 2 and C form a minimum sliding angle which is reached when the bar is completely elongated and points A and P 1 as well as B and P 2 coincide respectively (see Figure 15). Considering this minimum angle and the tensile force in the central cables, it is possible to dene the minimum force the bar has to apply to the point B to balance it. To this end, we write the equilibrium of the point B located at the same place as the point C. Considering the notations of the Figure 15, where the force x applied by the bar to the point B and the angle α 3 are the two unknown parameters of the system:

       T 1 cos(α 1 ) + (T 1 + x) cos(α 2 ) = 2T 1 cos(α 3 ) T 1 sin(α 1 ) + (T 1 + x) sin(α 2 ) = 2T 1 sin(α 3 ) (67)
Considering this nal geometry, we nd that the equilibrium of the point B can be ensured by the force x = 3, 09 kN and the angle α 3 = 42˚. At this point, two remarks should be made. First, this force is indicative and in practice the variations of the bar stiness change the equilibrium geometry presented in Figure 15. Second the geometry of Figure 15 is the exact solution of the problem corresponding to the limit case (for which the stiness of the bar is null). Due to the nature of the sliding cable which ensures the same tension between each couple of nodes and because in practice two points can not perfectly coincide, this limit case can not be reached numerically.

P 1 P 2 α 2 α 1 T 1 + x B α 3 C 2T 1 T 1
Figure 15: Limit equilibrium with a very low stiness of the bar.

Conclusion

Through these two comparative studies, we have tested on simple cases the accuracy and the stability of our numerical models of sliding cables. For both comparisons with analytic models (4.1) or with computations carried out with another algorithm (4.2), the results of the numerical simulations are in perfect agreements. This allows us to condently model a complete barrier.

5. A case study: numerical simulation of a full-scale experiment

Structure description and presentation of the experimental conditions

In this section, we present the numerical simulation of one of the experiments carried out in the framework of the French national project C2ROP (www.c2rop.fr) which brings together many public and private partners around the topics of landslide risk and protection devices. The numerical tool, initially developed in the laboratoire Navier for elastic gridshells [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF], works within the framework of Rhinoceros 3D and its plugin Grasshopper.

Barrier architecture

The scheme of the tested barrier is presented in Figure 16 and has been detailed in [START_REF] Coulibaly | Sliding cable modeling: An attempt at a unied formulation[END_REF].

A three-modules barrier is anchored perpendicularly to a vertical cli. The whole fence is 2.75m high and 15m long (5m for each modulus). This fence consists of 10 support cables (4 on each edge parallel to the cli and 1 on each lateral edge). With this layout, the sliding of the net is possible along each edge. Because of the large span of such structures, the wiring along the edge is generally complex. It is presented in Figure 17. It is symmetric between the upstream and downstream edges. The two lateral cables (in green in Figure 17) hold the small edges of the net. The two external cables link the lateral anchor with the head of the central post (cables A and A bis in blue) and two central cables connect the head of the rst and third posts (cable B and B bis in red). Each cable holds partially the net and is free near the post's head to avoid stress concentration. With this assembly, the net is hold by two cables in the centre of the modulus and by a single cable near posts. The set formed by the net and the support cables is hold to the cli by 4 steel posts (modelled by beam elements). Finally, the head of each post is maintained by 3 cables (two upstream and one downstream), anchored to the cli (see Figure 16). Brake elements (so called snake brakes in which a cable winds rubbing along a set of screws) are connected in series with each support cable. Unlike the case studies carried out previously on the sliding cable, the ring cable model is calibrated thanks to a comparison with a plane tensile experiment. This experiment was also carried out in the framework of the french national project C2ROP. The experimental setup consists of a square ring net (three by three) which is pulled in one direction and maintained in the perpendicular one (cf. Figure 19). The test is driven in displacement and load cells measure reaction forces in the two directions. The results are given according to the axial displacement of the moving edge (see Figure 20).

Upstream Downstream B sup bis B inf bis H 2 H 1 B inf B sup H lat

Numerical simulation and description of the ring model

The basic pattern of a single ring is described in Figure 8. By assuming that the ring deformation is uniform in the square net, the tensile behaviour of a single ring can be identied from the experimental results in the loading direction. Indeed, by referring to Figure 20, the equilibrium of the point A gives the tensile force in the ring, while its strain is calculated from the axial displacement. The constitutive law is then continuously identied, thanks to the least squares method, by a ve order polynomial:

K 0 ( ) = 3.56.10

where is the strain of the ring, which is dened according to d ring , the initial diameter of a ring, and ∆L, the overall axial elongation.

= d 2 ring + (d ring + ∆L 3 ) 2 - √ 2d ring √ 2d ring (69)
We note that a 5 degree polynomial is enough to correctly reproduce the curve behaviour.

However this polynomial ts the experimental behaviour only for the strains reached during the test. The behaviour modelled numerically is extended by the tangent at the maximum strain really allowed for larger strains (the rupture of the elements is not taken into account in the numerical computations). We remark then that the behaviour t on the axial load is also in accordance with the experimental results in the transverse direction. Indeed, for small axial displacements, it is obvious that the Poisson's ratio is close to one for both continuous real ring and discrete model (see Figure 20). Moreover, since for large displacements, the deformed geometry of the real ring and discrete model are the same, the Poisson's ratio observed in the simulations is always in agreement with the experiment. The material stiness of the ring is linearly linked to its dimensions so that the expression (68) of the constitutive law does not depend on the net density, which is practically the same as the one of the real net. All information relative to the elements stiness can now be regrouped in Table (

T 3 A N L 0 L 0 + ∆ l (a) (b)

Comparison with the numerical simulation

This section compares the experimental results and the numerical simulations for the overall response of the barrier (the resultant force applied to the winch) and also for the tensile load in two barrier components (support cables). In all this section we will refer to the drawing (17) which gives a nomenclature of the barrier components.

Overall barrier behaviour

Rather than comparing the resultant force applied to the block centre to the experimental force applied to the winch, we prefer here to consider the global stored elastic energy. Experimental results are thus integrated to the vertical displacement of the winch (see Figure 21). We observe that the numerical simulations are in good agreement with the experiment. For both, experimental and numerical results, the resultant force applied to the winch increases almost linearly according to the block's vertical displacement. This result is surprising. Indeed the large displacements in the structure and the sliding of the net along the support cables involves necessarily geometrical non-linearities. Moreover, the mechanical behaviour of the support cables is also non-linear. The complex overall behaviour of this structure, which is nevertheless well reproduced by the numerical model, requires more data to be more accurately described including the evolution of the local stress according to the actual geometry of the barrier. Detailed investigation of the behaviour is foreseen in a next paper.

Tensile stress in the central support cables

The tensile stress is recorded for the two (upstream and downstream) central cables called B in Figure 17 and located on the left side of the barrier. Because the experimental model is not completely symmetric there is a signicant dierence between the tensile stress in the cables B and B bis (probably because of some asymmetry in the prestress level). For this reason, we prefer to compare an average of the tensile stresses in the two upstream and downstream cables (see Figure 22). Further on, names of cables will refer to the Figure 18. We can observe that, for the two edges, the measures on the cables B and B bis give relatively remote results, particularly for the second stage of the behaviour during which the brakes scroll. The slight asymmetry of the structure has a limited impact on numerical results. We can also note that the eective threshold levels are dierent from the expected value of 25kN , in particular for the downstream support cables. This is probably due to existing prestress in the experimental structure prestress which is not taken into the simulation but which is implemented in the model and could in principle be introduced into the simulation). The comparison between the experiment and the simulations is hence carried out from the averaged results. For the downstream cables, the threshold of the brakes elements seems to happen for the same displacement of the block and the behaviours after it are similar. However, the rst stiness are dierent between both curves: in the experiment, during a rst stage the net is the only barrier component which is activated and the forces are then faster transferred to the support cable, whereas this transfer is gradual in the numerical simulations. Regarding now the tensile load in the upstream support cables, we observe that, this time, the results are close, even if the stiness before and after the thresholds of the brakes are slightly dierent.

The numerical results presented here are thus in good agreement with the experimental results. The dynamic relaxation method associated with frictionless cable reproduces well the behaviour of the whole structure and gives good results for the components considered separately. This validates a posteriori the hypothesis that, for the type of structures tested here, the friction does not have a signicant inuence on the equilibrium conguration and in the forces in the supporting members. It can thus be considered as a valuable alternative to nd very quickly equilibrium shapes of rockfall barriers in an early design phase and to test the behaviour of a large variety of alternative designs.

Furthermore, we recall here that this result is also a consequence of the existence of brakes attached at the end of the sliding cables that cut o dynamic peaks in the supporting structure. As said in the introduction, the net is not protected by the brakes and thus see the dynamic eects of the loading. Therefore, in early phase design, we advise that the engineer complete the quasi-static simulations by a verication of the bullet eect following for example [START_REF] Hambleton | Perforation of exible rockfall barriers by normal block impact[END_REF]. Then, for detailed design of the barrier, only full dynamic simulations with friction along the sliding cables will give the required accuracy for a reasonnably optimised sizing of the structural members. At this late design phase, there sould no more be any need for testing of alternative structures and the computational cost of the complete dynamic simulations should not be a problem.

Conclusion

In this paper we proposed a continuous model for frictionless sliding cables and its discretization. It was used for a simple estimation of the sliding cable mass parameters ensuring the stability of the nite dierence scheme in the framework of the dynamic relaxation method. This estimation was determined by evaluating an upper bound of the spectral radius of the tangent stiness matrix for two specic cases: the n-nodes sliding cable and the ring cable which is formed by closing a sliding cable. The numerical performance of these two elements have been tested and the numerical simulation of a quasi-static loading test on a complete barrier prototype was shown. Although quasi-static experiments on whole structure are rare in the literature, they are a necessary step to carry out dynamic simulations, because they allow to verify the assembly of the barrier components model, which are frequently calibrated by means of quasi-static tests. Hence in these case-studies, the simulations have shown very good agreement with reference problems.

In further work, dynamic simulations will be carried out by extending the proposed stability conditions to estimate the maximum time step ensuring the numerical stability. The extension is straight forward if lumped mass matrices are assumed, but requires some developments in the general case. The nal aim of the study would be to conduct a comparative study between the quasi-static and the dynamic calculations to further investigate the similarities and dierences of the two approaches and to produce equivalent acceleration models for the design of these structures.

Figure 1 :

 1 Figure 1: Top view of the rockfall barrier prototype developed in the framework of the french national project C2ROP (Photo courtesy of NP C2ROP).

Figure 2 :

 2 Figure 2: Pattern of a ASM4 ring net.

Figure 3 :

 3 Figure 3: Comparison of quasi-static /dynamic experimental global behaviour of the C2ROP fence (by courtesy PN C2ROP).

Figure 4 :

 4 Figure 4: "Curtain eect" due to block impact (by courtesy PN C2ROP).

Figure 5 :

 5 Figure 5: Equilibrium of a sliding cable under concentrated forces (assumingt i = -t(s - i+1 ) = t(s + i ))

Figure 6 :

 6 Figure 6: Limit of the frictionless sliding model under parallel forces.

Figure 7 :

 7 Figure 7: Sliding cable with continuously increasing rotation angle θ.

Figure 9 :

 9 Figure 9: Analytic model before and after loading.

Figure 10 : 4 . 2 .

 1042 Figure 10: First case study "simple cable": comparison between analytic and numerical results with and without prestress.

Figure 11 :

 11 Figure 11: Second case study: Multiple sliding cables model.

Figure 12 :

 12 Figure 12: Comparison of the DR method with a non-linear algorithm implemented in matlab.

Figure 13 :

 13 Figure 13: Cables load path for increasing displacements in the multiple sliding cables model.

4. 2 . 4 .

 24 Parametric study of the bar stiness In a second simulations series, we propose to study the impact of the parameter k b (the stiness of the bar) on the behaviour of the structure. The vertical positions of the point C is xed at the altitude Y C = 5 m and the stiness parameter changes in the range 10 3 N < k b < 10 8 N . Once again the computations are carried out with the two algorithms.

Figure 14 :

 14 Figure 14: Inuence of the bar stiness on the behaviour of the structure.

Figure 16 :

 16 Figure 16: Technical drawing of the tested barrier (Courtesy of PN C2ROP).

Figure 18 :

 18 Figure 18: Front view of the quasi-static experiment on the rockfall barrier prototype (Courtesy of PN C2ROP).

Figure 19 :

 19 Figure 19: Experimental setup of the plane tensile test on a ASM4 ring net (Courtesy of PN C2ROP).

Figure 20 :

 20 Figure 20: Experimental results of the plane tensile test.

Figure 21 :

 21 Figure 21: Experiment test vs. numerical simulation.

Figure 22 :

 22 Figure 22: Central support cables behaviour: experiment test vs. numerical simulation. (a) : Resultant force applied to the block, (b) : Elastic strain energy of the whole structure.

Table 1 :

 1 Comparative convergence study with DRM.

	m)

  9 5 -6.80.10 8 4 + 3.95.10 7 3 -4.07.10 5 2 + 3.55.10 3

Table 2 :

 2 2). Synthese of all barrier components stiness.

	Component	behaviour		Characteristics
	sliding cable + brake	bi-linear	before threshold threshold level after threshold
	(upstream/downstream)	in traction	ES = 1140kN	25kN	ES = 35kN
	sliding cable + brake	bi-linear	before threshold threshold level after threshold
	(lateral cables)	in traction	ES = 1190kN	25kN	ES = 35kN
	12 anchoring cables	traction only		ES = 1260kN
	rings	polynomial	5 th degree, t on the plane tensile test (68)
	bars of the block	compression only		ES = 10 7 kN

Acknowledgements

The experiments presented in this paper were carried out in the framework of the French national project C2ROP. The authors wish to thank all the members of the national project and particularly the members of the working group on soft barriers for their invaluable support, namely M-A. Chanut and J. Coulibaly from CEREMA, P. Robit, I. Olmedo from GTS, C. Galandrin from CAN and D. Bertrand from INSA Lyon.

By introducing l min the length of the smallest segment, the following obvious upper bound can be established:

Finally by grouping the upper bounds of the two spectral radii, we nd the following inequality:

According to (37), the optimal nodal masses of a sliding cable in the dynamic relaxation scheme are thus dened by:

Note that, conversely, the optimal time step of a sliding cable with n equal sliding masses in a centred explicit dynamic scheme would be dened by:

The ring cable model

A second application of the sliding cable model is proposed in this section. The ASM4 ring net is indeed composed of interlaced rings whose behaviour can be reproduced with a discrete model of a looped sliding cable (the square drawn in red in Figure 8). The discrete model used here is inspired from those of A. Volkwein and D. Coulibaly [START_REF] Volkwein | Numerische Simulation von exiblen Steinschlagschutzsystemen[END_REF][START_REF] Coulibaly | Nonlinear Discrete Mechanical Model of Steel Rings[END_REF] in which the ring is composed by the four contact points with the neighbour nodes.

The behaviour of one of these rings undergoing a two points tensile test is presented in Figure 8 taken from [START_REF] Nicot | Etude du comportement mécanique des ouvrages souples fr protection contre les éboulements rocheux[END_REF]. During the rst stage its perimeter remains the same and only by the following expression:

Then by injecting this expression of K max in (37), we can determine the critical ctitious mass in the case of a 4-nodes ring cable:

From these two relations on the spectral radius of the sliding cable and the ring cable, we are able to ensure the stability of these two elements implemented inside a dynamic relaxation code. The modelling and the computations of a whole barrier does not present other major diculties. In the next section, we hence propose two simple case studies to verify the results

given by the sliding cable model and to assess the accuracy and eciency of the method.

Then, a more elaborated example will be shown in the fth section. The structure studied here was already used as elementary problem by A. Volkwein and L. Ghoussoub in their PhD [START_REF] Volkwein | Numerische Simulation von exiblen Steinschlagschutzsystemen[END_REF][START_REF] Ghoussoub | Analyse de quelques éléments du comportement des écrans de lets pare-blocs[END_REF]. It is formed by a four points sliding cable and a module transporting a load. This module is xed to the two middle points of the cable and represents the loading undergone by a net during a block impact. The whole structure is hold by the two extremities of the sliding cable as shown in Figure 9. The module is formed by two bars

and injecting this expression in equation ( 61), we nd the relation between the weight of the module and the geometric deformation of the cable:

Comparison between analytic and numerical results

The structure presented in previous paragraph is modelled with its loading inside the dynamic relaxation tool. Simulations are carried out for dierent values of the load P to observe the evolution of the system. The initial length of the cable is l 0 = 40m and its stiness is K = 1250N.m -1 , the dimensions of the module are h = 11m and d = 14m. In the rst simulations series the sliding cable is not prestressed (T 0 = 0). Figure 10 shows the evolution of the P load according to the θ angle for the numerical and analytic results. We remark that for a non-prestressed cable the numerical results perfectly correspond to the analytic results, with an almost 0.3% error on the whole load path.

We then carry out the same simulations adding prestress to the cable with a tension of T 0 = 2000N . Once again, we observe that the numerical results are in good agreement with the analytic model, and this time, the average error is less than 0.1%. The stiening eect of the prestress is well reproduced by the model. For this case study (with or without prestress) the computational tool reproduces successfully the geometrical non-linearities due to the large vertical deformation of the cable.

Loading

The position of every post foot is xed but the rotation of the post is free. The positions of all cable anchors on the cli are also xed. The links between cables and posts or cables and brakes are hinges. The sliding of the net on the surrounding cables is free.

A quasi-static loading is then applied to the structure. A 740kg normalised polyhedralshaped concrete block is maintained by a winch and slowly placed in the middle of the net until the barrier reaches equilibrium under dead-weight and block load. Then, another winch is hooked on the bottom part of the block and pulled orthogonal to the net (see Figure 18).

Its vertical displacement is controlled and step by step increased so that the loading can be considered as quasi-static. Several load cells are put on the fence components enabling the stresses in all the support cables and in some anchoring ropes to be recorded (see black ags in Figure 17). The resultant load applied to the winch is also recorded.

Superstructure and loading modellings

The experimental loading being monotonic, the behaviour of all the elements is considered reversible (elastic but not necessarily linear). The relaxed length of each barrier component is given by the initial geometry of the structure (no prestress is taken into account). The stiness are identied from preliminary experiments on the components (brake, cable and net) separately. The behaviour of the ring net is studied in detail in section 5.3.2. Sliding cables and brakes being set in series in the experiment, they are merged in the numerical model into a unique bi-linear sliding cable element for simplicity reasons. The rst stiness corresponds to the sti behaviour of the brake (before its scrolling), the second stiness is low in comparison and corresponds to the brake's scrolling. The behaviour of anchoring cables is simply modelled thanks to a linear law giving the stress according to the relative displacement of the two cable endpoints. The steel posts are modelled as truss elements with a similar linear law.

The block shape used for the numerical simulations is not polyhedral but spherical. Its diameter is 0.85m and corresponds to the sphere tangent to the edges of the polyhedron.

Contact is modelled with many bar elements, which work only in compression and whose relaxed lengths correspond to the block radius. The tensile stiness is null and the behaviour is linear elastic in compression. The stiness is determined according to the maximum tensile force applied to the structure and the tolerance of the interpenetration level. With a maximum force of almost 210kN and a tolerance of 1cm, the stiness of the block bar is set to: K block = 5.10 6 kN/m. As in the experiment, the simulations are carried out by imposing the block displacement (the convergence is quicker than by imposing a force on the block).


T Laigre
email: t.laigre@brgm.fr

Y Balouin
email: y.balouin@brgm.fr

D Villarroel-Lamb
email: deborah.villarroel-lamb@sta.uwi.edu.

M Moisan

Primary Topic: Coastal observations and monitoring Secondary Topic: Foreshore, beach and dune morphodynamics OBSERVATION OF SHORELINE AND UPPERBEACH VEGETATION EVOLUTION BY VIDEO: FIELD ASSESSEMENTS IN THE CARIBBEAN

Introduction

Caribbean coastal ecosystems such as coral reefs, mangroves, seagrasses or upper beach vegetation acts as a physical barrier against waves and can play a major role in reducing the impact of coastal hazards for small islands communities. However they are highly impacted by human activities and global changes and their coverage is dropping in the Caribbean as well as in the rest of the world. Consequently, the protection effect is declining as the ecosystems deteriorates or dies. Therefore, there is a growing interest for management interventions that may restore or enhance the coastal protection value of these ecosystems based on research in the last decades related to physical processes that govern these systems. Nevertheless, upper beach vegetation is not as wellinvestigated and a very few studies examine their effect on erosion and run-up dampening, especially in a tropical context. The aim of this paper is to assess upperbeach vegetation/shoreline evolution using high frequency video imagery.

Field site

Anse-Maurice is a narrow beach located on the eastern shore of Grande-Terre Island in Guadeloupe Archipelago (France) in the Lesser Antilles. The beach consist of a small lagoon bounded by a narrow chaotic fringing reef. The reef is characterized by an assembly of mainly dead colonies of Acropora Palmata with an algae cover, though the bottom subtract remain very complex. The upper beach vegetation is composed by an association of indigene species (i.e. Ipomoea pes capre, Coccoloba uvifera) and hexogen species (mainly coconut trees). This ecosystem is impacted by human trampling and goat grazing which motivated the ONF (The French Forestry Agency) to implement a nature-based solution by establishing revegetation enclosures. The site is exposed to strong Atlantic swells during the winter season (from December to March) and episodically during the cyclonic season (from July to November).

Instruments and methodology

Two low-cost Solarcam® videos monitoring systems had been installed to provide a full range of view from the incident ocean waves to the upper beach. Those devices take a picture every 10 minutes. "Timestacks" were generated, which represent the evolution over time of a specific pixel line on a picture. In our case, after a georectification process (see [START_REF] Holland | Practical use of video imagery in nearshore oceanographic field studies[END_REF] at mid/long term scales, this technique allows for shoreline and upperbeach vegetation monitoring. This is possible because of the sharp colorimetric difference between water and sand on the one hand, and sand and vegetation in the other. Wave conditions for the period Hs, Tp and Dp are extracted from spectral wave model (WW3) simulations provided through the MARC platform (https://marc.ifremer.fr/). From Hs and Tp, wave power is calculated and presented in figure 2a. Sea level conditions are extracted from the tide gauge of Pointe à Pitre (https://data.shom.fr/).

Results and discussion

The figure 2.c) and d) represent 591 days of monitoring. Preliminary results show several periods of vegetation progression and regression as well as shoreline change. However, while the vegetation limit is stable overtime, the shoreline is in net retreat. Excluding the anthropic vegetation removal on the 13/11/19, rapid vegetation and shoreline changes are correlated with extreme swells events. Whereas progressive changes (over several months) seems to be partially linked to the sea level meso-cycle. Indeed, Furthermore extreme swell events have a greater impact when associated with high sea level (on the 05/08/ or 18/09/20 for example) while events with relatively low sea level (on the 13 and 20/01/20 for example) do not induce important changes. Analyses are ongoing to evaluate the role of seasonal trends in sea levels and wave direction on shoreline evolution, and assess the role of Sargassum beaching in these tendencies.

Figure 1

 1 Figure 1. a) Geographic location of the field site (red square), swell (green point) and sea level (yellow star) data extraction points, b) overall map of the field site with revegetation enclosure position in red polygons and camera positions in black points. Black and yellow lines indicate the position of presented cross-shore profiles (resp. fig 2 c and d).

Figure 2

 2 Figure 2. a) Significate wave height evolution, b) sea level evolution, c) timestack evolution of the black line in fig.1.b) d) ) timestack evolution of the yellow line in fig.1.b).

Primary Topic: Coastal observations and monitoring Secondary Topic: Foreshore, beach and dune morphodynamics

Acknowledgements

This work is part of the EU INTERREG Project CARIB-COAST. TL acknowledges financial support from BRGM, UWI, Guadeloupe Region, ODE Guadeloupe and Martinique, DEAL Martinique and EU through a PhD grant and technical support from M. Norden and T. Delahaye. The video system used in this work was funded by DEAL Guadeloupe, ADEME and BRGM.


Thibault Laigre

Yann Balouin

Manuel Moisan

Manuel Moisan Measurements

Measurements and modelling of the effect of a fringing reef and upperbeach vegetation on runup behaviour -A Caribbean Context

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Measurements and modelling of the effect of a fringing reef and upperbeach vegetation on runup behaviour -A Caribbean Context

T. Laigre, Y. Balouin, D. Villarroel-Lamb, M. Moisan

Site description

Pocket beach with a small lagoon and a fringing reef Numerical modelling

• Runup analysis on a specific storm event

Observations

• Daily extraction of maximum swash on images from camera over mid/long term (768 days)

Runup comparison on two profiles

Winter swells Cyclonic season

Winter swells

  Atlantic winter storms (December-March) • Episodically exposed to Cyclone (July-November) Revegetation enclosure from the French forestry agency Subject to sargassum beaching since low-cost video systems : • Short term (Storm scale) : • Identification of swash events • Identification overall agitation offshore and on the lagoon • Mid/long-term : • Upper beach vegetation limit monitoring • Swash maximum excursion monitoring Pre and post-storm RTK-GPS monitoring : • Ground truth for extractions from cameras •

Conclusions & Perspectives

Numerical modelling during storm :

• Numerical modelling sensitivity test with varying reef and vegetation properties in order to quantify the benefits (loss) of ecosystem state of health changes. • Modelling probabilistic events : 5,10,100… years return storms +sea level rise scenarios (IPCC)

Perspectives :

Vegetation enclosure at Anse-Maurice


Thibault Laigre
email: t.laigre@brgm.fr

Yann Balouin

Nico Valentini

Deborah Villarroel- Lamb

Manuel Moisan

Effects of a fringing reef and upperbeach vegetation on coastal flooding in a Caribbean context

Keywords: Caribbean, Numerical modelling, Coral reef, Upperbeach vegetation, Nature-Based Solutions

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Effects of a fringing reef and upperbeach vegetation on coastal flooding in a Caribbean context

Thibault Laigre, Yann Balouin, Nico Valentini, Deborah Villarroel-Lamb,

Manuel Moisan

To cite this version: The proposed presentation deals with the assessment of coral reef and upperbeach vegetation impact on coastal flooding at Anse Maurice, a reef-fringed pocket beach located at Guadeloupe Island, in the Caribbean region. Several field-monitoring campaigns have been conducted with the aim of characterizing the beach evolution. Upperbeach vegetation evolution and marine inundation was assessed using a fixed video system. In order to highlight the processes induced by the reef and upperbeach vegetation from the front reef to the nearshore focusing on runup levels, a numerical modelling experiment using XBeach in its non-hydrostatic mode was undertaken. Data from coastal ecosystems and morphology were used as inputs. Observations acquired from a hydrodynamic field campaign in September/October 2020 and the above-mentioned video system were used as calibration.

The latter is particularly useful to the analyses of swash events occurrence and to identify the overall agitation in the lagoon. Analyses of processes and runup attenuation by the reef and upperbeach vegetation during the distant Hurricane Teddy (September 2020) pointed out the role of reef and upperbeach vegetation ecosystems. Modelling results confirm observations, and clearly indicate that coastal inundation is mostly driven by storm events. During the storm events, the reef acts as a protection by dissipating most of the short wave energy, but also induces the generation of infragravity waves that appear to be the most important factor driving the berm overtopping. In that case, the vegetation acts as another barrier to the flow.

To assess the effect of ecosystems health, several simulations on the same storm event were undertaken with varying ecosystem properties. It was concluded that wave runup at the coast is firstly dependent on the reef properties, but the upperbeach vegetation also clearly influences the maximum swash excursion at the beach and the resulting coastal flooding.