The SAVE method provides an empirical procedure to quickly describe the overall seismic response of a building through the assignment of the seismic vulnerability classes. The 1998 European Macroseismic Scale, EMS'98 (Grünthal, 1998) considers five vulnerability classes (see Figure 2), denoted with letters ranging from A to E with decreasing levels of fragility, and assigns to each building the most likely vulnerability class as a function of the vertical structure typology, neglecting the influence of the other typological-structural characteristics (age, floor number, horizontal structure, etc).

The SAVE procedure reduces the level of uncertainty on the attribution of vulnerability classes considering the effects on the building seismic response of additional typological features, defined modifiers, such as the horizontal structural typology, the number of floors, the presence of ties and other features. The method is empirical and built on the extensive PLINIVS database in which typological and damage information of about 250,000 masonry buildings are collected. In the first step, the correlation between the vertical typologies and the level of damage is estimated. The vertical structures defined in the SAVE method are three: generic masonry, denoted with V 0 , characterized by the absence of information on the quality of the wall structure; V 1 , weak and irregular masonry; V 2 , regular and good quality masonry. The levels of damage follow the EMS'98 scale and are: D0, that is, no damage; D1, light, non-structural damage; D2, light structural damage; D3, high structural damage; D4, partial collapse; and D5, global collapse. Referring to the vertical structure V i , S P D V i is the Synthetic Parameter of Damage index, and it is computed as the barycentric abscissa of the damage distribution (Zuccaro eCacace, 2015). As an outcome of this step, three S P D ranges representing the vertical structures are defined (Table 1): class A includes most of the buildings with weak and irregular masonry; class B includes most of the buildings with generic masonry; and class C includes most of the buildings characterized by good quality masonry.

In the second step, the influence of the parameters affecting the average seismic behaviour of the building is estimated through the Synthetic Parameter of Damage (SPD) variation. To this purpose, buildings with vertical structure V i , typology k of the j parameter are classified with the number S P D V i P j k , where j refers to the horizontal structure and k to the floor deformability. The influence of the modifier k on the parameter P j on the vertical structure V i is defined as the difference between the S P D V i P j k and S P D V i . It is essential to point out that some typological parameters that characterise the seismic response of the building could be mutually dependent. For example, the presence of ties within a horizontal structure may be recurrent in the case of deformable floors. When calculating the influence of the presence of tie rods and deformable floors, it is, therefore, possible that the S P D V i P j k value is calculated on the same sample for both parameters. To avoid overestimation when determining the influence of a parameter, it is necessary to define its correlation with the remaining dependent parameters. For parameters independent of the others (e.g., the location of the building in the aggregate), it is not necessary to evaluate correlation factors. After that, once the typological characteristic of a building is known, the corresponding SPD can be calculated as follows: S P D = S P D v + ∑ s = 1 n q s + ∑ j = 1 m ∑ i = 1 m p i + p j ∙ c i j 2 m , (1) where q is the influence of each independent parameter, p is the influence of the dependent parameter, n is the total number of independent parameters, m is the number of dependent parameters, and c i j is the correlation coefficient between the classes of parameters p i and p j .

This Section illustrates the boundary value problem (BVP) at the base of the PRD method referring to a no-tension Heyman material. Specific mathematical restrictions on the stress and latent strain are introduced to model the masonry response. A 2D masonry structure is modelled as a continuum occupying the region Ω of the Euclidean space E 2 . T denotes the stress, u the displacement of the material points x ∈ Ω , while E the infinitesimal strain field assuming small displacements. The Heyman model is enforced through the so-called Normal Rigid No-Tension (NRNT) material restrictions: T ∈ S y m − , E ∈ S y m + , T ∙ E = 0 , (2) with S y m − , S y m + the mutually polar cones of negative and positive semidefinite symmetric tensors. Equations 2 are equivalent to the well-known normality rule: T ∈ S y m − , T − T * ∙ E ≥ 0 , ∀ T * ∈ S y m − , (3)

representing the necessary conditions to rigorously apply the classic theorem of limit analysis to unilateral masonry materials. For more information, the reader is referred to (Kooharian, 1952; Livesley, 1978; Giaquinta and Giusti, 1985; Fortunato et al., 2014; Fortunato et al., 2016; Chiozzi et al., 2017). Therefore the BVP on domains made up of NRNT material reads E = s y m ∇ u , E ∈ S y m + , u = u ¯   o n   ∂ Ω D , div T + b = 0 , T ∈ S y m − , T n = s ¯   o n   ∂ Ω N , T ∙ E = 0 , (4) with n is the unit outward normal to the boundary ∂ Ω partitioned into its constrained ∂ Ω D and loaed ∂ Ω N parts. A BVP solution is a triplet u , E u , T satisfying Equations 4. Moreover, because of Equation 4, u and E are linked each other. For NRNT materials, strain and stress are bounded measures and can be decomposed in the sum of regular . r and a singular . s parts: E = E r + E s , T = T r + T s , (5)

