Dario Benedetto
email: benedetto@mat.uniroma1.it

AND Thierry Paul
email: thierry.paul@upmc.fr

Stefano Rossi
email: stefano.rossi@math.uzh.ch

PROPAGATION OF CHAOS AND HYDRODYNAMIC DESCRIPTION FOR TOPOLOGICAL MODELS

Keywords: January 8, 2024. MSC2020 subject classifications. 35Q92, 35Q83, 82B40, propagation of chaos, monokinetic solutions, topological interaction, Cucker-Smale model

In this work, we study the deterministic Cucker-Smale model with topological interaction. Focusing on the solutions of the corresponding Liouville equation, we show that propagation of chaos holds. Moreover, considering monokinetic solutions, we also obtain a rigorous derivation of the hydrodynamic description given by a pressureless Euler-type system.

Introduction

In recent decades, physics of complex systems has increasingly dealt with the description of groups of animals exhibiting collective behavior, such as flocks of birds, fish schools, locust swarms, and migrating cells ( [START_REF] Emlen | Flocking behaviour in birds[END_REF], [START_REF] Cucker | Emergent behavior in flocks[END_REF], [START_REF] Lopez | Theraulaz From behavioural analyses to models of collective motion in fish schools Interface focus[END_REF], [START_REF] Bazazi | Collective motion and cannibalism in locust migratory bands[END_REF], [START_REF] Giniūnaitė | Modelling collective cell migration: neural crest as a model paradigm[END_REF]).

From a modeling point of view, these systems offer new challenges and various models have been proposed to describe their interaction (see for example [START_REF] Vicsek | Shochet Novel type of phase transition in a system of self-driven particles[END_REF], [START_REF] Bertin | Grégoire Boltzmann and hydrodynamic description for selfpropelled particles[END_REF], [START_REF] Vicsek | Zafeiris Collective motion[END_REF]). Among the first to be introduced, the Cucker-Smale model ( [START_REF] Cucker | On the mathematics of emergence Japanese[END_REF]) describes a bird as a self-propelling particle interacting with its neighbors. In this case, the interaction is such that neighboring birds tend to align their velocities and the strength of the interaction is described through weights which depend on the metric distance between the agents.

Around 2008, a new type of interaction between agents called "topological interaction" was introduced. In [START_REF] Ballerini | Interaction ruling animal collective behavior depends on topological rather than metric distance: evidence from a field study[END_REF][START_REF] Cavagna | Tavarone From empirical data to inter-individual interactions: unveiling the rules of collective animal behavior Math[END_REF], the CoBBS group in Rome, after collecting 3D observational data for flocks of starlings, observed that regardless of the density of the flock, each agent interacts on average with its first 6-7 neighbors. This suggests that the strength of the interaction between agents does not depend on the metric distance between them, but rather on the "topological" distance that takes into account the proximity rank of the latter with respect to the former (see also [START_REF] Bode | Wood Limited interactions in flocks: relating model simulations to empirical data[END_REF][START_REF] Camperi | Spatially balanced topological interaction grants optimal cohesion in flocking models[END_REF][START_REF] Ginelli | Chaté Relevance of metric-free interactions in flocking phenomena[END_REF][START_REF] Niizato | Emergence of the scale-invariant proportion in a flock from the metric-topological interaction[END_REF][START_REF] Shang | Bouffanais Consensus reaching in swarms ruled by a hybrid metrictopological distance[END_REF][START_REF] Martin | Multi-agent flocking under topological interactions[END_REF][START_REF] Shang | Bouffanais Influence of the number of topologically interacting neighbors on swarm dynamics[END_REF][START_REF] Wang | Chen Synchronization of multi-agent systems with metric-topological interactions[END_REF]).

In the following years, deterministic and stochastic models with topological interactions were introduced in the mathematical physics literature. In [START_REF] Blanchet | Degond Topological interactions in a Boltzmann-type framework[END_REF][START_REF] Blanchet | Degond Kinetic models for topological nearest-neighbor interactions[END_REF][START_REF] Degond | Propagation of chaos for topological interactions[END_REF][START_REF] Degond | Propagation of chaos for topological interactions by a coupling technique Atti Accad[END_REF] kinetic models of Boltzmann type are derived for topological interaction models based on jump processes. As far as deterministic models are concerned, in [START_REF]Haskovec Flocking dynamics and mean-field limit in the Cucker-Smale-type model with topological interactions[END_REF] a Cucker-Smale model is introduced where the interaction, instead of being metric as usual, is topological. In [START_REF]Haskovec Flocking dynamics and mean-field limit in the Cucker-Smale-type model with topological interactions[END_REF], the kinetic mean-field and hydrodynamic equations of Euler type are also written and derived in the case of a smoothed version of the model. Indeed, from a mathematical point of view, topological interactions fall outside the case of two-body interaction and present various problems in the derivation of mean-field and hydrodynamic equations. Specifically, the continuity estimate à la Dobrushin (see [START_REF] Dobrušin | Vlasov equations Funktsional[END_REF]), valid in the metrical case of regular two-body interactions, does not work here since the solutions are not weakly continuous with respect to the initial data, as shown in [START_REF] Benedetto | Mean-field limit for particle systems with topological interactions Math[END_REF]. In [START_REF] Benedetto | Mean-field limit for particle systems with topological interactions Math[END_REF] the existence of the dynamics and the mean-field limit have been rigorously proved for this same model, considering solutions of the limit equation with bounded density and proving that, for positive times t, the Wasserstein distance between the limit solution and the empirical measure at time t tends to zero as the number of particles increases, if this holds at the initial time.

In the present work, we focus on the problems of propagation of chaos and derivation of the hydrodynamic equations for this model. The starting point is the Liouville equation verified by the N -particle system and the aim is to show that the marginals of the N -body distribution function converge to tensorial powers of solutions of a suitable kinetic equation of Vlasov type. With the same approach, we show the validity of the hydrodynamic equations which describe the evolution of monokinetic initial data. The analysis will follow the approach in [START_REF] Golse | On the Mean Field and Classical Limits of Quantum Mechanics[END_REF] (see also [START_REF] Natalini | On the mean field limit for Cucker-Smale models Discrete Contin[END_REF], [START_REF] Natalini | The Mean-Field Limit for Hybrid Models of Collective Motions with Chemotaxis[END_REF], [START_REF] Menci | Microscopic, kinetic and hydrodynamic hybrid models of collective motions with chemotaxis: a numerical study Math[END_REF] for numerical considerations and [START_REF] Paul | Trélat From microscopic to macroscopic scale equations: mean field, hydrodynamic and graph limits[END_REF] for the case of agent systems).

Due to the aforementioned difficulties, the topological nature of the problem requires different techniques for the proof of these two results. Regarding propagation of chaos, the analysis makes use of another distance between probability measures, called discrepancy distance, in addition to the Wasserstein one. Furthermore, being the dynamics well-defined only for almost all initial data, we need to consider a regularization of the monokinetic initial datum in the proof of the hydrodynamic derivation, obtaining the result for any limiting point with respect to the regularization.

The outline of the paper is as follows: in Section 2 we recall the topological Cucker-Smale model and introduce the associated Liouville equation, as well as recalling the results obtained in [START_REF] Benedetto | Mean-field limit for particle systems with topological interactions Math[END_REF] which will be useful later on. In Section 3 we provide a proof of the propagation of chaos, which will be a direct consequence of the validity of the law of large numbers. In Section 4 we focus on the derivation of the hydrodynamic description given by Euler-type equations studying the so-called monokinetic solutions.

Model and general framework

A Cucker-Smale type model for the motion of N agents, in the mean-field scaling, is the system

$ ' & ' % 9 x i ptq " v i ptq 9 v i ptq " 1 N N ÿ j"1 p ij pv j ptq ´vi ptqq, i " 1, . . . , N (2.1)
where px i , v i q P R d ˆRd (d " 1, 2, 3, . . . ) and the "communication weights" tp ij u N i,j"1 are positive functions that take into account the interactions between agents. In classical models, p ij depends only on the euclidean distance |x i ´xj | between the agents. In topological models the weights depend on the positions of the agents through their rank:

p ij -K `M px i , |x i ´xj |q ˘, (2.2)
where K : r0, 1s Ñ R `is a positive Lipschitz continuous non-increasing function, and, for r ą 0, the rank function

M px i , rq - 1 N N ÿ k"1 X t|x k ´xi | ď ru (2.3)
counts the number of agents at a distance less than or equal to r from x i , normalized with N . Note that in this case p ij is a piecewise-constant function of the positions of all the agents.

We indicate by PpR k q the space of probability measures on R k . In the mean-field limit N Ñ `8, the one-agent distribution function f t " f pt, x, vq P PpR 2d q is expected to verify the equation

B t f pt, x, vq `v ¨∇x f pt, x, vq `∇v ¨pW rSf t , f t spx, vqf pt, x, vqq " 0, (2.4)
where Sf t pxq -ş f t px, vq dv P PpR d q denotes the spatial distribution and, given a probability measure f P PpR d ˆRd q and a probability measure ρ P PpR d q, W rρ, f s is the mean-field interaction given by (2.8)

W rρ, f spx, vq - ż R d ˆRd K pM rρspx,
In other words, f t is the push-forward of f 0 along the flow generated by the one-particle system (2.8), where the force W depends on f t itself.

Given Z N " px 1 , . . . , x N , v 1 , . . . , v N q P R dN ˆRdN , we define the empirical measure on R d ˆRd as

µ Z N - 1 N N ÿ i"1 δ x i b δ v i . ( 2
9 v i ptq " W rSµ Z N ptq , µ Z N ptq spx i ptq, v i ptqq, i " 1, . . . , N. (2.10)
We indicate by Z N pt, Z N q " pX N pt, Z N q, V N pt, Z N qq the solution of this system with initial datum Z N P R dN ˆRdN .

In [START_REF] Benedetto | Mean-field limit for particle systems with topological interactions Math[END_REF], in the framework of the mean-field theory, the rigorous derivation of (2.4) starting from (2.10) is obtained. More precisely, the following theorem is proved.

Theorem 0. [3, Theorems 3.4, 4.3 and 5.2] Let K : r0, 1s Ñ R `in (2.5) be a positive Lipschitz continuous non-increasing function. It holds that: i) except for a set of initial data Z N P R dN ˆRdN with Lebesgue measure zero, there exists a unique global in-time solution

