Evaluating the Seismic Capacity of Dry ‐ Joint Masonry Arch Structures via the Combined Finite ‐ Discrete Element Method

: The behaviour of dry ‐ joint masonry arch structures is highly nonlinear and discontinuous since they are composed of individual discrete blocks. These structures are vulnerable to seismic excitations. It is difficult for traditional methods like the standard finite element method (FEM) to simulate masonry failure due to their intrinsic limitations. An advanced computational approach, i.e., the combined finite ‐ discrete element method (FDEM), was employed in this study to examine the first ‐ order seismic capacity of masonry arches and buttressed arches with different shapes sub ‐ jected to gravity and constant horizontal acceleration. Within the framework of the FDEM, masonry blocks are discretised into discrete elements. A finite element formulation is implemented into each discrete element, providing accurate predictions of the deformation of each block and contact inter ‐ actions between blocks. Numerical examples are presented and validated with results from the ex ‐ isting literature, demonstrating that the FDEM is capable of capturing the seismic capacities and hinge locations of masonry arch structures. Further simulations on geometric parameters and fric ‐ tion coefficient of masonry buttressed arches were conducted, and their influences on the seismic capacities are revealed.


Introduction
Masonry arches have been the traditional forms of architectural heritages and historical constructions for a long time. They are assemblages of voussoirs made of stones or bricks, which are either bonded with low-strength mortar or completely mortarless. In this regard, a dry-joint assumption can be reasonably held. Moreover, masonry arch structures are vulnerable to seismic excitations. Due to the intrinsic discontinuities, their responses under seismic excitations are highly nonlinear. Thus, research into the seismic capacity of dry-joint masonry arch structures is required.
If fracture and sliding are not considered, the dynamic failure of masonry arches can be determined by the hinge failure mechanism following Heyman's hypothesis [1] that masonry voussoirs have zero tensile strength and infinite compressive strength, as well as large enough friction between them. Some pioneering work has been proposed based on analytical limit analysis. Oppenheim [2] investigated the dynamic response of an unreinforced masonry arch by a four-link mechanism, in which the arch was modelled as three rigid bars with predefined hinges. Subsequently, the same approach was employed to study the performance of circular stone arches subjected to base excitations [3,4]. Baratta et al. [5] examined the seismic responses of a simple masonry arch portal subjected to order seismic capacity of masonry arch structures under the combination of gravity and constant horizontal acceleration. Under such loading, the failure of masonry arch structures is a dynamic process. Therefore, the FDEM is an appropriate simulation tool, and the 2D FDEM program 'Y' [45] is employed in this study. The layout for the rest of the paper is as follows. Fundamentals on the FDEM modelling are addressed in Section 2. Numerical examples, including a variety of arch structures subjected to both gravitational and constant horizontal accelerations, are presented in Section 3. A parametric investigation regarding geometric parameters and friction coefficient are conducted in Section 4. Finally, the conclusions are reached in Section 5.

Voussoir Discretisation
Dry-joint masonry arches are composed of voussoirs. Within the framework of the FDEM, each voussoir is an individual discrete block and is meshed with numbers of threenode constant strain triangular (CST) elements. The discretisation of two adjacent voussoirs, A and B, are illustrated in Figure 1. In this study, the fracture behaviour of voussoir is not included. Hence, no element interface (i.e., joint element) needs to be defined between adjacent elements within a single block, implying voussoirs A and B cannot be further divided. Furthermore, the dry-joint assumption suggests that only contact and friction need to be considered between voussoirs A and B.
In this study, a No Binary Search (NBS) contact detection algorithm, called Munjiza-NBS contact detection algorithm [28], was implemented in the FDEM program 'Y'. One of the advantages of using this algorithm is that the total computational time to detect all the contact couples is proportional to the number of discrete elements, and further details can be found in Munjiza and Andrews [25]. Another important aspect of the FDEM is the contact interaction law defining the contact force, and it will be introduced in Section 2.3.