In the present formulation, only singular strains and stresses are considered, meaning that the displacement u and the stress vector s can admit jump discontinuities. The sets K and H of the admissible displacements and stresses are defined as follows: K = u ∈ S   /   E = s y m ∇ u ∈ S y m +   &   u = u ¯   o n   ∂ Ω D , H = T ∈ S ′   /   div T + b = 0 , T ∈ S y m − , T m = s ¯   o n   ∂ Ω N , (6) with S , S ′ two suitable function spaces. For more mathematical details about the regularity of these spaces, the interested reader is referred to (Giaquinta and Giusti, 1985). Thereafter, only potential jumps whose supports are lines lying within the domain are considered. Therefore, stress and strain fields are everywhere zero except for those curves, where they are modelled with Dirac delta lines. The BVP can be solved using two possible variational strategies, i.e., either solving the equilibrium problem through the minimisation of the complementary energy or the kinematic problem by minimising the total potential energy (TPE). The PRD method uses the latter through a straightforward displacement approach, which results in the search for a displacement u ∈ K for which there exist a stress T ∈ H such that T ∙ E u = 0 . The TPE for an NRNT material reads: ℘ u = − ∫ ∂ Ω N s ¯ ∙ u   d s − ∫ Ω b ∙ u   d a , (7) with u ∈ K , i.e. the space of kinematically admissible displacements. This function depends exclusively on the displacement u , and for the case at hand, reduces to the potential energy of the external loads only. Notably, the minimiser u ° , i.e., the element (if unique) or the elements of K on which the potential energy attains the minimum value (if bounded) ℘ u ° = min u ∈ K ℘ u , (8)

implicitly guarantees the equilibrium of the loads imposed on the structure. From a mathematical standpoint, the minimiser u ° can also be not unique, which means that multiple energy solutions to the same BVP can coexist. For detailed information, the interested reader can refer to (Giaquinta and Giusti, 1985) and (Anzellotti, 1985), where the existence of the minimum is proven for Elastic Normal No-Tension materials assuming as kinematically admissible space the set B D Ω and an additional restriction on the given loads, i.e., the safe load condition. For a discussion about the safe load condition and its numerical applications, the reader is referred to the CDF method (Iannuzzo et al., 2021). To the scope of the present paper, it must be noted that the TPE ℘ u is always bounded from below as long as the load is compatible, i.e., the set H is not void. This concept provides a strategy to apply the TPE to find collapse multipliers. As it will be shown, it will also reconnect its use to classic lower-bound approaches.

Indeed, the main idea is to describe the external loads as an affine distribution so they remain proportional to the original self-weight through a scalar parameter. From a physical standpoint, many experimental tests are commonly performed through a tilting test to look at possible collapse mechanisms (DeJong, 2009; Ochsendorf, 2022). The numerical modelling of those tests is straightforward and can be performed through two different strategies. The first one is a pure geometrical strategy, and it consists of directly and fictitiously rotating the original angle geometry of an angle α until the structure collapses. An equivalent strategy to model horizontal actions is to consider an additional distribution of horizontal pseudo-static forces b ⊥ such that the total volume forces are: b λ λ = b g + λ   b g ⊥ , (9) with ∙ ⊥ the linear operator that clockwise rotates a vector of π / 2 , and λ is the scalar load multiplier. This strategy is equivalent to classic pushover analysis. In this regard, the value of the horizontal scalar multiplier represents the gravity loads’ percentage applied on a horizontal plane that drives the structure to its last stable equilibrated configuration. This value λ c is called a pseudo-static horizontal collapse multiplier. The previous formulation refers only to body (mass) loads b . However, it must be extended whenever additional loads or tractions s ¯ depending on the mass are also involved. Because of (9), TPE reads: ℘ u , λ = − ∫ ∂ Ω N s ¯ g ∙ u   d s − ∫ Ω b g ∙ u   d a − λ ( ∫ ∂ Ω N s ¯ g ⊥ ∙ u   d s + ∫ Ω b g ⊥ ∙ u   d a )   . (10)

