Limit analysis of masonry arches and domes with finite strength: funicular analysis versus stability area method

semi-analytical Durand-Claye’s method. For benchmark case studies, such as symmetric masonry arches and domes with specific stereotomies subjected to axi-symmetrical load conditions, the set of statically admissible solutions compatible with equilibrium and strength requirements is graphically determined in terms of the horizontal thrust and its eccentricity at the crown, examining the shape of the stability area. Assuming an infinite value for the friction coefficient, the collapse condition is reached when the stability area shrinks to a single point. The results obtained from both of these methods are in excellent agreement. The influence of compressive strength on the bearing capacity of the structures is also discussed.


Introduction
The research proposed here aims at studying the collapse behavior of masonry arches and domes, with particular attention to the influence of the compressive strength.In the framework of the lower bound theorem of limit analysis, the maximum static collapse load is determined by solving an optimization problem using an algorithm developed specifically based on the concept of force density [3,53].The results are validated through a modern semi-analytical Abstract This study, framed within the context of the lower bound theorem of limit analysis, aims to assess the anti-funicular equilibrium of masonry arches and domes using a computational approach based on the constrained force density method.In contrast to the commonly adopted classical Heyman's assumptions, the approach proposed here considers the effects of finite compressive strength in the material.Assuming a fixed plan projection for a network with independent sets of branches, a suitable set of local constraints is enforced at each joint to account for the limit bending moment resulting from the material's assumptions, including limited compressive strength and zero tensile strength.Additionally, the stereotomy of the voussoirs is considered by assigning a geometric law to the joint inclination.The collapse load is determined by formulating a multi-constrained maximization problem.The method is validated using a modern version of the formulation of the graphical Durand-Claye's method [4,6,[20][21][22].
The objectives of this research arise from a somewhat paradoxical situation, which will be better illustrated in the following Sections.While the complexity of the mechanical response of masonry poses undoubted challenges in the modeling of masonry elements and structures, there is an increasingly current demand for fast methods for the limit analysis of such structures, particularly arches and vaults.These methods prove to be useful tools in the preservation of architectural heritage.
Faced with this framework, this contribution is configured as an attempt to consider more refined hypotheses compared to Heyman's assumptions [33] in order to take into account the limited compressive strength of masonry.Moreover, both approaches are capable of taking stereotomy (that is, the art of cutting stones) into account.

A challenging topic: modelling the masonry material
Structural analysis of masonry arches and domes is a challenging subject due to the complex behavior of masonry, a heterogeneous material, and the geometry of the voussoirs that compose these structures.As noted by [11], defining the concept of 'masonry' is not a straightforward task, as it encompasses various materials depending on the components used (brickwork, stonework, etc.), block types, assembly method, the presence or absence of mortar, and the influence of workmanship.An in-depth examination of the approaches used in the scientific literature to describe masonry behavior is beyond the scope of this work.Interested readers can refer to [37], which highlights the challenges in conducting advanced experimental tests to determine masonry properties, given the unique characteristics of this material.
While not an exhaustive list, several studies, including those by [42,51,54], have aimed to model the compressive strength of masonry.The first contribution, focusing on unreinforced masonry, derives compressive strength from component and mortar properties, employing a phenomenological approach within a discrete framework.It uses a fictitious microstructure composed of linear elastic particles and non-linear interface elements.[54] proposes simple homogenized models for characterizing masonry behavior, utilizing a representative volume element and a finite element micromodelling approach.Finally, the work by [42] introduces various techniques for determining the homogenized strength domain of running or header bond masonry walls, either subjected to in-plane or transverse loads.The contributions cited above aim to develop models close to reality, capable of accurately capturing the structural response of masonry structures (or masonry elements).However, a drawback is the high uncertainty on the properties of the material, which makes it difficult to choose the appropriate mechanical parameters.Another difficulty lies in the computational effort required when adopting complex constitutive relations.

Classical limit analysis: Heyman's hypotheses
Faced with the issues described in the previous Section-together with the observation that one of the most notable characteristics of masonry is its uncertain tensile strength, which is often significantly lower than its compressive strength-a simplified approach can be followed, by assuming that masonry transfers only compressive stresses.As is well known, the simplest hypotheses on masonry are those proposed by Heyman [33], according to whom masonry can be modeled considering zero tensile strength, infinite compressive strength and an infinite value of the friction coefficient between the blocks.In the theoretical context of limit analysis, these hypotheses are recognized as an acceptable description of the behavior of masonry, such as to allow an adequate assessment of the stability of masonry structures, as well as the interpretation of crack patterns.As regards curved masonry structures, the safe theorem [33], deriving from this framework and exploiting the lower bound theorem of limit analysis, claims that stability is guaranteed if a thrust line can be found that is in equilibrium with the external loads and is located inside the vault.
It should be noted, however, that assuming an infinite compressive strength for masonry may lead to an overestimation of the collapse load, since the mechanical properties of the material are not adequately taken into consideration.

