A Tool for the Rapid Seismic Assessment of Historic Masonry Structures Based on Limit Analysis Optimisation and Rocking Dynamics

Featured Application

Rapid seismic assessment of historic masonry structures.

Abstract

This paper presents a user-friendly, CAD-interfaced methodology for the rapid seismic assessment of historic masonry structures. The proposed multi-level procedure consists of a two-step analysis that combines upper bound limit analysis with non-linear dynamic (rocking) analysis to solve for seismic collapse in a computationally-efficient manner. In the first step, the failure mechanisms are defined by means of parameterization of the failure surfaces. Hence, the upper bound limit theorem of the limit analysis, coupled with a heuristic solver, is subsequently adopted to search for the load multiplier’s minimum value and the macro-block geometry. In the second step, the kinematic constants defining the rocking equation of motion are automatically computed for the refined macro-block model, which can be solved for representative time-histories. The proposed methodology has been entirely integrated in the user-friendly visual programming environment offered by Rhinoceros3D + Grasshopper, allowing it to be used by students, researchers and practicing structural engineers. Unlike time-consuming advanced methods of analysis, the proposed method allows users to perform a seismic assessment of masonry buildings in a rapid and computationally-efficient manner. Such an approach is particularly useful for territorial scale vulnerability analysis (e.g., for risk assessment and mitigation historic city centres) or as post-seismic event response (when the safety and stability of a large number of buildings need to be assessed with limited resources). The capabilities of the tool are demonstrated by comparing its predictions with those arising from the literature as well as from code-based assessment methods for three case studies.

1. Introduction

In the last decade, impressive advancements have been made with respect to the preservation of the structural integrity of historic masonry structures [1].

In the literature, there are several advanced methods of analysis that can be implemented in both numerical methods or analytical tools [2]. The majority of such studies address the structural analysis and safety problem of masonry structures by using numerical models based either on the Finite Element Method (FEM) [3,4,5] or Discrete Element Method (DEM) [6].

In the case of FE models, the masonry is typically modelled following either a continuum approach (designated as macro-modelling) [7] or discretizing explicitly both units and joints (designated as micro-modelling). In the case of DE models, the masonry is typically modelled as the assemblage of rigid blocks inter-connected through interfaces of a given stiffness, which makes it possible to capture their large displacement response and the opening and closing of joints [6]. It is worth noting that both of these methods require a large set of input parameters [8] and can be time-consuming and computationally expensive in the case of modelling failure, which becomes even more relevant in the case of dynamic analysis.

Hence, powerful analytical tools based on the theorems of limit analysis can be an appropriate alternative. The main advantage of limit analysis approaches is that they require less knowledge of the mechanical behavior of materials and are, therefore, more practical.

Specifically, the Upper Bound theorem of limit analysis is a useful tool for the rapid assessment of historic masonry structures [9]. However, the computation of the load multiplier depends on the macroblocks’ geometry that strongly influences the assessment of their structural integrity. Therefore, multiple (theoretically infinite) failure mechanisms need to be considered to evaluate the minimum of the kinematically compatible load multipliers [10]. Some researchers proposed using optimization routines that solve the minimization problem that is constrained under specific hypotheses [11].

Casapulla et al. [12,13] proposed a macro-block model coupled with a simplified procedure for predicting the collapse load and the failure mechanism of in-plane loaded masonry walls with non-associative frictional contact interfaces. The same authors [14] have recently revisited the previous macro-block approach implementing the frictional resistance computation.

Funari et al. [14] implemented a novel two-step analysis method in a visual programming environment, which can manage the data arising from both non-linear static or dynamic analysis to detect the most likely collapse mechanism through the Control Surface Method (CSM). In this framework, the macro-block geometry’s parameterization allows for the exploration of the domain of possible solutions using the upper bound method of limit analysis coupled with an optimization tool based on Genetic Algorithms.

Due to their relatively low computational cost, upper bound limit-analysis procedures have been implemented in building codes [15,16] in order to provide rapid verifications of the stability of masonry structures subjected to earthquakes. Such approaches include both strength and displacement based procedures, which, while rapid, tend to incorporate limited dynamic effects and thus yield predictions that are generally on the more conservative side [17] and that can consequently result in expensive, and occasionally unnecessary, retrofitting solutions. As an alternative, non-linear dynamic analysis, whereby the equations of motion of the macro-blocks are directly integrated, could be used to more accurately model the dynamic behaviour of these masonry structural assemblies. Specifically, by modelling the macro-block as a SDOF system undergoing rocking motion, better agreement has been found with experimental results, with substantially less scatter presented as compared to the corresponding code-based predictions [17].

Equations of the motion for single rigid prismatic blocks rocking on rigid interfaces (assuming no bouncing or sliding) were first derived by Housner [18]. Since then, equations of motion have been derived for more complex geometries such as masonry spires [19] and for more complex mechanisms such as facades subjected to additional masses due to floors and roofs, thrusts from vaults, and the restraining influence of tie-rods [20,21,22,23]. Equations have also been derived for mechanisms with two and three elements in the kinematic chain, including cracked wall sections [23], arches [22], and symmetric [23,24] and asymmetric [25] portal frames. Furthermore, through linearization of the equations of motion about the point of unstable equilibrium, local dynamic equivalence can be derived between more complex rocking systems and the single rocking block [25]. Consequently, the dynamic response of the different mechanisms (single, two and multi-block) can be approximated using the same general linearized equation of motion.