The last stable configuration, along with the collapse multiplier, can be defined through the following min-max problem: ℘ u ° , λ c = max λ ∈ R + min u ∈ K ℘ u , λ , (11)

Because of the uniqueness of the multiplier due to the normality rule, this min-max problem can be solved following two strategies: the first one translates this in classic limit-analysis approaches, while the second one solves the problem through a sequence of minimisation problems until a solution can still be found. The present contribution follows this last procedure as it directly considers different typological features in the objective function as specified in Section 4.

Following the SAVE method, the most relevant parameters in the vulnerability class assignment are the vertical and the horizontal structures, the presence of ties, the number of floors and the construction age. The first four parameters are used to generate the sample of virtual façades. Conversely, the last one is related to the exposure and cannot be directly considered to construct the analytical-fragility model. The construction age is a fundamental parameter when approaching the problem from an observational standpoint, as it often provides sufficient information about construction technology. When adopting an analytical approach, the construction technology is explicitly modelled, and the importance of including the construction age in the model vanishes. The four SAVE parameters used to scan the variability of the building stock under consideration:

1) vertical structure, which is defined in terms of hollow bricks, tuff, filled brick, irregular stonework, and regular stonework;

From a technological perspective, not all horizontal structures are compatible with tie roads, and thus, parameters 2 and 3 are collected in a unique variable called horizontal technology. Therefore, these four parameters allow the definition of the following six horizontal technology classes: steel, RC, wooden with ties, wooden without tie rods, vaults with tie rods and vaults without tie rods. In addition to the SAVE parameters, other geometrical variables are considered, such as the façade baseline length and thickness, inter-storey height, and openings, whose size and numbers are assumed as variables. The inter-storey height has been fixed in all the analyses at 3.5 m. The façade baseline length ranges from 4.00 m to 7.00 m, which statistically includes most of the 2D masonry geometries catalogued in the PLINIVS database. The possible intermediate heights range from 5.00 m to 6.00 m. The dimensions of the openings have been fixed at 1.2  m x 2.5  m, while the number of openings per floor is:

- at least one opening per floor when the façade baseline length is 4.00 m;

Combining this information, the parameter “length and openings” can be labelled as:

- L4O1: Length 4.00 m Opening 1;

The wall’s thickness has been considered as a function of the number of floors, and its variation over the height goes from 0.60 m to 0.45 m, as detailed in Section 4.1. Table 2 collects all the typological parameters used to define the building stocks. As the reader can note, combining all of them, the number of façades under consideration is 750.

For each of these typologies, the following information has to be defined: the vulnerability class on the base of the SAVE parameters, the collapse multiplier, which will be found through numerical simulations, and the frequency of occurrence in a typological database.

Looking at the vulnerability class assignment, the exploitable parameters for the SAVE method are the vertical structure, the horizontal technology, and the number of floors. Based on these parameters, Table 3; Table 4; Table 5 report the vulnerability class assignments for buildings with 1–2, 3-4 and 5 floors, respectively. These values are used to develop the fragility curves as reported in Section 4.

The approximate solution to the minimum problem (11) is obtained by restricting the search of the minimum in the class K pr of piecewise rigid displacements. This infinite-dimensional space is discretised by considering a finite partition Ω i i ∈ 1,2 , . . , M , (12) of Ω (see Figure 3) into a number M of rigid pieces, such that ∑ i = 1 M P Ω i < ∞ , (13)

P Ω i being the perimeter of Ω i . In particular, restricting to convex polygonal elements, the boundary ∂ Ω i of the n-polygon Ω i , is composed of n segments Γ , of length l , whose unit normal and tangent vectors are denoted n , t , while the endpoints with 0 and 1. The edges shared by two adjacent elements or lying on the constrained boundary are called interfaces and are collected in the set Γ . K pr M represents then the finite-dimensional approximation of K pr generated by this partition. The discretised version of the minimum problem reads: ℘ u ^ = min u ∈ K pr M ℘ u   . (14)

A generic piecewise rigid displacement u ∈ K pr M is in a one-to-one correspondence with the vector U of 3 M components representing the 3 M rigid body parameters. These parameters are restricted by the assumption that the strain must be positive semidefinite. For piecewise rigid displacements, the strain coincides with its singular part, namely: E = E s = v   δ Γ   n ⊗ n + 1 2 w   δ Γ   t ⊗ n + n ⊗ t , (15)

E s being concentrated along the interfaces among blocks, that is, within the present approximation, over the set Γ . From (4), it follows that sliding is not allowed on any interface and only detachment is possible: v = u ∙ n ≥ 0 , w = u ∙ t = 0 , (16)

