A discrete‐element approach accounting for P‐Delta effects

This paper presents a discrete macro‐element accounting for P‐Delta effects to describe the rocking response of masonry walls subjected to out‐of‐plane (OOP) loadings. Both constitutive and geometric nonlinearities are considered within the discrete macro‐element method (DMEM), which is a modeling approach characterized by a very low computational cost compared to traditional distinct element method (DEM) and detailed finite element (FEM) strategies. OOP failure mechanisms are the main cause of severe damage for unreinforced masonry (URM) buildings, subjected to seismic actions. These mechanisms are generally activated by low seismic excitation and displacements. However, after their activation, they can evolve towards large displacements related to rigid‐block‐like kinematics that strongly affects the nonlinear mechanical response. Therefore, geometric nonlinearities, often ignored, should be included in the analyses. A new simplified, still accurate P‐Delta formulation is presented, according to which the global equilibrium is imposed by referring to the undeformed system configuration, avoiding assembling the geometric stiffness matrix. Namely, the system load vector is updated at each step of the analysis, accounting for the additional moments generated by the in‐plane compression forces acting on the macro‐elements in the deformed configuration. The proposed model is validated against closed‐form analytical solutions of rigid‐block benchmarks in large displacements and the results of experimental tests already available in the literature. In addition, extensive parametric analyses are performed to investigate the role of different mechanical and geometric parameters characterizing the ultimate non‐linear response of masonry walls subjected to horizontal forces. The results show how the proposed model, including P‐Delta effects, accurately predicts the non‐linear rocking response of masonry walls until the attainment of the unstable configuration.


