# NSCD discrete element method for modelling masonry structures

Brahim Chetouane, Frédéric Dubois, Marc Vinches, Claude Bohatier

► To cite this version:

Brahim Chetouane, Frédéric Dubois, Marc Vinches, Claude Bohatier. NSCD discrete element method for modelling masonry structures. International Journal for Numerical Methods in Engineering, 2005, 64, pp.65-94. 10.1002/nme.1358. hal-01476972

HAL Id: hal-01476972

https://hal.science/hal-01476972

Submitted on 26 Feb 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
NSCD discrete element method for modelling
masonry structures

B. Chetouane¹, F. Dubois², M. Vinches¹ and C. Bohatier²

¹Alès School of Mines, CMGD, 6 Avenue de Clavières, 30319 Alès, France

²University of Montpellier II LMGC UMR 5508, Place Eugène Bataillon, 34095 Montpellier, Cedex 05, cc048, France

The aim of this paper is to present a discrete element approach for the study of stonework. In the present work, a masonry structure is considered as a collection of rigid or deformable blocks, interacting together by contact laws. In this paper, we use the non-smooth contact dynamics (NSCD) resolution method mainly used for the modelling of granular media. In the considered masonry structures we define, on an elementary cell, average local strain and stress tensors. These definitions are valid under dynamical loading of the structure, taking into account rotations. We present their use on academic and on real masonry structures.

KEY WORDS: discrete element method; NSCD; strain tensor; stress tensor; masonry

# 1. INTRODUCTION

The growing concern about the preservation of the built heritage leads us to study in more detail the mechanical behaviour of masonry structures. In most cases, the masonry is an heterogeneous material composed of blocks, with various shapes and sizes, separated by joints. Numerous models have been proposed in the literature in order to study the mechanical behaviour of masonry and particular interest in this context has been dedicated to the micro-mechanical models that enable investigations on the behaviour of masonry starting from the knowledge of its single constituents. Each component of the masonry structure can have its own behaviour, which can be complex. Discrete element methods are widely used to model granular materials [1] (sand, rocks, grains in silos, etc.), condensed media (colloid), biological media (bone cells) and some damaged or micro-fractured continuous media. They can also be used to model mechanisms, tensegrity structures or robotics systems.
In the following section, we present all the components of a particular discrete element method based on non-smooth contact dynamics (NSCD): space discretization of the bodies, dynamical behaviour, time integration of bodies movement, contact detection, interaction laws (impact laws, behaviour laws, etc.) and resolution methods. Then, we propose a definition of the stress and strain tensors in granular materials, in order to obtain an average information on this granular medium.

We have done a comparison of the results, for a wall, of an entirely rigid modelling, an entirely deformable one and a finite element continuous one. For the discrete models, the calculations were performed using the LMGC90 code, with the NSCD solver.

## 2. NSCD DISCRETE ELEMENT METHOD

### 2.1. *Dynamic equations in presence of contact*

Dynamic equations governing a body kinematics in the presence of contact (but without impact) could be written as follows:

$$ M\dot{u} = F(q, u, t) + P(t) + \mathbb{R} \qquad (1) $$

boundary conditions + initial conditions

where $q$ indicates a Lagrange variable (for example translations or rotations of the gravity centre of a rigid body or node mesh displacements of a deformable body), $u$ the velocity function, $\mathbb{R}$ the contact efforts contribution to this equation (torque or nodal reaction) and $F(q, u, t) + P(t)$ represents the sum of internal and external forces. $F(q, u, t)$ comprises certain standard terms and also centrifugal and gyroscopic terms, $M$ is the inertia matrix made up of mass and inertia components. For a deformable body, the dynamic equation can be obtained by a finite element discretization of the continuous problem. For a rigid body, this equation means the conservation equations of momentum. Using the Newton–Euler equations form we could write:

$$ \begin{aligned} m\dot{u} &= p(t) + r \\ I\dot{\omega} &= -\omega \wedge (I\omega) + M_{p(t)} + M_r \end{aligned} \qquad (2) $$

where $u$ is the velocity function, $\omega$ the angular velocity around the body centre of inertia, $p(t)$ and $M_{p(t)}$, respectively, the force and torque due to the applied loads, $r$ and $M_r$, respectively, the force and torque due to contact reactions, $m$ and $I$, respectively, the mass and inertia matrices. The second equation of the system may be non-linear due to the first term on the right-hand side. This term disappears in a 2D formulation and when the body is geometrically isotropic.

Due to the presence of impacts, the velocity function is non-smooth. Assuming it admits locally bounded variations, one may rewrite the problem in terms of differential measures [2, 3]. The measure formalism is particularly relevant because it allows one to describe continuous motion (Lebesgue measure) and impacts (Dirac measure) in the same way.
## 2.2. Time integrator

Whatever the integration scheme, we are looking for a discrete solution of differential equations, i.e. for instance: $t_0 < \cdots < t_n$. For the smooth case, we can use any available integrator. Those schemes have different properties in terms of precision, stability, energy conservation. In non-smooth cases, two approaches are possible:

* Event-driven approach where the temporal discretization is adjusted to non-smooth events. We try to separate the regular part of the motion (free flight, established contact and separation) from the non-smooth part of the motion (impacts and changing in the frictional status).

During smooth motion this method can use any high order integrator [4], but detection and treatment of non-smooth events need particular and careful formulation. In the presence of a great number of simultaneous contacts and finite accumulation of contacts, this kind of approach can prove itself useless because the time step between two events becomes too short. It is, however, widely used in the robotic field, where this situation does not occur.

* Time-stepping approach, where the temporal discretization is arbitrarily fixed and not necessarily constant. In this case, the integrator order is necessarily low. On the other hand, we can deal with a great number of discontinuities during one time step. Among the DEM using time-stepping approach, we could find two types of integrators: explicit schemes for methods based on so-called molecular dynamics [5] and implicit schemes for methods based on contact dynamics [3].

## 2.3. Interaction definition

We describe here the main ingredients for the definition of interactions between two bodies:

* the geometrical description of contactors;

* the interactions detection algorithm;

* the modelling of the interaction behaviour (variables, associated forces, etc.).

The geometrical description of the contactor could be done in several ways: analytical, polynomial, piece-wise linear, interpolated, smoothed, etc. The control of the evolution of the contact during the movement is necessary. The choice of the geometrical description is an important task, which influences the local computations of the variables as the relative gaps and the velocities between bodies.

Moreover, when studying a large collection of bodies, this description will influence the contact detection and the computation times. In our case, we will use an analytical description of contactors. However, in the case of faceted rigid bodies (polygons and polyhedrons) or deformable bodies, it is necessary to adopt an interpolated (or smoothed) description, which can be memory and computing time consuming.

The general principle used for the detection of the interactions is the following:

* Detection of a rough proximity of bodies (pre-detection). Starting from rough geometrical informations and alert distance values, we identify potential pairs of bodies in contact.

* Visibility test of bodies. A visibility table allows an arbitrary selection of potential pairs of bodies in contact.

* Coarse contact computation: local frame, gap, relative velocity, etc.
One can notice that the configuration of the contact is implicitly dependent on the problem solution. The fact of having small time steps leads us to an explicit contact detection. A leap-frog technique used here presents many advantages [6].

