Discrete Element Modeling of the Seismic Behavior of Masonry Construction

Discrete element models are a powerful tool for the analysis of masonry, given their ability to represent the discontinuous nature of these structures, and to simulate the most common deformation and failure modes. In particular, discrete elements allow the assessment of the seismic behavior of masonry construction, using either pushover analysis or time domain dynamic analysis. The fundamental concepts of discrete elements are concisely presented, stressing the issues related to masonry modeling. Methods for generation of block models are discussed, with some examples for the case of irregular stone masonry walls. A discrete element analysis of a shaking table test performed on a traditional stone masonry house is discussed, as a demonstration of the capabilities of these models. Practical application issues are examined, namely the computational requirements for dynamic analysis.


Introduction
The safety assessment of masonry structures under seismic loads demands numerical models capable of representing appropriately the response and the failure modes observed during earthquakes.Laboratory testing, namely using shaking tables, offers a controlled environment in which the main variables can be more accurately measured, providing important data for the validation of numerical models.A wide range of models is available nowadays, each one with its specific strengths and a preferred range of applications [1].Discrete element (DE) models, involving the representation of masonry by a system of interacting distinct blocks, are one of the numerical tools particularly suited for the simulation of phenomena such as sliding and separation along joints, which may induce progressive structural damage and collapse [2].
Many applications of DE models to masonry structures under seismic loading have been reported in the literature.Early works addressed mainly the study of historical monuments involving systems with a relatively small number of blocks [3].Presently, it is possible to address much more complex structures, composed of a large number of blocks, and assess their seismic capacity by means of either pushover methods [4] or dynamic analysis [5].DE models using these alternative analysis techniques, based on static or dynamic representations of the seismic action, have been applied in the simulation of lab experiments of masonry walls under out-of-plane loading [6,7].The failure modes observed in lab tests involving in-plane cyclic loading of masonry panels have also been effectively reproduced by these models [8].
In the literature, several discontinuum modeling techniques have been developed for masonry analysis, which share many of the fundamental concepts of DE models.For example, the rigid body and spring model [9] involves the representation of masonry by a system of rigid blocks, which do not necessarily correspond to the masonry units, as a homogenization procedure is invoked to obtain the contact properties.The applied element method [10,11] also concentrates the deformation at the interfaces between rigid blocks, which are discretized into a fine mesh of contact springs with nonlinear behavior.The fiber contact element method [12] is a related approach, involving the definition of contact properties from a system of fibers that connect adjacent blocks.Within a finite element framework, models based on rigid elements and nonlinear joints have also been advanced for masonry structures, namely resorting to limit analysis concepts [13,14].
In the present article, the fundamental concepts of discrete element models are concisely discussed in the next section, including computational efficiency issues.The generation of models to represent traditional stone wall structures is then addressed.Two different methods of creating irregular assemblies are described, one based on a Voronoi tessellation and the other on overlaying irregular courses.The responses of models generated with different random geometric parameters are compared.The issue of validation of the numerical models is examined, with particular reference to the results of a shaking table test of a physical model of a traditional stone masonry house, for which a simplified block representation was created.Finally, some concluding remarks are presented.

Block Representation
Discrete element models provide a discontinuous idealization of masonry (commonly designated as micro-modeling), in which the units and the joints are explicitly represented [1].There are multiple DE formulations and codes available nowadays, using different terminologies, but they all share common underlying concepts [2].The units are typically represented by polygonal or polyhedral shapes.The simplest models assume the blocks to be rigid bodies, with all the system deformation placed in the joints.This is a good assumption for hard stone blocks, but it also simplifies the analysis of large structures.Deformable blocks provide a more versatile framework.The block deformation is typically modelled by discretizing it internally into a finite element mesh.The element material is assumed elastic in most cases, but nonlinear constitutive models may also be assigned to the elements.For many problems in which the inter-block movements are dominant, the use of rigid blocks provides results practically equivalent to coarse-mesh deformable blocks.A comparative study for out-of-plane failure of walls is presented in [4].

Contact between Blocks
Most discrete element models adopt a simplified representation of the contact between blocks based on sets of point contacts.Therefore, no joint or interface elements are defined, contrary to finite element micro-models.The point contact approach is more versatile and facilitates the analysis in the large displacement range, as the system connectivity changes, and the type, location, and orientation of contacts need to be periodically updated.The drawback is that more contact points are required to achieve an accurate contact stress distribution.
The point contact approach assumes that the stress at a contact point is a function of the relative displacement of the two blocks at that point.The contact points are typically placed at vertex-face or edge-edge interactions.For deformable blocks, additional contact points are created for the nodal points of the internal element mesh that fall on block boundaries.Block faces are typically discretized into triangular facets, allowing contact areas to be assigned to each contact point, and therefore to relate contact stresses and forces.
An example of contact representation in the code 3DEC [15] is shown in Figure 1.Two types of elementary or point contacts may exist: vertex-to-face (VF) and edge-to-edge (EE).A contact point of VF type is created where a vertex touches a face, while a contact point of EE type is located where two edges intersect.Vertex-to-edge or vertex-to-vertex interactions are considered particular cases of the VF type.The common case of face-to-face interaction between two blocks is represented by several elementary contacts (Figure 1a).In this case, the sum of the areas of the point contacts equals the total area of contact between the blocks.On the other hand, when the blocks interact by two intersecting edges, a case of a true point contact, then a single elementary contact exists.For practical purposes, a true point interaction obviously has a finite stiffness, so a minimum contact area has to be defined, typically a small fraction of the average contact area.purposes, a true point interaction obviously has a finite stiffness, so a minimum contact area has to be defined, typically a small fraction of the average contact area.

