1 Introduction

Masonry arches are widely employed in architectural heritages and historical bridges due to their appealing appearances. Typically, they are assembled with bricks or stone voussoirs. In some aged masonry arches, the strength of mortar is too low to bond individual voussoirs reliably together, while in others, mortar is absent, and therefore, loads are undertaken by interactive and frictional forces between arch voussoirs. These facts lead to a widely held dry-joint assumption. In reality, masonry arches are frequently subjected to support settlements during their service life. Excessive and uneven support settlement is very hazardous, and it can result in the collapse of arches. Such behavior is highly nonlinear, and research on the sensitivity and vulnerability of dry-joint masonry arches induced by spreading supports is warranted.

Although masonry arch structures have been studied by many researchers in past decades, most of them focus on seismic safety [1, 2]. Dimitri et al. [3] investigated the dynamic behavior and resistance of multi-drum columns and masonry arches on buttresses subjected to step and harmonic impulses, and a set of parametric study with different geometries and loading conditions was performed. Later, Dimitri and Tornabene [4] studied the stability of unreinforced masonry arches and portals subjected to constant horizontal ground acceleration combined with the existence of gravity. In addition to the examination on seismic capacity, some work focuses on the maintenance. For example, the effect of bonding an FRP sheet to the intrados of a circular arch on the minimum value of the lateral thrust was evaluated by De Lorenzis et al. [5].

Research on collapse of masonry arches subjected to support settlement is still limited and has not been explored in detail yet [6]. Pioneering work can be traced back to 1960s when limit analysis was employed by Heyman [7]. In Heyman [7], masonry arch voussoirs were modeled with no-tension material and friction between voussoirs is considered sufficient to prevent sliding. Following the work of Heyman [7, 8], arch voussoirs are considered as rigid blocks [9,10,11,12] in most studies.

Experimentation is the most direct approach to investigate the failure of masonry arches with spreading supports. Ochsendorf [13, 14] examined two small-scale circular masonry arch models on horizontally spreading supports. Voussoirs were made from cast concrete, and a number of collapse characteristics, e.g., span increase ratio, initial and failure hinge angle, dip of crown, etc., were discussed. Later, an extensive experimental program about Gothic masonry arches of various shapes with both horizontal and vertical support movements was presented, as well as the comparisons with theoretical predictions [15]. Smars [16] investigated scaled masonry arch models with both horizontal and vertical support movements at a rate of 0.4 mm/s. With the assistance of high-speed camera, collapse mechanism of arch models was recorded well.

In spite of the attempts aforesaid, it is worth mentioning that most experiments were conducted on the reduced scale models. Furthermore, structure tests are not easy to perform and results may be influenced by a number of factors like assembly error or construction imperfection. To overcome such shortcomings, computational analysis becomes popular. Coccia et al. [17] proposed a computational approach to calculate the magnitude of displacements causing the collapse of circular arches. The approach is based on Heyman’s assumptions [7], and it was employed to study the failure behavior of arches with moving abutments by Ochsendorf [13]. Zampieri et al. [18] applied the principle of virtual work together with thrust line analysis to predict failure mechanism of masonry arches with different support movements, and the results were in good agreement with test data of masonry arches with mortar. Galassi et al. [19] presented a novel numerical procedure based on combinatorial analysis to identify collapse and limit settlements of masonry voussoir arches. Three hinges were captured, and failure characteristics such as hinge locations and span increase were verified. For more sophisticated structures, masonry buttressed arches with spreading supports were also examined by Zampieri et al. [20] using the same approach in Zampieri et al. [18]. Horizontal, vertical and inclined support settlements were investigated, and it was shown that the buttressed arch is more vulnerable to the horizontal displacement than the inclined or vertical settlements.