In Equation 16, u is the retive displacement between two element nodes of the discretisation. Notice that conditions (16), derived from the normality assumption, represent a condition of unilateral contact with no sliding among blocks. With N the number of the interfaces Γ , and v 0 , v 1 , w 0 , w 1 the normal and tangential components of the relative displacements on the endpoints 0, 1 of any segments belonging to Γ , restrictions (16) are equivalent to the 2 N inequalities v 0 ≥ 0 , v 1 ≥ 0 , (17) and the 2 N equalities w 0 = 0 , w 1 = 0   . (18)

Assuming homogeneous boundary conditions, restrictions (16) can be expressed in matrix forms as a function of the vector U collecting the unknown Lagrangian parameters and can be rewritten: A   U ≥ 0 , B   U = 0   . (19) Finally, the discretised version of the minimum problem (14), which approximates the minimum continuum problem (8), reads: ℘ U 0 = min U ∈ K M − b + λ   b ⊥   U , (20) where the objective function representing the potential energy of the external body loads b is minimised within the set: K M = U ∈ R 3 M   /   A   U ≥ 0 , B   U = 0   . (21)

The limit analysis problem can be then solved through a sequence of LP problems as:t ℘ U 0 , λ c = max λ ∈ R + min U ∈ K M ℘ U , λ

The minimisation problem (20) transforms the original minimisation problem (8) for a continuum, into a minimisation problem for a structure composed of rigid parts, acted on by given loads and given settlements and subject to unilateral contact conditions along the interfaces. Problem (20) is a standard linear finite-dimensional minimisation problem, as both objective function and constraints are linear. The existence of the solution of this approximate problem is trivially guaranteed if the original problem is bounded from below. For small problems, it can be solved exactly with the simplex method (Dantzig et al., 1955), while for large problems, interior-point algorithms represent efficient and fast alternatives (Mehrotra, 1992; Vanderbei, 2015; Dantzig, 2016).

The present Section gives an overview of the procedure adopted to parametrically model geometries and mechanical features of the buildings’ stock extracted from the PLINIVS database. The 2D analysis of a building is carried out by referring to its main façade. For example, looking at three-story buildings, Figures 4A, B describe the geometries and loading conditions of two facades having one and two openings per floor, respectively.

The inter-story height h is assumed to equal 3.5 m for all the façades independently of the floor number, while the openings’ dimensions, as mentioned above, are fixed at 1.2 m x 2.5 m. The façade’s depth d w varies over the building’s height H b according to the following relation d w H b = 0.60   m i f   0 ≤ H b ≤ 2 h 0.50   m i f   2 h ≤ H b ≤ 4 h 0.45   m i f   4 h ≤ H b ≤ 5 h   . (22)

Table 6 reports geometric data needed to generate a masonry building of the sample referring to the symbols introduced in Figure 4. It is worth emphasising that these geometrical features as well as the openings number, do not depend on the building height but only on the building baseline length. Indeed, structures with baseline lengths equal to 4 m show only one opening per floor; the ones whose length is equal to 6 m or 6.6 m, show two openings per floor; and façades with a baseline length equal to 5 m can have one or two openings per floor as illustrated in Table 6.

Regarding the loads, mass density values depend on the masonry typology as reported in Table 7.

The presence of slabs and the corresponding transmitted live loads are modelled through additional linear loads placed as in Figure 4, and whose gross values are reported in Table 8. The liner loads adopted in the numerical analyses can be derived by multiplying these values by the slab’s net orthogonal length, assumed to equal 5 m.

Lastly, the presence of tie roads depends on the number of floors and, according to the PLINIVS database, whenever they are present (i.e., for wooden and vaulted slabs only), one assumes:

- One-, two- and three-story façade has only one tie-rod at the uppermost level;

Buildings with reinforced concrete and steel slabs are assumed not to show any tire rod systems. The tie-rod strength ( T r ), i.e., the maximum stress it can withstand depends on the strength of the cable ( T c ) and of the anchor plate ( T a ) as well as on the compressive resistance of the masonry ( C a ), as from the following relation: T r = min T c , T a , C a   . (23)

It is assumed that the square anchor plate edge length is 0.30 m, the cable has a net diameter of 30 mm, and it is made up of steel with a yield tension of 240 MPa. Moreover, the calculation assumes a confidence factor of 1.35 according to the Italian national code (Ministero delle Infrastrutture e dei Trasporti, 2018) and the masonry compressive strength values as detailed in Table 9.