Solution methods and computational efficiency
The need to represent large displacement problems, demanding a progressive update of block positions and contact point locations, has prompted many DE codes to avoid matrix solutions, opting for explicit solution methods.For time domain analysis, the integration of the equations of motion, of either the rigid blocks or the deformable block nodes, using a central-difference algorithm, is a common choice.A feature of many DE codes is to use the same algorithm for static problems by means of dynamic relaxation, which resorts to artificial damping to obtain convergence to the equilibrium solution, or to a failure mechanism.Adaptive damping and mass scaling techniques may be employed to improve the convergence rate of the relaxation procedure.
These explicit time stepping algorithms are quite effective for quasi-static analysis, namely in the case of pushover methods for seismic assessment.For dynamic analysis, however, since mass scaling cannot be applied, time steps are often small due to numerical stability requirements [9].These limitations may be particularly severe when the stiffness-proportional component of Rayleigh damping is used.Therefore, run times can be demanding for large or complex models.The use of rigid blocks often allows larger time steps, so they are frequently applied in dynamic analysis, as done in the study described in section 4.

Material properties
The calibration of numerical models with experimental data is a key issue addressed by many authors.For example, a consistent method has been developed to identify the parameters of a DE model of a brick panel from laboratory tests [16].More often, the calibration is based on empirical trial-and-error procedures until an acceptable fit of key indicators is reached.Parametric studies are invaluable to find the potential influence of material properties that display a significant uncertainty.
A model for dynamic analysis of existing structures under earthquake loading has to provide realistic natural frequencies.Given the uncertainty in assigning the deformability to many types of masonry, particularly for historical structures, in situ measurements by ambient vibration are important to calibrate the model.For a rigid block representation, the model deformability is given by the joint stiffness parameters, which can be estimated directly given the measured frequencies.For example, even for simple stone masonry structures, the estimation of the displacements caused by an earthquake has been shown to be much improved if the structural frequencies are correctly matched [17].

Modeling the masonry morphology
Devising a block system to represent modern brick masonry poses no difficulty, as only regular shaped blocks are needed, and the procedure can be easily automated.In the case of existing stone

Solution Methods and Computational Efficiency
The need to represent large displacement problems, demanding a progressive update of block positions and contact point locations, has prompted many DE codes to avoid matrix solutions, opting for explicit solution methods.For time domain analysis, the integration of the equations of motion, of either the rigid blocks or the deformable block nodes, using a central-difference algorithm, is a common choice.A feature of many DE codes is to use the same algorithm for static problems by means of dynamic relaxation, which resorts to artificial damping to obtain convergence to the equilibrium solution, or to a failure mechanism.Adaptive damping and mass scaling techniques may be employed to improve the convergence rate of the relaxation procedure.
These explicit time stepping algorithms are quite effective for quasi-static analysis, namely in the case of pushover methods for seismic assessment.For dynamic analysis, however, since mass scaling cannot be applied, time steps are often small due to numerical stability requirements [9].These limitations may be particularly severe when the stiffness-proportional component of Rayleigh damping is used.Therefore, run times can be demanding for large or complex models.The use of rigid blocks often allows larger time steps, so they are frequently applied in dynamic analysis, as done in the study described in Section 4.

Material Properties
The calibration of numerical models with experimental data is a key issue addressed by many authors.For example, a consistent method has been developed to identify the parameters of a DE model of a brick panel from laboratory tests [16].More often, the calibration is based on empirical trial-and-error procedures until an acceptable fit of key indicators is reached.Parametric studies are invaluable to find the potential influence of material properties that display a significant uncertainty.
A model for dynamic analysis of existing structures under earthquake loading has to provide realistic natural frequencies.Given the uncertainty in assigning the deformability to many types of masonry, particularly for historical structures, in situ measurements by ambient vibration are important to calibrate the model.For a rigid block representation, the model deformability is given by the joint stiffness parameters, which can be estimated directly given the measured frequencies.For example, even for simple stone masonry structures, the estimation of the displacements caused by an earthquake has been shown to be much improved if the structural frequencies are correctly matched [17].