pX N pt, Z N q, V N pt, Z N qq P C 1 pR `, R 2dN q ˆCpR `, R 2dN q of (2.10). Moreover, if |x i | ď R x and |v i | ď R v for any i P t1, . . . N u, we have that |x i ptq| ď R x `tR v , |v i ptq| ď R v , i P t1, . . . , N u. (2.11)
ii) Let f 0 px, vq P L 8 pR d ˆRd q be a probability density such that supppf 0 q Ă BRx ˆB Rv . Given T ą 0, there exists a unique solution, in the weak sense defined by (2.7) and (2.8), of the kinetic equation (2.4) with initial datum f 0 , such that f P C w `r0, T s; L 8 pR d ˆRd q ˘, i.e. f t is a probability measure with bounded density for any t ě 0, weakly continuous in t.

Moreover, for any t ě 0, supppf t q Ă BRx`tRv ˆB Rv and }f t } 8 ď }f 0 } 8 e dγt , (2.12)

where γ " ş 1 0 Kpzq dz. iii) If f 0 is as in item ii), for any t ą 0 the empirical measure µ Z N ptq associated to the system (2.10) weakly converges to f t in the limit N Ñ `8, provided this is true at time zero and the initial data Z N are chosen such that the dynamics exists, as in item i).

In the present work, we focus on the statistical description of the dynamical system (2.10), considering, at time 0, N independently and identically distributed particles with law F N p0, Z N q " f bN 0 . At time t the particles are distributed with the law F N pt, Z N q, weak solution of the N -body Liouville equation

B t F N pt, Z N q `N ÿ i"1 v i ¨∇x i F N pt, Z N q `N ÿ i"1 ∇ v i ¨´1 N N ÿ j"1 K ´M rSµ Z N spx i , |x i ´xj |q ¯pv j ´vi qF N pt, Z N q ¯" 0,
(2.13) in the sense that F N pt, Z N q is the push-forward of F N p0, Z N q along the flow Z N pt, Z N q. Note that F N pt, Z N q, for t ą 0, is symmetric in the exchange of particles.

In the next section we show that, in the limit N Ñ `8, for any s ě 1 the s-particles marginals of F N pt, Z N q factorize, and the limit is described by the solution of the mean-field equation f t with initial datum f 0 . To quantify the convergence, we use the Wasserstein distance W 1 defined as follows. Let µ, ν P PpR k q be two probability measures with finite first moments,

W 1 pµ, νq :" inf πPCpµ,νq ż R k ˆRk |x ´y| dπpx, yq, (2.14)
where Cpµ, νq is the set of all couplings of µ and ν, i.e. probability measures on the product space with marginals µ and ν in the first and second variable, respectively.

Convergence of the marginals

In this section we prove the propagation of chaos for solutions of the N -body Liouville equation (2.13). We briefly explain the meaning of this expression and we refer to [START_REF] Chaintron | Diez Propagation of chaos: A review of models, methods and applications. I. Models and methods[END_REF][START_REF] Chaintron | Diez Propagation of chaos: A review of models, methods and applications[END_REF] for a review, and to [START_REF] Nguyen | Shvydkoy Propagation of chaos for the Cucker-Smale systems under heavy tail communication[END_REF] for a propagation of chaos result in the case of non-topological Cucker-Smale models. We introduce the s-particles marginals as follows, with an abuse of notation regarding the order of the space and velocity variables in F N pZ N q:

F N :s pZ s q " ż F N pZ s , z s`1 , . . . , z N qdz s`1 . . . dz N , s " 1, 2, . . . , N,
where z i " px i , v i q P R 2d . We expect that, if N is large, the details of the individual mutual interactions are negligible, and the description given by F N :s pt, Z s q is similar to the one given by f bs t , where f t weakly solves (2.4).

Note that f bs t pZ s q weakly solves

B t f bs t pZ s q `s ÿ i"1 v i ¨∇x i f bs t pZ s q `s ÿ i"1 ∇ v i ¨´W rSf t , f t spx i , v i qf bs t pZ s q ¯" 0.
This last equation describes the law of the system when each particle evolves independently from the others, with the interaction given by the mean-field force defined in (2.5). Then propagation of chaos holds if F N :s pt, Z s q converges to f bs t pZ s q for any s ě 1.

We will prove the following result.

Theorem 1 (Propagation of chaos for the topological CS model). Assume that the interaction function K is as in Theorem 0, namely a positive Lipschitz continuous non-increasing function. Let f P C w pr0, T s; L 8 pR 2d qq be a weak solution of the kinetic equation (2.4) with initial datum f 0 px, vq P L 8 pR d ˆRd q such that supppf 0 q Ă BRx ˆB Rv . Consider F N ptq weak solution of (2.13) such that F N p0q " f bN 0 . For any integer s P t1, . . . , N u, it holds that sup tPr0,T s W 1 pF N :s ptq, f bs ptqq ď se λpT qKp1`}f 0 }8e dγT q a C d pN q, where γ :" ş 1 0 Kpzq dz, λpT q is a constant depending on d, R x , R v and T , and K :" maxp1, LippKq, }K} 8 q, where LippKq is the Lipschitz constant of K, while

C d pN q :" $ ' & ' % N ´1{2 if d " 1 N ´1{2 logpN q if d " 2 N ´1{d if d ą 2.
(3.1)

Remark 1. In the proof we will need a quantitative version of the law of large numbers, i.e. an estimate of ErW 1 pρ, µ X N qs when X N " px 1 , . . . , x N q are N independently and identically distributed R d -valued random variables with law ρ P L 8 pR d q. This is a widely studied problem in probability, in connection with topics of optimal transport, random matching and combinatorial optimization (see, for instance, the bibliographical notes to Chapter 6 of [START_REF] Villani | Optimal transport: old and new Grundlehren der mathematischen Wissenschaften 338[END_REF]). Here we use the explicit bounds reported by Fournier and Guillin in [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure Probab[END_REF], which imply that ErW 1 pρ, µ X N qs ď cC d pN q, where C d pN q is defined in (3.1) and c is a constant depending polynomially on the diameter of supppρq and on d.

Remark 2. The time dependence of the constant λpT q in Theorem 1 comes from the growth of the support of the distribution, and therefore grows polynomially in time. In particular, it is also possible to prove the validity of the theorem for times T slowly growing with N .

In the proof of Theorem 1 we will employ as a technical tool the discrepancy distance, defined as follows:

Dpµ, νq :" sup

where c is a constant that depends on the dimension d. ii) Given X N " px 1 , . . . , x N q, Y N " py 1 , . . . y N q, with |x i ´yi | ď δ for some δ and any i " 1, . . . N , consider the two empirical measures µ X N and µ Y N . Then, for any probability measure ρ P L 8 pR d q supported on a ball BR ,

Dpµ X N , µ Y N q ď cR d´1 }ρ} 8 δ `cDpµ Y N , ρq. (3.4)
In the following lemma, we summarize the technical details of the proofs concerning the Lipschitz estimates of the interaction function W defined in (2.5).

Lemma 1. Given r x , r v ą 0, let f 1 , f 2 P PpR d ˆRd q be two probability measures supported in Brx ˆB rv and ρ 1 , ρ 2 P PpR d q two probability measures supported in Brx and such that }ρ 2 } 8 ă `8. For any ξ 1 , ξ 2 P Brx and η 1 , η 2 P Brv it holds:

|W rρ 1 , f 1 spξ 1 , η 1 q ´W rρ 2 , f 2 spξ 2 , η 2 q| ď 2r v LippKqDpρ 1 , ρ 2 qp cLippKq}ρ 2 } 8 r d´1 x r v `}K} 8 q pW 1 pf 1 , f 2 q `|ξ 1 ´ξ2 | `|η 1 ´η2
|q , where c is a constant that depends only on d.

Proof of Lemma 1. By the triangle inequality, |W rρ 1 , f 1 spξ 1 , η 1 q´W rρ 2 , f 2 spξ 2 , η 2 q| can be bounded by the sum of

A 1 :" |W rρ 1 , f 1 spξ 1 , η 1 q ´W rρ 2 , f 1 spξ 1 , η 1 q|, A 2 :" |W rρ 2 , f 1 spξ 1 , η 1 q ´W rρ 2 , f 2 spξ 1 , η 1 q|, A 3 :" |W rρ 2 , f 2 spξ 1 , η 1 q ´W rρ 2 , f 2 spξ 2 , η 2 q|.
Estimate of A 1 : by the definition of M in (2.6) and the definition of discrepancy distance in (3.2), we have for any r ě 0

|M rρ 1 spξ 1 , rq ´M rρ 2 spξ 1 , rq| ď Dpρ 1 , ρ 2 q.
By the definition of W in (2.5) and the hypothesis on the supports, we get that A 1 is bounded by

A 1 ď 2LippKqr v Dpρ 1 , ρ 2 q.
Estimates of A 2 and A 3 : it is easy to prove that, given x P R d and r 1 , r 2 ą 0, |M rρ 2 spx, r 1 q ´M rρ 2 spx, r 2 q| ď c}ρ 2 } 8 ˇˇr d

A 2 ď pcLippKq}ρ 2 } 8 r d´1 x r v `}K} 8 qW 1 pf 1 , f 2 q, and A 3 ď pcLippKq}ρ 2 } 8 r d´1 x r v `}K} 8 qp|ξ 1 ´ξ2 | `|η 1 ´η2 |q.
Proof of Theorem 1. For any Σ N " py 1 , . . . , y N , w 1 , . . . , w N q P R dN ˆRdN , we consider Σ N pt, Σ N q " py 1 ptq, . . . , y N ptq, w 1 ptq, . . . , w N ptqq P R dN ˆRdN ,

where, for i P t1, . . . , N u, py i ptq, w i ptqq has initial datum py i , w i q and evolves independently with the mean-field interaction: " 9 y i ptq " w i ptq 9 w i ptq " W rSf t , f t spy i ptq, w i ptqq, i " 1, . . . , N.

We associate to such Σ N pt, Σ N q the empirical measure µ Σ N ptq as in (2.9).

We define a coupling π N t between F N ptq and f bN t in the following way: at time t " 0, it is given by π N 0 pZ N , Σ N q :" f bN 0 pZ N qδpZ N ´ΣN q. For positive times, π N t is given by the push-forward of π N 0 along the product of the flows given by (2.10) and (3.7), i.e. for any ϕ P C b pR 2dN ˆR2dN q ż ϕpZ N , Σ N q dπ N t pZ N , Σ N q " ż ϕpZ N pt, Z N q, Σ N pt, Z N qq f bN 0 pZ N q dZ N .

Next, given i P t1, . . . , N u, we introduce the quantity D N ptq :"

ż ´|x i ´yi | `|v i ´wi | ¯dπ N t pZ N , Σ N q,
which does not depend on i, thanks to the symmetry of the law. We prove the weak convergence of the s-marginals F N :s to f bs by showing that D N ptq Ñ 0: namely, using the symmetry of π N t ,

W 1 pF N :s ptq, f bs t q ď s ÿ i"1 ż ´|x i ´yi | `|v i ´wi | ¯dπ N t pZ N , Σ N q ď sD N ptq,
where we used that |Z s ´Σs | ď