INTRODUCTION
Masonry structures constitute a significant percentage of the existing buildings worldwide, including historical and monumental constructions, representing an architectural and cultural heritage. Historical masonry buildings are generally characterized by masonry walls weakly connected to each other and to the other structural elements leading to complex structural behavior under earthquakes. The seismic response of masonry monumental structures is generally governed by the coupled in-plane and out-of-plane (OOP) responses of masonry walls. 1,2 In-situ post-earthquake observations clearly showed that for unreinforced masonry (URM) constructions, the OOP behavior of masonry walls is the most likely cause of structural damage or collapse. 3,4 In addition, once activated, OOP failure mechanisms can evolve towards large displacements, which are generally assumed to be represented by rigid-block-like kinematics, 5 strongly influencing the structural response. Hence, geometric nonlinearities assume great importance when performing a non-linear seismic analysis on a URM masonry wall undergoing OOP rocking mechanisms. 6 This entails a significant increase in complexity in the study of the non-linear behavior of URM structures already characterized by a highly non-linear constitutive behavior even at a slight intensity of the seismic load. On the other hand, the accurate numerical simulations of the non-linear incremental static or dynamic response of masonry walls subjected to seismic excitation, considering both constitutive and geometric nonlinearities, represent a very complex computational issue that has been the subject of extensive research over the last decades. Accurate analyses employing detailed finite element (FEM) or distinct element (DEM) models, accounting for both geometric and constitutive nonlinearities in large displacements, requiring the adoption of advanced numerical strategies to update the system configuration and the geometric stiffness matrix during the analyses. Consequently, these approaches present the disadvantage of being more complex, time-consuming, and computationally expensive. However, geometric nonlinearities can be taken into account in a simplified way by considering the so-called second-order (P-Delta) effects. Despite this, the number of P-Delta formulations available in the literature for assessing masonry structures is limited and generally restricted to the use of FEM simulations or oversimplified single-degree-of-freedom (SDOF) equivalent descriptions.
Numerous non-linear numerical strategies, characterized by different levels of accuracy and efficiency, have already been proposed to assess the non-linear response of URM walls when subjected to OOP loading. Among these, effective tools are based on the macro-block limit analysis approach, which describes the wall as a kinematism of rigid blocks and are used to evaluate the lateral load amplitude activating the OOP mechanism. 7-10 These methods are employed for the seismic assessment within the general framework of the force-based approach (FBA) as implemented in international structural codes (EC8 11 ; Italian Code NTC18 12 ). Within this framework, advanced limit-analysis formulations have been proposed to account for the stabilizing effects of retaining walls, [13][14][15][16] and the interaction with horizontal diaphragms. 17 In addition, more detailed mesoscale limit-analyses formulations have been proposed to investigate the torsion-shearbending moment strength domains in OOP mechanisms. [18][19][20][21] However, force-based limit-analysis approaches cannot provide information concerning the evolution of the mechanism after its activation and the effective displacement capacity of the system before reaching the overturning. 22,23 To overcome this limit, displacement-based limit-analysis approaches have gained popularity for the seismic evaluation of existing masonry structures. 5,[24][25][26][27] These methods aim at assessing, in addition to the ultimate lateral force, the displacement capacity of the system. The latter is generally assumed to coincide with the configuration when the wall's self-weight becomes unstabilizing. Then, the displacement capacity of the system is compared to the corresponding seismic demand evaluated by dynamic time-history analyses or spectrum analyses conducted on an equivalent SDOF system. 28 These strategies are straightforward and take into account P-Delta effects, and for these reasons, they are widely used in practice and have been implemented in some structural codes, such as the Italian code (NTC18). However, they neglect the stabilizing effects of inertia forces and require the preventive assumption of the failure mechanism. Another effective method largely used to assess OOP rocking walls is the rigid-block model, which describes the entire wall by a single rigid block overturning around its base vertexes and subjected to ground acceleration time histories. This model response is generally obtained by a direct integration of the equation of the motion in large displacements. This approximate approach, since the wall is assumed rigid, allows for a simplified dynamic analyses with a low computational cost. 9,[29][30][31] However, the rigid'body model does not account for masonry deformability, threedimensional wall's boundary conditions, and requires specific hypotheses on the energy dissipation model which may be unsuitable for describing masonry walls. Finally, the rigid-block hypothesis introduces further approximations when adopted for describing complex 3D failure mechanisms, as those characterizing the ultimate response of large URM walls interacting with retaining walls or other structural elements.
On the other hand, FEM and discrete-element approaches, including the distinct element method (DEM), can explicitly account for 3D OOP boundary conditions allowing for geometric nonlinearities, that is, second-order effects of largedisplacement kinematics. Among them, numerical methods characterized by different orders of complexity have been proposed and applied. Mesoscale FEM analyses 32 explicitly accounting for the actual masonry bond through the use of zero-thickness interfaces, or continuum FEM analyses through homogeneous equivalent media at the macro-scale 21 have been performed adopting co-rotational approaches, 33 which represent an accurate and relatively efficient tool under the hypothesis of small deformations and large displacements. DEM approaches are well suited for masonry structures with dry-and mortared joints to perform non-linear static and dynamic analyses accounting for large displacements and the separation of structural parts. In these approaches, the position and orientation of elements and contacts are updated at each step of the analysis, 35 making them particularly suitable to describe local failures, and characterizing irregular masonries. 35,36,66 However, FEM or DEM strategies need complex 3D constitutive laws and cumbersome calibration procedures 21,37 and require high computational costs, making them generally unsuitable for engineering applications. For the above reasons, reliable numerical models effectively used in research are rarely compatible with practical seismic assessments 38 and are difficult to be implemented consistently to technical code prescriptions and within commercial engineering-oriented software packages. 39 For these reasons simplified engineering-oriented structural macro-elements are often adopted, including the equivalent frame models 40,41 or the discrete macro-element method (DMEM). These strategies significantly reduce the computational cost of seismic analysis, making them suitable to be employed by practitioners. However, these approaches often neglect the OOP response of masonry walls, and their use is restricted to assessing the global response of masonry buildings. Despite this, recently, an enhanced frame-like macro-element incorporating the OOP response also accounting for second-order effects has been proposed in the literature. 42,43 The present study aims to cover this research gap by upgrading the discrete macro-element method (DMEM) introduced by Caliò et al. 44 including geometric nonlinearities related to P-Delta effects. The DMEM was originally proposed for the in-plane response of masonry walls 44 and subsequently extended to 3D kinematics to account for the coupled in-plane and OOP behavior of masonry walls when subjected to earthquake loadings. 45,46 This strategy discretises the masonry walls by shear-deformable macro-elements interacting with the adjacent ones by means of zero-thickness discrete or continuous 47 interfaces. The model presents many advantages, compared to previously proposed models for masonry structures, related to the geometrical consistency, the low computational cost, the possibility to couple discrete and finite elements, the straightforward model calibration, the possibility to be used at macro and meso-scale, 48 and other advantages duly expressed in the referenced papers. The model is also particularly suitable for infill frames and confined masonry structures 34,49 due to the accurate modeling of the beam-masonry interaction.
The novelty proposed in this paper consists of accounting for the P-Delta effects, in a computationally-efficient discrete macro-element (DMEM) strategy, by imposing the current equilibrium at each step of the analysis by modifying the global load vector taking into account the eccentricity of the external and internal loads with respect to the initial reference configuration. The model has been implemented within the HiStrA (Historical Structure Analysis) software 50 and has been validated by simulating parapet and vertical spanning brick-wall samples subjected to one-way OOP bending moments. The numerical predictions have also been compared with reference exact analytical rigid-body solutions. The influence of some model parameters, such as the cross-section discretisation employed in the interface elements and the wall deformability, have also been investigated through some parametric analyses. Finally, the model has been employed to simulate some experimental tests conducted on walls under OOP bending, 21,51-53 which has also been numerically investigated by other authors. 21,53 The comparisons between the numerical and experimental results demonstrate the accuracy of the proposed modeling strategy and its potential to be employed for practical assessments of URM buildings whose walls are subjected to OOP rocking mechanisms.