We use a node-to-node approach. However, in the case of rigid contactors where line-to-line contact occurs for polygons, we may use two simultaneous contact points in order to transmit torque. In fact, this assumes a slight concavity of the body shapes.

The behaviour laws of the contact between a candidate body (A) and an antagonist body (B) are described by point-to-point relations (see Figure 1). We suppose that we could define, at every moment, the two nearest points of (A) and (B) envelopes. Let (t, n) be a local frame linked to the antagonist body at the contact point; the normal vector n is toward the candidate body; t defines the tangent line of the antagonist body.

At any time, we can define for one contact:

* g the gap between two bodies; this distance could be seen as the normal component of **QP** vector in the local frame (Q, t, n),

* $U = (U_T, U_N)$ the relative velocity components in the local frame. It comes out from classical kinematic analysis that:

$$U = H^*(q)u \quad (3)$$

* $R = (R_T, R_N)$ the reaction components in the local frame. From duality considerations (equality of the power expressed with either local variables or global variables) it comes that the contribution $\mathbb{R}$ due to the local reaction $R$ satisfies the relation:

$$\mathbb{R} = H(q)R \quad (4)$$

The most used interaction laws are:

* Signorini's contact law: in its classical form, Signorini's law links the gap to the normal reaction by a relation, defined as 'complementarity' (Figure 2):

$$g \ge 0, \quad R_N \ge 0, \quad gR_N = 0 \quad (5)$$

In the case of contacts between rigid bodies, Moreau [2] shows that when the contact is established the previous relation is equivalent to the following one written in terms of

<figure><img src="image_2.png" /><figcaption>Figure 1. Interaction between two bodies.</figcaption></figure>
normal relative velocity:

$$U_N \ge 0, \quad R_N \ge 0, \quad U_N R_N = 0 \qquad (6)$$

* Coulomb's friction law: this threshold law is a relation between the tangential relative velocity and the tangential reaction (Figure 3):