ř s i"1 p|x i ´yi | d `|v i ´wi | d q.
From the definition of π t , we have, for any i P t1, . . . , N u,

D N ptq " ż ´|x i ptq ´yi ptq| `|v i ptq ´wi ptq| ¯df bN 0 pZ N q,
where px i ptq, v i ptqq for i " 1, . . . N solves (2.10) with initial datum Z N , and py i ptq, w i ptqq for i " 1, . . . N solves the decoupled system (3.7) with the same initial datum Z N . It follows that

D N ptq ď ż δpt, Z N q df bN 0 pZ N q, (3.8)
with δpt, Z N q :" max i"1,...,N p|x i ptq ´yi ptq| `|v i ptq ´wi ptq|q.

Since px i ptq, v i ptqq and py i ptq, w i ptqq have the same initial conditions, it holds that |x i ptq ´yi ptq| `|v i ptq ´wi ptq| ď ż t 0 |v i pτ q ´wi pτ q| dτ `ż t 0 ˇˇW rSµ Z N pτ q , µ Z N pτ q spx i pτ q, v i pτ qq ´W rSf τ , f τ spy i pτ q, w i pτ qq ˇˇdτ.

By Lemma 1 with r v " R v , r x " R x `τ R x , we bound the last integrand by

2R v LippKqDpSµ Z N pτ q , Sf τ qp cLippKq}Sf τ } 8 pR x `τ R v q d´1 R v `}K} 8 q `W1 pµ Z N pτ q , f τ q `δpτ, Z N q ˘.
By choosing the coupling π " 1 N ř N i"1 δ x i pτ q δ v i pτ q δ y i pτ q δ w i pτ q in definition (2.14), we can estimate the Wasserstein distance between the two empirical measures µ Z N pτ q and µ Σ N pτ q by δpτ, Z N q, so that W 1 pµ Z N pτ q , f τ q ď δpτ, Z N q `W1 pµ Σ N pτ q , f τ q.

By the triangle inequality

DpSµ Z N pτ q , Sf τ q ď DpSµ Z N pτ q , Sµ Σ N pτ q q `DpSµ Σ N pτ q , Sf τ q, and, by (3.4) with ρ " Sf τ , DpSµ Σ N pτ q , Sµ Z N pτ q q ď cpR x `τ R v q d´1 }Sf τ } 8 δpτ, Z N q `cDpSf τ , Sµ Σ N pτ q q, where, by (3.3),

DpSf τ , Sµ Σ N pτ q q ď c b pR x `τ R v q d´1 }Sf τ } 8 W 1 pSf τ , Sµ Σ N pτ q q.
In the sequel we indicate by λpτ q any positive, increasing polynomial function of τ , that depends on d, R v , R v , i.e. on the support of f τ and Sf τ , and with c any constant that depends at most on the dimension d.

Collecting the previous estimates, and using that x 1{2 ď p1 `xq{2 for x ě 0, we get, for a suitable λpτ q, δpt, Z N q ď ż t 0 Kλpτ qp1 `}Sf τ } 8 q ´δpτ, Z N q `W1 pµ Σ N pτ q , f τ q `bW 1 pSµ Σ N pτ q , Sf τ q ¯dτ, where we used that K " maxp1, LippKq, }K} 8 q. Note that }Sf τ } 8 ď cR d v }f τ } 8 and, by Theorem 0, }f τ } 8 ď e dγτ }f 0 } 8 . By the Grönwall's lemma and the fact that x ă e x , δpt, Z N q ď e KλpT qp1`}f 0 }8e dγT q ż t 0 ´W1 pµ Σ N pτ q , f τ q `bW 1 pSµ Σ N pτ q , Sf τ q ¯dτ.

(3.9) To estimate (3.8), we have to evaluate the expected value of the integrand in (3.9) w.r.t. the probability measure df bN 0 pZ N q. Since the empirical measure µ Σ N pτ q is supported in Σ N pt, Z N q, then by (2.7),

ż W 1 pµ Σ N pτ q , f τ q df bN 0 pZ N q " ż W 1 pµ Z N , f τ q df bN τ pZ N q ď λpτ qC 2d pN q,
where the rate C d pN q, as defined in (3.1), is obtained by using the Fournier and Guillin bound in [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure Probab[END_REF], and λpτ q grows polynomially in τ . By concavity we also obtain

ż b W 1 pSµ Σ N pτ q , Sf τ q df bN 0 pZ N q " ż a W 1 pSµ Z N , Sf τ q df bN τ pZ N q ď λpτ q a C d pN q.
We finally arrive at D N ptq ď e KλpT qp1`}f 0 }8e dγT q ´aC d pN q `C2d pN q ¯, from which the thesis follows.

Euler systems associated to monokinetic initial data

In this section we study the hydrodynamic description of the topological Cucker-Smale model, by considering the following pressureless Euler-type system for pρpt, xq, upt, xqq : r0, T s ˆRd Ñ R ˆRd :

$ ' ' ' ' ' & ' ' ' ' ' % B t ρpt, xq `∇x ¨pρpt, xqupt, xqq " 0 B t upt, xq `pupt, xq ¨∇x qupt, xq " ż R d K ´M rρptqspx, |x ´y|q ¯pupt, yq ´upt, xqqρpt, yq dy pρp0, xq, up0, xqq " pρ 0 pxq, u 0 pxqq, (4.1)
where KpM rρsq is defined as in (2.5) and (2.6), and pρ 0 pxq, u 0 pxqq is a regular compactly supported initial datum.

In the following, given T ą 0, we assume that (4.1) admits a unique solution pρpt, xq, upt, xqq which is regular and compactly supported for t P r0, T s. We leave this fact as a hypothesis, however see [START_REF] Ha | Kwon A hydrodynamic model for the interaction of Cucker-Smale particles and incompressible fluid Math[END_REF] for a proof in the case of the non-topological Cucker Smale model (see also [START_REF] Karper | Trivisa Hydrodynamic limit of the kinetic Cucker-Smale flocking model Math[END_REF][START_REF] Tadmor | Swarming: hydrodynamic alignment with pressure Bull[END_REF][START_REF] Lu | Tadmor Hydrodynamic alignment with pressure II. Multi-species Quart[END_REF]).

It is not difficult to show that, if pρpt, xq, upt, xqq is a regular solution of (4.1), then f pt, x, vq :" ρpt, xqδpv ´upt, xqq

is a weak solution of the topological Vlasov equation in the sense of (2.7), (2.8) with measure-valued initial datum f p0, x, vq " ρ 0 pxqδpv ´u0 pxqq: these are called monokinetic solutions. Note that in this case the field W rSf t , f t spx, vq is regular, so the flow in (2.8) is well-defined.

We will show how to obtain solutions of the Euler system (4.1) starting from solutions of the topological Liouville equation (2.13). The main obstacle is that, as stated in Theorem 0, the flow for the particle system is not defined for every initial datum, so in general it is not possible to consider a measurevalued solution of the Liouville equation (2.13). To overcome this problem we consider a regularization of the monokinetic initial datum:

f ε 0 px, vq " pρ 0 pxqδpv ´u0 pxqqq ˚x,v η ε px, vq, (4.3)
where η is a C 8 pR 2d q compactly supported mollifier and η ε pzq " ηpz{εq{ε 2d .

Our next result shows that after relaxing the regularity of the initial monokinetic measure as in (4.3), in the limit N Ñ `8, uniformly in ε, the marginals of solutions of the Liouville equation tend, in weak sense, to the tensorial powers of a monokinetic measure (4.2), built up out of the solution of the associated Euler system.

Theorem 2 (From particles to Euler). Let us consider a regular solution pρpt, xq, upt, xqq, t P r0, T s, of the Euler system (4.1) with regular compactly supported initial data pρ 0 px, vq, u 0 px, vqq. We suppose moreover that, for any t P r0, T s the supports of ρpt, xq and upt, xq are contained in BRx`tRv and that }upt,

xq} 8 ď R v , for some R x , R v ą 0. Let F ε
N ptq be a weak solution of the Liouville equation (2.13) with initial datum F ε N p0q :" pf ε 0 q bN , with f ε 0 as in (4.3).

Then, for any t P r0, T s and any s P t1, . . . , N u, setting V s " pv 1 , . . . , v s q, we have

W 1 ˆżR ds F ε N :s ptq dV s , ρptq bs ˙ď se KλpT q ş T 0 p1`}ρτ }8q dτ ´aC d pN q `4 ? ε ¯,
and, for any Lipschitz function Φ : R sd ˆRsd Ñ R, setting X s " px 1 , . . . , x s q, ˇˇˇˇż dX

s dV s ΦpX s , V s qF ε N :s pt, X s , V s q ´ż dX s s ź i"1 ρpt, x i qΦpX s , u 1 pt, x 1 q . . . u s pt, x s qq ˇˇˇď LippΦqse KλpT q ş T 0 p1`}ρτ }8q dτ ´aC d pN q `4 ? ε ¯,
where λpT q is a constant depending on d, R x , R v , T and C d pN q is defined in (3.1).

Remark 3. Note that the statement doesn't involve the Vlasov equation and that ε and N are independent, so the convergence holds also for N Ñ 8 and any sequence of limiting points as ε Ñ 0 of the distributions tF ε N ptqu N PN . In particular the rate of convergence is Op a C d pN qq along any sequence with ε " ε N " OpC d pN qq 2 . Note the algebraic behavior in N , in contrast with the logarithmic one, when the dynamics is mollified as, e.g., in [START_REF] Paul | Trélat From microscopic to macroscopic scale equations: mean field, hydrodynamic and graph limits[END_REF][START_REF] Bodineau | Saint-Raymond The Brownian motion as the limit of a deterministic system of hard-spheres Invent[END_REF].

Remark 4. Along the same lines, we can also generalize the content of Theorem 2, proving a convergence result for marginals F ε N :s of the Liouville equation with general initial conditions pg ε 0 q bN such that W 1 pg ε 0 , ρ 0 pxqδpv ´u0 pxqqq Ñ 0 as ε Ñ 0.

To prove Theorem 2, we need the following stability result for solutions of the Vlasov system. Proposition 1. Given T ą 0, for i P t1, 2u, let f i P C w pr0, T s, PpR d ˆRd qq be two weak solutions of the topological Vlasov equation (2.4) with initial data f i 0 P PpR d ˆRd q such that Sf i t are well-defined and belong to L 8 pR d q. Assume moreover that supppf i t q Ă BRx`Rvt ˆB Rv for t ą 0 and R x , R v ą 0. Then, for t P r0, T s,

W 1 `f 1 t , f 2 t ˘ď e KλpT q min i"1,2 ş T 0 p1`}Sf i τ }8q dτ max " W 1 pf 1 0 , f 2 0 q, b W 1 pf 1 0 , f 2 0 q * , (4.4)
where K :" maxp1, LippKq, }K} 8 q and λpT q is a constant depending on d, R x , R v and T .

