3D limit analysis of masonry pavilion domes on octagonal drum subjected to vertical loads

Within the framework of limit design applied to masonry structures, this paper aims at analyzing the different behavior of a pavilion dome according to the adopted construction and reinforcement technologies. By using the static theorem applied to the dome discretized in rigid macro-blocks of variable shape aligned along parallels and meridians, a mathematical model have been constructed in order to search for the load collapse multiplier, and thus to evaluate the degree of structural safety. Then, the associated failure mechanism is represented at the instant in which the collapse is reached. The program that implement the modeling is sufficiently versatile and, in addition to the mechanical characteristics, allows to define the intrados profile, the thickness variability, as well as to insert any window opening in the drum, the lantern at the top and the hoops at each level. The results shown here concern some numerical applications carried out on a theoretical dome, as well as those related to a first approach to the analysis of the dome of Santa Maria del Fiore in Florence by Brunelleschi.


INTRODUCTION
ithin the framework of limit design applied to masonry structures, this paper aims at analyzing the different behavior of a pavilion dome according to the adopted construction and reinforcement technologies; for example, if in correspondence of the ribs the bricks of the contiguous sails are well interlocked (as in the Santa W Maria del Fiore dome) the dome shows a much greater load-bearing capacity compared to the case of interlocking between the bricks completely missing (Santa Maria dell'Umiltà dome in Pistoia). In several previous works [1,2,3,4] the authors had already dealt with rotation domes subject to vertical and horizontal loads; by using the static theorem of the limit analysis applied to the dome discretized in rigid macro-blocks of variable shape aligned along the parallels and the meridians, several optimization models have been implemented in order to search for the load collapse multiplier, thus to evaluate the degree of structural safety. It should be noted that many different modelling approaches in the literature were used to investigate the behaviour of masonry constructions, and particularly of domes. Among those, it is worth mentioning the works by O'Dwyers on the funicular analysis of masonry vaults [5], and those using analytical and finite element models [6,7,8,10,11,12]. Rigid block models based on limit analysis and contact dynamics were also adopted [9,13]. In the previously papers, as in this, the yield domain conditions for the quadrilateral interfaces between the macro-blocks are expressed in terms of stress resultants, in observance of the following mechanical assumptions: inability to carry tension, unlimited or limited compressive strength (in the second case the domain is suitably linearized), frictional sliding with dilatancy. This last hypothesis, even if not realistic, has been assumed because in some cases the influence of failure for sliding is negligible in the kinematic mechanism of the domes; that allows us to base the formulation on simple optimization problems rather than more involved nonlinear analyses, with considerable reduction of computational costs. The yield domain conditions for the quadrilateral interfaces between the macro-blocks are expressed in terms of stress resultants, in observance of the following mechanical assumptions: inability to carry tension, unlimited or limited compressive strength (in the second case the domain is suitably linearized), frictional sliding with dilatancy. This last hypothesis, even if not realistic, has been assumed because in some cases the influence of sliding failure is negligible in the kinematic mechanism of the domes; that allows us to base the formulation on simple optimization problems rather than more involved nonlinear analyses, with considerable reduction of computational costs. It should be pointed out that the sliding failure was taken into account assuming a friction coefficient fc = 0.75. This value is perhaps excessive when compared with those taken by other authors in the most recent literature [14,15,16,17] where, for different materials, coefficients varying from 0.63 to 0.69 are used. However, even when a value included into that range is taken into account, the results -static and kinematic -are almost identical to those which are here presented (at least for the tested cases 7 and 11 of Tab. 1). Furthermore, analyzing the response even in the case of a further reduction in the friction coefficient, it was found that the result remains unchanged up to the value fc = 0.2 -0.22 for which the value of the collapse multiplier decreases and the kinematic mechanism shows an evident failure for sliding. Finally, for fc = 0.18 -0.20 there is a further reduction of the multiplier, but the collapse mechanism is qualitatively equivalent. Compared to previous tackled problems, in this work there are two new aspects: a different type of dome, because it is set on octagonal and not circular drum, and the role played by the interlocking between the bricks in correspondence of the ribs on the bearing capacity of the dome. As we consider vertical loads, we studied a dome segment -corresponding to one sixteenth -between two contiguous meridian symmetry planes (Fig. 1). In this context, we have developed a program through Excel and its Solver to solve the static optimization problem that provides the collapse multiplier. Once the collapse multiplier has been defined, again using Excel, the corresponding kinematic problem is implemented to represent the failure mechanism at the instant in which the collapse is reached. In order to simulate the lack of interlocking between the bricks at the ribs or its presence, an advantageous aspect in the present modeling is that the yield conditions on the related interfaces of the macro-blocks can be either introduced (the failure is possible) or not (the failure is not possible). Although for the considered application cases has always been used a discretization in only eight macro-elements (six for the dome and two for the drum), the program implemented in Excel appears sufficiently versatile and, in addition to the mechanical characteristics, allows to define the intrados profile, the thickness variability, as well as to insert any window opening in the drum, the lantern at the top and the hoops at each level. Some results of application cases are shown and, finally, a first approach was made to the analysis of the Santa Maria del Fiore dome in Florence, built by Brunelleschi. Of course, models with reduced block dimensions could be implemented to take into account failure mechanisms which were not considered in the present work. However, it should be noted that the configuration of the macro-blocks taken into consideration aimed at reproducing the main features of the crack patterns which are frequently observed on masonry domes, such as those investigated in this manuscript. As the authors pointed out in the manuscript, such crack patterns usually involve the ribs and the webs along radial and meridian interfaces which correspond to the ones considered in the model that was presented. However, it is clear that further analysis could be carried out to investigate the effect of interlocking and block size on the obtained failure mechanism.

PRELIMINARY AND GEOMETRIC FEATURES
ince only vertical loads are considered, a dome segment has been studied -corresponding to a sixteenth -between two contiguous meridian planes of symmetry, one passing through the center line of a sail, the other for a rib (Fig.1a). In the horizontal section of the dome, the angle between the two symmetry planes is π/8, and the segment can be studied -at least in this first approach -as fixed at the drum base and subject to dead loads due to their own weight and live loads corresponding to increasing overloads, as well as to the actions on the interfaces belonging to the symmetry planes, which are the redundant unknowns of the static problem. With regard to a dome segment of amplitude = /8 (Fig. 1b), a discretization in rigid blocks is made, each one bounded by two consecutive cross-sections i and i+1 defined by the angles θi and θi+1 respectively with the dome's Z-axis (Fig. 1c). symmetry plane of a sail symmetry plane passing for a rib The generic block like this defined have six faces: two of them -named radial -lying on the aforesaid radial cross-sections i and i+1, two other -named meridian -lying on meridian symmetry planes, and finally two other again lying on the extrados and intrados surfaces of dome respectively (Fig. 2).  Have been used two reference: the fixed one (O, X, Y, Z), to define the vertices coordinates of generic block (Fig.1), and a local one (G f , n, r, t) placed in G f centroid of each contact interface with a contiguous block (Fig.3a), to define the stress resultants interacting through this generic interface, as shown below. Figure 3: (a) Fixed and local references; (b) Own weight P acting in the G centroid of block and stress resultants applied on the i radial face and on j meridian face.

FORCES ACTING ON THE SINGLE BLOCK
ith reference to Fig. 3b, the generic block is subjected to the own weight P applied on its center of gravity G, and to the stress resultants applied on the centroid of each contact interface with a contiguous block, referred to the local reference n, r and t. The resultants on the radial interfaces are six, named normal force Nn (or simply N), shear forces Tt and Tr, twist moment M n , bending moments M t and M r , respectively applied on the interfaces i and i+1. Whereas, on the meridian interfaces j and j+1, lying on symmetry planes, the resultants are only three: the normal force Nn and the two bending moments Mt and Mr . With reference to i and j interfaces, all these forces can be expressed in matrix form by the S i and S j vector respectively: Besides, must be consider the overload acting on the dome -named OL -understood as a covering load applied in the centroid of extrados face (see Fig. 2), that in the present formulation increases through a collapse multiplier α (Fig. 4a). If a lantern is present at the top of dome, also a force PL corresponding to one sixteenth of the lantern own weight will act on the block top ( Fig. 4b), together with the OLL overload related to the small dome of the lantern (also increasing through the α multiplier). W Finally, if hoops are inserted, the possible action of them is represented by a Fh force applied on extrados face, lying in a plane parallel to the dome springing and oriented towards the inside. Having fixed a suitable N h pretension effort, the resultant force Fh = Nh· tan(/8) is positioned at half the height with respect to the height of the generic parallel ring with hoop (point C in Fig. 5).

EQUILIBRIUM CONDITIONS FOR THE SINGLE BLOCK
he equilibrium conditions for the generic block, in the global reference frame X, Y and Z are so formulated:

YIELD DOMAINS FOR THE GENERIC INTERFACE
ith reference to a generic quadrilateral interface (Fig.6), the six stress resultants on generic radial interface i+1 have to respect the yield conditions related to rocking and sliding domains; whereas, on generic interfaces j and j+1 lying on meridian symmetry planes, the stress resultants only are three (the normal force and the two bending moments) and just have to respect the conditions related to rocking domain. In the following, linearized yield domains circumscribing the actual nonlinear ones are shown. T W t G r Figure 6: Generic quadrilateral interface and local reference

Linear rocking yield domains
For the case of unlimited compression strength, a linear yield domain was proposed in [3,4] with reference to the quadrilateral section of Fig. 7; within a coherent kinematic mechanism, the rotation axis has to coincide with one of the four section sides, so we imposed four limit conditions: with d i (Fig.7a) distance between the G centroid and the generic side i (i=1, 2, 3, 4), and the Cartesian expression of the bending moment M (Fig.7b).  As the present paper only consider vertical loads, here we have assumed a limited compression strength of the masonry; in order to include a limit to the compressive stress into domain of Fig. 8, we have added four planes parallel to N-axis of the 3D-reference (Mt, Mr, N), and four other planes forming an opposite pyramid passing through the Q point at abscissa N=N 0 , being N 0 the limit normal force applied in the G centroid which leads the entire section to the collapse. This linear yield domain results circumscribed to the approximate nonlinear one proposed in [18], where an additional quadratic term on N was enclosed to take into account the limited compression strength. We note that the cross section of the two domains coincide at abscissa N=N 0 /2 ( Fig. 9, referred to the (M 1 , M 2 , N) 3D space, being M 1 and M 2 the components of M in the section principal reference). The twelve planes are summarized by the following equations: [5] domain in [18]  In particular, the yield domain normal force N -shear force T is defined by a cone (Fig.11a), with axis coinciding with the Naxis has been opportunely replaced, in this first approach, by a pyramid -circumscribed to the cone -having four faces. Therefore we impose four conditions making reference to the Cartesian components T t and T r of T shear force: Thus, in Fig. 10,  coincide with the angle of friction  0 . Instead, for the yield domain normal force N -twisting moment Mn reference was made to an "equivalent" circular section having radius R equal to the mean of the distances of G by the sides of quadrilateral section (Fig.11b): Therefore two conditions are imposed: being, in the yield domain by sliding of Fig. 10, tanφ = ⅔ R tanφ 0 . In matrix form, the conditions on the generic meridian interface and on radial one can be expressed respectively by:

GOVERNING EQUATIONS OF WHOLE STRUCTURE AND ASSESSMENT OF THE COLLAPSE MULTIPLIER 
amed m the number of balance equations for whole structure and n the number of the unknowns on the interfaces, in matrix form we have: where A is a (mxn) matrix of the coefficients of the unknowns X on the interfaces, listed in the (nx1) vector, while F is the (mx1) vector of the dead loads and of the possible actions of the hoops, and 1 F the (mx1) vector of the live load increasing by the unknown collapse multiplier α N Moreover, named d the number of the yield domain conditions, in matrix form we have: where Y is a (dx1) vector, D is a (dxn) matrix and finally TN is a (dx1) vector of known terms. Now the static theorem of limit analysis can be applied to evaluate the collapse multiplier αby implementing the following problem of mathematical optimization: maximize α subject to: being X and α the unknowns of the problem.

RECONSTRUCTION OF THE COLLAPSE MECHANISM
nce the unknowns α and X have been defined through the static theorem of the limit analysis (paragraph 5), is possible pursue the kinematic problem to identify the corresponding collapse mechanism. The unknowns of problem are: -(mx1) vector u which collects the degrees of freedom, six for every block of the discretized structure; -(nx1) vector Δ which collects the possible displacements between the interfaces (six for every radial interface and three for every meridian one); -(dx1) vector λ which collects the generalized strain rates associated to the yield conditions Y (Eq.12).
These unknowns are bounded by kinematic conditions: (14) and by the flow rule from which, by eliminating Δ, we obtain the compatibility conditions:

being T
A and T D the transposed matrices of A and D respectively (defined in paragraph 5).
the virtual work done by the forces through the associated displacements and by the stress resultants through the mutual displacements between the interfaces, the problem of mathematical optimization is set in the following way: Once solved the problem (17), the kinematic mechanism is represented at the instant in which the collapse is reached through suitable matrices of directional cosines. O

RESULTS
he above described discrete model has been implemented by using Excel and its Solver. Although for the considered application cases has always been used a discretization in only eight macro-elements (six for the dome and two for the drum), the implemented program appears sufficiently versatile and, in addition to the mechanical characteristics, allows to define the intrados profile, the thickness variability, as well as to insert any window opening in the drum, the lantern at the top and the hoops at each level. Once the collapse loads multiplier has been found, through the introduction of appropriate matrices of directional cosines and taking advantage of the symmetry, the program provide the collapse mechanism of a quarter of a dome, represented through axonometric and zenithal views, and a section. In order to test the validity of the procedure, was considered a single dome with the same geometric and mechanical characteristics, and was analyzed the effect on the load-bearing capacity of the following features:  Tab

Some compared cases with relative collapse mechanisms
We refer to compared cases shown in bold in Tab. 2. Only for practical display reasons the comparisons relating to the four features listed above with (a), (b), (c) and (d) will be illustrated below in order (c), (a), (b) and (d).
As previously mentioned, Figs. 13-17 show the collapse mechanism of a quarter of a dome, represented through axonometric and zenithal views, and a section.

Feature (c): absence/presence of window openings in the drum
With reference to the cases 1 and 2 shown in bold in Tab. 1, for the dome analyzed without failure at ribs and without windows (case 1 shown in Fig.13) the multiplier α is equal to about 24, while with a window the value decreases to about 14 (case 2 shown in Fig.14).    Finally, if instead of introducing the hoops (Fig.16), the thickness of the drum is increased from 1.5 to 2.5 meters, the α value increases from about 3 (case 7 shown in Fig.15) to about 13 (case 11 shown in Fig.17).

The case of S. Maria del Fiore
To validate the optimization program developed in Excel, also on a actual dome, a first approach was made to the analysis of the most famous dome on octagonal drum: the one built by Brunelleschi for the Cathedral of Florence between 1420 and 1436, of which numerous researchers have dealt [19][20][21][22].
Obviously -as has already been done by other authors who have carried out various types of analysis -it is assumed that the drum is a regular octagon, and with it also the dome; this allows us to study only one sixteenth of the dome, as for the theoretical one examined above. Actually, compared to the theoretical model assumed, the octagonal dome shows some irregularities concerning geometry and static behavior: -the sides have different lengths each-other, -the drum rests alternately on four solid walls (inclined at 45° with respect to the axis of the nave) and on four ogival arches (two parallel to the axis of the nave and two to that of the transept), with a consequent different crack pattern, also determined from the different structural elements that surround it. Moreover, the architectural reliefs carried out over the centuries present notable differences between them.
To overcome the uncertainties related to the dimensions of some elements that characterize the structure, we have taken the dimensions assumed by M. Como [19] using some rounding done by G. Conti [20].
In the Excel program, which for the dome (excluding the drum) provided a discretization in six elements, an further block was introduced on top, at the base of the lantern, in analogy to what was done by Como [19], while the other 30 voussoirs of its discretization have been incorporated by us, in groups of 5, into our 6 blocks. However, despite using the average value of the unit weight of the brick masonry assumed by him, we have taken into account the different unit weight of the portions of masonry made of sandstone present both in the part at the base of the dome and in the ring on top, named serraglio. Anyway, the weight of the dome (including the lantern) estimated by us is 29,272 t, quite close to the value reported by other scholars. The hypothesis that, despite the octagonal form, Brunelleschi have built a rotational dome (and therefore without the use of a supporting framework), remained secret for 500 years, but in recent times shared by all the most recent scholars of masonry structures, starting by Di Pasquale [21], supports the conviction that in correspondence of the ribs there is a continuity of masonry texture, with good interlocking between the bricks. This has allowed us to introduce the assumption that in the (j+1) faces (see Fig.2) of the blocks of the analyzed segment belonging to the meridian plane passing through the rib there is no failure. Furthermore, having the purpose of making a comparison with the pressure curve shown by Como [19], that starting from the top of the dome involves the whole upper part of the drum, 13 meters high, we had to make a second modeling. In fact, the limitation imposed by the program implemented in Excel -just two blocks to discretize the drum -led us to adopt a different modeling of the oculus (circular hole). However, as for the first modeling, the same horizontal maximum width of 5.8 meters assumed in [19] as average diameter was utilized: modeling 1 (later called Mod1). The oculus has been modeled with a hexagonal shape, but the drum is not considered in all its height, because the basic section of our calculation model was placed in the horizontal plane that cuts in half the aforementioned oculus. To this plan belong the faces (i+1) (Fig.2), i.e. those that have the smallest area and are therefore more prone to crushing failures. modeling 2 (later called Mod2). To be able to consider the whole upper part of the drum, the oculus has been modeled with a lozenge shape having vertical axis equal to the height of the drum (13 m), and the horizontal one equal to 5.8 m.
The mechanical and geometrical characteristics assumed for the dome are listed below in Tabs   With the program implemented in Excel, for Mod1 and Mod2, three different analyses were performed: (1) the classical one already applied to the theoretical dome (see previous sub-paragraph 7.1); (2) a similar analysis but applied to the dome in the current state, taking into account the cracks already present; (3) a comparison between the pressure curve obtained from M.Como and shown in [19] and those obtained with Mod2. More precisely, as regards points (1) and (2), the value of the collapse load multiplier obtained is the same, so this value can provide an indication, even if approximate, about the safety of the actual dome subject to vertical loads. On the other hand, the almost coincident values obtained with Mod2 are justified by the fact that, as for the Mod1, the failure occurs by crushing at the lower area section, corresponding to interface between the two blocks of the drum, belonging to the plane that cuts the oculus in half. Mod2 involves a slight increase in the volume of the oculus with a consequent slight decrease in the weight of the drum blocks, and this leads to the negligible growth of the α multiplier value. Finally, with regard to point (3), for a slice of dome equal to 1/8 of the entire dome, corresponding to a single sail of the dome subjected to its own weight and to that of the lantern on top, in [19] is shown the construction of the pressure curve that, starting from the upper end of the section of the voussoir at the crown, is tangent to the intrados of the dome near the haunches, and intercepts the center of the section at the drum base (Fig.20a). Instead, still using the program implemented in Excel, which analyzes a segment equal to 1/16 of the dome, two similar pressure curves were obtained using two different objectives in the optimization problem: for the first, it was imposed that the multiplier has value zero, in compliance with the condition that the curve intercepts the center of the section at the base of the drum (Fig.20b); for the second, it was imposed that the horizontal thrust transferred from the dome to the drum assumes the value of 2000 kN, in accordance with that of 400 t estimated by Como [19] (Fig.20c).
The results of the analyses referred to in points (1) and (2)    : Comparison between our pressure curves and that of Como [19]. (a) Curve in [19]. (b) Curve obtained by imposing that it intercepts the center of the section at the drum base. (c) Curve obtained by imposing the thrust assumes the value referred in [19].
Obviously, the simplicity of the proposed model, sufficient to successfully compare the answers obtained with those expected based on the experience, does not allow a comparison to be made between the solution obtained for Brunelleschi's dome and that provided by the model illustrated in [13]. The latter, having a more dense discretization, allows the insertion of the vertical elements of the herringbone and to obtain crises due to sliding, even if with very low values of the friction coefficient.

CONCLUSIONS
ithin the framework of limit analysis applied to masonry structures, this paper has aimed at analyzing the different behavior of a pavilion dome according to the adopted construction and reinforcement technologies. By using the static theorem applied to the dome discretized in rigid macro-blocks of variable shape aligned along parallels and meridians, a mathematical model has been constructed in order to search for the load collapse multiplier, so to evaluate the degree of structural safety. Then, the associated failure mechanism is represented at the instant in which the collapse is reached. The Excel program that implement the modeling is sufficiently versatile and, in addition to the mechanical characteristics, allows to define the intrados profile, the thickness variability, as well as to insert any window opening in the drum, the lantern weight on top and hoops at each level. The so far carried out applications, in the case of possible failure at the ribs, have shown the effectiveness of the hoops, that increases as the pre-tensioning stress increases or, alternatively, according to hoops number placed; however, for an reasonable overall force of pre-tensioning, the value of the multiplier obtained is anyway inferior to that computed in the case of a good interlocking between the bricks at the ribs. For the dome discretization into six blocks, the introduction of the hoops also seems to produce better effects if it's concentrated at the fifth ring starting from the top, rather than if distributed over several rings, for the same overall pre-tensioning force. The applications have also confirmed the effect of window openings in the drum: the value of the collapse multiplier decreases as the openings width increases. As further developments, in addition to using the present Excel program to investigate other aspects, such as the problem of the double-shell domes and the search for the minimum thickness of the pavilion domes, the authors are implementing a Matlab program, also eventually reducing the macroblocks size in order to refine the solutions already obtained and to study also the pavilion domes under horizontal loads. Obviously, the aim is to extending these new analyses also to the dome of Santa Maria del Fiore in Florence and compare to this dome the one of Santa Maria dell'Umiltà in Pistoia, apparently similar, in which, instead, cracks at the ribs appeared already during its construction. The comparison between the two different construction methods of the two domes (with and without interlocking between the bricks at ribs), could provide confirmation that, when it is assumed that there are no cracks at ribs, it follows that the two half-segments W in contact along a rib constitute a single block. A similar assumption has been considered also by Foraboschi [22] to analyse the behaviour of the Brunelleschi dome.