Stability analysis of dry stone lintels using combined finite-discrete element method

The use of the combined finite-discrete element method (FEM/DEM) in the analysis of seismic resistance of dry stone lintels is presented in this paper. Seven dry stone masonry walls with different types of stone lintels were exposed to linearly increasing horizontal base acceleration and to seismic excitation. Based on the analysis results, the seismic resistance and failure mechanisms were analysed for each type of lintel. Recommendations for the design of stone lintels, aimed at achieving the greatest possible seismic resistance, were also discussed and presented.


Introduction
Historic buildings, some of which are under protection of the United Nations Educational, Scientific and Cultural Organization (UNESCO), are a common sight in many bigger and smaller towns in Croatia.The structures are mainly scattered throughout coastal cities, towns and places that were once under Roman administration.One of the best-known cities containing such heritage is Dubrovnik, whose old town is under the UNESCO protection, and whose 1940-metre-long city walls, built between 1300s and 1700s to protect the city, are famous worldwide.No less known are the city of Split, where the Diocletian's palace dating back to 295 AD is located, and the Pula town centre that hosts the Triumphal Arch of the Sergii (also known as the Golden Gate) erected between 29 and 27 BC (Figure 1.a), the Twin Gates built towards the end of the second century (Figure 1.b), the Temple of Augustus erected between 2 BC and 14 AD (Figure 1.c), and a Roman amphitheatre, the so called Pula Arena, which is the sixth largest Roman amphitheatre in the world (Figure 1.d).A large variety of numerical methods have so far been developed for the analysis of masonry structures.They vary in the level of accuracy, amount of input data, results to be obtained, effects that may occur in structures for various excitations, and the time needed to perform calculations.The most widely used method for structural analysis, which is also utilised in the analysis of stone structures [1][2][3][4][5][6], is the finite element method.In this approach, the structure is simulated as an orthotropic continuum in which the mean stress and strain relation is obtained experimentally, or using homogenization techniques.One of the drawbacks of modelling stone structures with numerical models based on continuum, but without the use of more specific modelling techniques (for instance, without formulation of displacement field that would account for discontinuities), is the inability to take into account large discontinuities in the displacement field between stone blocks, and the inability to simulate mechanical interaction between multiple stone blocks, which is important for analysing structures exposed to intense seismic load and progressive failure of stone masonry structures.Numerical models based on the discrete element method have been developed for the analysis of problems involving mechanical interaction between several bodies that may have large rotations and displacements.This method, primarily aimed at simulating the sliding and detachment of connected rock masses along predefined crack paths [7], has found its application in the analysis of masonry structures [8][9][10].The main feature of the discrete element method, enabling its use in the analysis of masonry structures, is its ability to present a structure as a set of discrete elements linked together with contact elements.This kind of approach allows simulation of structural failure due to rotation, sliding, and impact load.A disadvantage of models based on the discrete element method is their inability to describe stress and strain within the discrete elements themselves, which is particularly significant when analysing the appearance and propagation of cracks.Many attempts have recently have been made to exploit the advantages offered by both the finite element method, and the discrete element method [11][12][13].One of such combined methods is the combined finite-discrete element method developed by Ante Munjiza [14,15].This approach is based on the simulation of a large number of discrete elements that can be in interaction.Each discrete element is discretized by its own mesh of finite elements, which enables its deformability.The material non-linearity, including the appearance and propagation of cracks and, ultimately, fragmentation of discrete elements, is enabled by means of a model with contact elements [16], implemented between the finite elements.To take into account all these effects, appropriate algorithms have been developed.These algortihms include, in every time interval, the contact detection and interaction, stress and strain monitoring in the final and contact element, accounting for the appearance and propagation of cracks, integration of time-dependent equations of motion that cover large displacements and rotations and, finally, visualisation of these effects.All these favourable properties make the combined finite-discrete element method highly suited for the analyses of stone structures [17,18].The use of the combined finite-discrete element method in the analysis of seismic resistance of dry stone masonry lintels is presented in the scope of this paper.Seven dry stone masonry walls with different types of stone lintels were exposed to the linearly increasing horizontal base acceleration and seismic load.The seismic resistance and failure mechanisms were analysed for each type of lintel on the basis of the results obtained.Furthermore, recommendations are discussed and propose for the design that would lead to greater seismic resistance of dry stone lintels.Stability analysis of dry stone lintels using combined finite-discrete element method