THE DISCRETE MACRO-ELEMENT METHOD
The DMEM has been proposed with the aim to introduce a low-computational-cost numerical strategy alternative to nonlinear FEM and DEM strategies. In its first proposal, the macro-model approach 44,54 was restricted to 2D kinematics to propose a simplified approach alternative to those based on equivalent frame models (EFM). Subsequently, the model was extended to 3D kinematics, introducing a discrete macro-element approach for an efficient simulation of in-plane and OOP response of masonry buildings. 45 The model has been implemented in the software HiStrA 50 and extensively validated against numerical and experimental results through non-linear static analyses 45,55 and time-history analyses simulating the experimental results of shaking table tests performed on URM prototypes. 46 Different applications on case studies are also available in literature. 56,57 According to this strategy, a masonry wall is divided into macro-portions, and each macro-portion is represented by shear-deformable spatial elements connected to each other through non-linear zero-thickness interfaces. Each macroelement inherits the geometry directly from the macro-portion of the structure that represents, and it is also identified by four vertexes (v 1 ,..,v 4 ), characterizing the middle plane of the element (π) ( Figure 1A) accounting for a generalized shear deformability described by a single Lagrangian parameter. The typical geometry of the general macro-element is sketched in Figure 1A where the main geometrical parameters are also reported.

The kinematics
A local cartesian reference system, fixed with respect to the absolute reference system, is associated within each macroelement with the origin coincident with the barycenter G of the element, where is also concentrated the competing mass. The local axes (i, j, k) are identified by the unit vectors i, j belonging to the plane of the element, π, with i oriented as the v 1 -v 2 direction, j oriented as the v 1 -v 4 direction, and k identifying the direction orthogonal to π ( Figure 1A). The element kinematics is described by seven Lagrangian parameters, collected in the vector d, expressed in Equation (1). The first six parameters, (U, V,W, Φ, Θ, Ψ), are related to the rigid spatial motion, three translations and three rotations of the barycenter G referred to the local reference system. A further Lagrangian parameter (γ) is related to the generalized macro-element shear deformation and is associated with the variation of the angle between the panel edges connecting the vertex v 1 to vertex v 2 and the vertex v 1 to vertex v 4 , respectively ( Figure 1B). The zero-thickness plane interfaces π i (i = 1,..4), ( Figure 1A), whose orientation is associated with unit vectors n i (i = 1,..4), cut and identify the element along its perimeter and rule the connections with the adjacent elements. Different thicknesses can be considered at each of the panel's vertexes (t 1 ,..,t 4 ), leading to direct and easy modeling of a linear variable wall thickness. Local reference systems e n1 , e n2 , e n3 (n = 1,..4) are associated with the perimeter interfaces, as shown in Figure 1C, for the description of the interfaces' kinematics. More details on DMEM kinematics can be found in ref. 45,47 It is worth noticing that no additional Lagrangian parameters is needed to describe the kinematics of interfaces. Therefore, the total number of degrees of freedom of the model is directly given by 7 × N being N the number of F I G U R E 2 (A) Spatial view of the discrete macro-element, (B) macro-element's middle plane, and (C) intrinsic reference system. macro-elements; this aspect strongly contributes to containing the model computational burden.
Considering a second-order linearized kinematics for the adopted macro-element, a direct expression between the displacement vector of the n-th (n = 1,. . . ,4) vertex of the irregular quadrilateral (u vn ), in the relative reference system, and the seven Lagrangian parameters Equation (1), can be established as follows: Being I 3x3 the 3 × 3 identity matrix and Ω vn a matrix related to the rigid rotation of the element, defined as follows: with: x vn , y vn , the coordinated of the vertex v n in the local reference system; (n = 1,..,4) a vector accounting for the contribution of the generalized shear deformation that, for each vertex, can be written as: Being α n (n = 1,. . . ,4) the corner angles, which are measured at the middle plane of the element ( Figure 1A). It is worth noticing that the generalized shear deformation is defined in the middle plane of the macro-element and is assumed independent of the z coordinate, referring to the local axis k. An intrinsic reference system (ξ, η) is introduced aiming at evaluating the displacement of any point P(x,y,z) belonging to the volume of the element (Figure 2A). The point P is projected on the plane of the element (π), identifying the point P π (x,y) in Figure 2B, which correspond to P π (ξ,η) in the intrinsic space ( Figure 2C) through the bi-linear polynomial transformations reported in Equation (5). According to an isoparametric formulation, the displacement of the point P π (ξ,η) in the local reference system can be obtained by interpolating the displacements of the vertexes, using the interpolating function in Equation (5), and the results is reported in Equation (6). Where: Equation (6) in compact form can be written as: Being ( , , ) a 3 × 7 matrix collecting seven shape functions governing the element kinematics Each shape function is a vector with three components that can be easily derived by Equation (6). It is worth noticing that the shape functions depend on the shear deformation parameter and are not affected by the local coordinate z. Consequently, each plan parallel to π presents the same shear deformation.