Based on these assumptions, once the masonry typology has been fixed, Eq. 23 yields the maximum tensile force that the tie-rod can sustain and that will be used in energy minimisation, as detailed in the next Section.

The minimum energy problem is formulated as an LP problem whose most general version reads: min − b g + λ   b g ⊥ ∙ U − B q s   q s ∙ U − λ   Q s ⊥ ∙ U − B F t   F t ∙ ∆ U (24) A s   U = 0 A n   U ≥ 0 A t   U = 0 with b g the vector modelling volume forces as lumped to the blocks’ centroid, q s the vector representing the live loads transmitted by the slabs and translated to the blocks’ centroid by the operator B q . It is worth pointing out that Q s ⊥ is the vector representing the seismic effects of the slabs’ loads, as detailed below; B F t F t r is the vector collecting the forces exerted by the tie rods and translated to the block’s centroids closest to the tie-rod endpoints, while ∆ U collects the unknown relative displacements between these blocks. The term B F t F t r is implemented only on digital models with deformable floors, such as wooden or vaulted floor systems, when tired rods are explicitly considered. It should be mentioned that the horizontal rotated component of the loads Q s ⊥ is modelled following two strategies depending on the slab’s typology Figure 5. For deformable floors (i.e., vaulted and wooden floor systems), it results: Q s ⊥ = B q s   q s ⊥ (25) while for rigid slabs (i.e., steel and RC): Q s ⊥ = B F q s ⊥   F q s ⊥ = B F q s ⊥   L q s ⊥   q s ⊥ (26) where F q s ⊥ = L q s ⊥   q s ⊥ is the force vector that lumps the linear horizontal load q s ⊥ through the operator L q s ⊥ and is concentrated on the blocks’ centroids placed on the left side of the building through the operator B F q s ⊥ . Looking at the constraints, the first enforces homogeneous boundary conditions, the second inequality the non-overlapping relation and the third one the no-sliding conditions.

Referring to three-story masonry buildings, Figure 6 shows five digital models and the corresponding discretisations for different façade baseline lengths. Particularly, the number of elements ranges from 6,194 for a baseline of 4 m, to 9,862 for a baseline of 6.6 m.

Table 10 shows the number of elements used to discretise the digital models as a function of the building’s baseline length and floor number. The robustness of the proposed discretisation against two different numerical approaches is validated here (Iannuzzo et al., 2021).

The computational time required to solve a single optimisation problem is 0.15 s using an agile laptop with an AMD Ryzen 7 5700U and using Mosek (Mosek, 2010) as a solver.

PRD results, obtained in terms of load multipliers, are translated into peak ground acceleration (PGA) according to the procedure detailed in (Ministero delle Infrastrutture e dei Trasporti, 2018), as P G A = λ   q S   g (27) with q the ductility factor, here assumed as q = 2 , S the subsoil factor, fixed to 1.25, and g = 9.81 m/s² the gravitational acceleration. Results are then depicted in Figure 7 in terms of Damage State Probability (DSP) as a function of the three vulnerability classes assigned by the SAVE procedure. The DSP index represents the cumulative percentage of all structures whose collapse is activated by a given PGA value.

The PRD results obtained are then approximated through lognormal curves F x = 1 2 π σ 2 ∫ 0 x 1 s e − ln ⁡ s − μ 2 2 σ 2 d s , (28) and the corresponding smooth curves are represented in Figure 8, again as a function of the three vulnerability classes. Equation 28 is the cumulative of a normal distribution with average value μ and standard deviation σ . In the present research, the three vulnerability classes are defined by the following values: − C l a s s   A : μ = 0.234 , σ = 0.631 − C l a s s   B : μ = 0.486 , σ = 0.349 − C l a s s   C : μ = 0.729 , σ = 0.235

In Figure 9, the PRD-based fragility curves are compared with the ones proposed in (Zuccaro et al., 2017), which evaluate the minimum acceleration triggering the first damage mechanism selected among predetermined mechanisms (including also out-of-plane ones) on a sample of buildings from the PLINIVS database. In this sense, it is worth noting that the PRD-based curves are associated with higher PGA levels since the curves evaluated in (Zuccaro et al., 2017) are strongly influenced by the out-of-plane mechanisms occuring at lower levels of PGA.

Restricting to vulnerability class B, shows the PRD-based fragility curve together with the ones proposed in (Cattari et al., 2014) for different damage levels. As expected, the numerical-based curve is in good agreement with the damage level D5, corresponding to the global collapse of the structure. Indeed, as applied, the PRD method can only detect the mechanism driving the structure to collapse without the possibility of accounting for lower damage levels.