Discretization of structure
In the proposed analysis, a dry stone masonry structure is considered as a set of discrete elements where every discrete element represents one stone block (Figure 1).Furthermore, each stone block is discretized by its own mesh of triangular finite elements.The material behaviour in finite elements is linear-elastic, while contact forces occurring between neighbouring blocks, as well as friction forces, are based on the principle of potential contact forces, as explained in more detail below.As the strength of stone is mostly much higher than the effective stress occurring in stone structures [17,19] and as, consequently, structural failure occurs mainly due to the loss of stability caused by relative displacements between blocks under the influence of a horizontal force, the possibility of stone cracking has not been considered in the analysis.

Deformability of finite elements
In the combined finite-discrete element method, the deformability of discrete elements is enabled through their discretisation via the finite element mesh.Given the need for an algorithm that would be as simple and fast as possible, and due to the fact that the calculation of contact forces is based on the same discretization scheme, the decision was made to adopt finite elements of the simplest possible geometry.For inplane problems, that is a three-node triangular finite element.The geometry of a 3-node triangular finite element is defined by global coordinates of each node (x, y), where (x i , y i ) represent the coordinates in their initial configuration, and (x c ,y c ) the coordinates in their current configuration (Figure 3).As discrete elements can change their position in space, their displacements can be divided into two distinct components: displacements of discrete elements as solid bodies by translation or rotation, and displacements causing deformations leading to changes in shape and volume.Displacements of a deformable body involving rotation and deformation in the vicinity of a certain point of a deformable body are defined by the deformation gradient F [14].Adopting the linear-elastic relationship between the stress and strain, the Cauchy stress tensor can be obtained from the following expression (3) where µ and λ are the Lamé constants, ε v is the volume deformation that is expressed by the following equation where the last element on the right-hand side of the expression (3) represents the contribution of the deformation velocity, and where µ is the damping coefficient and D is the rate of the deformation tensor.
The force belonging to each node and acting per unit of length on the side of a triangular element is obtained by direct integration of the stress tensor according to the following expression (5) where n x,c and n y,c are components of the unit external line normal to the side of the triangular finite element.

Contact detection and interaction
The objective of the contact detection algorithm is to locate pairs of neighbouring finite elements that are in contact, and to eliminate those that are too far away from one another and can no longer establish contact.Accordingly, the NBS (no binary Emili Zubović, Hrvoje Smoljanović, Boris Trogrlić search) algorithm for contact detection was implemented into the FEM/DEM model [20].The total time required for the detection of all contact pairs is proportional to the total number of discrete elements.Once the pairs of discrete elements are detected, the contact interaction algorithm [21] defines the contact forces between two discrete elements, one of which is then designated as the contactor and the other as the target (Figure 4).In the interaction algorithm, the distributed contact forces are defined by means of the penalty method, which is based on the principle of potential contact forces.When in contact, the contactor and the target overlap across the surface S bounded by the external edge Г βm∩βk .In such a case, the total contact differential force on the contactor df k is defined as (6) where P m and P k are the points in which the target and the contactor overlap, j is the corresponding potential function.By integrating (6) across the entire overlapping surface S, the total contact force is obtained as follows (7) which can also be written as (8) where n Г is the unit external line normal to the edge Г of the overlapping surface S.
In the framework of the contact forces algorithm, the Coulomb dry friction model for shear forces was implemented as follows where f t is the tangential elastic contact force, k t is the penalty coefficient for friction and d t the tangential vector of displacement between two elements [22].If f t is greater than the maximum friction force defined by the Coulomb law | f t |>µ| f n |, then elements slide along one another and the shear force between them is defined by means of the elastic normal force ft, according to the following expression where µ is the friction coefficient.