Although the computational procedures aforesaid are able to capture the failure of masonry arches subjected to support settlements, most of them followed the assumption in Heyman [7] that masonry arch voussoirs were modeled with no-tension rigid blocks of infinite compressive strength and sufficient friction without sliding. However, for large-size masonry arch structures, deformation and fracture of voussoirs cannot be ignored. Furthermore, the assumption of sufficient frictional force available results in the absence of sliding failure modes. The combined finite–discrete element method (FDEM) is a novel computational approach developed by Munjiza since 1990s [21,22,23,24]. It is a special branch of the DEM family. Within the framework of this method, structures are discretized using discrete elements, and a finite element mesh is incorporated into each discrete element. Thus, it can provide more accurate results, especially in the aspects of contact behavior and structural deformation. The details on the FDEM can be referred to the monographs [25,26,27]. The FDEM has been employed in discontinuous modeling of brittle and quasi-brittle applications such as glass impact fracture [28,29,30,31], shell fracture [32], mountain slope failure [33] and uniaxial rock test [34, 35]. Summary on structural applications with the FDEM was addressed by Munjiza et al. [36]. The FDEM was employed to analyze dry-stone masonry by Smoljanović et al. [37, 38] with a cohesive fracture model proposed by Zivaljic et al. [39]. A combined FDEM model has been developed to investigate masonry structures subjected to monotonic and cyclic loads by Smoljanović et al. [40, 41]. Recently, dynamic failure of masonry arch structures subjected to base impulses has been simulated using the FDEM with the consideration of block fracture [42]. The FDEM was also used in the analysis of the stability of dry-stone arches due to seismic action. The collapse mechanisms of 12 types of stone arch with different geometries subjected to three real horizontal and vertical ground motions were studied by Balic et al. [43].

Albeit there is some research work on masonry failure with the FDEM as mentioned above, very limited attention has been given to the failure characteristics of masonry arches induced by spreading supports. Since the collapse of masonry arches in such cases is kinematic driven, the FDEM is naturally more suitable as a simulation tool. Moreover, the deformation and fracture of arch voussoirs can be predicted well. Thus, research on the collapse of masonry arches due to spreading supports using the FDEM is still in demand.

In this paper, the collapse behavior of dry-joint masonry arches induced by spreading supports is simulated by using the two-dimensional FDEM program ‘Y’ [44]. Friction and voussoir fracture are taken into account to examine their influences. Lemos [45] showed that a two-dimensional analysis is reasonably accurate and acceptable for the in-plane failure of masonry arches. Following the masonry arch discretization and modeling, numerical examples are presented. A parametric study is conducted to investigate the influences of a set of parameters, e.g., embrace angle, number of voussoirs and ratio of thickness to radius. Friction thresholds between the pure hinge failure and the combined sliding–hinge failure are also provided in terms of these parameters. Finally, comparisons on collapse behavior of a masonry arch with and without voussoir fracture are conducted. In general, collapse simulation of dry-joint masonry arches induced by spreading supports with the FDEM is applicable and reliable.

2 Methods

Computational aspects regarding collapse simulation of masonry arches induced by spreading supports are introduced in this section. Masonry modeling and voussoir discretization are firstly introduced. The FDEM is a naturally suitable tool for kinematic problems, and it is capable of capturing solids transitioning from continua to discontinua. Key aspects of the FDEM with respect to masonry arch collapse are also briefly introduced in this section.

2.1 Masonry modeling

2.1.1 Basic geometry

A circular arch is illustrated in Fig. 1, in which L is the arch span, H is the arch clearance, β is the angle of embrace, h is the arch thickness, Ri and Re are the inner radius and the outer radius, respectively, and R is the mean radius, i.e., R = (Ri + Re)/2.

Fig. 1
figure 1

Geometry configuration of a circular arch

Full size image

With the movement of support(s), masonry arches fail mostly in the three-hinge mechanism. From the view of kinematics, the failure process can be described by four stages: (i) initiation of the three-hinge mechanism, (ii) stable development of the three-hinge mechanism, (iii) formation of additional hinge(s) and (iv) collapse. A typical three-hinge mechanism failure due to horizontal support movement is shown in Fig. 2. Figure 2a corresponds to the stage (i), i.e., one extrados hinge and two intrados hinges present immediately with the initiation of support movement. Ochsendorf [13] addressed that for symmetric configurations, extrados hinge forms at the crown (apex of the arch), and two intrados hinges appear with an angle of βo (measured vertically from each side of the crown in Fig. 2a) firstly due to the horizontal support movement. Then, with the increase in the support movement, the arch is still in a stable state following the three-hinge mechanism. Subsequently, an additional hinge presents (stage (ii)), the three-hinge mechanism deteriorates, and therefore, the masonry arch becomes unstable, i.e., in the stage (iii). After that, the arch falls quickly and a significant change in reaction forces can be observed. At this point, the arch is deemed collapsed (stage (iv)). Time elapsing between the stages (iii) and (iv) is very short. Referring to Fig. 2b, the angle between extrados and intrados hinges at the beginning of collapse is denoted as βu.