The proof, which requires a nontrivial extension of inequality (3.4), is given in the Appendix.

Proof of Theorem 2. To prove the thesis, it is sufficient to establish the following estimate for t P r0, T s:

W 1 pF ε N :s ptq, f bs t q ď se KλpT q ş T 0 p1`}ρpτ q}8q dτ ´aC d pN q `4 ? ε ¯, (4.5)
where f is the monokinetic solution (4.2) of the Vlasov equation associated to the solution pρptq, uptqq of the Euler system (4.1) and λpT q is a constant that depends on

d, R v , R x `T R v .
Let f ε ptq be the solution of the topological Vlasov equation (2.4) with initial datum f ε 0 . By the triangle inequality, we have

W 1 ´F ε N :s ptq, f bs t ¯ď W 1 ´F ε N :s ptq, pf ε t q bs ¯`sW 1 ´f ε t , f t ¯,
where we are using that W 1 ppf ε t q bs , f bs t q ď sW 1 pf ε t , f t q. The second term is managed by the stability estimate (4.4) in Proposition 1, which gives

W 1 ´f ε t , f t ¯ď e KλpT q ş T 0 p1`}ρpτ q}8q dτ max ! W 1 pf ε 0 , f 0 q, a W 1 pf ε 0 , f 0 q ) . (4.6)
By choosing πpdz, dz 1 q " ε ´2d ηppz ´z1 q{εqf 0 pzq dz dz 1 in the definition (2.14) of the Wasserstein distance, we can estimate W 1 pf ε 0 , f 0 q with ε ş |z|ηpzq dz. Then for ε ă 1, we estimate the maximum in (4.6) by c ? ε, where c depends only on η.

The estimate of the first term is similar to the one given in the proof of Theorem 1. In this case we have

W 1 ´F ε N :s ptq, pf ε t q bs ¯ď s ż R 2dN δ ε pt, Z N q dpf ε 0 q bN pZ N q,
where δ ε pt, Z N q :" max i"1,...,N p|x i ptq ´yε i ptq| `|v i ptq ´wε i ptq|q, with the independent flow Σ ε ptq :" py ε i ptq, w ε i ptqq " py ε i pt, Z N q, w ε i pt, Z N qq solving " 9 ˇˇW rSµ Z N pτ q , µ Z N pτ q spx i pτ q, v i pτ qq ´W rSf ε τ , f ε τ spy ε i pτ q, w ε i pτ qq ˇˇdτ.

y ε i ptq " w ε i ptq 9 w ε i ptq " W rSf ε t , f ε t spy ε i ptq, w ε i ptqq, i " 1, . . . ,
In order to avoid terms in }Sf ε τ } 8 , which could diverge for ε Ñ 0, we use carefully the triangle inequality in the last integrand, which we bound by |W rSµ Z N pτ q , µ Z N pτ q spx i pτ q, v i pτ qq ´W rρ τ , µ Σ ε N pτ q spx i pτ q, v i pτ qq|| W rρ τ , µ Σ ε N pτ q spx i pτ q, v i pτ qq ´W rSf ε τ , f ε τ spy ε i pτ q, w ε i pτ qq|.

By using that W 1 pµ Z N pτ q , µ Σ ε N pτ q q ď δ ε pτ, Z N q and K " maxt1, LippKq, }K} 8 u, by Lemma 1, we estimate the first term with 2R v KDpSµ Z N pτ q , ρ τ q `Kλpτ qδ ε pτ, Z N q, and the second one with

2R v Dpρ τ , Sf ε τ q `Kλpτ q `W1 pµ Σ ε N pτ q , f ε τ q `δε pτ, Z N q ˘. (4.7)
We now estimate the discrepancies in the previous two equations. By the triangle inequality DpSµ Z N pτ q , ρ τ q ď DpSµ Z N pτ q , Sµ Σ ε N pτ q q `DpSµ Σ ε N pτ q , ρ τ q, in which, by inequality (3.4), DpSµ Z N pτ q , Sµ Σ ε N pτ q q ď λpτ q}ρ τ } 8 δ ε pτ, Z N q `cDpSµ Σ ε N pτ q , ρ τ q. We bound the last term by using inequality (3.3):

DpSµ Σ ε N pτ q , ρ τ q ď λpτ qp1 `}ρ τ } 8 q `bW 1 pSµ Σ ε N pτ q , Sf ε τ q `aW 1 pSf ε τ , ρ τ q ¯.
Again by inequality (3.3), we estimate the discrepancy in (4.7) with Dpρ τ , Sf ε τ q ď λpτ qp1 `}ρ τ } 8 q a W 1 pSf ε τ , ρ τ q. Note that, by (4.6), W 1 pSf ε τ , ρ τ q is of order ? ε.

Collecting these estimates together, and using Grönwall's lemma, we arrive at

δ ε pt, Z N q ď e KλpT q ş T 0 p1`}ρτ }8q dτ ˆ´4 ? ε `ż T 0 dτ ż ´bW 1 pSf ε τ , Sµ Σ ε N pτ q q `W1 pµ Σ ε N pτ q , f ε τ q ¯df bN 0 pZ N q ¯.
Using again the Fournier and Guillin bound in (3.1) we get the estimate (4.5), from which the thesis follows.

Appendix: proof of Proposition 1

We consider, for t P r0, T s, two weak solutions f 1 t and f 2 t of (2.4) in the weak sense specified by (2.7) and (2.8), with support on BRx`tRv ˆB Rv . We indicate by Z i pt, zq, i " 1, 2, the corresponding flows and we assume that min i"1,2 ş T 0 }Sf i τ } 8 dτ is attained for i " 1. We define the intermediate dynamics p f pt, zq, Zpt, zqq, where f ptq is the push-forward of f 2 0 by the flow Zpt, zq " pxpt, x, vq, ỹpt, x, vqq which solves

$ ' & ' % 9 xpt, x, vq " ṽpt, x, vq 9 ṽpt, x, vq " W rSf 1 t , ft spxpt, x, vq, ṽpt, x, vqq xp0, x, vq " x, ṽp0, x, vq " v.
From the hypothesis, for i " 1, 2, the spatial density Sf i t is bounded, and then, by (3.5) and (3.6), the field K pM rSf i ptqspx, |x ´y|qq is locally Lipschitz in x and y. Hence pZ i pt, zqq, i " 1, 2, exists for any t P r0, T s, and also the couple p f pt, zq, Zpt, zqq is well-defined.

As in previous sections, we indicate by λpτ q any constant that depends on Concerning the second term in the inequality (4.8), we have

d, R v , R x `τ R v . We have W 1 `f 1 t , f 2 t ˘ď W 1 ´f 1 t , ft ¯`W 1 ´f t , f 2 t ¯. ( 4
W 1 p ft , f 2 t q ď ż R 2d
|Z 2 pt, zq´Zpt, zq| df 2 0 pzq ď δptqsup

zPsupppf 2 0 q
|Z 2 pt, zq´Zpt, zq|.

To estimate δptq, we note that |Z 2 pt, zq´Zpt, zq| ď ż t 0 ´δpτ q`ˇˇW rSf 2 τ , f 2 τ spZ 2 pτ, zqq ´W rSf 1 τ , fτ sp Zpτ, zqq ˇˇ¯dτ.

Using Lemma 1, we estimate the last term in the integrand by 2R v KDpSf 2 τ , Sf 1 τ q `Kλpτ qp1 `}Sf 1 τ } 8 qpW 1 pf 2 τ , fτ q `δpτ qq. (4.10)

Hence we need only to estimate DpSf 2 τ , Sf 1 τ q. By the triangle inequality DpSf 2 τ , Sf 1 τ q ď DpSf 2 τ , S fτ q `DpS fτ , Sf 1 τ q. The second term can be bounded by using inequality (3.3), obtaining that DpS fτ , Sf 1 τ q ď c b }Sf 1 τ } 8 pR x `τ R v q d´1 W 1 pS fτ , Sf 1 τ q. (4.11)

We now establish an estimate for DpSf 2 τ , S fτ q of the type of inequality (3.4). For any z 0 , z P R 2d , by definition of δptq, |z 0 ´Z2 pt, zq| ´δptq ď |z 0 ´Zpt, zq| ď |z 0 ´Z2 pt, zq| `δptq.

Let X be the set of functions φ P Cpr0, `8q, Rq, with first derivative continuous up to a finite number of jumps. Given φ P X and δ ą 0, let φ δ prqφ `pr `δq ´φ´p r ´δq, where, denoting by φ the function φprqş r 0 |φ 1 psq| ds, we have defined Since Φ δ ´Φ P X, the first term is bounded by cDpS fτ , Sf 1 τ q, while the second can be easily bounded by cpR x `τ R v q d´1 δpτ q}Sf 1 τ } 8 . We conclude that DpSf 2 τ , S fτ q ď λpτ qδpτ q}Sf 1 τ } 8 `cDpS fτ , Sf 1 τ q. (4.12)

φ ˘prq - $ ' & ' % 1 2 p φprq ˘φprqq, if r ě 0, ˘1 2
Collecting estimates (4.10), (4.11) and (4.12), we arrive at δptq ď ż t 0 λpτ qKp1 `}Sf 1 τ } 8 q ´δpτ q `bW 1 pSf 1 τ , S fτ q ¯dτ.

Using the Grönwall's lemma and that x ď e x , for t P r0, T s, The thesis follows from the last inequality, together with (4.8) and (4.9).

1 ´rd 2

 12 ˇˇ(3.5) and, given x 1 , x 2 P R d and r ą 0,|M rρ 2 spx 1 , rq ´M rρ 2 spx 2 , rq| ď cd}ρ 2 } 8 r d´1 |x 1 ´x2 |. (3.6) By (3.5) and (3.6), the Lipschitz constants of the function K ´M rρ 2 spx, |x ´ξ|q ¯pv ´ηq in the variables x, ξ P Brx and v, η P Brv are bounded by cLippKq}ρ 2 } 8 r d´1 x r v and }K} 8 , respectively. Then, by the definition (2.14) of Wasserstein distance,