The mechanical behavior
The element deformability is lumped in the zero-thickness interfaces and partially distributed in the generalized shear kinematics. 55 The mechanical interface calibration follows a direct fiber discretization approach, as qualitatively shown in Figure 3. In particular: a distribution of normal (or transversal) non-linear links ( Figure 3A) concentrates the axial and flexural deformability of masonry both in-plane and OOP; a single in-plane longitudinal link ( Figure 3B) governs the shear-sliding mechanism between masonry joints; two out-of-plane shear links control the OOP sliding mechanism as well as the torsional response ( Figure 3C). A further diagonal link ( Figure 3D) rules the generalized shear diagonal deformation. Despite the simplicity of the adopted mechanical scheme, such an approach accounts for the influence of the axial load on the flexural behavior and the contemporary occurrence of the bending moments along both the in-plane and OOP directions without the need to define complex bi or three-dimensional plasticity domains. The links are calibrated by imposing simple equivalences between the discrete model and a corresponding equivalent continuous anisotropic material according to different criteria depending on their purposes. 47 More specifically, transversal links can be calibrated according to different constitutive laws with fracture energy both in tension and compression; sliding behavior is calibrated according to Mohr-Coulomb criterion, accounting for the normal stress obtained by the transversal force of the interface; the diagonal behavior is limited and generally calibrated according to the Mohr-Coulomb or Turnšek and Čačovič 57 criteria or more specific constitutive laws for masonry structures that can also include the role of damage. 58,59

The external loading contribution
The vector ( , ) , expressed in Equation (7), provides the displacement components in the local reference system for each point of the macro-element, given the generalized shear deformation independent on the local coordinate z. Considering a force distribution f(x,y), belonging to the middle plane of macro-element (π), the corresponding vector of equivalent forces F 0 , associated with the macro-element degrees of freedom, can be evaluated through the external virtual work: Being A the load application area and̃( , ) a virtual displacements field that, in view of Equation (7), can be written in terms of intrinsic coordinates and the virtual displacements (̃T) corresponding to the adopted degrees of freedom, expressed in Equation (1).
In view of Equations (7) and (9), the external virtual work can therefore be expressed as: where J(ξ,η) is the Jacobian of the transformation. Equation (10) leads to the vector of equivalent forces F 0 of the macroelement for any load distribution applied on the middle plane of the element, as follow: If the force distribution is applied with a z eccentricity, the corresponding moment must be added to the vector F 0 , which represents the loading distribution applied in the middle plane of the element.

Evaluation of the P-Delta within the DMEM
Geometric nonlinearities can play a fundamental role in assessing the ultimate capacity of rocking masonry walls under earthquake actions. After the activation of the rocking mechanism, the wall response can evolve toward large displacements. In this latter case, due to the typical slender geometry of masonry walls, the stabilizing effects of vertical loads progressively reduce until the achievement of the critical condition in which vertical loads no longer play a stabilizing effect. A rigorous description of geometric nonlinearities within the context of DMEM would require the adoption of large displacement kinematics with the step-by-step update of the geometric stiffness matrix and current configuration of the system. As a consequence, it requires a more complex theoretical formulation and, above all, a significant increase in the associated computational cost. Aiming at maintaining the benefits of the DMEM, a simplified but robust and efficient P-Delta approach is proposed, characterized by a reasonable accuracy for engineering applications while maintaining a simple theoretical formulation and a very low computational effort. More specifically, the P-Delta effects are accounted for by updating the current positions of the external and alonginterface forces applied to the macro-element.
Let us consider a generic force distribution f(x,y) applied to a generic area A belonging to the middle plane of the macro-element. The corresponding vector moment ( − ) contribution with reference to the current configuration of the element, identified by the current degrees of freedom d, is given by Being ( , ) the displacement field of the loading application points. In view of Equation (7), Equation (12) can be written as: Considering j (j = 1,..,4) the resultants of the internal forces that the macro-element receives from the other elements through the interfaces, applied to the generic point c j ( , , ), the corresponding P-Delta moment results: Considering both the external and along-interface forces applied to the elements, the total P-Delta moment for each macro-element can be expressed as: Finally, the P-Delta (7 × 1) load vector, −Δ , can be written as:

MODEL VALIDATION
In this section, the novel P-Delta DMEM formulation is validated by performing non-linear static analyses on brickmasonry walls that have already been investigated in the literature. The walls are characterized by different geometrical layouts, materials, and boundary conditions. The numerical predictions obtained by the proposed DMEM model are compared with the analytical solutions obtained under the simplified hypothesis of rigid-block structure and with the results of some experimental tests on masonry-wall prototypes, available in literature. The analyses are performed using the HiStrA software, 50 where the model has been implemented. A modified iterative Newton-Raphson integration method with arch-length control is employed, according to which the system stiffness matrix is updated at the beginning of each step of the analysis. Finally, a convergence criterion based on the global equilibrium is employed. The convergence is achieved when the norm of the unbalance vector normalized by the norm of the total force vector (F 0 ), is less than a tolerance. The latter is fixed equal to 10 −3 in the analyses.

Numerical validation
In this section, the numerical predictions obtained by the proposed model are compared to the analytical results obtained by considering the rigid-block assumption. Under this hypothesis, the system does not show lateral displacements until the unstabilizing moment, due to the lateral forces, reaches the stabilizing value related to the wall self-weight and to the applied vertical loads. Once the incipient rocking condition is achieved, the rocking mechanism is activated and the monotonic load-displacement curve of the system follows a softening branch due to the effects of geometric nonlinearities until the critical displacement is reached. Nevertheless, the assumption of rigid-body behavior can be considered realistic for low values of vertical forces. Experimental and analytical studies 25,51,60 show that URM walls can deform significantly when subjected to high axial loads. The two main consequences of this are that (i) the pivot point has finite dimension so that the internal lever arm and the critical displacement (u cr ) are smaller than those corresponding to the rigid block, and (ii) the wall possess lateral deformability prior to incipient rocking ( Figure 4). As a consequence, the force-displacement capacity curve of the model with finite stiffness deviates from the ideal rigid behavior by assuming a curvilinear profile, as shown in Figure 4, where it can also be seen that the maximum force of the force-displacement curve is lower than the rigid threshold resistance F 0 . This latter behavior is encountered in the numerical simulations performed by the DMEM, being the interface links characterized by a finite stiffness. As a consequence, the rigid-body assumption has to be considered as a limit condition. As benchmark models, two isolated walls, whose mechanisms are characterized by the rotation around horizontal hinges, are considered. More specifically, a parapet wall (PW) of height h, thickness t and weight W, and a simplysupported wall (SSW), constituted by two rigid blocks of height h 1 and h 2 , weight W 1 and W 2 respectively and thickness t, spanning vertically between supports at ceiling/floor levels, are considered, as shown in Figure 5. A pre-compression force P is considered applied on the top section of the wall with an eccentricity e with respect to the middle section of the F I G U R E 4 Force-displacement relationships for rigid and deformable rocking walls.

F I G U R E 5 Rocking masonry walls: (A) parapet wall (PW) and (B) simply-supported wall (SSW).
TA B L E 1 Force capacity (F0) and critical displacement (ucr) associated with parapet wall and simply-supported wall mechanisms. wall. Only two one-way bending moment schemes are considered in order to avoid considering the influence of combined flexural-shear-torsion wall response. 61 The analytical expressions of the force-displacement (F-u) law is a rigid bi-linear curve and, therefore, can be described by evaluating the parameters F 0 and u cr , identifying the activation force and the ultimate displacement, respectively, reported in Table 1, for each of the two considered mechanisms. 62 In particular, u cr refers to the control points coincident with the geometric barycenter of the wall for the PW mechanism and the mid-height of the wall for the SSW mechanism. In the latter, for a sake of simplicity, it was assumed h 1 = h 2 and, as a consequence, W 1 = W 2 . Before undergoing non-linear rocking behavior, URM walls are characterized by a perfectly rigid response. Once the mechanism is triggered, the linear-softening branch becomes the reference curve until the displacement value equal to u cr is reached, corresponding to which the lateral force becomes zero. It can be observed that the displacement capacity is essentially a function of the wall thickness, whereas the strength capacity is significantly influenced by the wall boundary conditions. 51 In the simulations reported in the following, two different geometric layouts, characterized by different values of wall slenderness, are considered, as specified in Table 2, where l is the in-plane width of the wall.  The DMEM models are developed considering a single panel in the case of the PW and two panels for the SSW. In the latter case, it is considered h 1 = h 2 = h/2. The interface links of the macro-element are characterized by assuming a no-tension material, linear elastic in compression, with k n representing the stiffness of a unitary area, which is set equal to 5E + 08 N/m 3 , according to the assumption already adopted in the ref. 63 to approximate a quasi-rigid body. The results obtained by the DMEM for the two (PW and SSW) walls and the two different geometries are shown in Figure 6, compared to the corresponding analytical capacity curves of the rigid-block models. The results are normalized by the ultimate force (F 0 ) and critical displacement (u cr ) corresponding to the case of zero compression load P. In the analyses, the number rows of the interface transversal links, discretizing the wall thickness, was changed from 5 to 50 to investigate the influence of this discretization parameter on the wall response. In particular, the number of rows of links, being centered on the competent area, affects the approximation of the wall thickness and the location of the rotation point. Namely, using 50 rows, the effective model thickness (coincident with the distance between the two extreme links) is equal to 117.6 mm in the case of specimen 1 and 245 mm in the case of specimen 2, corresponding to an error of 2%. It is apparent that this error can be controlled by considering a different geometrical distribution of the non-linear links in the fiber discretization of the wall section. However, a classical fiber discretization has been adopted in the numerical simulations.