Time integration of equation of motion
In the combined finite-discrete element method, the shape and position of every discrete element in the plane is described by the current coordinates of finite element nodes.To calculate the current coordinates of nodes, it is necessary to take into account the mass of the corresponding system.In the combined finitediscrete element method, the system mass is concentrated in finite element nodes (Figure 5), which leads to a lumped-mass model.Time integration of the equations of motion in time was conducted explicitly for each corresponding node by means of the finite difference method [14], which is conditionally stable and whose stability and accuracy depend on the choice of a time step.The description of the update of variables can be written as (11) where x i , v i , f i , m i are the displacement vector, the velocity vector, the total mass vector, and the mass of each node, respectively, and Δt is the time step.

Use of FEM/DEM method for analysing seismic resistance of stone lintels
This section addresses the application of the FEM/DEM method in the analysis of seismic resistance of stone lintels.For this purpose, seven dry stone walls with different types of lintels were chosen.They were first exposed to a monotonically increasing horizontal base acceleration and, subsequently, to the base acceleration recorded during three real earthquakes.The geometry of walls and their respective lintels (Figures 6.a -12.a) was taken, with minor changes, from relevant literature [23].For the purposes of this research, the geometry of walls was partly modified with respect to the original pattern in that the dimensions of their stone blocks were taken to be identical up to Stability analysis of dry stone lintels using combined finite-discrete element method

Linearly increasing horizontal base acceleration
This section provides an overview of results for the analysed walls exposed to the linearly increasing horizontal base acceleration.The analyses were conducted for two different cases of base acceleration, and their changes over time are presented in Figure 13.These analysis were conducted to define the wall failure mechanism and to find a minimum horizontal base acceleration around the wall base line, at which a part of the wall starts to rotate around its end point due to the prevalence of the horizontal inertia force overturning moment over the gravity a certain height for each stone wall, and each wall was considered axisymmetric with respect to the vertical axis running through the centre of the wall.These alterations were made to facilitate comparison of different walls under the same load.The discretisation of walls by a finite element mesh was conducted in the Gmsh programme [24] and is presented in In the conducted analyses, the modulus of elasticity of stone blocks and the Poisson coefficient amounted to 48400 MPa and 0.2, respectively.The density of stone was assumed to be 2700 kg/m 3 , and the friction coefficient amounted to 0.6.The base of the walls was modelled as being absolutely stiff.The density of the final line of blocks was increased to take into account the influence of the supporting floor structure and to achieve the compressive stress of 0.2 MPa, which corresponds to the amount of stress most commonly occurring in stone structures [19].Horizontal base acceleration values at which the separation of stone blocks starts due to rotation of some parts of the walls are presented in Table 1.It can be observed that the minimum base acceleration that causes rotation of a part of a wall is almost the same for all types of walls and amounts to about 0.4 g, where g represents the adopted gravity constant of 9.81 m/ s 2 .Even if elastic properties of stone are not taken into account, it may reasonably be expected that any seismic excitation with peak acceleration lower than the one presented in Table 1 will not cause relative displacements between stone blocks.