Modeling the Masonry Morphology
Devising a block system to represent modern brick masonry poses no difficulty, as only regular shaped blocks are needed, and the procedure can be easily automated.In the case of existing stone masonry structures, the model generation tasks are much more time consuming.We should keep in mind that a numerical model is always an idealization of the real structure.In most cases, the analyst selects a simplified structure in such a way that the key features of the block structure are reproduced.Numerical experiments can be performed to evaluate the effect that a given simplification has on the results, and if this is acceptable given the purpose of the analysis.
In the case of monumental structures composed of large stone blocks, a photo survey may allow the definition of a numerical model that closely approximates the real morphology [18].However, in most DE models of stone masonry walls, an equivalent regular block system is employed, which respects the typical values of block sizes and the amount of interlocking.In the next two sections, the generation of block systems with irregular geometries is addressed for the case of masonry walls.In order to simplify the issue, only the non-regular pattern visible in the wall plane is modelled.In reality, wall cross-sections are often composed of two or three leaves of block and fill material.A 3D representation of such structures is only feasible in small models of masonry components.When 2D models are employed to study the out-of-plane stability of such walls, however, it is possible to achieve a much more detailed representation.De Felice [19] studied a three-leaf cross section of a traditional wall, simulating closely the observed irregular block geometries, including the inner leaf of rubble masonry.

Regular and Voronoi Block Patterns
In this study, the influence of the block pattern on the out-of-plane stability of stone masonry walls was analyzed with DE models.The regular block pattern is shown in Figure 2a [20].The wall dimensions were 20 × 10 m, with 0.80 m thickness.A Young's modulus of 2.5 GPa was assumed.As the blocks were assumed rigid, the joint stiffnesses were calculated to reproduce the global deformability: normal stiffness of 2.5 GPa/m, and shear stiffness of 1.0 GPa/m.For simplicity, a Coulomb friction model was adopted for the joints, with a friction angle of 35 • , and no cohesion or tensile strength.
masonry structures, the model generation tasks are much more time consuming.We should keep in mind that a numerical model is always an idealization of the real structure.In most cases, the analyst selects a simplified structure in such a way that the key features of the block structure are reproduced.Numerical experiments can be performed to evaluate the effect that a given simplification has on the results, and if this is acceptable given the purpose of the analysis.
In the case of monumental structures composed of large stone blocks, a photo survey may allow the definition of a numerical model that closely approximates the real morphology [18].However, in most DE models of stone masonry walls, an equivalent regular block system is employed, which respects the typical values of block sizes and the amount of interlocking.In the next two sections, the generation of block systems with irregular geometries is addressed for the case of masonry walls.In order to simplify the issue, only the non-regular pattern visible in the wall plane is modelled.In reality, wall cross-sections are often composed of two or three leaves of block and fill material.A 3D representation of such structures is only feasible in small models of masonry components.When 2D models are employed to study the out-of-plane stability of such walls, however, it is possible to achieve a much more detailed representation.De Felice [19] studied a three-leaf cross section of a traditional wall, simulating closely the observed irregular block geometries, including the inner leaf of rubble masonry.

Regular and Voronoi block patterns
In this study, the influence of the block pattern on the out-of-plane stability of stone masonry walls was analyzed with DE models.The regular block pattern is shown in Figure 2a [20].The wall dimensions were 20 x 10 m, with 0.80 m thickness.A Young's modulus of 2.5 GPa was assumed.As the blocks were assumed rigid, the joint stiffnesses were calculated to reproduce the global deformability: normal stiffness of 2.5 GPa/m, and shear stiffness of 1.0 GPa/m.For simplicity, a Coulomb friction model was adopted for the joints, with a friction angle of 35°, and no cohesion or tensile strength.
The rigid base block and the two vertical columns were fixed in all directions (Figure 2a).A distributed static load was applied to the wall in the out of plane direction, progressively increasing until failure ensued.Therefore, the wall was simply supported by the vertical column blocks, approximately simulating the effect of cross-walls, thus allowing the wall ends to rotate.
In addition to the regular system of Figure 2a, Voronoi polygon patterns were randomly created, with the same average block size as the regular mesh.One of these is shown in Figure 2b, at the stage in which the out-of-plane failure mode developed.Figure 3 compares the force-deformation curves and failure loads of the simulations with 3 different Voronoi block systems, with the same average size.It is interesting to observe that the three randomly generated patterns displayed fairly similar behaviors.In the same graphic, 2 cases of regular joints are shown.To be comparable with the Voronoi model, square blocks were used.The case of no offset had continuous vertical joints.The loss of interlocking led to a lower strength.An The rigid base block and the two vertical columns were fixed in all directions (Figure 2a).A distributed static load was applied to the wall in the out of plane direction, progressively increasing until failure ensued.Therefore, the wall was simply supported by the vertical column blocks, approximately simulating the effect of cross-walls, thus allowing the wall ends to rotate.
In addition to the regular system of Figure 2a, Voronoi polygon patterns were randomly created, with the same average block size as the regular mesh.One of these is shown in Figure 2b, at the stage in which the out-of-plane failure mode developed.
Figure 3 compares the force-deformation curves and failure loads of the simulations with 3 different Voronoi block systems, with the same average size.It is interesting to observe that the three randomly generated patterns displayed fairly similar behaviors.In the same graphic, 2 cases of regular joints are shown.To be comparable with the Voronoi model, square blocks were used.The case of no offset had continuous vertical joints.The loss of interlocking led to a lower strength.An offset of 0.1 m clearly increased the initial stiffness and failure load, which was only marginally lower than the 0.5 m offset case shown in Figure 2a.Cross-joint imbrication was critical, and the Voronoi pattern did not create discontinuous joints, so it tended to underestimate the block interlocking.