Element Motions
Referring to Munjiza [28], the translational and the rotational motions of an arbitrary discrete element i are formulated according to Newton's second law of motion, as where mi is the mass of discrete element i; a is the acceleration acts on element i with the exclusion of the gravitational acceleration; g is the gravitational acceleration; Ji is the moment of inertia about the centre of element i; ωi is the angular velocity about the element centre; and Fi and Ti are the resultant force and moment act on and about the centre of element i, respectively. According to Equations (1) and (2), the velocity and position of discrete element i can be explicitly obtained at each time step. An explicit central difference time integration scheme is used [28], and the nodal velocity of the next time step is given by Equation (3), as where a is the acceleration of the node, Δt is the time step, v is the velocity, and the subscripts 'next' and 'current' correspond to the next and the current time step.

Contact Forces
In the FDEM, the contact force f was evaluated and integrated through the overlapping area A. As shown in Figure 2, the infinitesimal contact force df due to the penetration of elemental area dA is defined in Equation (4), as df = df t + df c (4) in which dft and dfc are the force components of df from the target and the contactor, respectively, as where Pc and Pt are the points on the contactor and the target, respectively, and they share the same coordinate within A; φc and φt are predefined potentials; grad is the gradient; and Ep is the contact penalty and is usually 100 times of the Lame's first constant. Thus, the contact force f can be obtained by integrating dA over A, as Further details on the contact forces can be referred to Munjiza [28].

Numerical Examples
In this section, masonry arch structures subjected to the gravitational acceleration g and a constant horizontal acceleration λg were considered. λ is a scale factor, and the maximum λ that masonry arch structures can withstand is denoted as λ*, which indicates their first-order seismic capacity. The 2D analysis was conducted since it is acceptable, accurate enough and time-saving [46].

Arches without Buttresses
Dimitri and Tornabene [17] proposed an analytical model to evaluate the seismic capacity of masonry arches and buttressed arches based on the static theory of limit analysis, and their predictions of the model were also compared and verified with results from the DEM program UDEC [17]. In this section, two types of arches, i.e., the circular arch and the basket-handle arch, were investigated. Their geometric configurations are shown in Figure 3. The circular arch (Figure 3a) is simply characterised by the embrace angle β, the thickness t and the centreline radius Rc. Figure 3b shows the geometry of a basket-handle arch, which is composed of three circular arches: a larger arch with the centre O1, two  The seismic capacities of two circular arches and two basket-handle arches with the embrace angles β = 150° and β = 180° were investigated. Each arch is composed of 11 voussoirs. The geometric parameters are: (i) circular arches, Rc = 1.0 m and t/Rc = 0.2; (ii) baskethandle arches, Rc = 1.0 m, t/Rc = 0.2, d/Rc = 0.5 and r/Rc = 0.5. The geometric parameters were chosen such that direct comparison could be made with results in the literature [17]. The material properties are shown in Table 1. Mesh configurations of the arches used in the FDEM simulations are illustrated in Figure 4, in which the arches are yellow, the bases are blue and N is the total number of elements. The time step of the simulation was set to 1.0 × 10 −6 s, and a friction coefficient of 0.6 was defined.

Young's Modulus (GPa)
Poisson's Ratio Density (kg/m 3 ) The λ* was tested using every single value. For each λ, the simulation was performed for a period of time that was long enough so that the arch would collapse if the assigned λg was larger than its seismic capacity. Should the arch be still standing with little noticeable deformation at the end of the simulation, the structure was considered stable. Otherwise, it was unstable with the corresponding λg. Besides the current approach, a criterion based on the unbalanced force ratio/max velocity could have been a better way and shall be considered in future research. The explicit time integration in the FDEM was automatically terminated once the computed steps reached the predefined maximum step, or it could be terminated manually as long as the arch had collapsed under certain λg. According to Munjiza [28], viscous damping with a value of 2Δℎ was considered, where Δh is the characteristic size of the smallest element, E is the Young's modulus and ρ is the material density.
For λ > λ*, all of the four arches failed following the hinge mechanism. Take the instance of the circular arch with β = 180°; the entire failure process with λ = 0.3 was simulated by the FDEM and is presented in Figure 5. Initially, there were five hinges (see Figure  5a). Subsequently, one hinge closed, and the motion of the arch was dominated by the remaining four hinges until collapse, as shown in Figure 5b-d. Finally, the arch lost its stability and collapsed due to significant deformation. The dominant four hinges are denoted as hinge I, II, III and IV in Figure 5d. The seismic capacity factor λ* of the investigated arches are presented in Table 2, in which e is the relative difference between the simulation results and the results based on limit analysis. Clearly, the results from the FDEM simulation agree very well with the results in Dimitri and Tornabene [17], especially with the analytical results based on limit analysis. However, additional results on the circular and basket-handle arches with β = 120° show that the relative errors are higher, suggesting the limitation of the proposed model. The failure modes of circular and basket-handle arches are shown in Figures 6 and  7, respectively. It is observed that the failure modes predicted from the FDEM simulations are in excellent agreement with the results predicted by the UDEC [17], especially for the hinges and their locations.