Base acceleration in form of seismic excitation
A detailed presentation of the stability analysis conducted for lintels exposed to seismic excitation is provided below.To this end, the values of horizontal ground acceleration recorded during earthquakes in Petrovac (Montenegro), Erzincan (Turkey), and Selsund (South Iceland) in 1979, 1992 and 2000, respectively, were chosen.The records of the mentioned earthquakes were taken from the European Strong Motion Database, as shown in Figure 21.Since differences in the load bearing capacity of the analysed lintels were not clearly distinguishable in the real data found in the accelerograms, their amplitudes were gradually increased until structural failure.Due to a great number of results, only the response of the wall C is given in Figure 22.This wall was exposed to the values registered in the accelerogram of the Erzincan earthquake, scaled to the peak acceleration (a max ) of 0.6 g at various points in time.It can be concluded from these results that the duration of peak base accelerations higher than the minimum acceleration causing structural failure (Table 1) is too short for the part Stability analysis of dry stone lintels using combined finite-discrete element method of wall to overturn.However, alternating peaks in the base acceleration amplitude in both positive and negative direction lead to relative horizontal displacements of blocks, which is especially noticeable in the top part of the structure where the arch lintel transfers horizontal components of the reaction to the surrounding parts of the wall.
The structural response of the wall C to the same earthquake scaled to the peak acceleration of a max = 1.2 g is presented in Figure 23 for various points in time.The results show that the places where the cracking between blocks occurs during the seismic excitation actually correspond to the failure mechanism occurring during the linear base acceleration.Furthermore, it can be seen that the failure mechanism of the analysed wall was caused by a relative displacement of stone blocks during the seismic excitation.In other words, the structure collapses at the moment when the horizontal displacement of blocks becomes so large that the stone lintel loses its support.Furthermore, the presented results also reveal particularities of dry stone masonry walls when compared to monolith structures such as the reinforced concrete ones.More specifically, in the case of stone structures, any kind of seismic excitation having maximum peak acceleration greater than the minimum acceleration needed for a part of a structure to overturn, leads to relative displacement of stone blocks.A stone structure exposed to a series of earthquakes of relatively low amplitudes, where each subsequent tremor contributes to the displacement of stone blocks, would also collapse at a certain point.The determination of seismic resistance of stone structures on the basis of more than one earthquake, when the structure is always considered in its non-deformed configuration, can be misleading and deliver results pointing to a much greater seismic resistance than can realistically be expected.
The condition of all seven walls after their exposure to the Petrovac, Erzincan, and Selsund earthquakes is shown in Figures 24 to 26.Their amplitudes were increased up to a maximum peak acceleration of 1.20 g.The results presented point to the fact that the lintels capable of resisting larger horizontal displacements are more resistant to seismic action.
In the conducted analyses, these are the lintels type A, E, and F. In addition, one can observe that the lintel type G, which Emili Zubović, Hrvoje Smoljanović, Boris Trogrlić is composed of a greater number of smaller stone blocks, exhibits a rather poor seismic performance because even small horizontal displacements will lead to its collapse.Lower seismic resistance of the lintel type C when compared to the lintel type E can also be observed, although both may seem identical at first glance.The reason for this lies in the fact that the lintel type C has a much shallower arch so that, in addition to being more sensitive to horizontal displacement of supports, it also contributes to the detachment of the left and right supports due to horizontal reactions, which ultimately leads to structural failure.
The shape of stone blocks of an arch also plays an important role in its seismic performance.Comparing lintels D and E, where arch diameters equal the size of the opening they span, and where vertical reactions on the support are dominant, it can be observed that the lintel type E exhibits a higher seismic resistance.This can be attributed to the fact that the lintel E consists of longer stone blocks which, by virtue of their shape, prevent creation of hinges in the arch that lead to structural collapse, which is precisely the case with the lintel type D. Furthermore, it seems that the seismic resistance of flat lintels consisting of several stone blocks, such as in B, is poorer than that of arch lintels, given that the flat ones are more sensitive to horizontal displacements of stone blocks.Stability analysis of dry stone lintels using combined finite-discrete element method