Irregular coursed block patterns
The Voronoi pattern does not represent the typical coursed structure of traditional stone masonry, being essentially an expeditious way to create randomly shaped block systems.An alternative procedure was proposed to generate block patterns which displayed overlaid courses, but with irregular horizontal joints, course height, and cross-joint imbrication [20], as shown in Figure 4.The procedure developed to obtain this block pattern involved the sequential generation of the masonry courses using the five geometric parameters shown in Figure 5 [20].Each one of these parameters was defined in statistical terms by a mean value (m) and a deviation (d).The geometry of the bed joints was composed of a series of segments according to the following parameters: spacing (sm, sd); segment length (tm, td); and vertical deviation (hm, hd) from the mean trace.In the case shown, for simplicity, a uniform distribution was assumed.For instance, the spacing of bed joints was given by a random number in the interval [sm -sd, sm + sd].The cross joint were introduced sequentially, with each location defined by a spacing parameter (bm, bd), and an angle deviation from the vertical (am, ad).The system in Figure 4

Irregular Coursed Block Patterns
The Voronoi pattern does not represent the typical coursed structure of traditional stone masonry, being essentially an expeditious way to create randomly shaped block systems.An alternative procedure was proposed to generate block patterns which displayed overlaid courses, but with irregular horizontal joints, course height, and cross-joint imbrication [20], as shown in Figure 4.

Irregular coursed block patterns
The Voronoi pattern does not represent the typical coursed structure of traditional stone masonry, being essentially an expeditious way to create randomly shaped block systems.An alternative procedure was proposed to generate block patterns which displayed overlaid courses, but with irregular horizontal joints, course height, and cross-joint imbrication [20], as shown in Figure 4.The procedure developed to obtain this block pattern involved the sequential generation of the masonry courses using the five geometric parameters shown in Figure 5 [20].Each one of these parameters was defined in statistical terms by a mean value (m) and a deviation (d).The geometry of the bed joints was composed of a series of segments according to the following parameters: spacing (sm, sd); segment length (tm, td); and vertical deviation (hm, hd) from the mean trace.In the case shown, for simplicity, a uniform distribution was assumed.For instance, the spacing of bed joints was given by a random number in the interval [sm -sd, sm + sd].The cross joint were introduced sequentially, with each location defined by a spacing parameter (bm, bd), and an angle deviation from the vertical (am, ad).The system in Figure 4   The procedure developed to obtain this block pattern involved the sequential generation of the masonry courses using the five geometric parameters shown in Figure 5 [20].Each one of these parameters was defined in statistical terms by a mean value (m) and a deviation (d).The geometry of the bed joints was composed of a series of segments according to the following parameters: spacing (s m , s d ); segment length (t m , t d ); and vertical deviation (h m , h d ) from the mean trace.In the case shown, for simplicity, a uniform distribution was assumed.For instance, the spacing of bed joints was given by a random number in the interval The system in Figure 4 was one of three randomly created, and their force-displacement curves are compared in Figure 6 with those from regular jointed models with rectangular blocks (Figure 2a).The 3 irregular patterns displayed failure loads in the same range, which were also close to the case of the regular pattern with the smallest offset (0.1 m), and well above the curve for the continuous cross joints.It should be noted that the regular block failure loads in this figure, with rectangular shapes of ratio 2:1, were higher than the ones in Figure 3, which corresponded to square blocks.