W 1 0 b W 1

 101 ´f t , f 2 t ¯ď δptq ď e λpT qK ş T 0 p1`}Sf 1 τ }8q dτ ˆż T pSf 1 τ , S fτ q dτ.

  Br pxq denotes the closed ball in R d of center x and radius r. We also indicate by BR the closed ball BR p0q. Note that M rρspx, rq ď 1 for any x P R d and r ě 0.

	for any α P C b pR d ˆRd q, where f 0 is the initial probability measure and
	pxpt, x, vq, vpt, x, vqq is the solution of the following Cauchy problem
	$ '	9 xpt, x, vq " vpt, x, vq
	&
		9 vpt, x, vq " W rSf t , f t spxpt, x, vq, vpt, x, vqq
	' %	xp0, x, vq " x, vp0, x, vq " v.
			|x ´y|qq pw ´vqf py, wq dy dw, (2.5)
	with		ż
		M rρspx, rq -	dρ.	(2.6)
			Brpxq
	Here and after, A weak formulation of equation (2.4) is given requiring that the solution f t
	fulfills the identity
	ż	ż
	αpx, vq df t px, vq "	α pxpt, x, vq, vpt, x, vqq df 0 px, vq	(2.7)

  .9)It is easy to verify that if Z N ptq " px 1 ptq, . . . , x N ptq, v 1 ptq, . . . , v

		N ptqq solves
	(2.1), (2.2), (2.3), then µ Z N ptq is a weak solution of (2.4). Namely, M rSµ Z N spx, rq
	is exactly M px, rq defined in equation (2.3) above. Thus we can rewrite the
	agent evolution for the Cucker-Smale model (2.1) with topological interactions
	(2.2), (2.3) as
	"	9 x i ptq " v i ptq

  Proposition 0. [3, Propositions 2.4 and 2.5] i) Let ρ and ν be two probability measures on R d with support in a ball BR and such that ρ P L 8 pR d q. Then Dpν, ρq ď c a R d´1 }ρ} 8 W 1 pν, ρq,

	x,rą0	ˇˇż	Brpxq	dµ	´ż Brpxq	dν ˇˇ,	(3.2)

for µ, ν two probability measures on R d . We will need the following results concerning it.

  vq ´W rSf 1 τ , fτ spx, vq| ď Kλpτ qp1 `}Sf 1 τ } 8 qW 1 pf 1 τ , fτ q, from which |Z 1 pt, zq ´Zpt, zq| ď |z ´z|ż Kλpτ qp1 `}Sf 1 τ } 8 q `W1 pf 1 τ , fτ q `|Z 1 pτ, zq ´Zpτ, zq| ˘dτ,

					.8)
	The first term is under control since the continuity estimate à la Dobrushin
	holds. Namely, from Lemma 1,
	|W rSf 1 τ , f 1 τ spx, t
		0
	which allows us to obtain
	W 1	´f 1 t , ft ¯ď e λpT qK	ş t 0 p1`}Sf 1 τ }8q dτ W 1	`f 1 0 , f 2 0 ˘, t P r0, T s.	(4.9)

  It can be proved (see[3, Lemma 2.2]) that Dpρ 1 , ρ 2 q " sup φPX: }φ} X ď1 sup x ż φ `|x ´y| ˘`dρ 1 pyq ´dρ 2 pyq ˘, Fixed x 0 P R d , we define Φpxq " φp|x ´x0 |q and Φ δ in the same way. Then, it is not difficult to see that (see [3, Lemma 2.3]) Φ `X2 pt, zq ˘ď Φ δ ´Xpt, zq ¯. Xpτ, zqq ´Φp Xpτ, zqq ¯df 2 0 pzq " ż pΦ δ pxq ´Φpxqq dS fτ pxq " ż pΦ δ ´Φq ´dS fτ ´dSf 1

	where }φ} X -	ş `8 0	|φ 1 prq| dr.
	Hence
	ż	Φ ´dSf 2 τ ´dS fτ	¯" ż	Φpxq ´df 2 pzq ´d f pzq "
		ż	´ΦpX 2 pτ, zqq ´Φp Xpτ, zqq ¯df 2 0 pzq
	ď	ż	´Φδ p

φp0q, if r ă 0. τ ¯`ż pΦ δ ´Φq dSf 1 τ .


F Masi

I Stefanou
email: ioannis.stefanou@enpc.fr

P Vannucci
email: paolo.vannucci@uvsq.fr

V Maffi-Berthier
email: victor.maffi-berthier@ingerop.com

A DISCRETE ELEMENT METHOD APPROACH FOR THE PRESERVATION OF THE ARCHITECTURAL HERITAGE AGAINST EXPLOSIONS

Keywords: Blast actions, Masonry structures, Discrete Element Method, Architectural heritage

Most emblematic monumental and historical buildings are masonry structures, made of brick and stone. Assess their response and vulnerability against explosive threats is a crucial issue in the preservation of such architectural assets and their visitors. The investigation of the blast loading response of masonry structures requires a proper modelling approach of the material and the loads arising from an explosion.

Here we propose a Discrete Element Method approach to investigate the fast-dynamic response of planar masonry structures subjected to blast loading. The Discrete Element Method captures the detailed mechanical and geometrical characteristics of masonry as well as the inherent coupling between the in-and out-of-plane motions. Blast loads are applied through an external dynamic library, which accounts for the effect of surface rotation of masonry components due to the pressure load (incident angle effect).

The discrete model is validated through existing experimental tests. The numerical simulations are found to predict with high fidelity the dynamic response and failure of the specimens, corroborating our modelling assumptions. The proposed approach is central to understand the behaviour of masonry structures subjected to blast loading and the influence of the various micro-mechanical material parameters.

INTRODUCTION

In the present international context, there is a need for the scientific community and institutions to better assess the threat of explosions meant to destroy civil engineering structures. Emblematic historical and monumental buildings, which have been targeted repeatedly, are often made of brick and stone, for which the prediction of the response to blast loading is particularly challenging. These structures tend to be massive, providing a natural resistance to blast pressure, but also typically have very low tensile strength between bricks or stones, leaving their joints able to open, and rendering them susceptible to collapse.

The investigation of the mechanical behaviour of masonry structures attracts significant scientific research, mostly focusing on the mechanical behaviour under quasi-static and seismic loads. In the last decade, increasing interest was shown in the response and potential damage of masonry structures due to explosive events, such as deliberate and/or accidental blasts and vented explosions from gas leaks, see e.g. [START_REF] Gagnet | Assessment of resistance definitions used for blast analysis of unreinforced masonry walls[END_REF][START_REF] Keys | Experimental analysis of debris distribution of masonry panels subjected to long duration blast loading[END_REF][START_REF] Li | Experimental and numerical study of unreinforced clay brick masonry walls subjected to vented gas explosions[END_REF][START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF][START_REF] Masi | Response of monumental buildings to internal explosions[END_REF][START_REF] Masi | Rocking response and overturning of museum artefacts due to blast loading[END_REF].

Attention was focused on the response of planar masonry structures to blast loads. We refer here to the seminal work of Gabrielsen et al. [START_REF] Gabrielsen | Response of arching walls and debris from interior walls caused by blast loading[END_REF] who first experimentally investigated the response to blast loading of full-scale unreinforced masonry walls with and without development of arching actions. Several experimental campaigns investigated the blast performance of one-way arching masonry walls [START_REF] Gagnet | Assessment of resistance definitions used for blast analysis of unreinforced masonry walls[END_REF][START_REF] Keys | Experimental analysis of debris distribution of masonry panels subjected to long duration blast loading[END_REF][START_REF] Li | Experimental and numerical study of unreinforced clay brick masonry walls subjected to vented gas explosions[END_REF][START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF][START_REF] Varma | Damage to brick masonry panel walls under high explosive detonations, Pressure vessels and piping division[END_REF][START_REF] Dennis | Response of 1/4-scale concrete masonry unit (CMU) walls to blast[END_REF][START_REF] Abou-Zeid | Response of arching unreinforced concrete masonry walls to blast loading[END_REF]. Keys and Clubley [START_REF] Keys | Experimental analysis of debris distribution of masonry panels subjected to long duration blast loading[END_REF] investigated the debris distribution and failure pattern of masonry walls subjected to long duration blast loading. More recently, Michaloudis and Gebbeken [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF] analysed the response of unreinforced masonry walls constrained to rigid supports subjected to far-field and contact explosions through experimental and numerical investigations.

Several modelling strategies were developed to investigate the dynamic response of masonry structures to explosions. We refer herein to both meso-and macro-scale numerical approaches, see e.g. [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF][START_REF] Abou-Zeid | Response of arching unreinforced concrete masonry walls to blast loading[END_REF][START_REF] Wang | Prediction of fragment size and ejection distance of masonry wall under blast load using homogenized masonry material properties[END_REF][START_REF] Wei | Model validation and parametric study on the blast response of unreinforced brick masonry walls[END_REF][START_REF] Macorini | Nonlinear analysis of unreinforced masonry walls under blast loading using mesoscale partitioned modeling[END_REF][START_REF] Rafsanjani | Implementation and validation of a strain rate dependent anisotropic continuum model for masonry[END_REF][START_REF] Parisi | Blast resistance of tuff stone masonry walls[END_REF][START_REF] Silva | Rigid block and spring homogenized model (HRBSM) for masonry subjected to impact and blast loading[END_REF].

The Discrete Element Method (DEM) is herein proposed as an accurate and effective modelling strategy for masonry structures subjected to explosions. Blast loads are computed using an external dynamic library which accounts for the effect of surface rotation of building blocks as well as the evolution in time of their relative distance with respect to the impinging blast wave. Through appropriate validation of the proposed approach with experimental data, we prove that DEM is a reliable tool in predicting the response of masonry structures and their failure to far-field explosions. In such scenarios, the mechanical behaviour of the masonry joints plays a leading role in the overall strength and failure of the structure. Failure within the body of masonry bricks does not occur or is usually negligible, as corroborated by experimental tests [START_REF] Gagnet | Assessment of resistance definitions used for blast analysis of unreinforced masonry walls[END_REF][START_REF] Keys | Experimental analysis of debris distribution of masonry panels subjected to long duration blast loading[END_REF][START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF] as well as by theoretical investigations, see e.g. [START_REF] Stefanou | Three dimensional homogenization of masonry structures with building blocks of finite strength: A closed form strength domain[END_REF] where it is shown that failure of the building blocks is important only under high compression.

We further investigate the accuracy of a rigid block model under quasi-static and dynamic conditions. Nondeformable blocks may provide inaccurate predictions for the out-of-plane response due to an intrinsic overestimation of the flexural stiffness.

The proposed modelling approach is central to investigate in detail the response and failure of masonry structures subjected to far-field explosions as well as to identify the governing micro-mechanical material parameters. It can be used to conceive appropriate and effective protective devises to help in the preservation of the architectural heritage. The discrete model can further be used to perform full-scale numerical tests focusing on the strength of the structure at varying of the geometry, material parameters, as well as the intensity of the blast loading. Ultimately, our numerical model, relying on a detailed description of the system, offers the possibility to conceive in-scale experiments through the definition of appropriate scaling laws. Notice that performing blast experiments either in reduced or in full scale presents many difficulties [START_REF] Draganić | An overview of the methods for blast load testing and devices for pressure measurement[END_REF], due to the nature of the loading action (e.g. need of specialized personnel, risks, measuring devices, repeatability, uncertainties etc.)

