NSCD discrete element method for modelling masonry structures

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


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.

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: Mu = F(q, u, t) + P(t) + R boundary conditions + initial conditions (1) 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, 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: mu = p(t) + r I˙ = − ∧ (I ) + M p(t) + M r (2) where u is the velocity function, 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.

Time integrator
Whatever the integration scheme, we are looking for a discrete solution of differential equations, i.e. for instance: t 0 <· · ·<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 nonsmooth 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].

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: • 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 R due to the local reaction R satisfies the relation: 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): 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 1. Interaction between two bodies. normal relative velocity: • Coulomb's friction law: this threshold law is a relation between the tangential relative velocity and the tangential reaction ( Figure 3): There are other interaction laws. For example see References [7][8][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].

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]).

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 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: where our unknown of the problem. In the following we denote hR free = t f t i (F(q, u, s)+P(s)) ds. • Implicit time integrator. We use a classical method: • 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: • 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: Also considering the contact we obtain the block equation (normal and tangential part): 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 we fix all the contact reactions except those of the contact . It is quite simple to find a solution (U , R ) 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: 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].

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 5. Cell for the definition of stress and strain tensors in granular media.

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: (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: where grad V T k is the velocity gradient tensor inside the triangle T k and grad V t T k its transpose. The average velocity gradient, grad V T k is computed using the following interpolation: 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: where In order to study the deformation of the cells, we define two instantaneous invariants of the strain tensor as follows: •¯ 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.
•¯ d the cumulated deviatoric part of the strain rate tensor. It informs about sliding and shear zones in the structure.¯

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][25][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: wheref = (g − ) represents the volume forces including inertia forces; is the acceleration vector at any point of the body. Taking the tensorial product with an arbitrary vector x and integrating on a domain of volume V, one obtains an equivalent form: where j indicates the covariant differentiated with respect to x j component. This equation is equivalent to the equation below: where kj is the Kronecker symbol and n j is the j th component of the normal to the surface * . So the components of the average tensor on are obtained as follows: using T i = ij n j and ij kj = ik .
If now we use this expression with a rigid body (volume V), we can construct the average stress tensor on it as follows: where n c is the number of contact points for the considered body. An intrinsic form of the previous expression is 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 ( (G)) of the inertia centre and the angular velocity ( ) of the body: According to the expression of (M) given by (28), we can split in three contributions respectively 1 , 2 and 3 : By using the fact that for the inertia centre we have: OM d = mOG we can write, and if we decompose OM l as follows: OM l = OG + GM l and we use the fact that l=1,...,n c Finally, the 1 contribution may be written as One can observe that, if the terms 2 and 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. = 0): 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 2 we could first notice that: due to fact that: where we use the well-known property of the centre of mass GM d = 0. Now, introducing the anti-symmetric operator : with I 0 = (x 2 1 + x 2 2 + x 2 3 ) d and I qj = ((x 2 1 + x 2 2 + x 2 3 ) qj − x q x j ) d . ipq (dw p /dt) is the iq component of the anti-symmetric tensor j (dw/dt) defined as follows: It can be represented in the following matrix form: Assuming in the following that we work in the principal inertia frame, the inertia matrix I could be written as follows: ⎛ Indeed I 0 − I can be written as follows: ⎛ And finally 2 is equivalent to: Finally, to calculate 3 given by Equation (31) we could write where we use the decomposition GM = x j x j . Using the property of the anti-symmetric operator given as follows: So the matrix form of Equation (49) is Finally, collecting all the parts computed before, one can write the mean stress tensor of a rigid solid: If we consider the 2D case, w 1 = w 2 = 0 and w 3 =˙ , then We remark that the dynamic term of the average stress tensor given by Equation (56) is written with respect to angular velocity˙ , the angular acceleration¨ 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.

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 = 0.8. The block is supposed to be rigid with a density equal to 1820 kg m −3 .
The block is subjected to an horizontal movement with a velocity equal to V = 1 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 be the mean quasi-static stress tensor, is the relative gap between the non-diagonal terms of the tensor: where 2 = √ : is the quadratic norm of the tensor and 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 6. A block in a frictional movement along a rigid support. The relative gap of the tensor symmetry is equal to = 6.25 × 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).

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 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 kg m −3 . The friction coefficient stone/stone is equal to 0.80. The foundation is rigid with a density equal to 2500 kg m −3 .
We present, for two cases of excitations, the variation of the strain velocity tensor in the different cells of the wall.  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.

Macroscopic strain indicators.
We define a macroscopic strain indicator 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 expression is a time cumulated indicator taking into account the softening induced by this joint opening/closing mechanism at the stone/stone interface. is  given by the following expression: where i+1 is the macroscopic strain indicator at the end of a time step and i is its value at the beginning of the time step.  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 .
We can see that, in the dynamic case, the macroscopic strain indicator is bigger than in the quasi-static case. At the end of the dynamic excitation, 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.

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 kg m −3 . The foundation of this wall is supposed rigid, with a density equal to 2500 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.

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).

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 GPa and a Poisson's ratio equal to = 0.3. These blocks have a density equal to 1820 kg m −3 . For the Figure 15. Comparison between the two scenarii layer-by-layer: spherical component. Figure 16. Comparison between the two scenarii layer-by-layer: deviatoric component. discrete model, we suppose that the foundation is rigid, with a density equal to 2500 kg m −3 . 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. Figure 19 shows a comparison between the average stress invariants values for the two cases: deformable and rigid blocks.  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.

Comparison between deformable and rigid block calculations
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  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.  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).

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   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 GPa, a Poisson's ratio equal to = 0.25, a density = 1820 kg m −3 . In the two cases, the foundation is supposed rigid, with a density equal to 2500 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.

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 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 < R N in Coulomb's graph (see Figure 3). • Sliding contacts: They are the contacts that reach the threshold of sliding R T = 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.

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.