Experimental study
Shake table testing provides comprehensive data on the behavior of masonry structures under intense dynamic loading, which is invaluable for the validation of numerical codes.In the framework of the 9th International Masonry Conference, which took place in Guimarães, a comparative study of different numerical models was undertaken, involving blind predictions and a posteriori analyses of shake table tests on brick and stone masonry models performed at LNEC [21,22].The stone masonry structure, shown in Figure 7, consisted of a main gabled wall, with a door opening, and return walls on both ends.An opening was placed in one of the return walls, resulting in an asymmetry, inducing torsional movements.The walls, built of irregular stones, each had a thickness of 0.5 m.Unidirectional ground motion was applied in the out-of-plane direction of the main wall.Various finite element and discrete element models were employed in the comparative study [22].The main issues involved in the DE representation of the experiment, and the lesson learned, are discussed in The system in Figure 4 was one of three randomly created, and their force-displacement curves are compared in Figure 6 with those from regular jointed models with rectangular blocks (Figure 2a).The 3 irregular patterns displayed failure loads in the same range, which were also close to the case of the regular pattern with the smallest offset (0.1 m), and well above the curve for the continuous cross joints.It should be noted that the regular block failure loads in this figure, with rectangular shapes of ratio 2:1, were higher than the ones in Figure 3, which corresponded to square blocks.The system in Figure 4 was one of three randomly created, and their force-displacement curves are compared in Figure 6 with those from regular jointed models with rectangular blocks (Figure 2a).The 3 irregular patterns displayed failure loads in the same range, which were also close to the case of the regular pattern with the smallest offset (0.1 m), and well above the curve for the continuous cross joints.It should be noted that the regular block failure loads in this figure, with rectangular shapes of ratio 2:1, were higher than the ones in Figure 3, which corresponded to square blocks.

Experimental study
Shake table testing provides comprehensive data on the behavior of masonry structures under intense dynamic loading, which is invaluable for the validation of numerical codes.In the framework of the 9th International Masonry Conference, which took place in Guimarães, a comparative study of different numerical models was undertaken, involving blind predictions and a posteriori analyses of shake table tests on brick and stone masonry models performed at LNEC [21,22].The stone masonry structure, shown in Figure 7, consisted of a main gabled wall, with a door opening, and return walls on both ends.An opening was placed in one of the return walls, resulting in an asymmetry, inducing torsional movements.The walls, built of irregular stones, each had a thickness of 0.5 m.Unidirectional ground motion was applied in the out-of-plane direction of the main wall.Various finite element and discrete element models were employed in the comparative study [22].The main issues involved in the DE representation of the experiment, and the lesson learned, are discussed in

Experimental study
Shake table testing provides comprehensive data on the behavior of masonry structures under intense dynamic loading, which is invaluable for the validation of numerical codes.In the framework of the 9th International Masonry Conference, which took place in Guimarães, a comparative study of different numerical models was undertaken, involving blind predictions and a posteriori analyses of shake table tests on brick and stone masonry models performed at LNEC [21,22].The stone masonry structure, shown in Figure 7, consisted of a main gabled wall, with a door opening, and return walls on both ends.An opening was placed in one of the return walls, resulting in an asymmetry, inducing torsional movements.The walls, built of irregular stones, each had a thickness of 0.5 m.Unidirectional ground motion was applied in the out-of-plane direction of the main wall.Various finite element and discrete element models were employed in the comparative study [22].The main issues involved in the DE representation of the experiment, and the lesson learned, are discussed in the following section, with particular reference to the 3DEC rigid block analyses [23].Two other types of DE models were also included in the study: a combined finite-discrete element approach [24]; and a macro-element formulation [25].
of DE models were also included in the study: a combined finite-discrete element approach [24]; and a macro-element formulation [25].

DE representation of the stone masonry test model
DE allows various levels of detail to be included in the representation of a masonry structure.It would be possible to devise a model in which each real stone would be simulated, adopting a more or less complex geometry.However, in the analyses of the stone masonry specimen performed with 3DEC, a simplified representation was chosen [23].The aim was to follow a model generation methodology applicable to the practical study of a real structure, where the real block shapes are known.Only the main pattern of the stone arrangement was attempted, respecting the average unit dimensions, but assuming horizontal joints, and vertical joint offsets of the order of the real ones, as shown in Figure 8. Furthermore, a single block was placed across the wall thickness, unlike the twoleaf walls of the test structure.The model had about 200 brick-shaped blocks, which were assumed to behave as rigid bodies, each with 6 degrees-of-freedom.

Material properties
Consistent with the aim of applying the simplest modeling assumptions, the joints were assigned the standard Mohr-Coulomb model with brittle behavior available in 3DEC [15].Joint deformability is characterized by normal and shear stiffnesses, and joint strength by the friction angle, cohesion, and tensile strength.If either tensile or shear failure take place, the cohesion and tensile strengths are both set to zero, while the friction angle is unchanged.
In a rigid block model, all the deformation is concentrated at the joints, so the joint stiffness parameters have to be selected to reproduce the masonry elastic moduli.For the blind prediction before the test, a Young's modulus of 2077 MPa given by wallette tests was the available information.The estimation of the normal and shear joint stiffnesses was based on an average spacing of 0.3 m for

DE Representation of the Stone Masonry Test Model
DE allows various levels of detail to be included in the representation of a masonry structure.It would be possible to devise a model in which each real stone would be simulated, adopting a more or less complex geometry.However, in the analyses of the stone masonry specimen performed with 3DEC, a simplified representation was chosen [23].The aim was to follow a model generation methodology applicable to the practical study of a real structure, where the real block shapes are known.Only the main pattern of the stone arrangement was attempted, respecting the average unit dimensions, but assuming horizontal joints, and vertical joint offsets of the order of the real ones, as shown in Figure 8. Furthermore, a single block was placed across the wall thickness, unlike the two-leaf walls of the test structure.The model had about 200 brick-shaped blocks, which were assumed to behave as rigid bodies, each with 6 degrees-of-freedom.the following section, with particular reference to the 3DEC rigid block analyses [23].Two other types of DE models were also included in the study: a combined finite-discrete element approach [24]; and a macro-element formulation [25].