$$\|R_T\| \le \mu R_N \quad \left\{ \begin{array}{l} \|R_T\| < \mu R_N \Rightarrow U_T = 0 \\ \|R_T\| = \mu R_N \Rightarrow U_T = -\lambda \frac{R_T}{\|R_T\|} \quad \lambda > 0 \end{array} \right. \qquad (7)$$

There are other interaction laws. For example see References [7–9]. Moreover, if we study the contact between rigid bodies, it is more appropriate to add an impact law to the frictional contact laws. That law could track complex physical phenomena that occur in the short time of the contact. Those phenomena such as impacts and deformability are usually ignored, at the ordinary observation scale [10].

<figure><img src="image_3.png" /><figcaption>Figure 2. Signorini's law.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Coulomb's law.</figcaption></figure>
## 2.4. Strategies and resolution methods

As mentioned above, the choice and combination of the different components of a contact solver (type of space and temporal discretization, type of detection algorithm, type of interaction laws) give us a lot of contact resolution methods. The bibliography around existing approaches is very rich. We can find in Reference [11] a good survey about numerical simulation of dynamical non-smooth systems, and in Reference [12] a description of some strategies and resolution methods. In the following paragraph, we will present the latest two components of the NSCD solver:

* the way of formulating a contact problem;

* the resolution method.

We have summarized in Figure 4 the various ingredients of the problem.

At a global level, we have the dynamic system of the set of dynamical equations of the different bodies. At a local level, we have contact laws that link couples of points. To switch from one level to another, we have $H$ and $H^*$, matrices that depend on the form of contact zone, on the discretization of the body kinematic, etc.

The resolution approaches are the following:

* We can solve the problem simultaneously at local and global level. That is the approach used by Alart and Curnier [13].

* We can solve the problem at a global level. It is the elimination or substitution strategy. It can lead to the use of penalization or Lagrange multiplier method [14].

* We can solve the problem at a local level. It is the condensation strategy. It can lead to the use of LCP or SQP solvers [15].

Moreover, the fact that we have to solve simultaneously relations with normal components (contact) and tangential components (friction), leads also to various approaches:

* we can solve simultaneously the two types of relations (block approach);

* we can solve successively one problem first, and then the second one until we obtain the convergence (fixed point approach). In this case, we can treat the two types of relations at different levels. For example, Reference [16] have used a strategy to solve contact at a global level and friction at a local one.

After the choice of the resolution strategy, we determine the resolution method. Some methods are inspired from the optimization methods under constraints (for example Reference [17]).

## 2.5. Non-smooth contact dynamics (NSCD) method

One keypoint of this method is the use of a relevant mathematical formulation for non-smooth contact dynamic problems, i.e. able to cope with shocks (the velocity function only possess

<figure><img src="image_5.png" /><figcaption>Figure 4. Duality graph.</figcaption></figure>
bounded variation property and no continuity), with unilaterality (multi-valued contact law) and friction (threshold contact law). As stated before, one needs a dynamical equation written in terms of differential measures [2].

Following the presentation made before, the NSCD method consists in:

* Time-stepping method. Considering the integration property of a differential measure [2], one can perform the discretization of the dynamical equation. On a time interval $[t_i, t_f]$ of length $h$ we obtain:

$$ 
\begin{aligned}
\mathbb{M}(u_f - u_i) &= \int_{t_i}^{t_f} (\mathbb{F}(q, u, s) + \mathbb{P}(s)) \, ds + h \mathbb{R}_f \\
q_f &= q_i + \int_{t_i}^{t_f} u \, ds
\end{aligned}
\qquad (8) $$

where $h\mathbb{R}_f = (1/h) \int_{[t_i, t_f]} \mathbb{R} \, dv$ represents the total contact impulse over the time interval, our unknown of the problem. In the following we denote $h\mathbb{R}_{\text{free}} = \int_{t_i}^{t_f} (\mathbb{F}(q, u, s) + \mathbb{P}(s)) \, ds$.

* Implicit time integrator. We use a classical $\theta$ method:

$$ 
\begin{aligned}
\int_{t_i}^{t_f} (\mathbb{F}(q, u, t) + \mathbb{P}(t)) \, dt &= h((1 - \theta)(\mathbb{F}(q_i, u_i, t_i) + \mathbb{P}(t_i)) + \theta(\mathbb{F}(q_f, u_f, t_f) + \mathbb{P}(t_f))) \\
q_f &= q_i + h((1 - \theta)u_i + \theta u_f)
\end{aligned}
\qquad (9) $$

* Leap-frog explicit contact treatment. In fact $H$ and $H^*$ depend implicitly on the solution $(q_f)$. Rather than performing updates, we use a leap-frog approach, evaluating explicitly the contact condition in the $q_m$ configuration:

$$ q_m = q_i + h(1 - \theta)u_i \qquad (10) $$

* Non-smooth interactions. As presented before, we use Signori's and Coulomb's laws. Other laws may be used in order to take into account the mortar between blocks.

* Local method. In a local approach, we condensate (explicitly or numerically) the dynamic equations at contacts level, then we solve the local problem and finally we update the global system.

We condensate the global system under the following form:

$$ H^*\mathbb{M}^{-1}HhR - H^*\dot{q}_f = -H^*(\mathbb{M}^{-1}h\mathbb{R}_{\text{free}} + \dot{q}_i) \Rightarrow WhR - U = U_{\text{free}} \qquad (11) $$

Also considering the contact $\alpha$ we obtain the block equation (normal and tangential part):

$$ U^{\alpha} = U_{\text{free}}^{\alpha} + \sum_{\beta=1}^{n_c} W^{\alpha\beta} h R^{\beta} \qquad (12) $$

where $n_c$ is the contact points number.

The local solver is a non-linear Gauss Seidel by block method. It was developed by Jean and Moreau [3]. Considering the contact $\alpha$ we fix all the contact reactions except those of the contact $\alpha$. It is quite simple to find a solution $(U^{\alpha}, R^{\alpha})$ which verifies this linear equation and
the contact laws. In a 2D problem this solution is the intersection of two graphs. In 3D it requires the use of an iterative solver such as the 'Newton-like' one [13]. Once this contact is solved, we will solve iteratively all couples of equations, for each contact considering the other one as 'frozen'.

Once the contact reactions are obtained, the new global solution is obtained as follows:

$$ 
\begin{align} 
\mathbb{M}(u_f - u_i) &= h \mathbb{R}_{\text{free}} + h \mathbb{R}_f \tag{13} \\ 
q_f &= q_m + h \theta u_f \nonumber 
\end{align} 
$$

The NSCD method was used for the simulation of granular media made up of rigid bodies [2, 3] or of deformable bodies [18, 19], with complex interaction laws [13, 20]. It can even be paralleled for the cases with a very great number of bodies [21].

### 3. CONSTRUCTION OF MECHANICAL MACROSCOPIC TENSORS

In spite of considering the masonry as a collection of rigid blocks, the structure can be seen as a continuous deformable medium from a macroscopic point of view. The construction of strain and stress tensors can then appear as an issue.

The characterization of the deformation by the invariants of the strain and the strain rate tensors has great importance in the dynamic case. It gives precious information about the mechanism of dilation/contraction well known in the geotechnical field.

In the case of masonry structures, the equivalent mechanisms are the opening/closing of joints. This mechanism is responsible for the stiffness softening of the masonry structure [22].

The average stress tensor allows one to evaluate the strength of the structure with respect to a given failure criterion. This tensor is important to detect material instability such as shear zones.

In the following sections we introduce the definitions of these two macroscopic indicators in a representative volume composed of the considered grain and its neighbours (see Figure 5) and we verify their duality from an energy point of view.

<figure><img src="image_6.png" /><figcaption>Figure 5. Cell for the definition of stress and strain tensors in granular media.</figcaption></figure>
### 3.1. Construction of a strain tensor

The definition of a strain tensor remains a subtle notion in the case of a collection of rigid grains. The aim of this tensor is to describe the local repositioning of blocks. We can get information about the eventual sliding between grains, using the deviatoric part of this tensor. Its spherical part reflects the opening of joints and the extension of the sample. Bagi [23] chose a purely geometrical definition of this tensor, consisting in the definition of strain from the displacement of the centre of gravity of the grains. Satake [24] considered the same definition but adopted a Dirichlet's tessellation of the granular collection. The deformation cell (see Figure 5) is a patch composed of a central grain and its neighbours.

As shown in Figure 5, each vertex $i$ of the triangle ($i=1...3$) represents the centre of gravity of a block. Each block has its own velocity $u_i$ and the average velocity on the triangle is given by the following expression:

$$V = \sum_{i=1...3} N_i(x, y)u_i \quad (14)$$

where $N_i$ ($i=1...3$) are the shape functions usually defined in the finite elements technique. The strain rate inside the triangle $T_k$ is given by the following expression:

$$D^{T_k} = \frac{1}{2} [\text{grad } V_{T_k} + \text{grad } V_{T_k}^{\text{t}}] \quad (15)$$

where grad $V_{T_k}$ is the velocity gradient tensor inside the triangle $T_k$ and grad $V_{T_k}^{\text{t}}$ its transpose. The average velocity gradient, grad $V_{T_k}$ is computed using the following interpolation:

$$\text{grad } V_{T_k} = \sum_{i=1...3} \text{grad } N_i(x, y)u_i \quad (16)$$

Once the strain rate tensor of each triangle of the cell is determined, we determined the strain rate tensor averaging the tensors weighted by the respective surface of the triangles. The cell strain rate tensor is written as follows:

$$D = \frac{1}{V_{\text{cell}}} \sum_{T_k} D^{T_k} \cdot V_{T_k} \quad (17)$$

where $V_{\text{cell}} = \sum_{T_k} V_{T_k}$.

In order to study the deformation of the cells, we define two instantaneous invariants of the strain tensor as follows:

* $\bar{\epsilon}_s$ the cumulated spherical part of the strain rate tensor. It informs about the variation of the structure volume, i.e. opening/closing of the joints.

$$\bar{\epsilon}_s = \bar{\epsilon}_{s0} + \int_{t_0}^{t} \text{tr}(D_s) \, dt \quad (18)$$

* $\bar{\epsilon}_d$ the cumulated deviatoric part of the strain rate tensor. It informs about sliding and shear zones in the structure.

$$\bar{\epsilon}_d = \bar{\epsilon}_{d0} + \int_{t_0}^{t} \sqrt{D_d : D_d} \, dt \quad (19)$$
## 3.2. Construction of a stress tensor

Most authors who proposed a stress tensor definition in granular media have dealt with the quasi-static case [7, 24–26]. Even in this case, there is some disagreement about the symmetry of this tensor (see Reference [27]). But Cundall and Strack [5, 28], Christoffersen et al. [29] consider that the asymmetry of this tensor is negligible or absent.

Caillerie [30] finds that the asymmetry of the stress tensor depends on the way stresses are defined. The stresses defined in the static case are symmetric when there is no contact moment.

However, the stresses defined by virtual work may become asymmetric when there is no contact moment. Bardet et al. [27] studied the asymmetry of this tensor and showed that this asymmetry exists and decreases with the volume size. It is more detectable in elongated samples subjected to external moments in their boundary. They also showed that the asymmetry could be neglected in large samples subjected to small external moments far from the boundaries.

Moreau [31] obtains a similar expression of this tensor, considering it as the internal moment of contacts for each body and using a variational approach. He shows the symmetry of this tensor without rolling forces.

Recently, Fortin et al. [32, 33] proposed an expression of the stress tensor taking into account rolling and dynamical terms. They show that despite the little contribution of rotational inertia, this term is important to ensure the symmetry of the stress tensor in dynamical cases.

We proposed here to revisit this result using Moreau’s formalism.

In a continuous medium the stress tensor verifies:

$$ \operatorname{div} \sigma + \hat{f} = 0 \qquad (20) $$

where $\hat{f} = \rho(\mathbf{g} - \mathbf{\Gamma})$ represents the volume forces including inertia forces; $\mathbf{\Gamma}$ is the acceleration vector at any point of the body.

Taking the tensorial product with an arbitrary vector $\mathbf{x}$ and integrating on a domain $\Omega$ of volume $V$,

$$ \int_{\Omega} (\sigma_{ij, j} + \hat{f}_i) x_x d\Omega = 0 \qquad (21) $$

one obtains an equivalent form:

$$ \int_{\Omega} (\sigma_{ij} x_k)_j d\Omega - \int_{\Omega} \sigma_{ij} x_k, j d\Omega + \int_{\Omega} \hat{f}_i x_k d\Omega = 0 \qquad (22) $$

where $j$ indicates the covariant differentiated with respect to $x_j$ component. This equation is equivalent to the equation below:

$$ \int_{\partial \Omega} \sigma_{ij} x_k n_j d\Gamma - \int_{\Omega} \sigma_{ij} \delta_{kj} d\Omega + \int_{\Omega} \hat{f}_i x_k d\Omega = 0 \qquad (23) $$

where $\delta_{kj}$ is the Kronecker symbol and $n_j$ is the $j$th component of the normal to the surface $\partial\Omega$. So the components of the average tensor on $\Omega$ are obtained as follows:

$$ \langle \sigma_{ik} \rangle = \frac{1}{V} \int_{\Omega} \sigma_{ik} d\Omega = \frac{1}{V} \left( \int_{\partial \Omega} T_i x_k d\Gamma + \int_{\Omega} \hat{f}_i x_k d\Omega \right) \qquad (24) $$

using $T_i = \sigma_{ij} n_j$ and $\sigma_{ij} \delta_{kj} = \sigma_{ik}$.
If now we use this expression with a rigid body (volume V), we can construct the average stress tensor on it as follows:

$$ \langle \sigma_{ik} \rangle = \frac{1}{V} \left( \sum_{l=1,...,n_c} F_i^l x_k^l + \int_{\Omega} \hat{f}_i x_k \, d\Omega \right) \qquad (25) $$

where $n_c$ is the number of contact points for the considered body. $F_i^l$ is the ith component of the lth contact force applied to the body and considered, in this case, as an external force. $x_k^l$ is the kth component of vector $OM^l$ between the arbitrary point O and the contact point M^l.

An intrinsic form of the previous expression is

$$ \langle \sigma \rangle = \frac{1}{V} \left( \sum_{l=1,...,n_c} \mathbf{F}^l \otimes \mathbf{OM}^l + \int_{\Omega} \hat{f} \otimes \mathbf{OM} \, d\Omega \right) \qquad (26) $$

The expression of the tensor is independent from the choice of the point O as demonstrated below.

For a given point M in a rigid body, one can express its velocity and its acceleration in terms of the velocity ($V(G)$) and acceleration ($\Gamma(G)$) of the inertia centre and the angular velocity ($\Omega$) of the body:

$$ \mathbf{V}(M) = \mathbf{V}(G) + \mathbf{\Omega} \wedge \mathbf{GM} \qquad (27) $$

$$ \Gamma(M) = \frac{d\mathbf{V}(M)}{dt} = \Gamma(G) + \frac{d\mathbf{\Omega}}{dt} \wedge \mathbf{GM} + \mathbf{\Omega} \wedge (\mathbf{\Omega} \wedge \mathbf{GM}) \qquad (28) $$

According to the expression of $\Gamma(M)$ given by (28), we can split $\langle \sigma \rangle$ in three contributions respectively $\sigma_1$, $\sigma_2$ and $\sigma_3$:

$$ \sigma_1 = \frac{1}{V} \left[ \sum_{l=1,...,n_c} \mathbf{F}^l \otimes \mathbf{OM}^l + \int_{\Omega} \rho(\mathbf{g} - \Gamma(G)) \otimes \mathbf{OM} \, d\Omega \right] \qquad (29) $$

$$ \sigma_1 = -\frac{1}{V} \int_{\Omega} \rho \left( \frac{d\mathbf{\Omega}}{dt} \wedge \mathbf{GM} \right) \otimes \mathbf{OM} \, d\Omega \qquad (30) $$

$$ \sigma_3 = -\frac{1}{V} \int_{\Omega} \rho(\mathbf{\Omega} \wedge (\mathbf{\Omega} \wedge \mathbf{GM})) \otimes \mathbf{OM} \, d\Omega \qquad (31) $$

By using the fact that for the inertia centre we have: $\int_{\Omega} \rho \mathbf{OM} \, d\Omega = m \mathbf{OG}$ we can write,

$$ \sigma_1 = \frac{1}{V} \left[ \sum_{l=1,...,n_c} \mathbf{F}^l \otimes \mathbf{OM}^l + m\mathbf{g} \otimes \mathbf{OG} - \Gamma(G) \otimes \mathbf{OG} \right] \qquad (32) $$

and if we decompose $\mathbf{OM}^l$ as follows:

$$ \mathbf{OM}^l = \mathbf{OG} + \mathbf{GM}^l $$
and we use the fact that

$$ \left[ \sum_{l=1,...,n_c} \mathbf{F}^l + m\mathbf{g} - m\mathbf{\Gamma}(G) \right] \otimes \mathbf{OG} = \underline{\underline{0}} \qquad (33) $$

Finally, the $\sigma_1$ contribution may be written as

$$ \sigma_1 = \frac{1}{V} \left[ \sum_{l=1,...,n_c} \mathbf{F}^l \otimes \mathbf{GM}^l \right] \qquad (34) $$

One can observe that, if the terms $\sigma_2$ and $\sigma_3$ are equal to zero, and due to the equilibrium of torque, we obtain the classical form found in Reference [7] in a quasi-static case (i.e. $\mathbf{\Gamma} = \mathbf{0}$):

$$ \langle \sigma \rangle = \frac{1}{V} \sum_{l=1,...,n_c} \mathbf{GM}^l \otimes \mathbf{F}^l \qquad (35) $$

In fact, we have shown here that Equation (35) is valid as long as the rotational velocity and acceleration remain equal to zero, and not only in a quasi-static case.

Concerning the term $\sigma_2$ we could first notice that:

$$ \int_{\Omega} \rho \left( \frac{d\mathbf{\Omega}}{dt} \wedge \mathbf{GM} \right) \otimes \mathbf{OM} \, d\mathbf{\Omega} = \int_{\Omega} \rho \left( \frac{d\mathbf{\Omega}}{dt} \wedge \mathbf{GM} \right) \otimes \mathbf{GM} \, d\mathbf{\Omega} \qquad (36) $$

due to fact that:

$$ \int_{\Omega} \rho \left( \frac{d\mathbf{\Omega}}{dt} \wedge \mathbf{GM} \right) \otimes \mathbf{OG} \, d\mathbf{\Omega} = \frac{d\mathbf{\Omega}}{dt} \wedge \int_{\Omega} \rho \mathbf{GM} \, d\mathbf{\Omega} \otimes \mathbf{OG} = 0 \qquad (37) $$

where we use the well-known property of the centre of mass $\int_{\Omega} \rho \mathbf{GM} \, d\mathbf{\Omega} = 0$.

Now, introducing the anti-symmetric operator $\varepsilon$:

$$ (\mathbf{v} \wedge \mathbf{w})_i = \varepsilon_{ijk} v_j w_k \qquad (38) $$

one can write

$$ \left( \int_{\Omega} \rho \left( \frac{d\mathbf{\Omega}}{dt} \wedge \mathbf{GM} \right) \otimes \mathbf{GM} \, d\mathbf{\Omega} \right)_{ij} = \int_{\Omega} \varepsilon_{ipq} \frac{dw_p}{dt} x_q x_j \, d\mathbf{\Omega} \qquad (39) $$

$$ = \varepsilon_{ipq} \frac{dw_p}{dt} \int_{\Omega} \rho x_q x_j \, d\mathbf{\Omega} \qquad (40) $$

$$ = \varepsilon_{ipq} \frac{dw_p}{dt} [-I_{qj}(G) + I_0 \delta_{qj}] \qquad (41) $$

with $I_0 = \int_{\Omega} \rho(x_1^2 + x_2^2 + x_3^2) \, d\mathbf{\Omega}$ and $I_{qj} = \int_{\Omega} \rho((x_1^2 + x_2^2 + x_3^2) \delta_{qj} - x_q x_j) \, d\mathbf{\Omega}$.

$\varepsilon_{ipq}(dw_p/dt)$ is the $iq$ component of the anti-symmetric tensor $j(dw/dt)$ defined as follows:

$$ \varepsilon_{ipq} \frac{dw_p}{dt} = - \varepsilon_{ipq} \frac{dw_p}{dt} = - \left( j \left( \frac{dw}{dt} \right) \right)_{iq} \qquad (42) $$
It can be represented in the following matrix form:

$$j \left( \frac{dw}{dt} \right) = \begin{pmatrix} 0 & -\frac{dw_3}{dt} & \frac{dw_2}{dt} \\ \frac{dw_3}{dt} & 0 & -\frac{dw_1}{dt} \\ -\frac{dw_2}{dt} & \frac{dw_1}{dt} & 0 \end{pmatrix} \qquad (43)$$

Assuming in the following that we work in the principal inertia frame, the inertia matrix *I* could be written as follows:

$$\begin{pmatrix} \int_{\Omega} \rho(x_2^2 + x_3^2) \, d\Omega & 0 & 0 \\ 0 & \int_{\Omega} \rho(x_1^2 + x_3^2) \, d\Omega & 0 \\ 0 & 0 & \int_{\Omega} \rho(x_1^2 + x_2^2) \, d\Omega \end{pmatrix} \qquad (44)$$

Indeed $I_0\delta - I$ can be written as follows:

$$\begin{pmatrix} \int_{\Omega} \rho x_1^2 \, d\Omega & 0 & 0 \\ 0 & \int_{\Omega} \rho x_2^2 \, d\Omega & 0 \\ 0 & 0 & \int_{\Omega} \rho x_3^2 \, d\Omega \end{pmatrix} = \begin{pmatrix} A & 0 & 0 \\ 0 & B & 0 \\ 0 & 0 & C \end{pmatrix} \qquad (45)$$

And finally $\sigma_2$ is equivalent to:

$$\sigma_2 = -\frac{1}{V} \int_{\Omega} \rho \left( \frac{d\Omega}{dt} \wedge \mathbf{GM} \right) \otimes \mathbf{GM} \, d\Omega \qquad (46)$$

$$= \frac{1}{V} \left( j \left( \frac{dw}{dt} \right) \cdot (I_0 \underline{\delta} - \underline{I}) \right) \qquad (47)$$

$$= \begin{pmatrix} 0 & -B \frac{dw_3}{dt} & C \frac{dw_2}{dt} \\ A \frac{dw_3}{dt} & 0 & -C \frac{dw_1}{dt} \\ -A \frac{dw_2}{dt} & C \frac{dw_1}{dt} & 0 \end{pmatrix} \qquad (48)$$
Finally, to calculate $\sigma_3$ given by Equation (31) we could write

$$ \int_{\Omega} \rho(\Omega \wedge (\Omega \wedge \mathbf{GM})) \otimes \mathbf{OM} \, d\Omega = \int_{\Omega} \rho(\Omega \wedge (\Omega \wedge \mathbf{GM})) \otimes \mathbf{GM} \, d\Omega \quad (49) $$

$$ \left( \int_{\Omega} \rho(\Omega \wedge (\Omega \wedge \mathbf{GM})) \otimes \mathbf{GM} \, d\Omega \right)_{ij} = \int_{\Omega} \rho \varepsilon_{ikl} w_k \varepsilon_{lpq} w_p x_q x_j \, d\Omega \quad (50) $$

where we use the decomposition $\mathbf{GM} = x_j \mathbf{x}_j$.

Using the property of the anti-symmetric operator given as follows:

$$ \varepsilon_{ikl} \varepsilon_{ipq} = \delta_{kp} \delta_{lp} - \delta_{kq} \delta_{lp} \quad (51) $$

one can write

$$ \int_{\Omega} \rho \varepsilon_{ikl} \varepsilon_{lpq} w_k w_p x_q x_j \, d\Omega = \int_{\Omega} \rho w_p x_q w_i x_j \, d\Omega - \int_{\Omega} \rho w_k w_k x_i x_j \, d\Omega \quad (52) $$

$$ = w_q w_i [-I_{qj} + I_{0\delta_{qj}}] + w_k w_k [-I_{0\delta_{kj}} + I_{ik}] \quad (53) $$

So the matrix form of Equation (49) is

$$ \int_{\Omega} \rho(\Omega \wedge (\Omega \wedge \mathbf{GM})) \otimes \mathbf{OM} \, d\Omega = ([I_{0\delta} - \underline{I}] \cdot \Omega) \otimes \Omega - (\Omega \cdot \Omega)[I_{0\delta} - \underline{I}] \quad (54) $$

Finally, collecting all the parts computed before, one can write the mean stress tensor $\langle \sigma \rangle$ of a rigid solid:

$$ \begin{aligned} \langle \sigma \rangle = \frac{1}{V} \left[ \sum_{l=1,...,n_c} \mathbf{F}^l \otimes \mathbf{GM}^l \right] &+ \frac{(w_1^2 + w_2^2 + w_3^2)}{V} \begin{pmatrix} A & 0 & 0 \\ 0 & B & 0 \\ 0 & 0 & C \end{pmatrix} \\ - \frac{1}{V} \begin{pmatrix} A w_1^2 & A w_1 w_2 & A w_1 w_3 \\ B w_1 w_2 & B w_2^2 & B w_2 w_3 \\ C w_1 w_3 & C w_2 w_3 & C w_3^2 \end{pmatrix} \\ + \frac{1}{V} \begin{pmatrix} 0 & -B \frac{dw_3}{dt} & C \frac{dw_2}{dt} \\ A \frac{dw_3}{dt} & 0 & -C \frac{dw_1}{dt} \\ -A \frac{dw_2}{dt} & B \frac{dw_1}{dt} & 0 \end{pmatrix} \end{aligned} \quad (55) $$
If we consider the 2D case, $w_1 = w_2 = 0$ and $w_3 = \dot{\theta}$, then

$$ \langle \sigma \rangle = \frac{1}{V} \left[ \sum_{l=1,...,n_c} \mathbf{F}^l \otimes \mathbf{GM}^l \right] + \frac{\dot{\theta}^2}{V} \begin{pmatrix} A & 0 \\ 0 & B \end{pmatrix} + \frac{\ddot{\theta}}{V} \begin{pmatrix} 0 & -B \\ A & 0 \end{pmatrix} \quad (56) $$

We remark that the dynamic term of the average stress tensor given by Equation (56) is written with respect to angular velocity $\dot{\theta}$, the angular acceleration $\theta$ and the inertia matrix of the body. The term which will allow to obtain the symmetry of the average stress tensor is the anti-symmetric one. When the stress tensor is defined for a set of solids it is easy to show that the first term of Equation (55) takes into account external forces to this set and the inertia terms are piece-wise calculated. V represents the volume of the cell.

## 4. NUMERICAL EXAMPLES

### 4.1. A basic example

In this example, we illustrate the asymmetry of the stress tensor given by Equation (35). The studied problem is the motion of a cubic block (1 m × 1 m) in frictional contact with a rigid support; the friction coefficient is equal to $\mu = 0.8$. The block is supposed to be rigid with a density equal to 1820 kg m⁻³.

The block is subjected to an horizontal movement with a velocity equal to $V = 1 \text{ m s}^{-1}$.

The intensity of the contact force at the right corner is more important than at the left one due to the high friction coefficient between the support and the block. The block in movement has a tendency to lose contact with the support at the left corner (see Figure 6).

We follow the symmetry of the average stress tensor of the block versus time. Let $\langle \sigma \rangle$ be the mean quasi-static stress tensor, $\delta$ is the relative gap between the non-diagonal terms of the tensor:

$$ \delta = \frac{|\sigma_{12} - \sigma_{21}|}{\|\sigma\|_2} \quad (57) $$

where $\|\sigma\|_2 = \sqrt{\sigma : \sigma}$ is the quadratic norm of the tensor $\sigma$ and $\delta$ indicates the loss of symmetry of the stress tensor. In this case, this tensor loses its physical significance. Moreau [31] states that, without rotation, expression (35) gives symmetrical average stress tensors in granular media. In fact, expression (35) does not include the rotation inertia. Figure 7 shows that the symmetry of the tensor is quasi-perfect till the moment when the block loses contact with

<figure><img src="image_7.png" /><figcaption>Figure 6. A block in a frictional movement along a rigid support.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7. Evolution of the $\delta$ indicator.</figcaption></figure>

the horizontal rigid support. During the horizontal movement, the stress tensor is given by the following expression:

$$ \langle \sigma \rangle = \begin{pmatrix} 7843.1 & -9803.9 \\ -9810.0 & 12262.5 \end{pmatrix} \qquad (58) $$

The relative gap of the tensor symmetry is equal to $\delta = 6.25 \times 10^{-4}$. When the rotation velocity appears, the average stress tensor would lose its symmetry. The dynamic terms were not taken into account as stated in (55).

## 4.2. Strain tensor in an academic example

Here, we give an example of strain tensor in cells of a granular media. Every cell is composed of a central block and its neighbours.

The two calculations done correspond to two different types of excitations, quasi-static and dynamic, respectively. The first calculation is a settlement of the right part of the foundation (see Figure 8). The second calculation is a dynamic excitation equivalent to a sinusoidal velocity of amplitude equal to $1 \text{ m s}^{-1}$ and a frequency equal to 2 Hz during a period of 2 s (see Figure 9). The wall is composed by rigid stone blocks with a density equal to $1820 \text{ kg m}^{-3}$. The friction coefficient stone/stone is equal to 0.80. The foundation is rigid with a density equal to $2500 \text{ kg m}^{-3}$.

We present, for two cases of excitations, the variation of the strain velocity tensor in the different cells of the wall.
<figure><img src="image_9.png" /><figcaption>Figure 8. A 4 cm settlement of the foundation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. A sinusoidal velocity of the foundation.</figcaption></figure>

The cumulated dilation given in Figures 10 and 12 are the spherical invariants of the strain velocities tensors for the quasi-static and dynamic cases, respectively. Figures 11 and 13 give the deviatoric invariants of the strain velocities tensors for the quasi-static and dynamic cases, respectively.

4.2.1. *Macroscopic strain indicators.* We define a macroscopic strain indicator $\Delta$ as a relative surface ratio. Let $S_c$ be the surface of critical cells. That means the sum of surfaces of the cells where the dilation and the shear strain is greater than a threshold. Let $S_T$ be the sum of all cells surfaces. The $\Delta$ expression is a time cumulated indicator taking into account the softening induced by this joint opening/closing mechanism at the stone/stone interface. $\Delta$ is
<figure><img src="image_11.png" /><figcaption>Figure 10. Cumulated dilation for a settlement of the foundation of 4 cm.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Cumulated shear for a settlement of the foundation of 4 cm.</figcaption></figure>

given by the following expression:

$$ \Delta = \Delta_{i+1} = \Delta_i + \frac{S_c}{S_T} \qquad (59) $$

where $\Delta_{i+1}$ is the macroscopic strain indicator at the end of a time step and $\Delta_i$ is its value at the beginning of the time step.
<figure><img src="image_13.png" /><figcaption>Figure 12. Cumulated dilation at the end of the dynamic loading.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Cumulated shear at the end of the dynamic loading.</figcaption></figure>

In the following case, we consider the average strain invariants, along all the sample, as the values of interest.

Figure 14 shows the effect of the loading type on the evolution of $\Delta$.

We can see that, in the dynamic case, the macroscopic strain indicator $\Delta$ is bigger than in the quasi-static case. At the end of the dynamic excitation, $\Delta$ is equal to 12.50 but in the static case it is equal to 7.50. This indicator could provide good information for studying the discontinuous structures.
<figure><img src="image_15.png" /><figcaption>Figure 14. Macroscopic strain indicator evolution versus time.</figcaption></figure>

### 4.3. Influence of the model parameters on the stress tensor

In the following example, we study the equilibrium of a blocky wall composed of 248 blocks. At the beginning, the blocks are supposed rigid, with a density equal to $1820 \text{ kg m}^{-3}$. The foundation of this wall is supposed rigid, with a density equal to $2500 \text{ kg m}^{-3}$. The friction between the blocks is of Coulomb's dry friction type with a friction coefficient equal to 0.65. At first, we calculate the structure under gravity. We show that the results of this calculation are widely dependent on the construction history. These results are compared to those given by a discrete element approach with deformable blocks and eventually with a continuous finite element linear isotropic model.

#### 4.3.1. Influence of the construction history.
The calculation under gravity is done in two different ways. In the first case (scenario I), we directly calculate all the structure until the verification of a kinematic criterion of stability.‡ In the second case (scenario II), we activate the structure layer-by-layer and for every layer we do an intermediate calculation, taking into account the stored reactions of the previous calculation (Figures 15 and 16).

If we compare the two calculations block-by-block, on average, scenario II overestimates the stress tensor invariants by about 58%. Layer-by-layer, scenario II overestimates the stress invariants on every layer by about 7% for the spherical component and 16% for the deviatoric component (Figure 17).

#### 4.3.2. The influence of the joints in the masonry structures.
For the same example, we compare a discrete element model with elastic linear isotropic blocks and a continuous finite element model with the same properties. The blocks have a Young modulus equal to $E = 16 \text{ GPa}$ and a Poisson's ratio equal to $\nu = 0.3$. These blocks have a density equal to $1820 \text{ kg m}^{-3}$. For the

<footnote>‡The quadratic norm of maximal block velocity is lower than $10^{-6} \text{ m s}^{-1}$.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 15. Comparison between the two scenarii layer-by-layer: spherical component.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. Comparison between the two scenarii layer-by-layer: deviatoric component.</figcaption></figure>

discrete model, we suppose that the foundation is rigid, with a density equal to 2500 kg m⁻³. The friction between the blocks is of Coulomb's dry friction type with a friction coefficient equal to 0.65. For the second model, the calculation was done using ANSYS finite element method software. The comparison between the average stress layer-by-layer is given in Figure 18.

On average, the difference between the two calculations is 241%. So the calculation without taking into account the presence of the joints raises the loading level about 2.4 times. In fact, from a macroscopic point of view, a continuum FEM calculation overestimates the loading state of such a structure.

## 4.4. Comparison between deformable and rigid block calculations

Figure 19 shows a comparison between the average stress invariants values for the two cases: deformable and rigid blocks.
<figure><img src="image_18.png" /><figcaption>Figure 17. Histogram of the variation of scenario II compared with scenario I.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18. Comparison between a discrete element model and a continuous elastic linear isotropic one.</figcaption></figure>

We remark that there is only about 6% of the rigid blocks that have their average stress tensor invariants bigger than those of the deformable case. These values are usually less than twice these of the deformable ones. For the 94% of other blocks, the deformable blocks have their value at least equal to the rigid ones. By comparing the deformable and the rigid cases, the difference of the stress tensor corresponds to 36% for the spherical component of the tensor and 52% for the deviatoric component. An important observation is that the deformable case could present higher values of this stress tensor, compared to the rigid case. This detail is important for the following step of our work consisting in the substitution of rigid blocks by deformable ones.

In Figures 20 and 21, we compare the average stress invariants by layer, for the two cases, deformable and rigid blocks. Locally, the comparison between the two cases shows a difference of about 25% of the average stress tensor invariants. In the rigid case, the spherical component
<figure><img src="image_20.png" /><figcaption>Figure 19. Histogram of comparison between rigid and deformable blocks.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 20. Comparison between rigid and deformable blocks layer-by-layer: spherical component.</figcaption></figure>

by layer, on an average, is 25% bigger than in the deformable case and the deviatoric component is 25% smaller than in the deformable case. We remark that averaging the stress tensor in more than one block decreases the difference between the rigid approach, where the definition of stress tensor is done as a post-processing treatment and the deformable one where the notion of stress tensor is more realistic and attached to the intrinsic properties of the blocks. Table I gives macroscopic average stress invariants $p_{wall}$ and $q_{wall}$ consisting, on averaging, the stress invariants of the blocks or the layers. We remark that the difference between deformable and rigid cases is negligible. By checking the influence of this scale effect on results, on one hand, one can say that the stress tensor notion in rigid blocks, is more significant when it includes more than one block. On the other hand, according to the results of Table I, we could define two macroscopic variables $p_{wall}$ and $q_{wall}$ of a structure that could be defined either with a rigid or a deformable approach.
<figure><img src="image_22.png" /><figcaption>Figure 21. Comparison between rigid and deformable blocks layer-by-layer: deviatoric component.</figcaption></figure>

<table><caption>Table I. Averaging the stress invariants on the wall.</caption><thead><tr><th></th><th>Rigid case (scn.I)</th><th>Rigid case (scn.II)</th><th>Deformable case</th></tr></thead><tbody><tr><th>p<sub>wall</sub> [Pa]</th><td>1.57 × 10<sup>4</sup></td><td>1.59 × 10<sup>4</sup></td><td>1.55 × 10<sup>4</sup></td></tr><tr><th>q<sub>wall</sub> [Pa]</th><td>1.29 × 10<sup>4</sup></td><td>1.30 × 10<sup>4</sup></td><td>1.28 × 10<sup>4</sup></td></tr></tbody></table>

<table><caption>Table II. Contacts number and calculation time.</caption><thead><tr><th></th><th>Rigid case</th><th>Deformable case</th></tr></thead><tbody><tr><th>Contacts number</th><td>1391</td><td>1932</td></tr><tr><th>Calculation time</th><td>28 m n, 34.6 s</td><td>91 m n, 5.3 s</td></tr></tbody></table>

In spite of the variation of the results of the rigid approach, we usually choose the rigid configuration of blocks for the modelling of masonry structures instead of the deformable one, due to calculation time considerations (see Table II). In fact, in our case, a rigid calculation is three times faster than a deformable one. The rigid approach gives an average local information. The deformable approach is more realistic but also more expensive. Could we combine the two approaches, to get a more accurate calculation and an acceptable calculation time? The idea is to do a mixed rigid/deformable calculation: a rigid calculation in, for example, 80% of the structure and the other 20%, the most loaded blocks, by the introduction of their rheological behaviour (deformable blocks).
## 5. REAL CASE CALCULATION: PONT JULIEN (FRANCE 1ST CENTURY BC)

In this section, we present some calculations as an application of the notion of the average stress tensor on rigid blocks. We treat the example of a realistic modelling of a roman bridge called Pont Julien in Vaucluse (South of France). This bridge was constructed in the first century BC. It is composed of three arches. It has a length of 117.70 m, a width of 5.90 m. It is composed of three arches, respectively, of 10.20, 16.00 and 10.20 m spans. It has a maximal height of 9.00 m (see Figure 22).

In a first step we present a modelling of the bridge with the NSCD rigid discrete approach (see Figure 23). After that, a deformable discrete modelling with the same method (see Figure 24) allows us to compare the results: average stress tensors on the blocks, time calculation of the structure, Mohr–Coulomb failure criterion.

In this case, the structure is considered as a granular medium composed of 1118 blocks in

<figure><img src="image_23.png" /><figcaption>Figure 22. Pont Julien (France—1st century BC).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 23. Pont Julien discrete element modelling: rigid blocks.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 24. Pont Julien discrete element modelling: deformable blocks.</figcaption></figure>

contact with a rigid foundation. In the 2D deformable discrete approach, the blocks are meshed by triangular finite elements. They are supposed linear elastic with a Young modulus equal to $E = 16 \text{ GPa}$, a Poisson's ratio equal to $\nu = 0.25$, a density $\rho = 1820 \text{ kg m}^{-3}$. In the two cases, the foundation is supposed rigid, with a density equal to $2500 \text{ kg m}^{-3}$. The stone/stone friction is a dry Mohr-Coulomb type, with a friction coefficient equal to 0.80. The calculation is done till the static equilibrium is obtained. The accurate description of the block was performed from photography and measurements on site.

## 5.1. Calculation results

The stopping calculation criterion is of a kinematic type. It consists in the realization of the maximal velocity quadratic norm less than a given tolerance that we have choosen, in this case, equal to $10^{-6} \text{ m s}^{-1}$. The advantage of this type of criterion is, on one hand, to ensure that we have obtained a quasi-equilibrium state of the structure. On the other hand, we get a uniform stopping criterion of the two types of calculation in order to estimate correctly the calculation times. Figure 25, shows the total number of contacts between the blocks versus time. They are three types of contacts:

* Sticking contacts: They are the contacts which verify $R_T < \mu R_N$ in Coulomb's graph (see Figure 3).

* Sliding contacts: They are the contacts that reach the threshold of sliding $R_T = \mu R_N$ in Coulomb graph. For these contacts, we reach the allowed maximal value of shear effort.

* No contact state: It is the case where the gap is less than the alert distance introduced previously in NSCD method section, but the contact is not active.

Figures 26–28 give the network of weak, strong contacts and the principal stress directions, respectively. By weak contacts we mean those where the contact normal reactions have a value lower than the mean value on all the sample. Consequently, the strong contacts are those who have values bigger than the mean value in the sample.

The spherical invariant of the average stress tensor of the two calculations, at equilibrium, are given by Figures 29 and 30.
<figure><img src="image_26.png" /><figcaption>Figure 25. Contacts number and type in the structure.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 26. Strong contacts: normal reaction more than the average value.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 27. Weak contacts: normal reaction less than the average value.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure 28. Principal stress directions of the blocks.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 29. Spherical stress tensor invariant: rigid case.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 30. Spherical stress tensor invariant: deformable case.</figcaption></figure>
6. CONCLUSION

This work is an occasion to show the influence of the modelling approach on the results, when considering both static and dynamical loads. A rigid approach is faster than a deformable one but less realistic. The dispersion of the results of the average stress tensor can be reduced when averaging this tensor on more than one block. The stress tensor averaged on all the assembly of blocks has a constant value independent of the type of the approach but function of the geometry, the loads applied to the structure, the mechanical and physical properties of the blocks. In our future work, we will develop the mixed rigid/deformable blocks approach that seems important to optimize the calculation time and to give more realistic information on the behaviour of real structures. The automated process of substitution will represent a valuable improvement of the method, to allow its application to real engineering cases.

REFERENCES

1. Moreau JJ. Some numerical methods in multibody dynamics: application to granular materials. *European Journal of Mechanics - A/Solids* 1994; **13**(4-suppl.):93-114.

2. Moreau JJ. Unilateral contact and dry friction in finite dimension freedom dynamics. *CISM Courses and Lectures*, vol. 302. Springer: Berlin, 1988; 1-82.

3. Jean M, Moreau JJ. Unilaterality and dry friction in the dynamics of rigid body collections. *Proceedings of Contact Mechanics International Symposium* 1992; 31-48.

4. Abadie M. Simulation dynamique de mécanismes. Prise en compte du contact frottant. *Ph.D. Thesis*, Université de Montpellier 2, 1998.

5. Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. *Geotechnics* 1979; **29**(1): 47-65.

6. Pratt E, Lebon M, Raous M. Consistent time discretization of dynamical frictional contact problems. *European Revue of the Finite Elements* 1998; **7**(1-3):149-162.

7. Cambou B, Jean M. Micromechanics of granular materials. *Hermes Sciences*, 2001.

8. Jean M. The non-smooth contact dynamics method. *Computer Methods in Applied Mechanics and Engineering* 1999; **177**(3-4):235-257.

9. Dubois F. Simulation of multi-bodies contact problems. *Support of Course of the School Dynamic set of themes CNRS 'nonregular'*, 2002, http://maply.univ-lyon/dynimpact.

10. Stronge WJ. *Impact Mechanics*. Cambridge University Press: Cambridge, 2000.

11. Brogliato B, Ten Dam AA, Paoli L, Genot F, Abadie M. Numerical simulation of finite dimensional multibody nonsmooth mechanical systems. *Applied Mechanics Reviews* 2002; **55**(2):107-150.

12. Jean M. Digital simulation of the problems of contact with friction. *Days of the Tribological Company of France, Materials and Techniques* 1993; 1-3.

13. Alart P, Curnier A. Mixed formulation for frictional contact problems by Newton like solution methods. *Computer Methods in Applied Mechanics and Engineering* 1991; **92**:353-375.

14. Laursen T. *Computational Contact and Impact Mechanics*. Springer: Berlin, 2002.

15. Dubois F. Contact, frottement, grandes déformations élastoplastiques. Application à l'emboutissage. *Ph.D. Thesis*, Université de la méditerranée, 1994.

16. Horkay F, Jean M, Mehrez F. Unilateral contact and dry friction in numerical simulation of deep drawing. *Proceedings of NUMIFORM89*, 1989.

17. Ciarlet PG. *Introduction to the Matrix Numerical Analysis and Optimization*. Masson: Paris, NY, 1982.

18. Jean M. Frictional contact in collections of rigid and deformable bodies: numerical simulation of geomaterials. In *Mechanics of Geomaterials Interfaces*, Salvadurai M, Bolt JJ (eds). Elsevier: Amsterdam, 1995; 463-486.

19. Acary V. Contribution à la modélisation mécanique et numérique des édifices maçonnés. *Ph.D. Thesis*, Université de la Méditerranée, 2001.

20. Jean M, Acary V, Monerie Y. Non smooth contact dynamics approach of cohesive materials. *Proceedings of the Royal Society* (Special issue on non smooth mechanics) 2001.

21. Renouf M, Dubois F, Alart P. A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media. *Journal of Computational and Applied Mathematics* 2004; **168**:375-382.
22. Zucchini A, Lourenco PB. A micro-mechanical model for the homogenisation of masonry. *International Journal of Solids and Structures* 2002; **39**:3233-3255.

23. Bagi K. Geomaterial modelling of granular assemblies. *Acta Technica Academique Scientiarum* 1995; **107**(1-2):1-16

24. Satake M. Stress and strain in granular materials. In *Powders and Grains*, Kishino Y (ed.). Balkema AA: Rotterdam, 2001; 149-152.

25. Bagi K. Stress and strain in granular assemblies. *Mechanics of Materials* 1996; **22**:165-177.

26. Chetouane B, Dubois F, Vinches M, Bohatier C. Stress and strain tensors in granular medium application to masonry structures. *Proceedings of IMECE 03 International ASME Mechanical Engineering Congress and Exposure*, Washington, DC, 16-21 November 2003.

27. Bardet JP, Vardoulakis I. A discrete numerical model for granular assemblies. *International Journal of Solids and Structures* 2001; **38**:353-367.

28. Cundall PA, Strack ODL. The distinct element method as a tool for research in granular media, part I and II. *Report to National Science Foundation, Engineering, Department of Civil and Mineral Engineering*, University of Minnesota, Minneapolis, MN, 1978-1979, 76-20711.

29. Christoffersen J, Mehrabadi MM, Nemat-Nasser S. A Micromechanical description of granular material behavior. *Journal of Applied Mechanics* (ASME) 1981; **48**:339-344.

30. Caillerie D. Tenseur des contraintes dans un milieu granulaire. *Proceedings of GRECO Géomatériaux*. Aussois, Isere, 1991 (in French).

31. Moreau JJ. Numerical investigation of shear zones in granular materials. In *Proceedings of the HLRZ-Workshop One Friction, Arching, Dynamics Contact, Grassberger P, Wolf D (eds). World Scientific: Singapore, 1997; 233-247*.

32. Fortin J. Numerical simulation of multi-bodies systems dynamics applied to granular media. *Ph.D. Thesis*, Université de Lille I, 2000.

33. Fortin J, Millet O, De Saxcé G. Construction of the mean stress tensor for granular medium. *European Journal of Mechanics - A/Solids* 2003; **22**:567-582.

34. Moreau JJ. Contact et frottement en dynamique des systémes de corps rigides. *European Review of Finite Elements* 2000; **9**:9-28 (in French).

35. Salençon J. Calcul à la rupture et analyse limite. In *Presses de l'Ecole Nationale des Ponts et Chaussées*, Juin (ed.), 1983.