Mechanism
Observing the graphs in Figure 6, it can be noted that the proposed model well describes the overall response of the walls since the force-displacement curve tends to the theoretical rigid-block response as the number of links increases. Small differences between the numerical and the analytical responses are observed in the pre-peak and peak load stages due to the finite stiffness of the numerical model, as already qualitatively highlighted in Figure 4. These differences are more evident for specimen 1-PW, which is characterized by a higher value of slenderness. In all the investigated cases, the response is very close to the analytical limit case when more than twenty rows are employed in the discretization. In the following simulations, the interfaces are discretized using fifty rows of links to guarantee a good model accuracy, although for practical applications, ten rows could be sufficient to obtain satisfactory results. It is worth noticing that the number of links does not significantly affect the computational burden and computing time since the overall degrees of freedom are associated with the number of macro-elements only. Figure 7 shows force-displacement relationships varying F I G U R E 7 Force-displacement relationships varying k n . the interface stiffness k n from 5E + 7 N/m 3 to 2E + 10 N/m 3 . It can be observed that the numerical response is strongly affected by this parameter both in terms of peak force and ultimate displacement capacity. As expected, the greater the interface stiffness, the closer the numerical curve approaches the results predicted by the analytical response of the rigidblock, at the expense of a more significant computational burden. Higher values of normal stiffness are not consistent with the behavior of real masonry walls and reduce the model performances leading to an increase in the average number of iterations. The comparisons here reported aim to validate the proposed DMEM model and also to show the limits of the rigid-block behavior assumption. Slenderer wall (specimen 1) appears more sensitive to the variation of the interface stiffness, ranging the ultimate force from 40% to 60% of the corresponding theoretical peak force of the PW and SSW rigid block (Figure 7). In specimen 2, the corresponding peak-force reductions result in 60% and 80% (Figure 7). These values are consistent with the simplified SDOF approach proposed by, 25,51,64 showing a peak force close to 75% of the theoretical rigid-block peak force. A less influence of the interface stiffness on the critical displacement is observed. The largest difference on the ultimate displacement between the DMEM and the rigid-block models is approximately 20%, observed in the case of PW of specimen 1, while the interface stiffness does not significantly affect the ultimate displacement of the SSW specimens.
As a further investigation, the influence of an applied axial load is also evaluated by comparing numerical and analytical results in terms of force-displacement relationships (F-u) by varying the axial compression force applied at the top of the wall. It can be observed how the presence of a vertical pre-compression plays a fundamental role on the ultimate response of rocking walls as it leads to an increase of the lateral-load capacity but causes the decrease of the critical displacement of the wall ( Figure 8).
As already observed in previous results, the most significant differences correspond to the PW and SSW mechanisms of specimen 1. In particular, PW walls show a progressive inconsistency between the rigid-block curve and the DMEM prediction with the increase of the axial load, reaching a peak lateral force of 35% and 43% of the rigid-block body in case of specimen 1 and specimen 2, respectively, when the pre-compression is 0.15 MPa. These differences reduce in the corresponding specimens of the SSW walls, reaching 60% and 80% of the peak lateral load of the rigid body, respectively. A lower but still significant dependency of the compression load on the ultimate displacement is observed, with the maximum discrepancy between the numeric and analytical model of 20% observed in the case of specimen 1 of PW. These results are particularly significant regarding the assessment of historical constructions, when the PW mechanism is more likely to be activated and, as evidenced by the above results, the presence of a significant compression load can make the rigid-block predictions less accurate.

F I G U R E 8
Force-displacement relationships varying the wall pre-compression force.