DE representation of the stone masonry test model
DE allows various levels of detail to be included in the representation of a masonry structure.It would be possible to devise a model in which each real stone would be simulated, adopting a more or less complex geometry.However, in the analyses of the stone masonry specimen performed with 3DEC, a simplified representation was chosen [23].The aim was to follow a model generation methodology applicable to the practical study of a real structure, where the real block shapes are known.Only the main pattern of the stone arrangement was attempted, respecting the average unit dimensions, but assuming horizontal joints, and vertical joint offsets of the order of the real ones, as shown in Figure 8. Furthermore, a single block was placed across the wall thickness, unlike the twoleaf walls of the test structure.The model had about 200 brick-shaped blocks, which were assumed to behave as rigid bodies, each with 6 degrees-of-freedom.

Material properties
Consistent with the aim of applying the simplest modeling assumptions, the joints were assigned the standard Mohr-Coulomb model with brittle behavior available in 3DEC [15].Joint deformability is characterized by normal and shear stiffnesses, and joint strength by the friction angle, cohesion, and tensile strength.If either tensile or shear failure take place, the cohesion and tensile strengths are both set to zero, while the friction angle is unchanged.
In a rigid block model, all the deformation is concentrated at the joints, so the joint stiffness parameters have to be selected to reproduce the masonry elastic moduli.For the blind prediction before the test, a Young's modulus of 2077 MPa given by wallette tests was the available information.The estimation of the normal and shear joint stiffnesses was based on an average spacing of 0.3 m for

Material Properties
Consistent with the aim of applying the simplest modeling assumptions, the joints were assigned the standard Mohr-Coulomb model with brittle behavior available in 3DEC [15].Joint deformability is characterized by normal and shear stiffnesses, and joint strength by the friction angle, cohesion, and tensile strength.If either tensile or shear failure take place, the cohesion and tensile strengths are both set to zero, while the friction angle is unchanged.
In a rigid block model, all the deformation is concentrated at the joints, so the joint stiffness parameters have to be selected to reproduce the masonry elastic moduli.For the blind prediction before the test, a Young's modulus of 2077 MPa given by wallette tests was the available information.The estimation of the normal and shear joint stiffnesses was based on an average spacing of 0.3 m for horizontal joints and 0.5 m for the vertical ones [23].However, for the post-diction analyses, the natural frequencies measured experimentally were provided.In order to match these, the joint stiffnesses had to be reduced by a factor of 2, implying that the tested masonry walls were softer than the wallettes.
The initial joint strength parameters in the prediction runs were based on the diagonal compression of the wallettes.Assuming a friction angle of 30 • , a cohesion of 0.32 MPa was estimated, and the tensile strength was taken as half of the cohesion.In the post-test analyses, the same cohesive and tensile strengths were used, but the friction angle was increased to 35 • , in order to obtain a better fit of the experimental response.The higher friction possibly accounted for the joint irregularity joints in the physical model.For the base joint, between model and shaking table, where no slip occurred, a friction angle of 45 • was adopted.

Pushover and Dynamic Analyses
Prior to the dynamic runs, a quasi-static analysis was performed with the DE model to obtain the failure load in pushover-type tests.First, the gravity load was applied, with the base block fixed in all directions.To perform the pushover analysis, a static horizontal load, proportional to the block weight, was applied to all the blocks, in the out-of-plane direction of the façade.The horizontal load was increased in increments corresponding to a mass force of 0.05 g, until collapse.The failure mode, shown in Figure 9a, agreed quite well with the deformation pattern observed experimentally, involving the out-of-plane rotation of the façade wall, while the side wall with the opening broke into a four component mechanism (Figure 7c).The pushover collapse load was 0.65 g, well below the 1.05 g record that led the model to near-collapse in the shake table tests.In any case, the rather simplified geometry of the DE rigid block structure appeared to be able to represent the key features of the experimental behavior.
Buildings 2019, 9, 43 8 of 11 natural frequencies measured experimentally were provided.In order to match these, the joint stiffnesses had to be reduced by a factor of 2, implying that the tested masonry walls were softer than the wallettes.The initial joint strength parameters in the prediction runs were based on the diagonal compression of the wallettes.Assuming a friction angle of 30°, a cohesion of 0.32 MPa was estimated, and the tensile strength was taken as half of the cohesion.In the post-test analyses, the same cohesive and tensile strengths were used, but the friction angle was increased to 35°, in order to obtain a better fit of the experimental response.The higher friction possibly accounted for the joint irregularity joints in the physical model.For the base joint, between model and shaking table, where no slip occurred, a friction angle of 45° was adopted.