Still, derivation of these equations of motion can be challenging, especially for structures with complex geometries, or mechanisms involving more than one element in the kinematic chain. To that end, a CAD-interfaced analytical tool for the non-linear dynamic analysis of masonry collapse mechanisms was developed by Mehrotra and DeJong [26]. The tool can directly derive and solve the equation of the motion for any user-defined structural geometry. The broad applicability of this tool was demonstrated by the analysis of a variety of masonry structures such as churches [26] and monumental houses [27]. However, to perform the dynamic rocking analysis following this methodology, the exact collapse mechanism needs to be defined a priori, which requires in-field surveys/site inspections and ultimately, depends extensively on user experience and engineering judgment.

The main aim of this paper is to develop a novel, multi-level integrated modelling procedure that uses a combination of upper bound limit analysis [14,28] and non-linear dynamic (rocking) analysis [26] for the seismic collapse assessment of any user-defined structural configuration. In the first step, parametric modelling of the macro-block geometry is conducted, enabling exploration of the domain of possible solutions using the upper bound method of limit analysis. A heuristic solver based on the Nelder-Mead method [29] is then used to refine the geometry of the macro-blocks and search for the minimum value of the load multiplier. Once the macro-block (i.e., collapse mechanism) has been defined, the digital tool then computes the kinematic constants defining the corresponding (rocking) equation of motion, which can be solved for full time-histories. The response of the structure is provided both in terms of the full time-history response as well as in the form of the maximum predicted rotation. An overview of the proposed analysis procedure can be found in Figure 1.

The proposed methodology has been entirely integrated into a user-friendly visual programming environment that allows the user to easily connect data from different sources while keeping a clear understanding of the relationships between them thanks to the flowchart-like representation of the different components of the code. The procedure is made available using the environment offered by Rhinoceros3D + Grasshopper [30], which is well known by students, researchers, and practicing structural engineers. Unlike other more time-consuming advanced methods of analysis, the proposed method allows the users to perform a seismic assessment of masonry buildings in a rapid and computationally-efficient manner, while simultaneously providing more accurate predictions than simplified/code-based methods. Such an approach is particularly useful for territorial scale vulnerability analysis (e.g., for risk assessment and mitigation in historic city centres) or as post-seismic event response (when the safety and stability of a large number of buildings need to be assessed with limited resources).

The paper is organized as follows. The theoretical background of the proposed methodology is described in Section 2, while Section 3 reports its numerical implementation. The three benchmark case studies conducted using the proposed methodology, including the code-based assessment according to the Italian building code are presented in Section 4 and Section 5. The relevant conclusions are finally discussed in Section 6.

2. Methodology

In this section, the theoretical background of the proposed multi-level approach is described. First, the geometry of the collapsed macro-block is obtained by using a heuristic solver coupled with the Upper Bound theorem of limit analysis. Subsequently, the macro-block is considered as geometry data input to compute the kinematic constants defining the corresponding rocking equation of motion.

2.1. Upper-Bound Limit Analysis

The Upper bound limit analysis proposed is theoretically based on the mechanical model developed in Reference [31], which has been implemented into a visual environment in Reference [14]. According to Heyman’s no-tension model, masonry is considered to have an infinite compression strength and no capability to resist tension [9]. Load multipliers are computed, taking into account the friction at contact interfaces. The load multiplier $\lambda$, leading to the collapsed macro-brock, is obtained by using the principle of virtual work:

$$\lambda {\displaystyle \sum }_{i=1}^n{W}_i{\delta }_{\left(x,y\right),i}-{\displaystyle \sum }_{i=1}^n{W}_i{\delta }_{\left(z\right),i}-{\displaystyle \sum }_{s=1}^n{F}_{real,s}{\delta }_{\left(x,y\right),s}=0$$

(1)

where $W_i$ are the inertia forces arising from floors and roofs as well as the self-weight of the masonry walls, ${\delta }_{\left(x,y,z\right)}$ are the virtual displacements of the point in which the forces are applied, $F_{real,s}$ are the real frictional forces computed as a weighted value as a function of the inclinations of the crack line [32], i.e.,:

$$F_{real,s}=F_{max,s}\left(1-\frac{{\alpha }_c}{{\alpha }_b}\right)$$

(2)

where ${\alpha }_c$ is the actual crack inclination ${\alpha }_b$ is the crack inclination upper threshold (which depends by the geometry of the block [31]) and $F_{max,s}$ is the frictional force computed under the hypothesis of maximum frictional i.e.,:

$$F_{max,s}=W_{i,s}\left({\alpha }_b\right)f_c$$

(3)

where $W_{i,s}\left({\alpha }_b\right)$ is the weight of the macroblock OAB (Figure 2) and $f_c$ is the friction coefficient.