Buttressed Arches
In this section, two types of buttressed arches are investigated, as shown in Figure 8. The figure is drawn based on Dimitri and Tornabene [17]. The width of the buttress is B, the height from O to the base is H1, the height of the buttress is H, the height from the top of the buttress to the arch crown is h and other parameters are similar to the arches shown in Figure 3.  Table 1. The time step of the simulation was set to 1.0 × 10 −6 s, and a large friction coefficient of 0.9 was adopted to avoid potential sliding. Similarly, failure of all buttressed arches followed the four-hinge mechanism, too. Herein, the buttressed baskethandle arch with β = 180° was chosen and presented representatively, and the failure process with λ = 0.18 is shown in Figure 10. At the first instance (Figure 10a), four hinges initiated quickly. Three upper hinges (one extrados hinge and two intrados hinges) appeared within the arch, and a lower hinge formed at the bottom of the right buttress; therefore, the entire structure behaved as a linked mechanism. The hinges continued to develop till the collapse of the structure.
Basket-handle The seismic capacity factors λ* of the buttressed basket-handle arches are presented in Table 3, including the results from Dimitri and Tornabene [17]. Apparently, the FDEM simulation results agree well with the results based on the limit analysis and the UDEC simulation. The failure modes of buttressed basket-handle arches are illustrated in Figures  11 and 12. All the structures failed due to the instability of the four-hinge mechanism, and the simulation results by the FDEM are in excellent agreement with those by UDEC, especially the hinges and their locations. It is worth mentioning that the locations of hinges are not as same as the corresponding single arches in Section 3.1, i.e., one hinge formed at the bottom of the right buttress. With the ground acceleration at the base, the overturning moments of the right buttress induced by the reaction of the right arch springer and the inertial force of itself are all clockwise. When they surpassed the anticlockwise moment of self-weight, the lower hinge appeared at the exterior corner of the bottom of the buttress.  A sensitivity analysis regarding the influences of Young's modulus E on the seismic capacity factor λ* was performed, and the results are given in Figure 13. Buttressed and unbuttressed circular arches with different embrace angles were considered. It is shown that when E is small, i.e., E < 5 GPa, the corresponding λ* is also slightly small. Once E ≥ 10 GPa, λ* reaches a plateau and is kept constant afterwards. According to Figure 13, the conclusion can be reached that the influence of block stiffness on the seismic capacity is limited unless Young's modulus E is quite small, e.g., E < 0.1 GPa. Since Young's modulus of most ordinary masonry blocks is around 40 GPa (a lot far from the 'soft' region), a rigid block assumption is an entirely relevant and justified assumption.