The paper is structured as follows. In Section 2 we present the Discrete Element model and the main assumptions of our modelling approach. We briefly describe in Section 3 the effects of explosions on targets and the model used to apply the loading on the masonry structure. We perform in Section 4 numerical tests of the response of a masonry wall subjected to far-field explosions, assuming (elastic) deformable masonry blocks. The results are compared and validated with existing experimental tests. In Section 5 we investigate the validity and accuracy of a rigid block model to predict the out-of-plane response of a masonry wall.

DISCRETE ELEMENT METHOD AND MODELLING ASSUMPTIONS

The Discrete Element Method is used to investigate the dynamics of masonry structures subjected to blast loading. The model allows to directly access the various micro-mechanical parameters, such as the constitutive behaviour of the interfaces and bricks, as well as the geometry of the building blocks. Furthermore, a discrete approach allows to simulate the progressive failure of masonry and capture with increased fidelity the post-peak, softening, dynamic behaviour of a masonry structure with bricks undergoing large displacements and rotations.

The DEM simulations are herein carried out by means of the 3DEC software [19]. A central finite differences scheme is used in this code for integrating in time the equations of motion of each block. This code provides the means to apply the conceptual model of a masonry structure as a system of blocks, either rigid or deformable. Block deformability is considered through internal discretization of blocks into finite differences. Herein deformable blocks are used to achieve a detailed description of the structural response. The validity of a rigid blocks modelling assumption is further investigated in Section 5. A soft-contact algorithm is used to model the interactions between neighbouring blocks through interfaces/joints.

Elastic behaviour of the interfaces

A Discrete Element (DE) model with deformable blocks consists of deformable units which are subdivided into a finite-difference mesh of triangular elements, interacting along zero-thickness joints.

The elastic behaviour of the interfaces is defined through the following stress-displacement relationship [19]

� t n t t � = � k n 0 0 k t � � u n u t � or t = k u , (1)
where t and u are the vectors collecting the normal, t n , and tangential, t t , forces per joint's unit area and the joint normal, u n , and tangential, u t , displacement, respectively. The elastic stiffness matrix k collects the normal and tangential to the interface stiffnesses: k n and k t , respectively. Normal and tangential stiffnesses can be computed from the properties of the masonry components (bricks and joints), considering the expansion in dimension of the masonry blocks and the soft-contact assumption. The identification of the numerical values of the elastic parameters is obtained assuming a stack bond masonry and uniform stress distributions both in the block and mortar, see Figure 1. The joint normal stiffness k n is computed by requiring the analogy between a real masonry prism and the discrete element model, in terms of series of axial springs. The equivalent stiffness for the masonry prism in Figure 1 is

k n,eq -1 = � E b h b � -1 + � E m h m � -1 → k n,eq = E b E m E m h b + E b h m (2)
with E b , E m the Young's moduli of the brick and the joint, respectively, and h b , h m the height of the block and the joint (masonry prism).

Figure 1. Scheme to determine joint normal stiffness. The value of k n is found by requiring the equivalence between the axial stiffness of the real masonry prism (left) and the one of the DEM with expanded blocks and zerothickness interfaces (right).

The Young's and shear moduli of deformable blocks, E * b and G * b , are modified to account for the expansion in dimensions (with respect to the real dimensions of the masonry bricks)

E * b = E b �1 + h m h b � , G * b = G b �1 + h m h b � , (3)
where G b is the shear modulus of the masonry bricks. Equating expression (2) with the equivalent stiffness of the DE model along the normal direction, i.e.,

k n,DE = (h b +h m ) E * b ⁄
, we recover the joint normal stiffness k n = E m h m ⁄ , and similarly, along the tangential direction, the tangential stiffness is computed as k t = G m h m ⁄ .

Plastic and softening behaviour of the interfaces

In the absence of more detailed experimental data regarding the constitutive behaviour of masonry joints the Coulomb criterion seems to be a reasonable choice. Several experimental observations (e.g. [START_REF] Van Der Pluijm | Out-of-Plane bending of Masonry[END_REF][START_REF] Lourenço | Characterization of cyclic behavior of dry masonry joints[END_REF][START_REF] Vélez | Dry stone masonry walls in bending-part I: Static tests[END_REF]) justify its use up to moderate compression. Herein we assume an infinite compressive strength of the interfaces. Such hypothesis is verified a posteriori by monitoring the compressive stresses during each computation. The maximum shear (or tangential) force per joint's unit area t t (in N/m²) is limited by the Coulomb failure surface:

f 1 = t t -c -t n tan ' ≤ 0 ( 4
)
where c is the cohesion of the interface and ' the friction angle. Compression is here considered negative. In shear/tensional regime a tension cut-off is often used as shown in Figure 2. In other words, the maximum normal force per joint's unit area t n (in N/m²) is limited by the tensile strength according to:

f 2 = t n -f t ≤ 0 (5)
where f t is the tensile strength of the interface. The joint forces F t and F n are related to the corresponding displacements by Hooke's law

F t = A j t t , F n = A j t n (6)
where A j is the joint's area. The two inequalities (4, 5) define the elastic domain of a masonry interface. The evolution and contraction under combined shear and normal plastic deformation of the elastic domain can be considered in terms of the various micro-mechanisms related to progressive softening of the joint. According to experimental results on interfaces, a softening behaviour is observed as depicted in Figure 2. In this way cohesion, maximum tensile strength, and friction and dilatancy angles can evolve from their initial values c ; f t ; ' ; ª to some smaller residual values c res ; f t,res ; ' res ; ª res . All these values can be determined by typical experimental tests on interfaces.

Figure 2. Initial and residual strength surfaces used for modelling joints behaviour.

Regarding the plastic flow rule, this is given by the following potentials:

g 1 = t t -c -t n tan ª (7)
g 2 = t n -f t ( 8
)
where ª is the dilatancy angle. If ª = ' we say that the plastic flow rule is associative (normality condition), otherwise (ª < ') the material obeys a non-associative plastic flow rule. In both cases, the following general relation between the rate of change of forces and the rate of change of total displacements stands

t = k p u , (9)
with k p the plasticity matrix,

k p = k n •+ tan ' tan ª � tan ª 1 tan ' tan ª tan ' � , (10)
where

• = k n k t , (11)
t and u the vectors collecting the rate of change of the forces and of total displacement.

Constitutive behaviour of deformable blocks

An isotropic, elastic material behaviour is herein considered for the blocks. Infinite tensile and compressive strength are further assumed. This hypothesis may be strong in the case of relatively near-field explosions and especially in contact detonations, see e.g. [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]. Nevertheless, damage within the body of masonry bricks is generally negligible in far-field explosions, where the collapse of the masonry structure is governed by failure at the interfaces [START_REF] Gagnet | Assessment of resistance definitions used for blast analysis of unreinforced masonry walls[END_REF][START_REF] Keys | Experimental analysis of debris distribution of masonry panels subjected to long duration blast loading[END_REF][START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]. This is also corroborated by the fact that the tensile as well as the compressive strength are strain rate dependent. Extensive researches showed that at increasing loading rates the resistance of geomaterials, such as mortar, concrete, and marble, increases as well, mainly due to the finite growth rate of micro-cracks [START_REF] Freund | Crack propagation in an elastic solid subjected to general loading-i. constant rate of extension[END_REF][START_REF] Freund | Crack propagation in an elastic solid subjected to general loading-ii. non-uniform rate of extension[END_REF] and the viscosity of the material [START_REF] Weerheijm | Concrete under impact tensile loading and lateral compression[END_REF]. At increasing strain rates, an increase of the tensile strength among other parameters, is observed. The dynamic increase factor for tensile strength for geomaterials (such as mortar, tuff, granite, etc.) usually varies between 1 and approximately 7 in function of the involved strain rates, see e.g. [START_REF] Chen | Experimental study on dynamic tensile strength of cement mortar using split Hopkinson pressure bar technique[END_REF][START_REF] Ross | Effects of strain rate on concrete strength[END_REF].

In each computation, stresses as well as strain rates are monitored to verified that the related dynamic strength of the material is such that failure within blocks does not occur.

BLAST LOADS AND MODELLING

Explosion produces a blast wave of high-pressure accompanying high-temperature and supersonic expansion of gases. The abrupt increase of the pressure carried by a blast wave can produce severe structural damage. When the primary shock meets a target, it generates on it the so-called reflected overpressure, P r , which is the difference between the pressure determined by the explosion increased by the reflection at target's surface and the ambient one, P o . Figure 3 shows the schematic time variation of P r , which is determined by the arrival time of the shock wave, t A , the overpressure peak, P ro , the positive phase duration, t o , negative phase duration, t o-, and the underpressure peak, P ro-. These parameters are functions of the distance R and the explosive weight (conventionally expressed in TNT equivalent).

The pressure acting on a target due to blast loading is the algebraic sum of the hydrostatic overpressure and the dynamic pressure C D q ∶= ½u|u|/2, with C D the drag coefficient (function of the target shape and Mach and Reynold numbers), ½ the density, and u the velocity of gas particles.

The simulation of a blast can be conducted by using different approaches [START_REF] Remennikov | A review of methods for predicting bomb blast effects on buildings[END_REF][START_REF] Larcher | Explosions in Complex Geometries -A Comparison of Several Approaches[END_REF][START_REF] Masi | A study on the effects of an explosion in the Pantheon of Rome[END_REF]. Herein we refer to empirical models based on experimental results available in the existing literature. Figure 3. Time evolution of overpressure (i.e. the pressure measured relatively to the atmospheric one) due to an explosion acting on a target.

Blast model

We model blast actions following the work of Hyde [START_REF] Hyde | Conwep: Conventional weapons effects program[END_REF] with the empirical model ConWep, which relies on the best-ft interpolations of the experimental results from Kingery and Bulmash [START_REF] Kingery | Air blast parameters from TNT spherical air burst and hemispherical burst[END_REF]. The interpolations allow to determine the blast parameters and pressure loading from the knowledge of the trinitrotoluene (TNT) equivalent explosive weight, W, and the stand-off distance, R. The time evolution of the reflected pressure is modelled with the well-established modified Friedlander equation,

P r (t) = P ro �1 - t t o � exp �-d t t o � , ( 12
)
where d is the exponential decay coefficient. The impulse associated to the positive, i r+ , and the negative phase, i r-, reads

i r+ = ∫ P r (t) dt = �e -d + d -1� t o 0 P ro t o d 2 , ( 13
) i r-= ∫ P r (t) dt = ∞ t o P ro t o d 2 e -d , (14)
respectively. Equation 16allows to determine the exponential decay coefficient, d, by equating it with the best-fit interpolation of i r+ from experiments [START_REF] Kingery | Air blast parameters from TNT spherical air burst and hemispherical burst[END_REF].