F I G U R E 9
Force-displacement relationships considering (continue line) and neglecting (dashed line) P-effects for specimen 2.
Finally, in Figure 9, the role of P-Delta effects is quantified by comparing the predictions obtained by the previous and novel DMEM formulations, neglecting and considering the P-Delta effects, respectively. The comparisons highlight how neglecting the P-Delta effects not only does not allow the prediction of the ultimate displacement but also affects the accuracy in terms of peak load prediction. From Figure 9, it is also evident how the P-Delta effects become significant at a rather low displacement magnitude. Considering the cases investigated in this section, that threshold can be assumed to be 10% and 20% of the ultimate displacement for PW and SSW, respectively.

Experimental validation
In this section, the proposed modeling strategy is employed to simulate a series of experimental tests conducted on isolated walls subjected to a one-way OOP bending moment, carried out by ref. 51 and numerically investigated by different authors. 21,[51][52][53] The comparison among the numerical DMEM prediction, the experimental observations and the results of more refined modeling strategies enables evaluation of the accuracy of the proposed model in describing the OOP rocking behavior of URM walls.

Experimental tests
The experimental campaign has been performed at the Chapman Structural Testing laboratory of the University of Adelaide and the results are reported in ref. 51 and ref. 52. The tests have been carried out on simply-supported walls with and without pre-compression, representing load-bearing and non-load-bearing walls in URM buildings, which were subjected to monotonic static loads and dynamic excitations along the OOP direction ( Figure 10A). The specimens consisted in single-leaf brick-masonry walls, 110 mm thick, 1500 mm height, and 950 mm width. Standard extruded clay brick units with dimensions 230 × 110 × 76 mm 3 and a typical Australian mix of 1:1:6 (cement:lime:sand) mortar were used in the tests, and 10 mm mortar joints were adopted as well per standard construction practices in Australia. The average density of specimens was determined to be 1800 kg/m 3 .
The tests were conducted on uncracked and cracked specimens to investigate the influence of pre-existing structural damage on the OOP wall response. In the present study, the simulations are conducted considering only the uncracked specimens subjected to static loads. One specimen was subjected to an initial axial load designed to provide a uniform pre-compression σ v = 0.15 MPa to simulate the conditions of a real wall in a multistorey building. In all the static tests, the lateral load (F) was uniformly applied at the wall mid-height using a hand pump-driven hydraulic actuator. A rigid steel frame prevented the lateral OOP displacements at the base and top sections of the wall ( Figure 10A). In the specimen with zero axial load, both the vertical displacement and the rotation at the top section of the wall were allowed. In the specimen with initial axial load, the vertical displacement and rotation at the top section were partially restrained due to a spring, which was put in unilateral contact with the wall to transfer the initial vertical load and simulate the presence of a slab interacting with the wall. However, the elastic stiffness of that spring was not measured experimentally.

Numerical simulations and comparisons
The analyses have been performed by adopting, for the model, a refined mesh. Namely, a number of 18 macro-elements, corresponding to the number of brick layers of the specimen, and 19 interfaces, equally distributed along the wall height and discretized with seven links, have been considered. In the analyses, the top restraint of the model with precompression is modeled by considering a rigid element at the top section of the wall, connected with it by a no-tension interface (simulating the dry joint between the wall and the steel apparatus used in the test). The rigid element is then fully restrained against rotation and elastically restrained against vertical displacements ( Figure 10C), whose stiffness has been evaluated by fitting the experimental results, resulting in 1560 N/mm. The elastic and non-linear masonry material parameters required to calibrate the numerical model have been chosen according to ref. 53 and ref. 21 and are summarized in Table 3. In particular, the masonry tensile strength (σ t ) is assumed, according to ref. 65  and compression is used for the transversal links. Finally, an elasto-plastic constitutive law with a Mohr-Coulomb yield criterion is used for the shear sliding interface mechanisms. The results of the numerical analyses are reported in Figure  11 for the case of the specimen with zero axial load ( Figure 11A) and for the load-bearing wall ( Figure 11B), compared to the experimental results and to the predictions of more refined modelling strategies.
In both the analyzed specimens, the pre-peak branch is almost linear elastic with high stiffness, up to a maximum force that depends on the boundary and loading conditions. Once the maximum force is reached, which corresponds to the cracking of the mid-height interface element, the curve shows a sudden drop in resistance in the case of = 0. Then, the softening continues up to the residual strength of approximately 0.4 kN. Differently, the load-bearing wall recovers the lateral strength showing a much more ductile response due to the effect of the top restraint. The overall DMEM response agrees well with experimental results and FEM predictions for both the specimens investigated.
As a further investigation, parametric analyses have been performed on the wall with zero pre-tension to investigate the role of the mesh size and the number of interface transversal links. The simulations were conducted considering two different mesh discretisations. Namely, the fine mesh, already adopted in the simulation above, and a coarse mesh obtained by doubling the element size concerning the fine mesh. Finally, three different interface discretization, comprising 5-10-20 rows of transversal links along the wall thickness, have been considered. Figure 12A shows the DMEM responses obtained by varying the mesh size while keeping constant the number of rows of links equal to 20. Figure 12B shows the results obtained by changing the number of rows of links while adopting the fine mesh. In the graphs, the experimental results and numerical FEM predictions are reported for comparison. Observing Figure 12, it is possible to conclude that the investigated parameters do not significantly affect the accuracy of the proposed model. Moreover, the DMEM provided reliable predictions of the wall response, even adopting a coarse mesh and a low number of rows of links, with clear benefit in terms of the model's efficiency.
Furthermore, for the load-bearing wall, numerical results obtained with different values of the stiffness of the top spring (k spring ) are reported in Figure 13. In particular, two limit values are considered, corresponding to +50% and −50% of the reference value used to fit the experiments in Figure 11B. From Figure 13, it is apparent how this parameter affects the postpeak wall response and its ultimate displacement capacity. In the final analysis, the influence of masonry deformability on the rocking behavior of the specimen is investigated. Analyses are performed on the specimen with no pre-compression ( = 0) varying the masonry Young's modulus by doubling and halving the reference value reported in Table 3. The results are reported in Figure 14 concerning the refined mesh ( Figure 14A) and the coarse mesh ( Figure 14B) discretization. Comparing the graphs reported in Figure 14A and in Figure 14B, it can be concluded that the mesh discretization does not significantly affect the wall response, confirming the results reported in Figure 11A. Moreover, it is clear how the masonry deformability influences the pre-and peak-load stages of the response. More specifically, the lateral wall strength showed a variation of approximately 20% due to the Young's modulus variation.