A Multi-Span Pointed Arch Church
Nodargi and Bisegna [47] investigated the seismic capacity of a multi-span pointed arch church, as shown in Figure 14a, using the approach of limit analysis. It is a three-span buttressed pointed masonry arch structure. Similar to the arch structures in Sections 3.1 and 3.2, it was subjected to the gravitational acceleration g and a constant horizontal acceleration λg. Various potential failure sections were studied by Nodargi and Bisegna [47], e.g., 'horizontal' stereotomy in the buttresses and 'polar' stereotomy in the arches, as shown in Figure 14a. Accordingly, the FDEM discretisation is illustrated in Figure 14b, denoted by Scheme A. The buttresses were discretised by horizontal blocks with a height of 0.26 m. The pointed arches were discretised by polar voussoirs with an embrace angle of β = 2°, and a special block was defined at the crown of each arch because of its pointed shape. As a comparison, another discretisation scheme with monolithic buttresses and the same polar voussoirs was also modelled, as shown in Figure 14c, denoted by Scheme B. Similar to the arch structures in Sections 3.1 and 3.2, only the dry-joint contact and friction between adjacent blocks or voussoirs were considered in the FDEM simulations. The total number of elements is 786 for Scheme A, and 342 for Scheme B. Material properties are as same as in Table 1. In Nodargi and Bisegna [47], 'a sufficiently large friction angle is assumed to be available', and the friction coefficient was set to 0.6 in the FDEM simulation, as it is large enough to avoid sliding. The time steps for Scheme A and B were 1.0 × 10 −7 s and 1.0 × 10 −6 s, respectively.  Figure 15 shows the failure modes of the multi-span pointed arch church under gravity and a constant horizontal acceleration from both the FDEM simulations and the limit analysis [46]. It can be observed that they are in excellent agreement. Seven hinges were predicted by both the FDEM simulations and the limit analysis [47], and the locations of these hinges are the same. As shown in Figure 15, the right part of the church inclined to the right, while the left part still stood vertically without noticeable deformation. Three hinges appeared within the arch at the middle span, i.e., one extrados hinge and two intrados hinges; therefore, the arch turned into a three-hinge mechanism. Another two hinges formed in the arch at the right span, i.e., one extrados hinge and one intrados hinge, due to the deformation of the structure. The other two hinges formed at the bottom of two right buttresses due to their rotations under the horizontal acceleration λg.  Table 4 presents the seismic capacity factor λ* of the multi-span pointed arch church from the FDEM simulations and the limit analysis [47]. Apparently, the λ* predicted by the FDEM (both Scheme A and B) is as same as that from the limit analysis, and the difference is only 0.1%. To further examine the differences between Scheme A and B, an extrados hinge, denoted as hinge II in Figure 15b,c, was selected. The time histories of hinge II in the two schemes with λ = 0.0884 are shown in Figure 16. Each curve terminates at the instant when hinge II disappears since the two adjacent voussoirs at the hinge separated completely in the final and fell down.  The computations were run on a Desktop with a 3.70GHz Intel i7-8700K CPU. It was found that the computational time for Scheme A is 2.19 times that for Scheme B. The timehistory curves of Scheme A and B are almost coincident with each other before t = 4 s, as the angles of hinge II are very small at the beginning. Subsequently, both angles increased dramatically. The development of hinge II from Scheme B is always slower than that from Scheme A. The surviving time of hinge II from Scheme B is slightly longer (about 3%) than that from Scheme A. Since the difference is very small, and to reduce the computational time, monolithic buttresses can be assumed in the behaviour of this type of arches.

Parametric Investigation
Further simulations were conducted to investigate the effects of geometric characteristics (i.e., H/h and B/Rc) and the friction coefficient μ on the seismic capacity of buttressed masonry arches. Buttressed circular and basket-handle arches were selected as base cases, as shown in Figure 8. Basic geometric parameters of the buttressed arches are: (i) Rc = 1.0 m and t/Rc = 0.2 for both buttressed circular and basket-handle arches; (ii) d/Rc = 0.5 and r/Rc = 0.5 for buttressed basket-handle arches; and (iii) embrace angle β = 120°, 150° and 180° for both buttressed arches. The material properties are as same as in Table 1, and the time step is 1.0 × 10 −6 s.

H/h
Herein, the ratio of H/h (i.e., the buttress height over the arch height) was investigated. No sliding between voussoirs was considered. For both buttressed arches, the ratio of B/Rc was set to 0.8. The height of buttress H varied from 0.25 m to 3.0 m. The variations of λ* against H/h for both circular and basket-handle buttressed arches were plotted in Figure 17. Overall, the seismic capacity factor λ* declines along with the increasing H/h. A smaller H/h suggests shorter buttresses, corresponding to a lower centre of gravity from the ground; therefore, they are more stable under seismic motions. On the contrary, a larger H/h indicates higher buttresses, which are more apt to overturn under seismic motions. For the buttressed circular arch with β = 120°, a sudden drop of the slope was found between H/h = 1.0 and 1.3, giving rise to the discontinuity on the slope curvature. It is also observed from Figure 17 that the seismic capacity factor λ* deceases significantly along with an increase in embrace angle β. It is reasonable, because for a given H/h and since Rc is constant, a larger β leads to a larger h and, hence, a larger H, which will cause the structure to be more vulnerable to seismic motions, and thus a smaller seismic capacity factor λ* is obtained, and vice versa. It is also suggested that in order to obtain a relatively large λ*, H/h ≤ 1.0 is preferred.