COMPARISON WITH EXPERIMENTAL TESTS

The proposed numerical model is herein compared and validated with existing experimental tests. Among the experiments available in the literature, we select one of the most well-documented [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]. Notice that performing blast experiments either in reduced or in full scale presents many difficulties, due to the nature of the loading action which may result in large uncertainties of the obtained results.

Michaloudis and Gebbeken [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF] analyse the response of unreinforced masonry walls subjected to far-field and contact explosions through experimental and numerical investigations. Among four tests, two involve masonry walls subjected to the explosion of W 1 = 810 kg and W 2 = 1150 kg of TNT at R = 37 m from the targets. Due to the large stand-off distance, the blast wave impinges almost uniformly and instantaneously the target. However, no information is given concerning the evolution of the blast pressure in the experimental tests.

The brickwork consists of a running bond with bricks of nominal dimensions a × b × w = 80 × 240 × 120 mm, see Figure 4. The boundaries of the walls are constrained, through mortar interfaces, to stiff fixed supports.

In Test 1 (W 1 = 810 kg, R = 37 m), the observed maximum outward and inward deflection at the centre of the wall are 77 mm and 37 mm, respectively. In Test 2 (W 2 = 1150 kg, R = 37 m), a breach at the centre of the wall originates mostly only due to the damage of the joints. Failure within the body of individual bricks is not observed or is negligible [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]. The maximum dimensions of the breach are equal to 4 bricks along the length of the wall and to 13 bricks along the height. In the numerical DE model, a constant thickness of the mortar h m = 10 mm is assumed in the lack of more detailed information for the walls. Table 1 presents the numerical material parameters which have been selected from the literature, see e.g. [START_REF] Stefanou | Three dimensional homogenization of masonry structures with building blocks of finite strength: A closed form strength domain[END_REF][START_REF] Kingery | Air blast parameters from TNT spherical air burst and hemispherical burst[END_REF][START_REF] Wild | Pozzolanic properties of a variety of European clay bricks[END_REF] Superscripts h and b refers to head and bed joints, respectively. Subscript res identifies the residual (softening) value of the parameters, when differ from the initial ones (see Fig. 2). In the lack of information about the boundary conditions acting on the structure, the wall is considered simply supported at its four extremities, as complete rotational restraint or freedom may be difficult to achieve in practice. Blast loads are computed and applied with the dynamic library presented in Section 3.

Convergence analyses for contact and finite difference discretization are performed with the same criterion used in [START_REF] Petry | Cyclic test data of six unreinforced masonry walls with different boundary condition[END_REF]. Two different loads are applied in two stages. First, a gravity load is applied, and the system is left to reach equilibrium. Second, increasing out-of-plane forces, proportional to the block mass, are applied to each block, up to wall failure. Due to the inherent, dynamic character of most DEM codes, including 3DEC, care was taken in order to assure quasi-static conditions until failure. The loads are applied adequately slow in order to assure a small ratio of the kinetic energy over the total energy of the system. As far it concerns the specific DEM code used herein, 3DEC, the number of contact points situated across the wall thickness has a direct impact on the model response. This is especially important for out-of-plane loads as the one resulting from an explosion. From the convergence analyses, we find that at least 10 contact points are required to accurately modelling the out-of-plane deflection of the wall under study. The finite difference mesh for deformable blocks consists of tetrahedra with average characteristic length equal to 30 mm. As follows, the numerical simulations consist of two stages. A gravity load is first applied assuming a linear elastic behaviour, then, once the system has reached the equilibrium, blast loads are applied. In order to model the fast-dynamic response of the structure, no mechanical damping formulation is used. Table 2 shows the results obtained with the DEM model for Test 1 (W 1 = 810 kg, R = 37 m) and compare them with the test data and the numerical computations presented in [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF], which rely on similar modelling assumptions. The agreement of the results of the proposed model with the observed deflections is remarkably good (also compared to the numerical results in [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]). The maximum outward and inward displacements deviate from the experimentally observed one by 1.5 % and 3.2 %, respectively. We present in Figure 5 the evolution in time of the numerically measured deflection at the centre of the wall. In the free-oscillating response, the system gradually dissipate energy as a result of the slippage along interfaces, until it reaches an equilibrium. The centre of the wall displays a permanent deflection of approximately 7.1 mm along the outward direction.

Table 1: Material parameters used in the numerical simulations with deformable DE.

The assumption of infinite strength of the blocks in compression/tension and of the interfaces in compression is verified a posteriori. In both cases (blocks and interfaces), stresses are found to be 60 % lower than the corresponding material strength.

Table 2: Comparison of the maximum outward and inward deflection at the centre of the wall for Test 1 (W 1 = 810 kg, R = 37 m) between the experimental values and the numerical predictions using deformable blocks (DEM) and the values obtained in [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF].

Test 2

We present in Figure 6 the response of the wall from the numerical DEM simulations of Test 2 (W 2 = 1150 kg, R = 37 m). The out-of-plane dynamic response and the consequent formation of the breach can be clearly observed. 3 compares the dimensions of the breach from the experimental evidence and the numerical results. In Figure 7, we compare the final state of the masonry wall with a representative scheme obtained from the test results, see [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]. The DE model is found to well captures the form of failure, as well as the extension in height of the breach. However, less bricks are involved along the width direction. The number of bricks that travel from the wall is 22, while in the experiment is approximately equal to 40. The slight difference may be due to the complex fluidstructure interaction phenomena that take place, in the reality, as the breach originates. The blast pressure, acting continuously on the openings, may be the cause of a larger breach in the width, i.e., along the direction where the masonry brick are more likely to rotate. Furthermore, it is worth stressing that, in the test, some bricks involved in the breach, even if few in number, undergo damage, which is not considered in our model.

As for test 1, the assumptions of infinite compression strength of blocks and infnite tensile strength of blocks and interfaces are a posteriori verifed. The compression stresses are both in the masonry bricks and the interfaces below the material strength [START_REF] Ross | Effects of strain rate on concrete strength[END_REF]. Tensile stresses exceed in some bricks th the material strength. Nevertheless, the enhancement in tensile strength due to high loading rates requires to be accounted for [START_REF] Chen | Experimental study on dynamic tensile strength of cement mortar using split Hopkinson pressure bar technique[END_REF][START_REF] Ross | Effects of strain rate on concrete strength[END_REF]. The numerically measured strain rates within discrete elements are approximately 5-10 s -1 , which correspond to a dynamic tensile strength of approximately 4.5 -6 MPa [START_REF] Chen | Experimental study on dynamic tensile strength of cement mortar using split Hopkinson pressure bar technique[END_REF][START_REF] Ross | Effects of strain rate on concrete strength[END_REF]. Only some bricks involved in the formation of the breach are subjected to tensile stresses as high as 5 MPa. The proposed DE model shows a good agreement with the tests both for the dynamic response (Test 1) and the form of the failure (Test 2). This is a key aspect that is necessary to consider in every numerical strategy and approach which may be used to investigate the response and failure modes of masonry structures, and in particularly when the effects of blast loading are considered. A correct evaluation of the formation of breaches and weak planes may allow, through further investigations, the design of effective protections.

Effects of the (non-)associativity of the sliding behaviour of masonry joints

In the above calculations, a non-dilatant sliding behaviour was assumed for the interfaces, and namely à = 0°, as usually displayed by masonry joints, see e.g. [START_REF] Lourenço | Characterization of cyclic behavior of dry masonry joints[END_REF]. We investigate herein the effect of an associative behaviour of the interfaces, i.e. when à = ' = 35°, on the response of the wall for both tests (1 and 2). All material parameters, except dilatancy, are kept the same as in Table 1. Figure 8 displays the time evolution of the deflection at the centre of the wall for dilatant and non-dilatant joints. The associative case shows reduced out-of-plane response and increased strength, see also [START_REF] Godio | Effects of the dilatancy of joints and of the size of the building blocks on the mechanical behavior of masonry structures[END_REF]. In Test 1, the maximum deflection is 22.9 mm, namely 3.4 times smaller the one obtained with à = 0° (i.e. 78.2 mm). The maximum and minimum principal stress within masonry bricks are respectively 4 and 15 times larger than the values obtained with zero dilatancy. In Test 2, associativity of the interfaces results in a structure able to withstand the blast loading of 1150 kg of TNT. In this case, the maximum deflection is 80.5 mm, which is comparable to the results obtained with a non-dilatant behaviour, but under an explosive charge 340 kg smaller, i.e. 810 kg (Test 1).

IS A RIGID BLOCKS ASSUMPTION APPROPRIATE TO MODEL THE OUT-OF-PLANE RESPONSE?

Once the DE model has been validated with existing experimental results involving explosions, we investigate the simplified modelling assumption of infinitely rigid blocks which is often preferred in the literature, see e.g. [START_REF] Godio | Effects of the dilatancy of joints and of the size of the building blocks on the mechanical behavior of masonry structures[END_REF][START_REF] Bui | Masonry walls under membrane or bending loading cases : experiments and discrete element analysis[END_REF][START_REF] Çakt | Discrete element modeling of a scaled masonry structure and its validation[END_REF], for its reduced computational cost with respect to the more detailed model with deformable blocks. A rigid blocks assumption is usually reasonable under relatively low compressive loads, where the deformation is principally concentrated at the interfaces [START_REF] Stefanou | Three dimensional homogenization of masonry structures with building blocks of finite strength: A closed form strength domain[END_REF]. Moreover, from an energetic approach, the assumption of rigid blocks is supposed to provide safe estimations (lower bound) of the overall strength. Nevertheless, masonry structures subjected to out-of-plane loading, as it is the case with the experiments herein considered [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF] and as frequently occurs with explosive loads, the rigid blocks assumption can give unsatisfactory results, see e.g. [START_REF] Parisi | Blast resistance of tuff stone masonry walls[END_REF].

Notice that several contact points through the thickness of the masonry structure undergoing out-of-plane displacement are required both with rigid and deformable blocks formulations, for more we refer to [START_REF] Godio | Effects of the dilatancy of joints and of the size of the building blocks on the mechanical behavior of masonry structures[END_REF][START_REF] Lemos | Numerical issues in the representation of masonry structural dynamics with discrete elements[END_REF][START_REF] Lemos | Contact representation in rigid block models of masonry[END_REF]. However, in a rigid block model, the stress distribution at the interfaces follows a linear diagram, as the faces remain planar, thus an accurate discretization of contacts is crucial. For instance, it has been proved that in the frame of the DEM code herein used, 3DEC, at least 3 contact points along the thickness are required to obtain an accurate representation of the bending stiffness [START_REF] Lemos | Contact representation in rigid block models of masonry[END_REF].