Parametric analyses
In this section, the role of the main masonry nonlinear mechanical parameters on the wall response is investigated by considering the same specimens investigated in the previous section without applying any axial load. Figure 15 shows the first parametric investigation on the influence of tensile strength The analyses are performed considering an elasto-brittle constitutive law in tensile with three different values of tensile strength. It can be observed how the masonry strength significantly influences the peak load of the mechanisms. More specifically, the ultimate load corresponding to a tensile strength of 0.15 MPa is over double the ultimate load corresponding to the no-tension material model.
Finally, in Figure 16, three values of the tensile fracture energy are considered keeping constant the tensile strength equal to 0.15Mpa. Higher fracture energy values increase the peak force, ranging from 200% to 500% of the peak force corresponding to the rigid-block ultimate force. However, they cause a moderate increase in the displacement capacity that becomes negligible in the case of specimen 2. From Figure 15 and Figure 16, it can be concluded that tensile strenght and fracture energy affect the post-peak stage of the response up to 10% and 20% of the critical displacement in the case of the PW ans the SSW mechanism, respectively with more significal effects observed in the case of the specimen 1 and PW mechanism. This result confirms how the ultimate capacity of rocking walls is mainly governed by geometrical parameters rather than the material strength.

CONCLUSIONS AND FUTURE DEVELOPMENTS
The paper presents a further development of the DMEM by including an effective way of accounting for the geometric nonlinearities due to the applied in-plane compression forces, loads, which are particularly significant when rocking failure mechanisms of URM walls are activated, as typical of URM buildings and monumental structures. A simplified but robust P-Delta formulation is presented and implemented within a standard iterative Newton-Raphson method in the HiStrA software, which implements the DMEM approach. The proposed formulation does not require assembling and updating the geometric matrix of the system. Conversely, the geometric nonlinearities are addressed during the analyses by updating the global load vector according to the F I G U R E 1 6 Force-displacement relationship varying the tensile fracture energy (G t ).
macro-elements' current configurations. As a result, the proposed model does not require a significant increase of computational cost compared to DMEM strategy not accounting for the P-Delta effects.
The proposed model has been validated against analytical, numerical and experimental results available in the literature demonstrating its capacity to describe the non-linear response of rocking masonry walls subjected to different boundary and loading conditions. The results show good accuracy of the proposed model, which was able to describe the pre-peak and post-peak response as well as the ultimate lateral strength and the displacement capacity of masonry walls exhibiting common OOP rocking mechanisms. Namely, cantilever and vertical spanning mechanisms for different geometrical layouts and loading conditions confirmed the model's potential to be employed for real structural assessments of masonry structures whose non-linear response is characterized by the activation of rocking failure mechanisms.
Finally, parametric analyses have been conducted to investigate the role of (i) modeling parameters, like the mesh discretization and the interface link discretization; (ii) masonry deformability; (iii) axial load. Future developments will aim to extend the proposed model to the dynamic field to accurately simulate the effects of the geometric nonlinearities and the energy dissipation mechanisms on the rocking response and ultimate capacity of URM walls.

A C K N O W L E D G E M E N T S
The authors have nothing to report.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.