Simplified models considering the compressive strength
With respect to the approach described in Sect.1.2, it is observed that, in particular for some typologies of masonry structures (for example segmental arches and domes with depressed profile), considering the influence of a limited compressive strength is crucial from a mechanical and engineering point of view.
To overcome this drawback, some authors have conducted studies on the impact of limited compressive strength on the mechanical behavior of masonry structures, including masonry arches and vaults.They have proposed straightforward mechanical models for the material and simplified structural frameworks.
In the context of the lower bound theorem of limit analysis, the semi-analytical graphical procedure developed by Durand-Claye [21,22] for evaluating masonry arches and domes of revolution with symmetric geometry and loading conditions can be referenced.According to its modern version, as developed by [4,8], the bearing capacity of such structures can be rigorously determined by considering the limit bending moment at the joints, which depends on the limited compressive strength of the material.
Other approaches, utilizing the upper bound theorem of limit analysis, attempt to account, albeit approximately, for the effects of limited compressive strength through geometrical criteria.For instance, [32] discusses the mechanisms that arise in masonry arches by introducing a geometric safety factor, while [17] examines the collapse load of (reinforced) masonry arches with limited compressive strength.In the event of a collapse mechanism, crushing is modeled as interpenetration between the blocks, resulting in hinges formed at internal points of the joints.
Furthermore, there are contributions proposing non-linear elastic analyses.These analyses start by defining a constitutive relation for masonry-like materials with bounded compressive strength, and the non-linearity of this relation implies that loads are assigned incrementally in the applications.For certain case studies of masonry arches, represented as non-linear elastic beams, the solution is found either numerically [38], or explicitly [5].
In light of the complex context described above, the method presented here is theoretically grounded in the lower bound theorem of limit analysis.Its aim is to address masonry vaulted structures with intricate shapes, while considering both strength and geometric requirements.This approach conceptually aligns with funicular methods and thrust network analysis.

From graphic statics to the modern funicular approaches
Before delving deeper, a brief historical overview is crucial to recognize the roots of this approach in fundamental contributions centered around the thrust line concept, further developed in the context of graphical analysis.For more detailed insights, refer to [2].
In masonry arches, Robert Hooke [35] pioneered their structural assessment through graphical analysis, introducing the inverted hanging chain analogy.The thrust line concept was formalized by Gerstner [30], distinguishing the 'line of resistance' and 'line of pressure'.In the eighteenth century, this method extended to masonry domes, incorporating the effects of hoop forces, utilizing thrust line analysis, graphic statics, and membrane theory.Eddy [23] made the first graphical analysis contribution to masonry domes, studying the equilibrium of a hemispherical dome through graphical reinterpretation of membrane theory.Lévy [39] and Wolfe [57] followed, with recent contributions revisiting these methods, including [29,50,58].
This graphical approach laid the foundation for various computational procedures modeling thrust networks.Recent studies by [13,47], and [46] adhere to Heyman's hypotheses.Within this framework, [40] focuses on computational aspects, offering iterative procedures for both vertical and horizontal loads.Other static methods model masonry vaults as notension membranes [26,28,52] or derive truss-like stress paths from Airy stress functions [9,10,27,31,43], using mathematical tools for approximating stress in continuous bodies through discrete force networks.
Modern computerized static approaches facilitate the study of structures with complex geometry.Significant applications on spiral stairs by [12,48] showcase advanced computational strategies and consideration of non-trivial load conditions.The interest in methods based on the lower bound theorem of limit analysis is demonstrated by recent contributions like [44,45], performing limit analysis of masonry domes subjected to pseudo-static seismic forces, including graphical equilibrium methods such as the one provided by [49], basing its analysis on the safe theorem proposed by Heyman.

Collapse behavior: funicular versus stability area methods
Unlike the previously mentioned contributions, the numerical method presented here stands out for its ability to rigorously account for the limited compressive strength of masonry and the stereotomy of the voussoirs.In contrast to the advancements in this area proposed by [24,25], the current approach tackles this challenge through mathematical programming, enabling the management of networks with various topologies, building upon the procedure introduced by [14].A comprehensive description of this funicular method is available in [3], with a specific emphasis on identifying the funicular networks corresponding to the minimum thrust and exploring a range of statically admissible solutions.Unlike the paper just mentioned [3], the contribution proposed in the present work extends the funicular method to the limit analysis of vaulted masonry structures, providing a new specifically developed formulation.The aim is to determine the collapse load multiplier as the maximum of the static multipliers.To achieve this, an ad hoc multi-constrained maximization problem is implemented.
The investigation focuses on arches and domes with a predetermined stereotomy, subjected to selfweight and a vertical point load at the crown, with varying intensity.Equilibrium is represented by a three-dimensional network consisting of branches, which experience axial forces exclusively, and nodes situated along vertical straight lines passing through the centers of gravity of the voussoirs composing the vault.These nodes are constrained to have fixed horizontal coordinates, thereby ensuring a fixed plan projection, and preserving the lines of action of the vertical loads during the optimization process.Local constraints are applied at each joint to consider both the limited compressive strength of masonry and the stereotomy of the voussoirs.
As previously mentioned, this approach also considers the influence of the stereotomy of the voussoirs on the collapse load.As emphasized by [1], stereotomy significantly affects the behavior of a masonry vaulted structure when assuming finite or zero friction coefficients.However, the collapse condition is influenced by the stereotomy even under the assumption of an infinite friction coefficient, given the absence of tensile strength and the presence of finite compressive strength, as the cross-sectional area between the voussoirs varies depending on the stereotomy.
The validation of this numerical funicular method is performed using the aforementioned semi-analytical stability area method [21,22] in its modern version [4,8].This method enables the determination of all statically admissible solutions for symmetric masonry arches and axisymmetric domes of revolution.
The collapse load multiplier is determined by means of an iterative procedure.The analysis of several reference case studies is presented, demonstrating a strong agreement between the two approaches.