The minimum load factor capable of activating the failure mechanism is computed through an optimization routine that varies the inclination of the lines that simulate the cracked surfaces. As reported in Reference [31], the crack inclination, which represents the variable of the optimization procedure, must respect the conditions of geometric compatibility depending on the dimensions of the blocks. Hence, $\lambda$ is obtained by solving the following constrained minimization problem.

$$\{\min \lambda :\frac{l}{2nh}\le \tan {\alpha }_c\le \tan {\alpha }_c$$

(4)

where $l,h,n$ are the brick length, the brick height and the number of vertical brick rows, respectively.

2.2. Non-Linear Dynamic (Rocking) Analysis

Once the geometry of the collapsing macro-block is defined, the rocking dynamic analysis can be performed. Following the approach presented in Reference [20], the rocking equation of motion assumes the following general linearised form:

$$\tilde{I}\ddot{\varphi }-\tilde{K}\left(\varphi -{\varphi }_{cr}\right)=-\tilde{B}{\ddot{u}}_g+\tilde{M}$$

(5)

where $\tilde{I}$ is the moment of inertia of the macro-block computed about the axis of rotation (hinge), $\tilde{K}$ is the rotational stiffness of the structure, $\varphi$ is the rotation, while ${\varphi }_{cr}$ is the critical (unstable equilibrium) rotation of the macro-block, about which the equation of motion is linearised in order to obtain dynamic equivalence with the single rocking block. Additionally, $\tilde{M}$ is the moment due to external static forces (if any), while $\tilde{B}{\ddot{u}}_g$ is the moment provided by the ground motion.

Through the following transformation of variables (where g is the acceleration due to gravity):

$$\theta =\varphi \frac{\tilde{K}}{g\tilde{B}}$$

(6)

Equation (5) can be re-written as:

$$\ddot{\theta }=p_{eq}^2\left(\theta -\tilde{\lambda }-\frac{{\ddot{u}}_g}{g}\right)$$

(7)

where p_eq is the equivalent rocking frequency parameter and $\tilde{\lambda }$ is an approximation of the load multiplier that activates the mechanism (i.e., not explicitly taking into account frictional forces) as opposed to the minimum frictional load multiplier $\lambda$ as defined in the previous sub-section. In this work, rocking is defined to initiate when ${\ddot{u}}_g>g\tilde{\lambda }$.

Note that both p_eq and $\tilde{\lambda }$ are the equivalent rocking parameters defining the rocking equation of motion, and depend on the kinematic constants $\tilde{I}, \tilde{K}, \tilde{B}, \tilde{M}$ and ${\varphi }_{cr}$, which in turn are dependent on the geometry of the macro-block(s) and the type of mechanism taking place (i.e., single/two/multiple block). For more detailed (i.e., mechanism-specific) expressions of these terms, the reader is directed to Reference [26].

Furthermore, the rotation ${\theta }_{ov}$ upon the exceedance of which the macro-block will overturn and collapse also needs to be computed. As in the case of p_eq and $\tilde{\lambda }$, this rotation depends only on the kinematic constants:

$${\theta }_{ov}=\left({\varphi }_{cr}\tilde{K}-\tilde{M}\right)\frac{1}{g\tilde{B}}$$

(8)

Numerical integration of the equation of motion is then conducted using the explicit Runge-Kutta method of order 5(4), whereby the error is controlled assuming the accuracy of the fourth-order method, while steps are taken using the formula with fifth-order accuracy. The solution procedure is iterative: starting from a given set of initial conditions, the algorithm computes the rotation and angular velocity at each time-step, which are subsequently used as input (i.e., the initial conditions) for the following time-step. In the case of impact, the energy dissipated by the block(s), which results in a reduction in angular velocity, is accounted for through the coefficient of restitution η_COR, which depends both on the block geometry as well as on the type of rocking (i.e., one-sided or two-sided). More detailed expressions for this term can also be found in Reference [26].

3. Numerical Implementation

The multi-level analysis described above is implemented within the visual programming environment offered by Rhinoceros3D + Grasshopper. The proposed method is compartmentalised into clusters, where each group performs logical functions.

The first step of the visual script is devoted to importing the geometry of the structure that the user wants to investigate (BLOCK 1, Figure 3). This step of the procedure may be performed by following different approaches, i.e., parametric modelling of the structure, importing of the NURBS geometry previously modelled in Rhino3D or importing of the geometry directly using laser scanning. The same group of the visual script contains a series of components devoted to the setting of parameters required: friction coefficient, dimensions of the blocks, specific weight and the position of the rotational hinge (whereas its height with respect to the ground is assumed as a variable). In the same group, the user has to define the geometry of the openings, if present.

The second group is formed by a series of Grasshopper components able to split the structure by means of crack lines which define the geometry of the collapsing macro-block (BLOCK 2, Figure 3). This group contains a series of sliders that allow the visual program to cover a broad range of macro-block geometries according to the geometrical constraints defined in Equation (3). These sliders correspond to unknown variables in the optimization procedure.

The subsequent group (BLOCK 3, Figure 3) contains an ad-hoc programmed C# component, which formulates the principle of virtual work of the problem. The output of this cluster is the calculated load multiplier (BLOCK 4, Figure 3). In this work, to solve the minimisation problem reported in Equation (4), the Nelder–Mead method is adopted [33]. The Nelder-Mead method (NMM) is a commonly applied numerical method used to find the minimum or maximum of an objective function in a multidimensional space. The NMM component (BLOCK 5, Figure 3) is linked to the parameters that define the geometry of the collapsing macroblocks (inclination of the cracks) and the load multiplier that has to be optimised.

Once the geometry of the macro-block (i.e., the collapse mechanism) and its associated minimum load multiplier λ have been obtained using the optimisation procedure described, a GhPython script (the Python interpreter component for Grasshopper) follows (BLOCK 6, Figure 3). This component takes as input the macro-block geometry and hinge location, and automatically extracts the required geometric properties (such as volume, centroid and moment of inertia) defining the kinematic constants, which, in turn, are used to compute the parameters defining the rocking equation of motion. For a more detailed description of how the script works, the readers are directed to References [14,26].

The rocking equation of motion can then be solved for full time-histories, (BLOCK 7, Figure 3) for both harmonic and random, i.e., earthquake, ground motions. Note that in the case of the latter, the corresponding ground motion records would need to be provided as user-defined input (BLOCK 8, Figure 2).

The solution of the equation of motion requires additional Python packages such as NumPy [34] and SciPy [35], which are not available in IronPython, which is the implementation of Python used by Rhino and by extension, Grasshopper. Thus, to make the proposed modelling strategy work within the Grasshopper environment, a GH Python Remote component [36] was added to the cluster (BLOCK 7), as it provides a connection to an external instance of Python through which the required NumPy and SciPy functions can be executed.

4. Validation of the Upper-Bound Limit Analysis Procedure

In this section, the proposed methodology is evaluated using three case studies. The first case study, which involves a masonry wall without openings [31], is used to validate the upper bound limit analysis procedure implemented in the visual program-interface. The second and third case studies involve a masonry wall with openings and a church geometry inspired by that of San Nicolò di Capodimonte (Camogli, Genova, Italy).

4.1. Solid Masonry Wall Loaded in-Plane

A masonry wall loaded in its plane is first considered to validate the suitability of the proposed procedure. The geometry of the wall, as well as the loading configuration, are taken in agreement with those reported in Reference [31]. In the present example, two kinds of blocks are considered. Figure 4a,b show the masonry pattern for m = 1 and m = 0.5, respectively (where m is the unit aspect ratio of the blocks [31]). A friction coefficient equal to 0.75 is considered. Figure 5a shows the optimisation process, which is based on 20 iterations of the Nelder- Method-Optimisation component [37]. During the first generation of the solution, the solver starts randomly. At each succeeding generation, the result tends to converge to the value that minimises the load multiplier under the constraints imposed in Equation (3).

Considering a numerical approximation of two digits, the solver converges, after few iterations, in 2 s with an Intel^® Core™ i7–6700HQ processor. In the case of m = 1, the proposed solution scheme computes a load multiplier equal to 0.25, whereas the micro-model approach [38] computed a value of 0.27. When m = 0.5, the NMM solver converges to a value of the multiplier equal to 0.47, which is slightly lower than the one arising from [38] (0.49). It is worth noting that in both the analysed cases, the load multipliers computed using the proposed macro-block approach are slightly lower than those arising from the micro-block approach. These conservative results are in agreement with the concept of the macro-block approach, which allows users to get an estimation of the seismic safety of the structures while significantly reducing the computational cost.

The same considerations can be made from the point of view of the geometry of the failure mechanism, see Figure 5b. In particular, when m = 1, the collapsed macro-block is unable to mobilise any friction forces ($F_{real,s}\to 0$). Instead, when m = 0.5, the collapsed macro-block is characterised by an angle placed at an intermediate position among the three different cracks detected within the micro-block model [31,38].

4.2. Masonry Wall with Openings Loaded In-Plane

The second case study aims to investigate the capabilities of the proposed model on a masonry wall with openings. The geometry, the boundary conditions, and the mechanical parameters are taken in agreement with Reference [39], which integrated a computer procedure to determine the collapse using a non-standard limit analysis where the masonry was modelled as a discrete system of rigid blocks. The masonry wall is subjected to a gravity load and a horizontal uniform lateral load as a factor of the unit weight. The structural model, as well as the failure mechanism obtained by using the non-associative friction solution in Reference [39], are reported in Figure 6. As shown in Figure 6b, the crack pattern is mainly featured by the overturning of an assemblage of rigid blocks which assume the structural behavior of a macro-block.

Based on both engineering judgment and in-situ observations after earthquake events [40], cracks follow openings as they are the weakest part in a wall. Under seismic action, the stresses on the openings’ corners assume the highest values and the corners are usually attractors for crack initiation.

In literature, there are few research works, based on the macro-block approach [32,40], analysing cases in which there are openings in the masonry wall. Casapulla et al. [32] investigated the role of the openings in the rocking mechanism of masonry wall ends. Their model constrains the lower crack fixed a priori as strictly related to the width of the pier and height of the lower spandrel. For instance, the analytical model developed in Reference [32], considers only the inclination of the crack above the openings as a variable of the failure mechanism. Realistic collapse mechanisms include that the crack crossing the lower spandrel (below the opening) may change its inclination, as well as the possibility that the crack, concerning to the considered opening, involves the upper spandrel varying its starting point. Despite other analytical models [31,32], the proposed formulation is able to cover these cases by managing a larger set of variables and improving the research panorama in which the cracks could propagate. For that purpose, the proposed visual program is implemented by a C# component, which is adopted to take into account the constraints that the cracks inclinations must respect.

Figure 7a shows the optimisation process, which is based on 20 iterations of the Nelder- Method-Optimisation component [37]. The solver converges to a value of 0.21, after less than 10 iterations, in 3 s with an Intel^® Core™ i7-6700HQ processor. As shown in Figure 7b, the geometry of the obtained failure mechanism is in good agreement with those arising from the micro-block approach (see Figure 5b) [39]. The load multiplier obtained at the end of the optimisation procedure is 20% more conservative than that computed in Reference [39], which is reasonable given the level of simplification and the large size of the masonry units when compared to the size of the piers and spandrels.

To clarify the advancements brought by the proposed approach, a comparison between the limit analysis model (Model*) developed in Reference [20] and the one proposed here (Model) is reported in Figure 8. The main difference between Model* and Model are the constraints about the position of the lower and upper cracks, which, in the present case, are removed by allowing the optimization problem to cover a larger set of kinematically compatible solutions. The results reveal how the proposed Model can obtain a lower value of $\lambda$, whereas Model* computes a load multiplier which is not conservative if compared with that arising from the micro-block approach.

4.3. Church Loaded Out-of-Plane

The third case study is a masonry church characterised by a plan having the form of the Latin Cross. The 3D model of the church, and the mechanical parameters, are taken in agreement with those reported in Reference [41] (Figure 9a). In particular, the *.dwg 3D model has been downloaded from the following link https://data.mendeley.com/datasets/ycxvmj77x5/1.

In the present work, the failure mechanism generated by the seismic action activated along the direction perpendicular to the façade (i.e., out of plane) is considered. Figure 9b shows the reference failure mechanism obtained by using a non-associative friction solution and micro-block modelling [41].

Figure 10a shows the results of the optimisation process. Considering a numerical approximation of two digits, the solver converges to a value of 0.21, after less than 10 iterations, in 5 s with an Intel^® Core™ i7–6700HQ processor. It is worth noting how the geometry of the failure mechanism is in agreement with that obtained by adopting a micro-block model and a non-associative friction flow rule (Figure 10b). The proposed macro-block model also shows excellent prediction in terms of the load multiplier, which is 5% more conservative than the one obtained using the micro-block approach.

5. The Step towards Full Dynamic Seismic Assessment Tools

In this section, the results of the proposed upper-bound limit analysis tool for the masonry wall with openings and the church presented in the previous section are evaluated towards full dynamic seismic assessment procedures through comparisons with predictions of the equivalent static displacement-based procedure of the Italian building code [15,16], while the influence of ground motion characteristics (pulse vs. non-pulse-type) on dynamic response is also investigated.

5.1. Code Based Assessment According to the Italian Building Code

Once the load multipliers and mechanisms of the case studies have been computed, a force-based and equivalent static displacement-based assessment was conducted, using the procedure defined in the Italian building code [16,17]. The code-based assessments were performed for the so-called life-safety limit state (SLV) in the code, which typically corresponds to a return period of 475 years.

Regarding the force-based assessment, the code defines the acceleration capacity as the spectral acceleration of mechanism activation $a_0^{*}$, which is computed by using the following equation:

$$a_0^{*}=\frac{\tilde{\lambda }g}{e^{*}CF}$$

(9)

where $\tilde{\lambda }$ is the load multiplier that activates the rocking mechanism (Section 2.2), $g$ is the acceleration due to gravity, CF is the confidence factor (which is here assumed equal to 1.00) and $e^{*}$ is the participating mass fraction related to the effective participating mass M*. The latter two parameters are computed using the following equations:

$$M^{*}=\frac{{\left({\sum }_{i=1}^n{W}_i{\delta }_{x,i}\right)}^2}{g{\sum }_{i=1}^n{W}_i{\delta }_{x,i}^2}$$

(10)

$$e^{*}=\frac{g{M}^{*}}{{\sum }_{i=1}^{n+m}{W}_i}$$

(11)

where $W_i$ and ${\delta }_{x,i}$ are defined in Section 2.1.

The corresponding acceleration demand is equal to the PGA a_g multiplied by the soil factor S, in order to get the maximum acceleration expected on the site, divided by the behavior factor $q$, which is typically equal to 2 for masonry buildings. Here, the objective is to determine the PGA $a_{g,FB}$ for which this safety criterion is violated, which is found by setting the demand and capacity equal to each other, resulting in:

$$a_{g,FB}=\frac{q{a}_0^{*}}{S}=\frac{q\tilde{\lambda }g}{S{e}^{*}CF}$$

(12)

In the case of the displacement-based procedure, the ultimate displacement capacity $d_{SLV}{}^{*}$ is defined as 40% of the displacement $d_0{}^{*}$, which would cause the structure to overturn. Specifically, $d_0{}^{*}$ is computed as the spectral displacement, which corresponds to the null value of the spectral acceleration. Thus:

$$d_{SLV}{}^{*}=0.4d_0{}^{*}$$

(13)

The corresponding displacement demand ${\Delta }_d$ is then computed by evaluating the elastic spectral response acceleration $S_e$ at the equivalent secant period of the mechanism T_SLV with the latter calculated at $d_{SLV}{}^{*}$. The acceleration $a_{SLV}{}^{*}$, corresponding to this displacement, is obtained from the capacity curve generated for each mechanism (see Figure 11). The control points adopted to generate these curves were selected to coincide with the centroid of the macro-blocks. As Figure 11 illustrates, the curves are linear for both mechanisms and are defined by $a_0^{*}$ and $d_0^{*}$, where $d_0^{*}$ is the overturning displacement of the structure, as defined above. The displacement demand ${\Delta }_d$ can thus be calculated as:

$${\Delta }_d=\frac{T_s^2}{4{\pi }^2}g{S}_e\left(T_s\right)$$

(14)

where:

$$T_s=1.68\pi \sqrt{\frac{d_{SLV}{}^{*}}{a_{SLV}{}^{*}}}$$

(15)

In the present study, the displacement demand is computed by means of the 5% damped elastic response spectrum, as defined in Eurocode8 [42]. Note that Eurocode8 recommends two different design spectrums–Type 1 for high and moderate seismicity regions and Type 2 for low seismicity regions and near field earthquakes. The Type 1 spectrum was used for these analyses as the church is located in Italy, which is a zone more prone to moderate magnitude earthquakes at closer distances. The most common, intermediate, ground type B, was assumed (very dense sand or gravel or very stiff clay, with a soil factor S that increases the PGA by 20% with respect to the rock base). Moreover, as 5% damping is assumed, the damping correction factor η is set to 1. The parameters defining the spectrum are summarized in

Table 1. Parameters used to define the Eurocode8 response spectrum.

Parameter	S	TB (s)	TC (s)	TD (s)	η
Value	1.2	0.15	0.5	2.0	1.0

, while the equations describing the spectrum are given by:

$$S_e\left(T\right)=a_g\ast S\left[1+\frac{T}{T_B}\left(\eta \ast 2.5-1\right)\right]\hspace{1em}\hspace{1em}\hspace{1em}0\le T \le T_B$$

(16)

$$S_e\left(T\right)=a_g\ast S\ast \eta \ast 2.5\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{T}_B\le T \le T_C$$

(17)

$$S_e\left(T\right)=a_g\ast S\ast \eta \ast 2.5\ast \frac{T_C}{T} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{T}_C\le T \le T_D$$

(18)

$$S_e\left(T\right)=a_g\ast S\ast \eta \ast 2.5\ast \frac{T_C{T}_C}{T^2}\hspace{1em}\hspace{1em}\hspace{1em}{T}_D\le T$$

(19)

The PGA a_g,DB, for which the displacement-based safety criteria is violated, is then found by setting the displacement demand ${\Delta }_d$ equal to the displacement capacity $d_{SLV}{}^{*}$ and solving for $a_g$.

Table 2. Comparison of peak ground accelerations (PGAs) for which the Italian code safety check is violated using both the force-based and displacement-based procedures.

Mechanism	ag,FB (g)	ag,DB (g)
Masonry wall with openings	0.33	0.35
Church	0.27	0.48

reports the PGAs for which the safety check of the Italian building code is violated, using both linear kinematic analysis (force-based assessment) and non-linear kinematic analysis (displacement-based assessment). For the masonry wall with openings, the PGAs obtained from both the force-based and displacement-based methods compare reasonably well. However, for the church mechanism, the displacement-based procedure predicts a PGA which is almost double that obtained from the force-based approach. This is most likely due to the larger scale of the church, which makes it less vulnerable to collapse than the smaller scale masonry wall. This scale effect is partly accounted for by the displacement-based procedure, but not by the force-based approach, which instead relies only on the slenderness of the mechanism.

In order to further evaluate the building scale effects, a parametric analysis in terms of scale factor (SF) was performed for the church. To this end, the following values of the scale factors are considered: SF = 1; SF = 1/2; and SF = 1/3. From a comparison of the capacity curves, as illustrated by Figure 12, it can be seen that the displacement capacity is significantly influenced by the scale of the macro-block, whereas the value of $a_0^{*}$ remains constant for all the macro-block configurations. It is also noted that the displacement capacity is rather large when the scale increases, which raises some doubts over the capacity of the structure to withstand such a displacement level before disintegration.

Table 3. Comparison in terms of PGAs—SF for witch the Italian building code is violated using both the force-based and displacement-based procedures.

SF	ag,FB (g)	ag,DB (g)
1	0.27	0.48
1/2		0.24
1/3		0.17

reports the PGAs for which the safety check of the Italian building code is violated using the non-linear kinematic analysis for all the three idealized configurations. This result further illustrates how the seismic safety evaluation, performed following the displacement-based assessment, is influenced by the building scale, whereas the force-based method is independent of the scale and instead depends entirely on the slenderness of the macro-block.

5.2. Comparison of Rocking Analysis Results with Code-Based Assessment Procedures

In this section, dynamic analyses are conducted on the collapse mechanisms to provide a preliminary demonstration of the ability of the proposed tool to address safety and to compare the same with code-based procedures. To that end, full time-history analyses are performed for a set of both artificial and recorded accelerograms. In the case of the former, the program SIMQKE_GR [43], which provides a graphical interface to the program SIMQKE-1 [44], was used to generate a set of ten code-compatible artificial accelerograms (Figure 13), using as input the 5% damped Eurocode8 elastic response spectrum, as defined in the previous sub-section. Each record has a total duration of 20 s, while the average duration of the intense phase of the ground motion is 14 s. In the case of the recorded accelerograms, a set of 100 ground motions were obtained from the PEER NGA-West2 Ground Motion Database [45] (Figure 14), also scaled to the 5% damped response spectrum as defined above. The set of analysed ground motions comprise 50 different pulse-type and 50 different non-pulse type records, with an average close to the code spectrum and average durations of the intense phase of 13 s and 20 s, respectively. Additionally, as both mechanisms are assumed to undergo one-sided rocking, each accelerogram is run with reversed polarity as well, resulting in a total of 20 (artificial accelerogram) + 200 (recorded accelerogram) analyses being conducted for each mechanism.

In order to evaluate the predictions of the tool, an equivalent static displacement-based analysis is also conducted for each ground motion record, using, in this case, the actual response spectrum for each accelerogram. Additionally, to better compare the demand-capacity ratios obtained from the code-based results to those from the non-linear dynamic (rocking) analyses, the maximum rotations ${\theta }_{max}$ predicted by the tool were normalised by the equivalent code-mandated maximum allowable rotation ${\theta }_{ov,c}$, which is simply equal to 40% of ${\theta }_{ov}$ (where ${\theta }_{ov}$ is defined in Equation (8)).

5.2.1. Masonry Wall Loaded in Its Plane with Openings

In the case of the in-plane loaded masonry wall with openings, the kinematic constants of the mechanism are first calculated using the Python script in Grasshopper [30]. The kinematic constants are then used to determine the equivalent rocking parameters p_eq, $\tilde{\lambda }$ and η_COR defining the rocking equation of motion, which are listed in

Table 4. Masonry wall with openings. Equivalent rocking parameters calculated by the GhPython script.

peq (s−1)	$\tilde{\mathit{\lambda }}$ (rad)	ηCOR
1.55	0.20	−0.61

. Note that the coefficient of restitution η_COR is negative due to one-sided rocking being assumed for this mechanism.

In the first set of analyses, the rocking equation of motion was solved by scaling the suite of ground motions to a_g = 0.35 g for both the artificial and recorded accelerograms, as this was the PGA predicted to cause collapse by the displacement-based code approach. In the case of the artificial accelerograms, a comparison between the code-based predictions and those obtained from the rocking analyses revealed the code to consistently be the more conservative of the two (Figure 15, positive polarity indicated in black, negative in white). Specifically, using the actual response spectra for each accelerogram, the code predicted collapse in 70% of the cases (points lying above the solid horizontal line), while the maximum response predicted by the rocking model never exceeded 7.5% of ${\theta }_{ov,c}$ or 3% of ${\theta }_{ov}$.

Similar trends were also observed in the case of the recorded accelerograms. Specifically, in the case of both pulse and non-pulse type motions, the code predicted larger demands than the rocking model in 100% of cases (i.e., all points lying above the dashed line in Figure 16). In the case of pulse-type motions in particular, the code predicted collapse in 54% of the cases, while the maximum response predicted by the rocking model never exceeded 66% of ${\theta }_{ov,c}$ or 26% of ${\theta }_{ov}$. However, in the case of the non-pulse-type motions, the code only predicted collapse in 26% of the cases, while the maximum response predicted by the rocking model never exceeded 38% of ${\theta }_{ov,c}$ or 15% of ${\theta }_{ov}$. Nevertheless, in both cases, it can be seen that the code is much more conservative than the rocking model, and that the PGA of 0.35 g is insufficient to cause rocking collapse.

Thus, in a second set of analyses, the scaling of the ground motion was progressively increased until complete rocking collapse occurred, with collapse in this case being defined as θ = π/2. As the recorded ground motions are less conservative than the artificial ground motions, only the former are considered. From these analyses, it was found that for the masonry wall with openings, a minimum PGA of 0.50 g is required for at least one record to cause collapse, with the PGA having to be increased to 1.30 g for 50% of the recorded accelerograms (average of pulse and non-pulse type records) to cause complete overturning of the structure (Figure 17), which is 3.7 times higher than the code-based prediction of 0.35 g. These results clearly indicate the need for sounder code-base safety approaches such as the one proposed here. Additionally, the pulse-type ground motions were found to have more of a destructive effect on the rocking response (which was also illustrated by Figure 16), thus highlighting the sensitivity of the rocking model (as well as the code-based model to an extent) to ground motion characteristics.

The solution time for each of these rocking analyses was also recorded, and it was found that, when using an Intel^® Core™ i7-6700HQ processor, solutions are obtained in an average of 9 s.

5.2.2. Church Loaded Out-of-Plane

In the case of the church, the collapse mechanism found through the upper bound limit analysis was that of out-of-plane (OOP) failure of the façade with the partial overturning of the side walls as well. The equivalent rocking parameters computed for this mechanism by the GhPython script are listed in

Table 5. Church. Equivalent rocking parameters calculated by the GhPython script.

peq (s−1)	$\tilde{\mathit{\lambda }}$ (rad)	ηCOR
1.19	0.16	−0.46

. Note that the coefficient of restitution η_COR is again negative due to one-sided rocking being assumed for this mechanism as well.

In the first set of analyses, the rocking equation of motion was solved by scaling the suite of ground motions to a_g = 0.48 g for both the artificial and recorded accelerograms, which was the PGA predicted to cause collapse by the code’s displacement-based approach.

In the case of the artificial accelerograms, a comparison between the code-based predictions and those obtained from the rocking analyses again revealed the code to consistently be the more conservative of the two (Figure 18, positive polarity indicated in black, negative in white). Specifically, using the actual response spectra for each accelerogram, the code predicted collapse in 60% of the cases (points lying above the solid horizontal line), while the maximum response predicted by the rocking model never exceeded 21% of ${\theta }_{ov,c}$ or 8.4% of ${\theta }_{ov}$.

In the case of the recorded accelerograms, the rocking model actually predicted larger demands than the code in 4% of the cases for pulse-type ground motions (Figure 19), which could be due, in part, to the higher level of scaling employed as well as the increased slenderness of the church mechanism, as compared to the wall with openings. More specifically, from these analyses, it was found that the code predicted collapse in 64% of the cases, while the rocking model predicted collapse (according to the code and not according to equilibrium) in 2% of the cases (i.e., points lying to the left of the vertical line in Figure 19). However, the maximum rocking response never exceeded 51% of ${\theta }_{ov}$. In the case of non-pulse type motions, the code predicted larger demands than the rocking model in 100% of cases, but only predicted collapse in 44% of cases, while the maximum response predicted by the rocking model never exceeded 58% of ${\theta }_{ov,c}$ or 23% of ${\theta }_{ov}$. Nevertheless, for both types of ground motions (artificial and recorded, pulse and non-pulse type), it can be seen that the code is much more conservative than the rocking model, and that the PGA of 0.48 g is insufficient to cause rocking collapse.

Again, in a second set of analyses using recorded ground motion, the scaling of the ground motion was again progressively increased until complete rocking collapse occurred (i.e., θ = π/2). From these analyses, it was found that for the church, a minimum PGA of 0.63 g is required for at least one record to cause collapse, with the PGA having to be increased to 1.40 g for 50% of the recorded accelerograms (average of pulse and non-pulse type records) to cause the complete overturning of the structure (Figure 20), which is 2.9 times higher than the corresponding code-based prediction of 0.48 g. Again, the pulse-type ground motions were found to have more of a destructive effect on the rocking response (which was also illustrated by Figure 19), thus highlighting the sensitivity of the model (and to an extent the code-based approach) to ground motion characteristics. The results once more indicate the need for sounder code-base safety approaches such as the one proposed here.

The solution time for each rocking analysis was again recorded, and it was found that, when using an Intel^® Core™ i7-6700HQ processor, solutions are obtained in an average of 10 s at most.

6. Conclusions

This paper presents a new multi-level procedure that combines upper bound limit analysis with non-linear rocking dynamics to model seismic collapse of masonry structures in a computationally-efficient manner. The procedure has been entirely implemented within the visual programming environment offered by Rhinoceros3D + Grasshopper [30], which would enable it to be used by both academics and practicing engineers in a fast, simple and intuitive way.

The collapse mechanisms predicted by the tool are first validated using three case studies from the literature, comprising a masonry shear wall without openings, a masonry wall with openings and a masonry church. Good agreement was observed in terms of both the geometry of the failure mechanism and load multiplier. Subsequently, the collapse mechanisms of the wall with openings as well as the church were discussed with respect to seismic safety using code-based equivalent static approaches comprising both force and displacement-based procedures, as well as rocking dynamics.

From this study it was found that:

• The upper bound limit analysis coupled with the NMM tool is able to find the minimum load multiplier and the associate collapse mechanism in less than 5 s.
• The scale effect is not considered in the force-based approach, which can lead to considerably different results to those obtained from the displacement-based procedure.
• Both the displacement-based method and rocking dynamics were able to capture the sensitivity of the response to ground motion characteristics - specifically the more destructive influence of the pulse-type ground motions. However, the code-based procedure tends to predict larger displacement demands than the rocking model, which in practice may result in expensive and at times unnecessary interventions, while also not being able to reproduce motion evolution over time.
• The solution to the rocking equation of motion is obtained, on average, within a maximum of 10 s, thus demonstrating the ability of the proposed tool/procedure to conduct a large number of analyses in a time frame compatible with engineering practice.
• The definition of safety requires a full stochastic approach, which is the objective of the work currently being developed by the authors.