Fig. 2
figure 2

Arch failure in three-hinge mechanism: a initial hinge angle βo; b collapse hinge angle βu

Full size image

2.1.2 Voussoir discretization

Masonry arches are generally built with voussoirs, and each voussoir is a discrete block. Figure 3 shows the discretization of two adjacent voussoirs (i.e., A and B) from a dry-joint masonry arch within the framework of the FDEM. Each is meshed with three-node constant strain triangular (CST) elements. Element interfaces exist within each single voussoir in order to simulate the potential fracture. However, no element interface is implemented between the voussoirs A and B although the adjacent element edges coincide with each other, suggesting that voussoirs A and B can simply separate without obstacle of joint force, i.e., dry joint. Contact and friction can act due to the relative movement between such two adjacent voussoirs. Considering the surface roughness of masonry arch voussoirs [42], a friction coefficient of 0.6 is employed in the following FDEM simulations if it is not specifically stated.

Fig. 3
figure 3

Element discretization of masonry arch voussoirs

Full size image

2.2 FDEM

2.2.1 Element motions

In the FDEM, element motions follow the traditional DEM equations. Translational and rotational motions of a discrete element i are evaluated based on the Newton’s second law as follows:

$$ m_{i} \frac{{{\text{d}}^{2} }}{{{\text{d}}t^{2} }}{\mathbf{r}}_{{\mathbf{i}}} + m_{i} {\mathbf{g}} = {\mathbf{F}}_{{\mathbf{i}}} $$
(1)
$$ {\text{J}}_{i} \frac{\text{d}}{{{\text{d}}t}}{\varvec{\upomega}}_{{\mathbf{i}}} = {\mathbf{T}}_{{\mathbf{i}}} $$
(2)

where mi is the mass of discrete element i; ri is the position; g is the acceleration of gravity; Ji is the moment of inertia about the center of element i; ωi is the angular velocity about the element center; Fi is the resultant force acting on the element center; and Ti is the resultant moment about the center of element. Based on Eqs. (1) and (2), the velocity and position of an arbitrary discrete element can be explicitly obtained at each time step.

2.2.2 Contact force

Contact forces in the FDEM are determined from the overlapping area S (Fig. 4). The penetration of elemental area dA yields an infinitesimal contact force df,

$$ {\text{d}}{\mathbf{f}} = - {\text{d}}{\mathbf{f}}_{t} + {\text{d}}{\mathbf{f}}_{c} $$
(3)

where the subscripts ‘t’ and ‘c’ correspond to the target and the contactor, respectively. dft and dfc are defined in Eqs. (4) and (5), respectively, as follows:

$$ {\text{d}}{\mathbf{f}}_{t} = - E_{p} {\text{grad}}\varphi_{c} (P_{c} ){\text{d}}A $$
(4)
$$ {\text{d}}{\mathbf{f}}_{c} = - E_{p} {\text{grad}}\varphi_{t} (P_{t} ){\text{d}}A $$
(5)

where Pc and Pt are the points on contactor and target that share the same coordinate in S; φc and φt are pre-defined potential functions; grad is the gradient; and Ep is the contact penalty. Contact force f is obtained in Eq. (6) by integrating dA over S, as

Fig. 4
figure 4

Contact forces in the FDEM

Full size image
$$ {\mathbf{f}} = E_{p} \mathop \int \limits_{S} [{\text{grad}}\varphi_{c} (P_{c} ) - {\text{grad}}\varphi_{t} (P_{t} )]{\text{d}}A. $$
(6)

2.2.3 Fracture

Fracture of large-size voussoirs is possible during the collapse of masonry arches. In the FDEM, cracks are assumed to initiate and propagate along the element edges. As shown in Fig. 3, interfaces where fracture models are implemented in voussoirs A and B are created between each pair of adjacent discrete elements. The transition from continua to discontinua is determined by the interfacial deformation δ (Fig. 5). When δ reaches its critical value δc, adjacent discrete elements separate completely. Albeit elements can only separate along their boundaries, mesh bias can be reduced to the minimum if unstructured and small enough elements are employed [24].