Paper outline
The article is organized as follows.Section 2 outlines the procedure employed to formalize the limit analysis optimization problem using the force density method and mathematical programming; it provides an overview of the forces, eccentricities, and constraints at the joints.In Sect. 3 the Durand-Claye method is discussed, focusing on the determination of the collapse load multiplier for symmetric arches and domes.Section 4 presents the results of numerical simulations and a validation of the funicular numerical approach using the two methods; additionally, it explores the impact of a finite compressive strength on the mechanical behavior of arches and domes.Finally, Sect. 5 provides some concluding remarks, along with insights into ongoing research.

Force density method
The equilibrium of funicular networks is managed using the "force density method" [53].A threedimensional network is considered, consisting of m elements (subjected to axial forces only) and n s nodes.The network is situated within a Cartesian reference system Oxyz, z being the vertical axis.The arrays x s , y s , and z s collect the coordinates of the n s nodes.x , y , and z represent the subsets related to the n unre- strained nodes, which are those subjected to external loading.In contrast, x f , y f , and z f pertain to the n f restrained nodes.The topology of the spatial network is described by the connectivity matrix C s , where the submatrix C corresponds to the unrestrained nodes, while C f relates to the restrained ones.The arrays u = C s x s , v = C s y s , and w = C s z s collect the differ- ences in coordinates between the ends of each element along the x, y, and z axes, respectively.The length of each member of the network can be calculated as , and these values are collected in the array l , or equivalently in the matrix L = diag(l) .With s denoting the array that gathers the forces in the m elements, the force density vector, which stores the force-to-length ratios for each member of the spatial network, is expressed as q = L −1 s.
For networks with a fixed plan projection, see e.g.[13][14][15]36], the horizontal equilibrium of the unrestrained nodes can be written as: where x s0 and y s0 are the arrays collecting the pre- scribed x and y coordinates of the nodes, whereas p x and p y gather the nodal forces along the x and y axes, respectively.Equation (1) implies that m − r inde- pendent force densities exist, stored in , being r the rank of the coefficient matrix.The r dependent force densities in ̃ read: where B and d are matrices whose constant entries can be derived by performing Gauss-Jordan elimination on Eq. ( 1). (1) In the context of implementing limit analysis procedures, the array that collects the vertical nodal forces, denoted as p z , can be expressed as: where p zd (respectively, p zl ) represents dead loads (resp., live loads), being a load multiplier.Hence, upon introduction of Q = diag(q) , the equilibrium along the z axis reads: For any given set of , ̃ can be determined using Eqs.(2), and (4) can then be solved to calculate z f .

Limit condition at the joints
The i-th joint, which separates two adjacent blocks, is highlighted in bold in Fig. 1.This section is intersected by the line of action of the funicular force F i at point P i , while C i represents the centroid of the rectangular joint.In Sect.4, non-rectangular sections are addressed by considering the largest rectangle that can be inscribed within the original shape of the joint.The principal axes of inertia of the rectangular section are denoted as i and i , with n i as the normal vector.The dimensions of the rectangle are given by l i, × l i, .Using e x , e y , and e z as unit vectors aligned with the x, y, and z directions, respectively, the expression for the funicular force F i is as follows: The magnitude of the normal component of F i , referred to as N i , can be found as (5) Fig. 1 The i-th joint between two adjacent voussoirs of a dome Vol:. (1234567890) eccentricity with respect to i can be evaluated by computing the moment of N i about the same axis, M i, , and scaling by N i , i.e.: In the above expression, abs (⋅) is the absolute value of the scalar argument, while r i is the vector drawn from C i to any point belonging to the line of action of F i .Similarly, the eccentricity of N i with respect to i reads: see [3] for further details.
According to [4,24,38], when M i, = 0 , the limit state related to any i-th no-tension joint is characterized by the limiting value of the (positive) bending moment, M lim i, .This limit is expressed as a function of the axial force, N i ≤ 0 , the compressive strength,  c < 0 , and the section's dimensions, l i, × l i, : This expression will be used in Sects.3 and 4 for the collapse analysis of the benchmark case studies using Durand-Claye's method (i.e., symmetric masonry arches and domes of revolution loaded with axi-symmetric loads), characterized by unidirectional bending moment M i, and compressive axial force N i .
For the implementation of the funicular method, a more general criterion is required to account for cases of combined bending and compression, see the proposal by [41].According to this condition, a uniform distribution of compressive stress, equal to c , is assumed within a limited region of joint i, which resists the eccentric normal force N i ≤ 0 .This area, with dimensions l i, − 2e i, × l i, − 2e i, , is shown in Fig. 1 around P i .Thus, the limit condition is: which is fully equivalent to (8) when e i, = 0.