Pushover and dynamic analyses
Prior to the dynamic runs, a quasi-static analysis was performed with the DE model to obtain the failure load in pushover-type tests.First, the gravity load was applied, with the base block fixed in all directions.To perform the pushover analysis, a static horizontal load, proportional to the block weight, was applied to all the blocks, in the out-of-plane direction of the façade.The horizontal load was increased in increments corresponding to a mass force of 0.05 g, until collapse.The failure mode, shown in Figure 9a, agreed quite well with the deformation pattern observed experimentally, involving the out-of-plane rotation of the façade wall, while the side wall with the opening broke into a four component mechanism (Figure 7c).The pushover collapse load was 0.65 g, well below the 1.05 g record that led the model to near-collapse in the shake table tests.In any case, the rather simplified geometry of the DE rigid block structure appeared to be able to represent the key features of the experimental behavior.For the dynamic runs, the model was first brought to equilibrium under the gravity load.Then, a time domain dynamic analysis was performed.The model was loaded by applying at the base block the dynamic records measured in the shake table during the tests.Rayleigh damping was employed, including both the mass-and stiffness-proportional components, with a value of 2% of critical damping at the fundamental frequency of 10 Hz.This value was selected to improve the simulation of the final test with the strongest motion, which caused widespread damage.Thus, in the early test stages, with lower input, the numerical model was expected to display higher response peaks than the experiment.It should be noted that in the DE model, further energy dissipation took place by friction on the sliding joints.The typical time step in the dynamic runs, determined by the code for For the dynamic runs, the model was first brought to equilibrium under the gravity load.Then, a time domain dynamic analysis was performed.The model was loaded by applying at the base block the dynamic records measured in the shake table during the tests.Rayleigh damping was employed, including both the mass-and stiffness-proportional components, with a value of 2% of critical damping at the fundamental frequency of 10 Hz.This value was selected to improve the simulation of the final test with the strongest motion, which caused widespread damage.Thus, in the early test stages, with lower input, the numerical model was expected to display higher response peaks than the experiment.It should be noted that in the DE model, further energy dissipation took place by friction on the sliding joints.The typical time step in the dynamic runs, determined by the code for reasons of numerical stability, was in the order of 6 × 10 -6 s.The use of the stiffness-proportional component of Rayleigh damping was mostly responsible for such a small value of the time step.
The model configuration after the last dynamic stage, with a peak acceleration of 1.05 g, is depicted in Figure 9b, showing a clear resemblance with the observed behavior (Figure 7c).The peak displacements at the monitored points of the structure were also generally in reasonable agreement with the experiment [23].Naturally, we cannot expect to match the details of the time response curves.
In conclusion, the reported analysis was based on a rather simplified block model.It is certainly possible to improve the geometrical resemblance of the 3DEC representation, namely including the two-leaf wall structure, or the horizontal joint non-planarity, as discussed above.In addition, the simple joint constitutive model adopted does not allow the representation of the progressive accumulation of damage in the successive tests.It should be noted, however, that data to support those refinements is not usually available.In addition, it was found that the material deformability provided by the static tests of elementary block assemblages did not yield the natural frequencies measured experimentally.Therefore, the calibration of the numerical model with the actual structural frequencies observed was an essential step to improve its performance in the dynamic runs.

Discussion
Discrete element models are a flexible analysis tool for discontinuous media that can be exploited in many different ways.During model generation, various levels of detail can be introduced, so that a better resemblance to the real block structure can be progressively achieved by introducing finer detailing.The new technologies that record accurately the wall morphology make it possible nowadays.When these are not available, generation methods, such as discussed above, can be used to create numerical models that match the main patterns and features of the particular type of masonry under study.The results of the research project reported in the previous section indicate that a relatively simple block pattern, which did not attempt to reproduce each individual block, was able to provide a fairly good agreement with the experimental data.More data is required for a refined model, and higher computational costs are to be expected.Therefore, the need to simplify is obvious.More research effort needs to be devoted to evaluating the influence of model simplification on the numerical results, providing guidance for engineering practice.
Pushover analysis provides a powerful tool for seismic assessment.In DE models, the physically possible failure mechanisms can be viewed and evaluated.The effect of material properties can be checked by parametric studies, or by using a statistical variation of parameters, with reasonable computational costs.Time domain dynamic analysis is much more computationally intensive.The explicit algorithms employed in most DE codes allow a rigorous simulation of the evolving geometry and contact conditions, but typically require very small time-steps.A key requirement for a meaningful dynamic analysis of existing structures is a good simulation of the natural frequencies, as several studies have shown [17,23].In situ measurements appear, therefore, highly desirable to guarantee that the numerical model matches the dynamic characteristics of the real structure.
Shake table tests of masonry typically display a significant variability in the results.This has been observed even for rocking of a single stone block, for which successive tests with the same earthquake record sometimes varied significantly [26].For a column of marble drums, a relatively wide range of responses was also obtained [27].Large amplitude block rocking is a phenomenon known to be sensitive to the initial conditions and loading.The change in local contact conditions between blocks in successive tests, whenever slip takes place, can account for part of these variations.DE models with slightly different system geometry, properties, or input loads, have also been shown to lead to a clear dispersion of the results [3].From a practical point of view, this sensitivity implies that a single dynamic analysis is not sufficient, and thus multiple runs, varying the main input parameters of the model, are always necessary to reach reliable conclusions.