Fig. 5
figure 5

Element interfaces: a before deformation; b after deformation

Full size image

The fracture model used herein is similar to Hillerborg’s [46]. Masonry voussoirs are assumed to be quasi-brittle materials following the Mode-I fracture [47]. Two elements start to separate when the bonding stress σ reaches the tensile strength ft. Meanwhile, σ decreases with the increase in δ. When δ reaches δc, σ drops to zero. This process is described using a descending σ–δ curve named as ‘strain softening’. Once σ equals zero, these two elements separate completely and in-between joint force vanishes. The cohesive model for voussoir rupture was proposed by Zivaljic et al. [39]. Note δp is the onset separation when σ = ft, full relationship between σ and δ is given in Eq. (7) according to Munjiza [25],

$$ \sigma = \left\{ {\begin{array}{*{20}l} {\left[ {2\frac{\delta }{{\delta_{p} }} - \left( {\frac{\delta }{{\delta_{p} }}} \right)^{2} } \right]f_{t} } \hfill & {0 \le \delta \le \delta_{p} } \hfill \\ {f_{t} z } \hfill & {\delta_{p} < \delta \le \delta_{c} } \hfill \\ {2\frac{\delta }{{\delta_{p} }}f_{t} } \hfill & {\delta < 0 } \hfill \\ \end{array} } \right. $$
(7)

in which z is a heuristic parameter defined by Eq. (8) referring to Hordijk [48],

$$ z = \left[ {1 + (c_{1} D)^{3} } \right]e^{{ - c_{2} D}} - D(1 + c_{1}^{3} )e^{{ - c_{2} }} $$
(8)

where D is a fracture index within [0, 1]. If δp < δ ≤ δc, D = (δ − δp)/(δc − δp); if δ > δc, D = 1, indicating the complete separation. In Zivaljic et al. [39], \( {\text{c}}_{ 1} \) = 3.0 and \( {\text{c}}_{ 2} \) = 6.93 were adopted for masonry voussoirs. Based on the above formulations, a normalized strain softening curve is presented in Fig. 6, and the area below the descending curve equals the fracture energy Gf, as shown in Eq. (9).

Fig. 6
figure 6

Normalized bonding stress z versus fracture index D

Full size image
$$ G_{f} = \mathop \int \limits_{0}^{{\delta_{c} }} \sigma {\text{d}}\delta $$
(9)

3 Validations

In this section, the collapse behavior of dry-joint masonry arches induced by spreading supports is simulated with the FDEM. No fracture is considered. In Sect. 3.1, masonry arches under concentrated forces applied on the centers of voussoirs are examined. Section 3.2 investigates the collapse characteristics of symmetric masonry arch models. The gravity (i.e., g = 9.8 m/s2) is included in all examples. The simulation results are in good agreement with the data from existing experiments and other numerical simulations.

3.1 Arches with non-symmetric concentrated loads

Galassi et al. [19] conducted experimental and numerical predictions on the collapse layout and limit support movement of circular dry-joint masonry arches. In addition to the gravity, concentrated loads are applied on the centers of some voussoirs. Table 1 shows the geometry and loading conditions for the FDEM modeling. Cases are in abbreviation of A1 and A2, respectively. It is noted that data in brackets are from the tests in Galassi et al. [19]. Although some minor geometric differences can be found between the models for the FDEM simulation and the model in the tests due to manufacture imperfection and assembly errors, they are negligible. The voussoirs are assumed to be elastic, and the material parameters are given in Table 2. The meshes and voussoir ID in the FDEM simulations are illustrated in Fig. 7. There are 26 elements in total in Case A1, while there are 60 elements in Case A2. In both Cases A1 and A2, the left support moves outward horizontally (Fig. 7). Support movement rate is 1.0 mm/s so that quasi-static effect can be achieved. The time step in the FDEM simulations is set to 1.0 × 10−7s.

Table 1 Geometry and load conditions for Cases A1 and A2
Full size table
Table 2 Material properties of masonry voussoirs
Full size table
Fig. 7
figure 7

Configurations of the masonry arches in FDEM: a Case A1; b Case A2

Full size image

The deformed masonry arches and hinges (highlighted with circles) from the FDEM simulations are shown in Figs. 8 and 9, and the corresponding snapshots from tests are also presented for comparison. Figure 8 shows the kinematic failure sequences of the arch in Case A1 with different span increase d. It is apparent that the FDEM simulation accurately identifies the deformed three-hinge mechanism and the hinge locations. In Fig. 8a, the three-hinge mechanism initiates as long as the support moves. Here, the angle between two adjacent voussoir edges where a hinge appears is defined as the hinge angle. When the support moves further (Fig. 8b, c), the three-hinge mechanism develops and the hinge angles become larger and larger. Since the three-hinge mechanism is capable of withstanding the gravity and additional loads, the masonry arch is still in a stable state. However, the three-hinge mechanism will reach its limit state with the increasing support displacement, and subsequently, an additional hinge appears. Due to the instability of the four-hinge mechanism, the arch collapses promptly (Fig. 8d). The limit support settlement, i.e., the maximum settlement that can be withstood by the arch, is 42.26 mm in the FDEM simulation.

Fig. 8
figure 8

Failure sequences of Case A1

Full size image
Fig. 9
figure 9

Failure sequences of Case A2

Full size image

Sequences for Case A2 are presented in Fig. 9. Similar to Case A1, a three-hinge mechanism is observed too. Once the fourth hinge has formed, the masonry arch becomes unstable and it will collapse promptly. The limit settlement for Case A2 is 44.91 mm according to the FDEM simulation. Like in Case A1, failure of Case A2 also can be attributed to the instability of the four-hinge mechanism. Based on the collapse modes, it is apparent that the simulation results are in good qualitative agreement with the experimental observations.

The limit settlement ΔLc and the span increase ratio η for Cases A1 and A2 are presented in Table 3. The span increase ratio η is defined as

$$ \eta = \frac{{\Delta L_{c} }}{L} $$
(10)

where L is the span of arch defined in Fig. 1. Apparently, both the limit settlement and the span increase ratio from the simulations are close to those experimental and numerical data in Galassi et al. [19]. The relative difference e is defined in Eq. (11) and presented in Table 3, in which Ω is the data to be compared and Ωexp is the experimental data. In Case A1, e = 5.4% for both the FDEM and the combinatorial analysis approach [19]. However, e for FDEM is 10.3% lower than that from the numerical data by Galassi et al. [19] in Case A2.

Table 3 Comparisons on the limit settlement and span increase ratio
Full size table
$$ e = \left| {\frac{{\varOmega - \varOmega_{\exp } }}{{\varOmega_{\exp } }}} \right| \times 100\%. $$
(11)

3.2 Arches with symmetric configurations

Two small-scale masonry arch models with symmetric configurations were investigated by Ochsendorf [13], and the geometries of the two arch models (denoted by Arches O1 and O2, respectively) are shown in Table 4. Each arch model is composed of 16 uniform voussoirs, which are made from cast concrete with the material properties in Table 5. The voussoir ID and the mesh adopted in the FDEM simulation are illustrated in Fig. 10. There are 64 elements in total for both Arches O1 and O2, and the characteristic element size is 0.05 m. In the FDEM simulation, the right support of both arches moves outward (Fig. 10) at a rate of 1.0 mm/s. The time step in the FDEM simulation is 1.0 × 10−7 s. No loads are applied on the arches except the gravity.

Table 4 Geometries of Arches O1 and O2
Full size table
Table 5 Material properties of voussoirs
Full size table
Fig. 10
figure 10

Mesh of arches with symmetric configurations: a Arch O1; b Arch O2

Full size image

Figure 11 shows the deformation of Arch O1 at different span increase d. A three-hinge mechanism is clearly shown in Fig. 11a, and the arch is stable until d = 64.33 mm. Afterward, two additional hinges present simultaneously between the end voussoirs and the supports as shown in Fig. 11b. Due to the instability of the five-hinge mechanism, the masonry arch tends to collapse instantly.

Fig. 11
figure 11

Failure process of Arch O1 (horizontal displacement contour, unit: m): a d = 0.06433 m; b d = 0.06462 m

Full size image

Failure process of Arch O2 at different d is presented in Fig. 12. The three-hinge mechanism can be observed at the beginning, as shown in Fig. 12a, and it is stable up to around d = 50 mm. Unlike Arch O1, the hinges initially present at βo = 60° and then move up to the next dry-joint when Arch O2 is about to collapse as illustrated in Fig. 12c. Due to the instability of the five-hinge mechanism, Arch O2 finally collapses around d = 51.18 mm with βu = 50°.

Fig. 12
figure 12

Failure process of Arch O2 (vertical displacement contour, unit: m): a d = 0.02 m; b d = 0.05 m; c d = 0.05112 m; d d = 0.05125 m

Full size image

Figure 13 shows the deformed arches at their moment of collapse in the experiments [13]. Generally, the failure modes predicted by the FDEM are in excellent agreement with the experimental observations. The collapse of Arches O1 and O2 was also simulated using the program ArchSpread [13]. Thus, comparisons are made on the collapse characteristics in Table 6. The relative differences e on span increase ratio η are 7.1% and 9.2% for FDEM simulation for Arches O1 and O2, respectively. And their counterparts from numerical analysis by Ochsendorf [13] are 9.7% and 12.8%, respectively. Apparently, the results based on the FDEM have a good accuracy against the experimental results, and the relative difference on span increase ratio η is better than that from the numerical approach in [13]. There are also some little discrepancies, which may be attributed to the imperfect casting of arch voussoirs in the experiments. For instance, the embraced angle for each voussoir of Arch O1 is 11.25° in the FDEM modeling, and hence, βo = βu = 56.25°. However, such precise concrete casting is highly difficult in the experiments and hence the embraced angle for each concrete voussoir is approximately 10° in Ochsendorf [13].

Fig. 13
figure 13

Deformed masonry arches from experiments [13]: a Arch O1; b Arch O2

Full size image
Table 6 Comparison on collapse characteristics of Arches O1 and O2
Full size table

3.3 Summary

Numerical simulations on the collapse responses of masonry arches subjected to moving supports are presented and validated with the experiments and simulations conducted by other researchers. In Sect. 3.1, masonry arches with non-symmetric concentrated loads are investigated. The collapse of masonry arches with symmetric configurations is examined in Sect. 3.2. The failure modes are qualitatively compared with the experimental snapshots, and a very good agreement is achieved. More detailed comparisons are made on span increase ratio, limit settlement, initial and failure hinge angles, etc. Apparently, the numerical results based on the FDEM agree well with the results from known sources. Based on above analysis, some characteristics for the collapse of masonry arches induced by spreading supports are summarized as follows:

  1. 1.

    For non-symmetric configurations (either non-symmetric masonry arch or symmetric masonry arch with non-symmetric loads or boundary conditions), the responses are non-symmetric. A three-hinge mechanism dominates from the onset of support movement up to the formation of an addition hinge. Then, the masonry arch collapse shortly due to the instability.

  2. 2.

    For symmetric configurations, the three-hinge mechanism also governs as long as the support moves. However, when it develops to a critical level, two additional hinges present simultaneously. Then, the arch collapses quickly because of the instability of the five-hinge mechanism.

4 Parametric study

A parametric study is performed in this section to show the FDEM results and provide comparison examples for future research. A base case is firstly proposed with geometric parameters as follows: angle of embrace β = 120°, mean radius R = 10 m, h/R = 0.15. The arch is made up of 12 identical voussoirs, and the material properties for voussoirs are same as those in Table 5. No fracture is considered and no loads are applied to the arch except the gravity. A movement at a rate of 0.05 m/s is applied to the right support, and a time step of 1.0 × 10−6s is used in the simulations. The following parametric study is based on this example with the consideration of horizontal outward support movement. During the analysis, one parameter is changed at a time, while others are kept unaltered.

4.1 Embrace angle

The embrace angle β determines the shape of a masonry arch. In this study, β is investigated within the range from 60° to 180° at a constant step of 20°, and the results are shown in Fig. 14. It is observed that with the increase in β, the span increase ratio η increases monotonically and slightly for 60° ≤ β ≤ 100°, while it decreases monotonically and dramatically for 100° ≤ β ≤ 180°. The maximum η (η = 17.29%) is reached at β = 100°. For β = 180°, the corresponding η is only 7.55%, suggesting masonry arches with such an embrace angle are vulnerable to support movement. Therefore, an optimal span increase ratio can be achieved for masonry arches when the embrace angles are in the range of 80°–100°.

Fig. 14
figure 14

Variation of span increase ratio η against embrace angle β

Full size image

4.2 Number of voussoirs

The variation of span increase ratio η against the number of voussoirs n is plotted in Fig. 15. It is found that η drops significantly from n = 4 to n = 12. However, η becomes stable with further increase in n, i.e., it varies slightly around 12.5%. This trend suggests that when the number of voussoirs n reaches certain value (e.g., n = 12 in this specific configuration), the span increase ratio η is hardly affected. It is worth mentioning that although the span increase ratio η is relatively larger when the number of voussoirs n is small, an excessively small n (e.g., n < 6) is rarely adopted in the design of masonry arches. Thus, for a masonry arch built with a reasonable number of voussoirs, the span increase ratio would be stable and within a narrow range.

Fig. 15
figure 15

Variation of span increase ratio η against the total number of arch voussoirs n

Full size image

4.3 Ratio of thickness to radius

With a constant radius R, a lager h/R means that the arch is built with larger voussoirs, which leads to a larger self-weight. Meanwhile, it also results in larger contact areas between the voussoirs, i.e., larger frictions. Therefore, the failure mode of arch becomes more complicated. The simulations show that it changes from the pure hinge failure to a combined failure mode including sliding and hinges with the increase in h/R, as shown in Fig. 16. For a friction coefficient 0.6, the threshold between these two types of failure mode is h/R = 0.35.

Fig. 16
figure 16

Two failure modes with different h/R: a pure hinge failure with h/R = 0.3; b sliding-hinge failure with h/R = 0.5

Full size image

The variation of the span increase ratio η against the h/R within the pure hinge failure domain is given in Fig. 17. Apparently, the span increase ratio increases with the increasing h/R, exhibiting almost linearly. It may be attributed to the increase in voussoir size, as larger size voussoirs are able to resist larger support displacement.

Fig. 17
figure 17

Variation of span increase ratio η against h/R

Full size image

4.4 Friction coefficient

As mentioned in Sect. 4.3, the friction plays an important role in the failure mode of masonry arches. In this section, an in-depth investigation has been conducted. Through numerous simulations, the variations of the friction coefficient µ against the ratio of thickness to radius h/R, the number of voussoirs n and the embrace angle β are plotted in Fig. 18, showing the thresholds of friction coefficient µ which distinguishes the two failure modes: the pure hinge failure and the combined sliding–hinge failure.

Fig. 18
figure 18

Thresholds of friction coefficient µ between different failure modes against other parameters: a thickness-to-radius ratio h/R; b number of voussoirs n; c embrace angle β

Full size image

Figure 18a provides the variation of threshold µ against h/R. When h/R increases, the size of voussoirs also increases, and hence, the weight of voussoirs increases. As a result, a larger friction coefficient µ is required to prevent voussoirs from sliding. Figure 18b shows the threshold µ against the total number of voussoirs n. Since n does not change the geometry of arch, the minimum µ for pure hinge failure mode is almost at a constant level, i.e., around µ = 0.23.

Figure 18c presents the variation of the threshold µ against the embrace angel β. The minimum µ is in descending trend in general. Specifically, the curve can be divided into two parts: (1) 40° ≤ β ≤ 100°, the threshold deceases dramatically; (2) 100° ≤ β ≤ 180°, the corresponding µ is varies very slightly around µ = 0.24. This can be attributed to the change of thrust in masonry arch. When β is small, little thrust is provided by masonry arches. Accordingly, the compression between voussoirs is also very small. This requires a relatively larger friction coefficient µ to compensate the reduction in friction force to avoid the combined sliding–hinge failure mode. On the contrary, when the embrace angle β becomes larger, a larger thrust can be established in the arch, so do the compression forces between voussoirs. Consequently, only a smaller µ is needed to assure that the pure hinge failure dominates.

5 Voussoir fracture

In the previous sections, fracture of voussoirs is not considered. However, fracture and fragmentation of voussoirs are widely observed in masonry arches subjected to support movements, especially for the masonry built with large-size voussoirs. In Sect. 2.2.3, the cohesive fracture model is introduced. Herein, it is employed to account for the rupture of voussoirs, and the example in Sect. 4 is adopted as a base case to investigate in further.

A finer mesh with the characteristic size of 0.2 m is employed for voussoirs in the FDEM simulation, as shown in Fig. 19. The material properties about the fracture of voussoirs are based on Smoljanović et al. [40, 41]: Gf = 150 J/m2 and ft = 2.0 MPa, which have been proven to be acceptable and reliable [42]. Other parameters (e.g., elastic material parameters, time step, support movement rate, etc.) are kept as same as those in Sect. 4. Munjiza and John [24] suggest that the characteristic element size should be much smaller than the length of plastic zone Δ when rupture is considered, and the lower bound of Δ is given by

$$ \Delta = \frac{E}{{4f_{t} }}\delta_{c} $$
(12)

where E is the elastic modulus, ft is the tensile strength and δc is the critical separation. For arch voussoirs, Eq. (12) yields Δ ≈ 0.75 m, suggesting the plastic zone may spread over 3–4 elements, and the employed mesh is able to capture the fracture behavior of voussoirs well.

Fig. 19
figure 19

Configurations of the base model considering voussoir rupture: a geometry; b mesh

Full size image

The failure process of the base arch example with the consideration of voussoir fracture is schematically shown in Fig. 20. It can be observed from Fig. 20a that similar to the non-fracture simulation, a symmetrically located three-hinge mechanism has formed initially. When the right support moves to some critical value, two additional hinges appear and the arch becomes unstable. The end voussoirs start to rotate toward the adjacent supports. Once they impact on the supports, fracture occurs (Fig. 20b). The masonry arch has failed at this moment although it has not collapsed completely. When the support moves further a bit, plenty of fractures are identified and they spread to adjacent voussoirs toward the crown, as shown in Fig. 20c. Fragments fall and the stability of the arch is lost, leading to the complete collapse of the masonry arch.

Fig. 20
figure 20

Failure process of the base case with fracture (vertical velocity contour, unit: m/s): a d = 2.0 m; b d = 2.1425 m; c d = 2.1575 m

Full size image

The angle between the left support and the adjacent end voussoir is denoted as α, and Fig. 21 plots the variation of α against the normalized span increase d/L until the end voussoirs impacts on the supports. The results based on non-fracture voussoirs are also presented for comparison. Figure 21b illustrates that the two curves are highly similar. Initially, α increases with the increase in d, and this trend is kept until α reaches its maximum. Afterward, α drops significantly since the end voussoirs rotating toward the supports. The two curves are nearly identical in the ascending stage. However, the ultimate span increase ratio η of the masonry arch with voussoir fracture is slightly larger than that from the non-fracture masonry arch. The possible reasons for this phenomenon are listed as follows:

Fig. 21
figure 21

Variation of α against the normalized span increase (d/L): a definition of α; b rotation curve

Full size image
  1. 1.

    Fracture of voussoirs consumes some kinetic energy, and hence the motions of voussoirs are eased;

  2. 2.

    Fragments generated after the fracture impede the hinge rotation and postpone final collapse to some extent.

6 Concluding remarks

The combined finite–discrete element method (FDEM) is employed in this paper for simulating the collapse of dry-joint masonry arches induced by spreading supports. The results are well compared and validated with those from experimental and computational sources, and a good agreement is achieved.

  1. 1.

    In Sect. 3.1, deformed masonry arches in the form of three-hinge mechanism simulated by FDEM are validated with tests. The results validate the reliability of the FDEM simulation not only qualitatively but also quantitatively. For non-symmetric configurations, a four-hinge mechanism presents before the collapse.

  2. 2.

    Section 3.2 investigates the collapse of concrete masonry arch models under gravity. Since configurations are symmetric, a five-hinge mechanism is observed before failure. The failure characteristics from the FDEM simulation are in good agreement with the experimental and numerical data.

A parametric study is performed on the embrace angle, number of voussoirs, thickness-to-radius ratio and friction coefficient. It is found that the optimal embrace angle to resist a horizontal support movement is around 80°–100°. Regarding the thickness-to-radius ratio, a higher span increase ratio would achieve when h/R is raised. The thresholds of friction between the pure hinge failure and the combined sliding–hinge failure are also presented. A base case with the consideration of voussoir rupture is simulated. It is found that the fracture and fragments may limit the rigid motions of voussoirs to some extent.

The simulations show that the FDEM is highly applicable and reliable to investigate the collapse of masonry arches induced by spreading supports. Although only dry-joint circular arches are investigated in this paper, the applications can be extended to the masonry arches with mortar, as well as in other shapes.