Optimization problem
In [14], the funicular analysis of networks with a fixed plan projection was addressed, encompassing general topology and loading.This was achieved by formulating a minimization problem in terms of the independent set of force densities and the vertical coordinates of the restrained nodes.The thrust was used as objective function, whereas constraints on the elevation of the unrestrained nodes were formulated.An analytical investigation into the equilibrium of an arch under vertical loads was preliminary performed using the proposed approach.It was observed that the equations have the same form as those found when optimizing elastic structures.Notably, the coordinates of the unrestrained nodes were found to depend linearly on z f or the inverse of .There exist methods of sequential convex programming whose approximations were especially conceived for such kind of problems [18].This is the case of the Method of Moving Asymptotes (MMA) [55], which will be used for the numerical simulations in Sect.4, see also [16].The proposed minimum thrust problem was endowed with constraints addressing limit conditions at the joints in [14].
Within the framework of the lower bound theorem of limit analysis, a multi-constrained optimization problem is formulated in this contribution to maximize the load multiplier that scales the vector of the nodal live loads p zl , see Eq. (3), as follows: In the above statement, the unknowns include: (i) the load multiplier ; (ii) any reduced set of independent force densities ; (iii) the vertical coordinates of the restrained nodes z f .The arrays z min f and z max f store the lower and upper bounds for the vertical coordinates of the n f restrained nodes, respectively.Equa- tion (10b) facilitates the computation of the set of dependent force densities ̃ from the independent ones.Equation (10c) ensures the equilibrium of the unrestrained nodes in the vertical direction.The occurrence of any positive independent force density is prevented by the side enforcements on ̃ , combined with the local constraints in Eq. (10d).Additionally, Eq. (10e) ensures that crushing is prevented at the i-th no-tension joint, as detailed in Sect.2.2.By setting c → −∞ in Eq. (10e), the optimiza- tion problem in Eq. ( 10) yields funicular networks adhering to Heyman's assumptions.Choosing finite values for c relaxes these assumptions, recognizing that strength is not infinite at all joints.It is worth noting that Eq. (10e) implicitly incorporates the hypothesis of zero tensile strength at the joints.It is finally remarked that Eq. ( 10) can be applied to general types of networks and loading, including horizontal loads, provided that a solution for Eq. ( 1) exists, see [14].However, the application is rigorously valid only for vertical loads, as equilibrium in this case remains unaffected by changes in the elevation of the nodes.The simulations presented next pertain to the case where p x = p y = 0.