With rigid blocks, the normal and tangential stiffnesses of the interfaces are appropriately modified with respect to the expressions previously derived. In this case, the stiffness of the entire system is reduced to the stiffness of the zero-thickness interfaces, k n , and namely

k n = E b E m E m h b + E b h m (15)
A similar expression can be derived for joint tangential stiffness k t :

k t = G b G m G m h b + G b h m (16)
with G b and G m the shear moduli of the masonry bricks and the mortar, respectively. The material parameters that define the plasticity behaviour and the corresponding softening remained unchanged, see Sect. 2. We present in Table 4 the material parameters used in the numerical DE model, derived from those given in Table 1 (for deformable DE) and considering the non-deformability of the blocks, see Eq.s [START_REF] Parisi | Blast resistance of tuff stone masonry walls[END_REF] and [START_REF] Silva | Rigid block and spring homogenized model (HRBSM) for masonry subjected to impact and blast loading[END_REF]. The discretization of the contact interfaces is investigated as presented in Section 4. From the convergence analyses, a discretization consisting of 20 contact points along the wall thickness is selected. Boundary conditions are kept the same as for the deformable DE.

As follows, we investigate the effects of the non-deformability of blocks in the numerical DE model. This is first accomplished in quasi-static conditions, applying a constant out-of-plane pressure to a localized region of the structure. Second, the dynamic response to blast loading is studied by numerical simulations of tests 1 and 2. Table 4: Material parameters used in the numerical simulations with rigid DE. Superscripts b and h refer to bed and head joints, respectively. The joints stiffnesses are differentiated along the two directions in order to consider the different height of masonry bricks, see Eq.s [START_REF] Parisi | Blast resistance of tuff stone masonry walls[END_REF] and [START_REF] Silva | Rigid block and spring homogenized model (HRBSM) for masonry subjected to impact and blast loading[END_REF]. f res (kPa) 0

Response under quasi-static loading

A constant uniform pressure, P = 100 kPa, is applied to a centred region of the wall, 7 courses high and 3 1 2 wide. We assume a linear elastic behaviour of the interfaces. We present in Figure 8 the deformed shaped obtained with rigid DE and deformable ones. Table 5 compares the strain energies and maximum deflection obtained with the two models. Rigid blocks are found to underestimate the overall out-of-plane deflection with a maximum displacement equal to 1.55 mm with respect to the one predicted with deformable blocks, i.e., 2.26 mm. The reason lies on an intrinsic rotational locking (term we select in analogy with the well-known shear locking phenomenon observed in Finite Element Method formulations). The rotational locking happens independently from the fineness of the contact discretization, as the rigid block model consists of 20 contact points along the wall thickness. This is clearly shown in Figure 9 by the different slopes at the boundaries (x = ±1.2 m and y = ±1.24 m) for rigid and deformable DE.

Referring to Table 5, we notice that the overall (material) strain energy, blocks plus interfaces, is almost equal for both models, however, the contributions of the joints strain energy highly differ. In particular, the compressive and tensile strain energies with rigid DE are found to be approximately two times larger than with deformable blocks. This is a direct consequence of the inadequacy of a rigid body behaviour for the case where relatively high compressive forces develop. In the flexural response of the wall, the (deformable) blocks can deform significantly, phenomenon that is not capture by a rigid body model. Moreover, it should be noticed that in the specific DEM code used herein, the modelling assumption of interfaces represented by only two elastic parameters, namely normal and tangential stiffnesses, may cause the enhancement of the observed rotational locking with rigid blocks. Figure 9. Comparison between rigid and deformable block models in terms of deflection along vertical (y axis) direction (left) and horizontal (x axis) direction (center). The curves are obtained by interpolation of the out-ofplane displacement of the centroid of blocks (marks). The constant and uniform pressure load is applied to the region highlighted in gray (right).

Table 5: Comparison of the results from the models with rigid and deformable blocks under quasi-static conditions. Material strain energy refers to the algebraic sum of the joints and blocks strain energies.

Response under fast-dynamic excitations

The dynamic response of the same model used in Section 4 and derived from the experimental tests in [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF] is herein investigated using the hypothesis of infinitely rigid blocks.

We compare in Table 5 the numerical predictions for Test 1 (W 1 = 810 kg, R = 37 m) using rigid DE with the deformable blocks model. The rigid DE model is found to underestimate the outward deflection by 37.1 % the experimentally measured value. A rather good agreement is found instead for the inward deflection. We present in Figure 10 the time evolution of the out-of-plane displacement at the centre of the structure as computed with both DE models. It is worth noticing the difference of the mean value of the free-oscillating In the case of Test 2 (W 2 = 1150 kg, R = 37 m), the rigid body model overall predicts the dimensions of the breach with good agreement with respect to the test results, see Figure 11. In particular, the width (4 courses of bricks) is accurately predicted by the numerical simulations, while the extension in height is overestimated, see Table 7. In this case, still due to the observed rotational locking, the rigid block model is found to overestimate the breach dimension, see Fig. 11.

CONCLUDING REMARKS

We investigated the dynamic response and failure of masonry structures subjected to blast loads. The Discrete Element Method (DEM) was herein proposed as an accurate and effective modelling strategy as it captures the progressive failure and response of a masonry structure with bricks undergoing large displacements and rotations.

The proposed model relies on simplifying assumptions as it considers elastic deformable blocks and the non-linear material response as being limited to the mortar interfaces. Indeed, the mechanical behaviour of the joints plays a leading role in the overall strength and failure of the masonry structure under far-field explosions, while failure within the body of masonry bricks does not occur or is usually negligible [START_REF] Gagnet | Assessment of resistance definitions used for blast analysis of unreinforced masonry walls[END_REF][START_REF] Keys | Experimental analysis of debris distribution of masonry panels subjected to long duration blast loading[END_REF][START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF][START_REF] Stefanou | Three dimensional homogenization of masonry structures with building blocks of finite strength: A closed form strength domain[END_REF].

The accuracy of the DEM model and corresponding modelling assumptions were investigated and validated through comparisons with existing, well-documented experimental tests [START_REF] Michaloudis | Modeling masonry walls under far-field and contact detonations[END_REF]. We proved that our numerical model can predict with high fidelity the dynamic response and failure under far-field explosions of masonry structures.

Once the numerical model was validated, it was used to assess the accuracy of a rigid block model in predicting the out-of-plane response under quasi-static and fast-dynamic conditions. The hypothesis of infinitely rigid blocks was found to overestimate the flexural stiffness of the system. As already discussed in [START_REF] Godio | Effects of the dilatancy of joints and of the size of the building blocks on the mechanical behavior of masonry structures[END_REF][START_REF] Lemos | Numerical issues in the representation of masonry structural dynamics with discrete elements[END_REF][START_REF] Lemos | Contact representation in rigid block models of masonry[END_REF], the contact discretization of the (rigid) Discrete Element (DE) model is a critical issue and has a direct impact on the model response. Nevertheless, the additional investigations we presented showed that, even in the case of a fine discretization of contacts, a rigid DE model displays an intrinsic rotational locking between adjacent blocks. This may result in a rather significant underestimation of the out-of-plane displacement of the structure, which was found to be around 70 % of the actual one (from experimental tests and a deformable block model) for the cases herein considered.

It is worth noticing however that a rigid block model has a reduced computational cost with respect to the more detailed description offered by deformable DE. In the simulations herein performed, the computing time for rigid blocks was approximately 50-60 % faster than with deformable ones.

Differently from the static and dynamic out-of-plane response, we found that the rigid block model gives sufficiently accurate predictions for the failure (strength) of the masonry structure subjected to blast loads.

The results obtained in this paper can be useful for improving the understanding on the behaviour of masonry structures under out-of-plane deformations and in particularly under blast loads. Moreover, the proposed model can be used to perform reliable numerical tests of those systems for which the experimental set-up may be arduous and complex due to the fluid-structure interaction blast phenomena as well as the non-standard geometry of studied masonry structures [START_REF] Masi | Masonry vaults under explosive loads[END_REF], e.g. arches, vaults, domes, etc. This is a central point in the preservation of our architectural heritage against explosions as the numerical predictions of the failure modes can be used to design appropriate and effective protective devices. The proposed Discrete Element Method approach can be also a useful tool in investigating the response and strength of large masonry structures to blast loading with numerical methods based on continuum mechanics (upscaling/homogenization).

Breach

Figure 4 .

 4 Figure 4. Geometrical model for test 1 and 2 in [4]. The masonry wall has thickness w = 120 mm and coursing consisting of a running brick bond. Bricks have nominal size a × b × w = 80 × 240 × 120 mm.

Figure 5 .

 5 Figure 5. Evolution of the deflection at the centre of the wall for Test 1 (W 1 = 810 kg, R = 37 m) computed by the proposed DE model and from the computations in [4].

Figure 6 .Figure 7 .

 67 Figure 6. Response of the DE model for Test 2 (W 2 = 1150 kg, R = 37 m), a-e), and final state of the masonry wall, f).

Figure 8 .

 8 Figure 8. Influence of the dilatant behaviour of masonry joints in terms of the deflection at the centre of the wall for Test 1 (W 1 = 810 kg, R = 37 m), left, and Test 2 (W 2 = 1150 kg, R = 37 m), right.

Figure 10 .

 10 Figure 10. Evolution of the deflection at the center of the wall obtained from the numerical simulations with rigid and deformable DE for Test 1 (W 1 = 810 kg, R = 37 m).

Figure 11 .

 11 Figure 11. Comparison between the numerical results of the rigid and deformable block models for Test 2 (W 2 = 1150 kg, R = 37 m). The experimental breach extension is schematically represented by the dashed black line.

Table 6 :

 6 Deformable blocks display a permanent deflection of approximately 7.1 mm in the outward direction. The rigid body model, instead, oscillates around a mean value equal to zero. The reason lies again on the rotational locking that a rigid DE description presents. The underestimated out-of-plane response results in a lesser amount of plastic energy dissipation at the interfaces. Comparison of the maximum outward and inward deflection at the centre of the wall for Test 1 (W 1 = 810 kg, R = 37 m) between the rigid and deformable block models and the experimental test.

	y
	x
	P = 100 kPa

Table 3 :

 3 Dimensions of the breach for Test 2 (W 2 = 1150 kg, R = 37 m). Comparison between the observed values and the numerical predictions using rigid and deformable blocks (DE).

	Maximum deflection	Experimental	Numerical predictions	Relative error
			rigid DE deformable DE	rigid DE deformable DE
	Outward (mm)	77.0	48.4	78.2	37.1 %	1.5 %
	Inward (mm)	37.0	39.2	38.2	5.6 %	3.2 %