B/Rc
The influence of the ratio B/Rc (i.e., the buttress width over the radius of the reference circle) on the seismic capacity factor λ* was studied. For both circular and basket-handle buttressed arches, the ratio of H/h was set to 3.0. The width of buttress B up to 1.0m with a constant step of 0.2 m was studied. Like in Section 4.1, a large friction coefficient μ was set to avoid sliding. The variations of λ* against B/Rc are shown in Figure 18. Pre-examinations on the minimum B/Rc to guarantee the arches are statically stable (i.e., they can stand under gravity) were conducted, and the values are tabulated in Table 5. Any buttress with a smaller B will lead to the collapse of the structure due to the horizontal thrust of the arch induced by gravity. Thus, all the curves in Figure 18 start from the minimum B/Rc that the arch structures are statically stable under gravity.  For both circular and basket-handle buttressed arches, λ* increases monotonically with an increase in B/Rc. Since Rc is constant, the wider the width B, the larger the λ*. For buttressed circular arches, the angle of embrace β has nearly no influence on λ*, as the three curves are quite close to each other. For buttressed basket-handle arches, the same applies except the λ* with β = 180° and B/Rc = 1.0, as a much higher seismic capacity is obtained. The reason is that with the combination of β = 180° and B/Rc = 1.0, the horizontal thrust of the arch is too small to push the buttress down, and the failure of the structure is dominated by hinge rotations of arch voussoirs only, resulting in a much higher λ*.

Friction Coefficient
The influence of friction coefficient μ on the seismic capacity factor λ* was also examined. For both circular and basket-handle buttressed arches, B/Rc was set to 0.8 and H was set to 3.0 m. The friction coefficient μ varied in a range from 0.1 to 1.0, covering most practical engineering materials. The variations of λ* against the friction coefficient μ are illustrated in Figure 19. As it is shown, for small friction, i.e., μ ≤ 0.5, λ* increases significantly with the increase of μ. However, after μ reaches some thresholds, λ* becomes constant, i.e., it is barely affected by the friction coefficient. This can be attributed to the change of failure modes. If μ is below the threshold, the sliding failure will occur, as shown in Figure 20. Once μ is beyond the threshold, no sliding occurs, and the hinge failure mode dominates. Therefore, careful choice of the material regarding the friction angle needs to be conducted with respect to the choice of the geometry of the arch.

Concluding Remarks
In this paper, the combined finite-discrete element method (FDEM) was employed to study the first-order seismic capacity of masonry arch structures subjected to both gravity and a constant horizontal acceleration. The simulation results are validated with data from the existing literature.
In Sections 3.1 and 3.2, the seismic capacity factors for circular and basket-handle arch structures were predicted by the FDEM, and the results were compared with those from the limit analysis and the UDEC. The FDEM results show that the failures are dominated by a four-hinge failure mechanism should sliding be avoided. Both the seismic capacity factor λ* and the hinge locations were well verified. Supporting the high buttresses, the buttressed arches are much more vulnerable to seismic motions than the corresponding unbuttressed arches. In Section 3.3, the seismic capacity of a multi-span buttressed pointed masonry arch structure was studied with the FDEM, and the obtained results were almost identical to that from the limit analysis, illustrating the reliability and robustness of the FDEM in analysing arch structures with sophisticated shapes.
A parametric investigation on H/h, B/Rc and the friction coefficient μ was conducted in Section 4 to evaluate their effects on λ*. It was assumed that sliding was avoided when H/h and B/Rc varied, and thus the failure was dominated by pure hinging. It was found that λ* decreases with the increase of H/h, while it increases along with an increase in B/Rc. Moreover, a small μ can result in pure sliding or combined sliding-hinging failure, leading to a lower λ*. If sliding occurs, λ* increases along with an increase in μ. When μ reaches some threshold, pure hinging failure dominates and λ* becomes almost constant, indicating that sliding is unlikely to happen any longer. In reality, the load multiplier λ is usually no higher than one. Though only simple investigations have been performed over limited parameters, a more comprehensive parametric study and general conclusions regarding the friction, shape of the arches etc., can be expected in future work.
In general, the FDEM has been proven as a useful and reliable tool in predicting the first-order seismic capacity of masonry arch structures subjected to both gravity and constant horizontal acceleration. Thus, this method can be applied to evaluate the mechanical behaviour of masonry arch structures subjected to more sophisticated loads [48], or with steel-ties between buttresses [49], in future.