Limit analysis using the Durand-Claye's method
Durand-Claye's method [21,22] is employed here to validate the funicular numerical algorithm.
In this paper, the modern version of the stability area method presented in [4,7] is adopted, which explores the stability of masonry arches and domes under the influence of their self-weight.Reference is also made to the further analysis carried out by [8], where the collapse of domes with oculus and lantern is dealt with.Compared to the contributions just mentioned, the current focus is on determining the collapse value of the load multiplier, , through an iterative procedure.This is done by considering more complex stereotomies of the voussoirs and assuming that the inclination of the joints is analytically defined.
In the following sections, symmetric masonry arches and domes of revolution subjected to both their self-weight and the weight of a vertical point load acting downward at the crown are considered.This point load represents a simple case of live load and corresponds to p zl in Eq. (3).For the purposes of the Durand-Claye method applications, this point load is denoted as simply p zl (see Fig. 2).
As for the analysis of masonry domes, the 'slicing technique' is employed.This procedure is theoretically based on historical contributions to the structural assessment of masonry domes, aiming to find statically admissible solutions through membrane theory.Some of these contributions are mentioned in Sect.1.4, such as [23,39,57].Due to the low tensile strength, the dome cracks along its meridians, causing the hooping action along its parallels to vanish.The most significant consequence of meridional cracking is the occurrence of a horizontal thrust at the base of the dome.This thrust results in dividing the cracked dome into a series of 'slices' delimited by vertical planes of symmetry, hereinafter referred to as 'lunes,' each with an amplitude of Δ (as shown in Fig. 3).Each 'lune' behaves as an independent half-arch with a variable width, as described in [34] and [19].For the applications to domes, an absolute Cartesian coordinate system Oxyz is introduced, with O located on the axis of revolution of the dome.The (x, z) plane corresponds to the vertical plane of symmetry of the 'lune' (as illustrated in Fig. 3).Any cross-section (joint i in Fig. 3) is represented as a rectangle with an area of l i, × x C i Δ , where C i is the midpoint of seg- ment D i E i .
For a masonry arch (or a dome's 'lune'), equilibrium conditions are imposed on the portion of the structure between the ideal vertical cross-section and any joint i (as seen in Figs. 2 and 4, right).Symmetry conditions ensure that an eccentric horizontal thrust, denoted as f, acts in correspondence with this vertical section.If the center of pressure is at point P 0 (as illustrated in Fig. 4, right), the eccentricity of f, referred to as e 0 , is measured with respect to the cen- troid of this vertical section, C 0 .
The graphical procedure initially introduced by Durand-Claye enables the identification of the complete set of statically admissible solutions for the entire arch (or dome's 'lune') by visualizing the socalled 'stability area' (Fig. 4, left).This area represents the locus of points in the plane (f , e 0 ) that cor- respond to statically admissible solutions.These solutions are determined by considering both equilibrium and strength capacity of any cross-section.
In the equilibrium equations related to the portion of the structure under consideration, a vertical point load acting at the crown is included, corresponding to p zl ∕n , where n = 2 for the (half) arch and n = 2 ∕Δ for the dome's 'lune'.
Reinterpreting the graphical constructions originally introduced by Durand-Claye in terms of stress resultants (Fig. 4, right) allows us to firmly place this approach within the theoretical framework of limit analysis.The stability area can be determined by considering the formal expressions of the bending moment, M i, , and normal force, N i , expressed as functions of f , e 0 , and i , while fixing a value for the load multiplier , as shown in Eq. ( 11): Then, in order to consider the mechanical properties of masonry, the limit bending moment, M lim i, , at each joint i, according to Eq. ( 8) is imposed as a bound for M i, (Fig. 4, right).Focusing on the arbitrary i-th joint, the inequalities implicitly define the region A i of the plane (f , e 0 ) (see the light yellow region in Fig. 4, left), which represents the statically admissible solutions for the portion of arch (or dome's 'lune') comprised between the vertical crown section and joint i.The procedure must be repeated for all the joints i; the intersection between all the areas A i defines the area of stability A for the entire arch (or dome's 'lune') (see the yellow solid region in Fig. 4, left).
(11) M i, = M i, (f , e 0 , i , ) and For a masonry arch, static admissibility is ensured if the stability area forms a non-vanishing region, such as the yellow solid area in Fig. 4, left.This occurs when the corresponding load multiplier is lower than the maximum of the statically admissible load multipliers.In this case, the values of the minimum and maximum thrusts, denoted as f min and f max , along with their eccentricity at the crown, e 0 , can be determined by the coordinates (f , e 0 ) of the two black points in Fig. 4, left [3].
By progressively increasing the value of , the stability area diminishes until it disappears.Under the assumptions made for describing masonry behavior, this condition is met when the stability area reduces to a single point, resulting in only one statically admissible solution (f , e 0 ) .From a kinematical per- spective, this limit solution corresponds to a collapse mechanism characterized by mutual rotation between the voussoirs or crushing at the critical joints.
Regarding masonry domes, modifications are necessary for Durand-Claye's method.Initially, the stability area is constructed for a single 'lune' treated as an arch of variable width, and the load multiplier is incrementally increased until the stability area vanishes.To determine if this condition signifies the collapse of the dome, the associated mechanism must be analyzed.If a kinematically admissible mechanism is identified not only for the single 'lune' but also for the entire dome, the collapse load multiplier, m,DC , is determined.Otherwise, must be further increased until a kinematically admissible mechanism emerges, as discussed in [7].In Sect. 4 situations where crushing is present will be shown and discussed, providing additional details on these issues.

Numerical examples
Symmetric masonry arches and domes subjected to their self-weight (dead load) and a vertical point load at the crown (live load) are considered.A material with a unit weight of m = 15 kN/m 3 is assumed throughout this section.The point load is specified as × 1 kN, acting in a downward direction.
The optimization approach will be validated using the stability area method, involving an arch (see Sect. 4.1), and two domes (see Sects.4.2 and 4.3), initially neglecting the contribution of hoop forces.Further investigations will be conducted using networks composed of both meridians and parallels.
Collapse load multipliers will be provided along with the computed funicular polygons/networks.In the relevant representations associated with the numerical funicular method, the symbols • and + indicate joint sections where crossing branches activate a strength constraint.The former symbol is used when the branch intersects the cross-section of the joint below the centroid, while the latter symbol is used otherwise.Additionally, the joint numbering ( i = 1, 2, 3, … ) will commence from the joint nearest to the vertical crown section ( i = 1).
It should be noted that, while the funicular method is based on the resolution of an optimization problem aimed at maximizing the statically admissible load multiplier, , to directly obtain the collapse load multiplier, an iterative application of Durand-Claye's method is needed to find the collapse load multiplier: is progressively increased in order to identify the collapse condition, corresponding to the reduction of the stability area to a single point.As regards Durand-Claye's method, it is emphasized that the collapse load can be obtained with the desired precision, since the stability area is delimited by analytical curves in the (f , e 0 ) plane, by resulting in an effective tool for the validation of the numerical approach.
The numerical analyses have been performed by means of appositely developed algorithms implemented in MATLAB [56] (funicular method) and in Mathematica [59] (Durand-Claye's method).

Example 1
The segmental arch with non-conventional stereotomy represented in Fig. 5a is considered.The intrados lies along a circle centered at C in = (0, 0, 0.50) m, of radius r in = 3.50 m.The extra- dos lies along a circle centered at C ex = (0, 0, 0) m, with radius r ex = 4.50 m.Half the angle of embrace- ment, , is 30 • .The out-of-plane thickness of the arch is th = 0.50 m.The arch is made of thirteen voussoirs of amplitude 4.62 • , whose stereotomy is defined by radial lines originating from point C st = (0, 0, −1) m.The polygons used for the funic- ular analysis consist of fourteen branches and fifteen nodes, as depicted in Fig. 5b.According to the Gauss-Jordan elimination performed on the system of equations that govern the horizontal equilibrium, only one independent force density exists.Any of the members of the polygon can be selected in this regard.The numerical simulations use the one depicted in red, i.e. the branch selected by the relevant function providing the reduced row-echelon form of Eq. ( 1).It is remarked that the optimization problem is formulated in four unknowns: the load multiplier, the independent force density, and the vertical coordinates of the two outer nodes.
The analyses are conducted by searching for the collapse load multiplier while assuming various hypotheses regarding compressive strength.This is done to gain a better understanding of the influence of this mechanical parameter on collapse behavior.The values of c considered here range from c = − 1000 MPa, representing a virtually unlimited compressive strength, to c = − 0.5MPa, corresponding to a neg- ligible strength.Within this range, values compatible with the actual compressive strength of the masonry material are considered.
Table 1 presents the collapse load multipliers ( m,DC computed using the stability area method; m obtained through the funicular analysis approach) for different values of c .The graph in Fig. 6 illustrates the relationship between the collapse load multiplier and compressive strength.The results obtained using both methods are in excellent agreement and demonstrate a linear relationship between compressive strength and the collapse load multiplier.
In Fig. 7, two different funicular polygons at incipient collapse are shown, corresponding to c = − 1000MPa (Fig. 7a) and c = − 10MPa (Fig. 7b).The values of the collapse load multiplier are m = 120217.56and m = 1196.93 ,respectively.These results highlight the significant influence of compressive strength on the bearing capacity of the  arch.Conversely, the shape of the funicular polygons remains unchanged as c varies and identifies three critical joints within half of the arch: i = 1 and 7 (symbol + ), and i = 4 (symbol • ).The forces in the branches, however, are scaled based on the collapse multiplier.These observations apply to all the results reported in Table 1 and in Fig. 6.
To validate the funicular approach, the stability area corresponding to the collapse condition for the arch in Example 1 is determined.For instance, when setting c = − 10MPa, the stability area reduces to a single point with a collapse load multiplier of m,DC = 1198.86(as shown in the left part of Fig. 8).The right part of Fig. 8 displays the thrust line at the collapse condition.
A strong agreement is observed between the results obtained using the funicular method and the numerical value of the collapse load multiplier, as well as the positions of critical joints.Critical joints are locations where the bending moment reaches its limit, with joints i = 1 and i = 7 exhibiting a limit positive bending moment, and joint i = 4 exhibiting a negative value.The centers of pressure at these critical joints are marked with + for positive limit bending moments and • for negative limit bending moments.Notably, the influence of compressive strength affects the positions of the centers of pressure at the critical joints, placing them inside the joint instead of precisely at the extrados/intrados boundary (as shown in the right part of Fig. 8).
Practically, the critical joints are identified by analyzing the stability area and the intersections between curves in the (f , e 0 ) plane corresponding to +M lim i, (red curves) or −M lim i, (blue curves).

Example 2
The second benchmark study involves the dome of revolution shown in Fig. 9a.The dome is characterized by a thin depth.The intrados follows a circular path with its center at C in = (0, 0, 0) and a radius of Meridian-only networks for the funicular analysis are represented in Fig. 9b, while general networks comprising both meridians and parallels are illustrated in Fig. 9c.According to the Gauss-Jordan elimination performed on the equations that govern the horizontal equilibrium of the network, the number of independent branches equals the number of meridians minus two, with one additional independent member per parallel.In the above figures, the independent set of branches used in the simulations-retrieved by the function processing Eq. ( 1)-is marked in red.
Table 2 displays the collapse load multipliers obtained using both Durand-Claye's method ( m,DC ) and the funicular method employing meridian-only networks ( m ).Similar to Example 1, various values of compressive strength are considered.The results from both methods show a strong agreement, with only marginal discrepancies.These discrepancies are due to different simplifications of the geometry of blocks and joints, as well as variations in the software used for implementing the methods.The collapse load multiplier values are very close, with a percentage error ranging from 1.4% to 1.5%.
It is worth noting that in applications involving Durand-Claye's method, the cross-section is simplified as described in Sect.3.This involved assuming O ≡ C st as the origin of the absolute Cartesian system (Oxyz), as illustrated in Fig. 3.
Funicular networks corresponding to an incipient collapse condition are obtained for all the compressive strength values listed in Table 2.For instance, we plot the meridian-only funicular networks at incipient  collapse in Fig. 10 for c = − 1000MPa (Fig. 10a) and c = − 2.5MPa (Fig. 10b).The corresponding col- lapse load multipliers are m = 13.90 and m = 13.50, respectively.
Unlike the previous case, the influence of compressive strength is less pronounced due to the distinctive geometry of the dome, characterized by a rather thin profile.For values of c = − 5 , − 10 , − 15 , − 20 , − 1000MPa, the shape of the funicular networks remains consistent and identifies two critical joints where strength constraints are fully activated at convergence: i = 1 (indicated with the symbol + ) and i = 5 (indicated with the symbol • ; as shown in the example in Fig. 10a).From a kinematic perspective, for these compressive strength values, the collapse load multipliers do not identify a collapse mechanism for the dome.However, the constraints in Equation (10e) are nearly fully active for the branches at the base of the dome, as evidenced by the location of the meridians at joint i = 9.
For c = − 0.5, − 1, − 2.5MPa the funicular net- works are, again, identical.The activation of a strength constraint occurs at three joints, i = 1 , i = 9 (symbol + ), and i = 5 (symbol • ), which corresponds, from a kinematical perspective, to a three-hinge collapse mechanism for the dome (see the example of Fig. 10b).
Once again, the forces in the branches for both types of funicular networks described earlier (as seen in Fig. 10a, b) are scaled in accordance with the collapse multiplier.
The positions of the three critical joints identified by the funicular method for c = − 0.5, − 1, − 2.5 MPa coincide with those identified by Durand-Claye's method for all the values of c listed in Table 2. Specifically, at joints i = 1 and i = 9 , the positive bending moment reaches its limit value, +M lim i, , while at joint i = 5 , the limit negative bend- ing moment −M lim i, is attained.As an example, in Fig. 11, the thrust line obtained by setting c = − 10 MPa, corresponding to m,DC = 13.99 , is displayed.
In Fig. 12, the collapse load multipliers, m,DC and m , obtained using Durand-Claye's method and the funicular method, respectively, are plotted against the compressive strength.In this case, the relationship between compressive strength and the collapse load multiplier is nonlinear: the slope of the graphs approaches zero as the compressive strength value increases.This trend in Fig. 12 is notably different from what was observed in Example 1, as shown in Fig. 6.For the dome in Example 2, the collapse load multiplier values are very close to each other for c ≥ −10MPa, which corresponds to values compatible with the actual strength of masonry.For this type of dome, Heyman's approach may be considered suitable.Further confirmation of this observation comes from the results obtained when considering both meridians and parallels: the collapse load multiplier values remain unchanged compared to those associated with the meridian-only funicular networks, as presented in Table 2.At incipient collapse, the parallels are not active.These results validate the 'slicing technique' for analyzing the collapse behavior of this dome type.The vanishing of the stability area identifies a three-hinge collapse mechanism, which is kinematically admissible not only for the single 'lune' but also for the entire dome.Funicular networks employed in this study have the same features of those already introduced when commenting on Fig. 9b, c.The search for the collapse load multiplier initially focuses on networks with branches exclusively along the meridians.Table 3 presents the results for m,DC and m under varying assumptions regarding compressive strength.The results obtained through the funicular method, based on meridian-only networks, consistently align well with those provided by the stability area method.As observed in Example 2, minor discrepancies can be attributed to slightly different geometric assumptions for the dome's layout, resulting in a percentage error ranging from 1.5% to 2 % .Figure 14 illustrates the lin- ear relationship between the collapse load multipliers, m,DC and m , and compressive strength.
It is interesting to note that the shape of the networks remains consistent across different values of compressive strength, specifically, for c = − 15, − 20, − 1000MPa.The forces in the branches scale proportionally to the collapse load multiplier.Figure 15a displays one of these networks with c = − 1000MPa ( m = 91848.22 ).In this case, three critical joints are identified: i = 1 and i = 7 (indicated by + for strength constraint activation), and i = 3 (indicated by • ).These findings are highly con- sistent with the results obtained using Durand-Claye's method.Not only do they match in terms of the numerical value of the collapse load multiplier, but also in the positioning of the critical joints.At joints i = 1 and i = 7 , the positive bending moment reaches its limit value, +M lim i, , while joint i = 3 corresponds to the limit negative bending moment, −M lim i, (refer to Fig. 16, where the thrust line is depicted with c = − 1000MPa and m,DC = 93723.88).
A similar situation arises with c = − 10MPa ( m = 937.31 ),as shown in Fig. 15b.In this case, the shape of the funicular network remains the same as before, and three critical joints are identified.The strength constraint at the joint near the crown, i = 1 , is now denoted by the symbol • : as will be explained below, this corresponds to a condition that is very close to pure crushing at joint i = 1 , where e i, = e i, = 0.
In Fig. 15c, the meridian-only funicular network at the limit condition is presented for c = − 5MPa.Although the shape of the network is the same as that shown in Fig. 15a, b, only two critical joints are present, i = 1 (symbol • ) and i = 3 (symbol •).
With further decreases in the compressive strength, the shape of the meridian-only funicular network at incipient collapse also changes.For example, Fig. 15d displays the network obtained with c = − 0.5MPa, which appears much more arcuated compared to the networks in Fig. 15a-c.A strength constraint is activated at convergence at two joints: i = 1 (symbol • ) and i = 4 (symbol •).
To validate the above results, Durand-Claye's method confirms the identification of two critical joints, both for c = − 5 and c = − 0.5MPa.Moreo- ver, at these low values of compressive strength, crushing is observed at the first joint, i = 1 , where the center of pressure is found to be centroidal.As an example, the thrust line at the limit condition corresponding to c = − 0.5MPa is illustrated in Fig. 17 ( m,DC = 43.01 ).As observed with the funicular net- work, the shape is more arcuated than that obtained with larger compressive strength values.Pure crushing is indicated with the symbol • with a + inside.
In order to better understand the collapse behavior of the dome and the influence of the compressive strength, a second set of analyses is performed, by considering networks with both meridians and  The shape of the funicular networks with meridians and parallels remains consistent for various values of c , including −5 , −10 , −15 , −20 , and −1000MPa, with the forces in the branches being scaled accordingly.As an illustration, in Fig. 18a, one of these networks is presented (specifically, for c = − 10MPa and mp = 985.99 ).At the onset of collapse, the first parallel, located near the crown section, becomes activated.A different behavior is observed when dealing with a very low value of the compressive strength, such as c = − 0.5MPa (resulting in mp = 50.90).In this scenario, all par- allels are activated, causing the dome to undergo compression in both the meridional and hoop directions.This behavior is notably distinct from that exhibited by the meridian-only network.
A summary of the investigations for which funicular polygons and networks have been discussed in Sects.4.1, 4.2, and 4.3 is given in Table 4.

Conclusions
The numerical method proposed in this study, based on the lower bound theorem of limit analysis, seeks to maximize the statically admissible collapse load multiplier by solving a multi-constrained optimization problem.This problem accounts for the finite compressive strength of the material and the stereotomy of the voussoirs by interpreting them as constraints applied at the joints, each inclined according to a geometric criterion.
This study examines an arch and two domes with different geometries under the influence of gravity loads and a vertical point load applied at the crown.When dealing with masonry domes, a first set of analyses focus on meridian-only funicular networks.These structures are modeled as assemblies of 'lunes' using the 'slicing technique'.For this set of analyses, as well as for those related to the arch, the method is validated using the modern version of Durand-Claye's approach, i.e., a semi-analytical graphical technique that visualizes the set of all statically admissible solutions by drawing the stability area.Unlike the funicular method, this approach cannot directly find the collapse load multiplier through maximization, but involves progressively increasing the load multiplier value until the stability area vanishes.The results obtained using these two different methods show strong agreement, not only in terms of the numerical value of the collapse load multiplier but also in identifying the locations of critical joints.
A second set of analyses focuses on modeling domes with networks composed of both meridians and parallels.The results demonstrate the suitability of the 'slicing technique' for the first dome, characterized by a thin profile, and the advantageous impact of hoop forces for the second dome, which features a lowered profile with a non-conventional stereotomy.
The effect of limited compressive strength on the collapse load multiplier and the configuration of the funicular network is examined for each of the structures under investigation.The presence of either a linear or non-linear relationship between the collapse load multiplier and compressive stress is discussed in the context of arch/dome typology and the applicability of Heyman's hypotheses.Of particular interest are the conditions that lead to pure crushing at specific joints.
To better understand the application scope of the two methods in this study-the constrained force density method and the Durand-Claye method-a comparison is insightful, focusing on aspects like model setup time, calculation speed in different scenarios, and the type and quality of results.
The Durand-Claye method, a semi-analytical approach, offers solutions for symmetric masonry arches and domes with symmetrical loads, allowing for analytical shape determination.The model setup requires updates for each case study, typically quick for two-dimensional geometries.In contrast, the  optimization-coupled force density method allows the study of various masonry structures, including unsymmetrical vaults, with generic shapes and load conditions.Capturing three-dimensional behavior increases model setup time, yet parametric geometry can alleviate this.The Durand-Claye method assesses equilibrium through stability area analysis, while the force density method is a numerical approach.Depending on the required output, one method may be more convenient than the other.The Durand-Claye algorithm is efficient for stability assessment, while the force density method, using maximization for collapse load determination, provides a oneshot solution.The fixed plan projection assumption in funicular networks ensures a numerically efficient optimization problem.Computational costs for problems herein are minimal, solved in minutes on a standard desktop.Despite the dependence on network size, the computational burden is expected to be lower than similar structural optimization applications.Notably, sequential convex programming may not guarantee global optimum convergence; however, it yields a safe solution, with optimality checked by the stability area provided by the Durand-Claye method.
Future developments of this research aim to incorporate a limited friction coefficient into the limit analysis of masonry arches and domes, and to consider more complex geometries and load conditions.In regard to the latter point, peculiar attention will be paid to horizontal loads, with the aim of investigating the capacity of arches and domes in earthquake-prone areas.This will be implemented by considering suitable static loads that are equivalent to the seismic action.

Fig. 2
Fig. 2 Scheme of a masonry arch (or dome) subjected to a vertical point load p zl at the crown

Fig. 3 Fig. 4
Fig. 3 Dome of revolution: geometry of one of the 'lunes' of amplitude Δ according to the 'slicing technique'; simplified geometry of the cross section

Fig. 5
Fig. 5 Example 1. Reference section of an arch (a), with nodes and branches of the polygon used for the funicular analysis (b).The independent branch is marked in red.(Color figure online)

Fig. 9
Fig. 9 Example 2. Reference profile of a dome of revolution with angle of embrace = 80 • (a), with nodes and branches of the networks used for the funicular analysis with meridians only (b) either with meridians and parallels (c).The independent branches are marked in red.(Color figure online)

Fig. 10
Fig. 10 Example 2. Meridian-only funicular networks at incipient collapse for a dome with: a c = − 1000 MPa ( m = 13.90 ); b c = − 2.5 MPa ( m = 13.50 ).Colors refer to the magnitude of the branch forces, in kN

4. 3 3 A
Example dome of revolution is addressed, characterized by a flattened profile and a non-conventional stereotomy, as shown in Fig. 13.The intrados follows the curvature of a circle centered at C in = (0, 0, 0.50) m and a radius of r in = 3.50 m.The extrados conforms to a cir- cle with its center at C ex = (0, 0, 0) m and a radius of r ex = 4.25 m.The angle of embracement is = 30 • , which encompasses half of the dome.Excluding the keystone, each of the twenty-four 'lunes' ( Δ = 15 • ) consists of six voussoirs, each with an amplitude of 4.62 • , and their stereotomy is defined by radial lines originating from the point C st = (0, 0, −1)m.

Table 1
Example 1Collapse load multipliers for different assumptions on the strength in compression: m,DC is computed using the stability area method; m refers to the funicular analysis approach c (MPa) Fig. Example 1. Collapse load multiplier versus compressive strength (MPa): funicular analysis and stability area method

Table 2
Example 2Collapse load multipliers for different assumptions on the compressive strength: m,DC is computed using the stability area method, m refers to the funicular analysis approach c (MPa)Vol.: (0123456789)

Table 3
Example 3Collapse load multipliers for different assumptions on the compressive strength: m,DC is computed using the stability area method, m and mp refer to the funicular analysis approach, considering meridian-only and meridian-and-parallel networks, respectively; Δ p = mp ∕ m − 1