Conclusion
The use of the combined finite-discrete element method in the analysis of structural stability of stone lintels is presented in this paper.This method is considered suitable for the analysis of the problem treated because it allows modelling of each stone block as an individual discrete element.Furthermore, it permits modelling of the interaction at the contact between the blocks, including dry friction and energy loss due to collision of two stone blocks.It also permits modelling deformability of stone blocks and, finally, it takes into consideration large displacements and rotations.In addition, this method enables accounting for the splitting of stone blocks due to a force that exceeds the strength of material in tension and shear, and this via contact elements implemented between the mesh of finite elements of each block.However, such a scenario was not considered in the analysis as the stresses in the structure were relatively low, and the calculation time was considerably reduced.Seven different types of stone lintels were exposed to the horizontal monotonically increasing base acceleration, as well as to base acceleration recorded during three real earthquakes, the amplitude of which was gradually increased until structural collapse.The analyses point to the way in which structural collapse takes place in these seven stone structures, to the nature of the failure mechanisms for each wall, and to the way in which different stone block patterns in lintels influence seismic performance of their respective structures.Ultimately, the type of stone lintel having the highest seismic resistance is discussed.As a measure relevant for seismic resistance of structures, the minimum amplitude of the horizontal base acceleration leading to structural collapse was adopted.The results of the numerical analyses conducted in the paper allow us to conclude that the studied structures have significant seismic resistance.The minimum peak acceleration causing relative displacements of stone blocks was determined based on the linearly increasing horizontal base acceleration.The failure mechanism obtained on the basis of the linearly increasing horizontal base acceleration corresponds to the pattern of cracking between the stone blocks exposed to seismic excitation.Furthermore, it can be seen that the resistance of walls to seismic action is by far greater than their resistance to the linearly increasing base acceleration.This is so because the peak acceleration taking place during the seismic excitation does not last long enough to cause a part of a structure to overturn.Instead, the peak acceleration acts alternately on both sides, which consequently causes the opening and closing of cracks between stone blocks.During every repeated cycle of this kind, permanent cracks appear between vertical dry joints and, finally, the collapse occurs when the horizontal displacement of blocks is so large that the lintel can no longer transfer the vertical component of the reaction.For that reason, the lintels that can withstand greater horizontal displacements have a higher seismic resistance.These are the lintels A, E, and F given in the analysis presented in the paper.It is suggested that the type G lintel, i.e. the lintel consisting of small stone blocks, should be avoided in design, as even small displacement of blocks can lead to structural collapse.As for the arch lintels, the horizontal reaction of the arch on the surrounding wall should be as low as possible, or should not be present at all, just like in the case of arches whose diameter is equal to the length of the opening.Such arches are more resistant to horizontal displacements of stone blocks.It is precisely for that reason that the type E lintel is more resistant than the type C lintel.The shape of stone blocks used in an arch also plays an important role in seismic resistance.The very shape of stone blocks should prevent creation of hinges in the arch, i.e. it should impede relative rotation between the stone blocks as that would lead to the collapse of the arch (the scenario with lintel D).This can be accomplished by using more elongated stone blocks, as in the case of the type E lintel.

Figure 1 .
Figure 1.Ancient stone structures in Pula: a) Triumphal Arch of the Sergii (Golden Gate); b) Twin Gates; c) Temple of Augustus; d) Pula Arena

Figure 2 .
Figure 2. Discretization of dry stone masonry structure

Figure 3 .
Figure 3. Initial and current configuration of triangular finite element

Figure 4 .
Figure 4. Contact differential force in the vicinity of points P m and P k

Figure 14 .Figure 20 .
Figure 14.Failure mechanism for type A lintel: a) positive base acceleration; b) negative base acceleration

Figure 22 .
Figure 22.Dry stone frame of type C during Erzincan earthquake scaled to peak acceleration of ag = 0.60 in time periods: a) t = 0.00 s; b) t = 1.94 s; c) t = 2.86 s; d) t = 4.07 s; e) t = 7.70 s; f) t = 18.49s

Figure 23 .
Figure 23.Type C frame during Erzincan earthquake scaled to peak acceleration of 1.20 g in time periods: a) t = 0.00 s; b) appearance of the first crack in the time t = 3.58 s; c) t = 3.90 s; d) t = 4.38 s; e) t = 4.68; f) t = 19.49s

Figure 24 .Figure 26 .
Figure 24.Condition after Petrovac earthquake scaled to maximum peak acceleration of 1.2 g for: a) wall A; b) wall B; c) wall C; d) wall D; e) wall E; f) wall F; g) wall G