Concluding Remarks
Discrete element models are an important tool in the analysis of masonry, given their ability to reproduce the observed patterns of behavior, particularly failure modes defined by the block structure.Analysis of masonry structures under intense earthquake loading has been one of the main areas of application of these models.For a successful application in engineering practice, any numerical tool needs to prove accurate and robust.Validation studies involving comparisons with shake table tests or the observed response during earthquakes are fundamental to provide confidence in a numerical method for seismic assessment studies.Complex models, however, often require data that is not readily available for existing masonry constructions.Model simplification strategies, whether in the representation of the block geometry or in the constitutive assumptions, are a key to providing meaningful results with the existing data in practical situations.The improvement of modeling procedures remains a major research goal.

Figure 1 .
Figure 1.Representation of block interactions by elementary point contacts of vertex-to-face type (VF) and edge-to-edge type (EE): (a) Face-to-face block contact; (b) Edge-to-edge block contact.

Figure 1 .
Figure 1.Representation of block interactions by elementary point contacts of vertex-to-face type (VF) and edge-to-edge type (EE): (a) Face-to-face block contact; (b) Edge-to-edge block contact.

Figure 2 .
Figure 2. Masonry wall models.(a) Regular block pattern; (b) Voronoi block pattern, showing failure mode for out-of-plane loading.

Figure 2 .
Figure 2. Masonry wall models.(a) Regular block pattern; (b) Voronoi block pattern, showing failure mode for out-of-plane loading.
Buildings 2019, 9, 43 5 of 11 offset of 0.1 m clearly increased the initial stiffness and failure load, which was only marginally lower than the 0.5 m offset case shown in Figure2a.Cross-joint imbrication was critical, and the Voronoi pattern did not create discontinuous joints, so it tended to underestimate the block interlocking.

Figure 3 .
Figure 3. Failure load for wall models: regular model of square blocks with cross-joint offsets of 0 and 0.1 m, and three random Voronoi patterns.

Figure 4 .
Figure 4. Masonry wall model with irregular coursed structure.

Figure 3 .
Figure 3. Failure load for wall models: regular model of square blocks with cross-joint offsets of 0 and 0.1 m, and three random Voronoi patterns.
clearly increased the initial stiffness and failure load, which was only marginally lower than the 0.5 m offset case shown in Figure2a.Cross-joint imbrication was critical, and the Voronoi pattern did not create discontinuous joints, so it tended to underestimate the block interlocking.

Figure 3 .
Figure 3. Failure load for wall models: regular model of square blocks with cross-joint offsets of 0 and 0.1 m, and three random Voronoi patterns.

Figure 4 .
Figure 4. Masonry wall model with irregular coursed structure.

Figure 4 .
Figure 4. Masonry wall model with irregular coursed structure.

11 Figure 5 .
Figure 5. Definition of geometric parameters in block generation procedure.

Figure 6 .
Figure 6.Failure load for wall models: regular model of rectangular blocks with cross-joint offsets of 0 and 0.1 m, and three irregular coursed patterns.

Figure 5 .
Figure 5. Definition of geometric parameters in block generation procedure.

Figure 6 .
Figure 6.Failure load for wall models: regular model of rectangular blocks with cross-joint offsets of 0 and 0.1 m, and three irregular coursed patterns.

Figure 6 .
Figure 6.Failure load for wall models: regular model of rectangular blocks with cross-joint offsets of 0 and 0.1 m, and three irregular coursed patterns.

Figure 7 .
Figure 7. Stone structure: (a) general view and (b) return wall with opening; (c) damage of the return wall after the test (adapted from [23]).

Figure 8 .
Figure 8. Simplified 3DEC model of the stone structure.

Figure 7 .
Figure 7. Stone structure: (a) general view and (b) return wall with opening; (c) damage of the return wall after the test (adapted from [23]).

Figure 7 .
Figure 7. Stone structure: (a) general view and (b) return wall with opening; (c) damage of the return wall after the test (adapted from [23]).

Figure 8 .
Figure 8. Simplified 3DEC model of the stone structure.

Figure 8 .
Figure 8. Simplified 3DEC model of the stone structure.

Figure 9 .
Figure 9. Numerical model configuration: (a) pushover failure mechanism; (b) final configuration at end of dynamic analysis.

Figure 9 .
Figure 9. Numerical model configuration: (a) pushover failure mechanism; (b) final configuration at end of dynamic analysis.