Distinct element modeling of the in‐plane response of a steel‐framed retrofit solution for URM structures

This paper presents the results of a numerical study aimed at assessing the effectiveness of an innovative steel modular retrofit system on improving the in‐plane lateral capacity of unreinforced masonry (URM) load‐bearing structures. The investigated retrofit solution consists of modular steel frames fastened to the external surface of masonry piers through chemical anchors. Distinct element models are built and calibrated against experimental data coming from a series of in‐plane quasi‐static tests performed on piers in both bare and retrofitted conditions. Two different masonry typologies are considered herein: one made up with solid clay bricks and lime mortar, assembled in a header bond pattern, and the other with typical Italian hollow clay units and cement‐lime mortar, arranged in a Flemish bond pattern. A simplified microscale approach, combined with a mesh‐refinement strategy, is adopted to decrease computational demand without losing numerical accuracy. All in‐plane failure mechanisms, such as mortar and unit tensile failures, as well as masonry crushing for high compressive stresses, are included in the developed numerical models. Moreover, a methodology to explicitly account for the contribution of the investigated retrofit in the distinct element method (DEM) framework is presented and discussed. The proposed modeling strategy yielded a good prediction in terms of lateral stiffness, lateral strength, hysteretic response, and crack pattern for unreinforced and strengthened specimens. A novel set of metrics is proposed to quantitatively assess the retrofit performance benefits.


INTRODUCTION
In recent years, the increasing interest in the combined improvement of seismic performance and energy efficiency of existing structures, has led to the development of innovative retrofit solutions able to fulfill structural, energy, and sustainability requirements. Among different structural typologies, unreinforced masonry (URM) construction is a widespread structural system that deserves special attention since a major part of the existing building stock consists of non-engineered structures with poor seismic detailing. [1][2][3][4] In this framework, Progetto Sisma s.r.l. has developed a retrofit solution called "Resisto 5.9," a modular steel system integrated with a homogeneous thermal insulation coating (e.g., polystyrene panels) also mitigating the differential thermal expansion between steel and masonry. Steel modules are connected to the external surface of the masonry through chemical anchors and among them by means of steel bolts, with the aim of improving the in-plane capacity of masonry piers. The retrofit system was also conceived in order to enhance the connection between orthogonal walls and among walls and horizontal elements, constraining the out-of-plane overturning of walls. An experimental campaign, to assess the effectiveness of this retrofit solution, was conducted in 2020−2021 at the EUCENTRE and University of Pavia facilities, Italy. Specifically, in-plane quasi-static tests were performed on unreinforced and retrofitted piers belonging to different masonry typologies, representative of existing Italian constructions. Experimental results indicated a promising improvement of the in-plane response of tested masonry piers providing a noticeable increase of deformation capacity, especially in cases of a predominant shear failure. 5,6 However, the obtained results inevitably refer to the specific features of tested masonry piers and to the specific details of the tested retrofit system. Indeed, laterally loaded URM members may exhibit different in-plane failure mechanisms depending on several external and internal factors, such as the aspect ratio, boundary conditions, vertical load, masonry bond pattern and mechanical properties. Failure modes influence the in-plane strength, displacement, energy dissipation capacities of piers, but can also modify the interaction between the installed retrofit and the masonry, affecting the associated performance improvement. Reliable numerical models calibrated and validated against the available tests, could enable the extension of the experimental outcomes, through simulations of a large variety of masonry structures which could enable the optimization of retrofit system detailing to further enhance seismic performance. In this regard, to simulate the response of strengthened masonry members, a reliable numerical model should be able, first, to account for the complex behavior of masonry subjected to lateral forces (e.g., mortar joints dislocation, crack localization, compressive, and tensile failure of masonry units, etc.), and in second place to explicitly account for the retrofit contribution.
Numerical modeling approaches usually employed in literature to simulate the response of masonry structures can be classified in three main categories with an increasing level of complexity: simplified, continuum, and discontinuum models. Several applications have demonstrated that simplified numerical techniques are able to satisfactorily simulate the seismic response of large-scale unreinforced and reinforced masonry structures with a relatively low computational cost and ease of calibration. [7][8][9][10][11][12] In continuum-based models (or macro-models), there is no distinction between units and mortar and the masonry is defined as an anisotropic continuous medium through plasticity models or macro-scale constitutive laws. 13 The damage within a masonry structure is then smeared through the continuous medium, simplifying the effects of crack localization, joint opening and frictional phenomena. 14,15 The limitations of simplified and macromodels can be overcome by using discontinuum models (or micro-models), which are able to simulate masonry material behavior, by explicitly representing units, mortar, and unit-mortar interfaces. 16 Despite the high level of accuracy, micromodeling cannot be easily employed to simulate the response of large structures because of the significant computational effort and large amount of input parameters, required for each masonry constituent. 17 In this study, a simplified micro-modeling approach based on the distinct element method (DEM) 18 was employed to replicate the in-plane response of masonry structures. In the DEM framework, adjacent units, that can be rigid or deformable, are ideally connected by zero-thickness spring layers and mechanically interact with each other along their boundaries according to contact constitutive laws representative of unit-mortar interface behavior. Although various literature works have demonstrated how DEM-based models are able to satisfactorily predict the response of URM structures, [19][20][21][22][23][24][25][26] the use of DEM to simulate potential retrofit solutions still remains a challenging topic, which has not yet been extensively explored.
This paper discusses an innovative numerical strategy to explicitly include the contribution of the investigated steelframed retrofit solution on the in-plane cyclic lateral response of masonry piers, using a DEM commercial software developed by Itasca, 3DEC. 27 Finite elements (FEs) and nonlinear links, available in the software environment, are employed to model the retrofit frames and their connections to the masonry leaf. The proposed modeling strategy is NOVELTY • Distinct element modeling of the in-plane quasi-static response of unreinforced and strengthened masonry elements. • Proposal of a novel mesh-refinement strategy to reduce computational effort and analysis time without losing the accuracy of the numerical outcomes. • Use of a new procedure to explicitly include the contribution of framed retrofit solutions in the distinct element method (DEM) framework. • Definition of a novel set of metrics to quantitatively assess the retrofit system performance.
validated by simulating the experimental in-plane response of unreinforced and retrofitted tested specimens, belonging to two different masonry typologies: one made up with solid clay bricks and lime mortar, arranged in a header bond pattern, and the other with typical Italian "Doppio UNI" hollow clay units and cement-lime mortar, assembled in a Flemish bond pattern. 5,6 A detailed discussion on the DEM idealization adopted for the two masonry typologies is provided, with a particular focus on the selection of the masonry unit discretization (i.e., mesh size), which sensibly affected pier numerical damage mechanism and strength degradation, to reduce computational effort without losing accuracy on final outcomes. Finally, retrofit performance benefits are assessed through a novel set of comprehensive metrics, selected to account for the level of damage in both mortar joints and units.

DEM FRAMEWORK FOR MASONRY STRUCTURES
In the DEM framework, masonry is represented as an aggregate of distinct brick-shaped blocks. When two blocks are detected to be in contact, subcontacts are generated along the contact surface. In this research, masonry units are modeled as deformable blocks, and divided into multiple finite-difference (FD) regions, consisting in constant-strain tetrahedral elements with three independent degrees of freedom. In the case of deformable blocks, the number of generated subcontacts for each contact depends on the FD mesh of interacting blocks, as shown in Figure 1. Zero-thickness contact springs are employed in both the normal and shear contact directions, with stiffnesses k n and k s , respectively. In the elastic range, contact stresses in the normal (σ) and shear (τ) directions are calculated as follows, depending on the related displacement increments ∆u n and ∆u s : Regarding the computational procedure, 3DEC employs a dynamic time-integration algorithm that solves the equations of motion for each block in the system using an explicit FD method. 27 At each timestep, the integration of the equation of motion provides the new block positions, and the contact displacement increments. The subcontact constitutive laws are then used to calculate the new subcontact forces, which are applied to the blocks in the next timestep. The explicit nature F I G U R E 1 Representation of masonry according to the distinct element method.

F I G U R E 2
Constitutive models adopted to describe the behavior of masonry units and joints. of the implemented algorithm can make DEM analysis computationally-expensive when models with numerous blocks are required. As discussed in detail in Sections 2.2 and 2.3, a variety of modeling techniques were employed to reduce the computational effort without losing accuracy of numerical outcomes.

Contact and block constitutive models
The shear behavior of detected contacts was described by a Coulomb-slip model, which requires the definition of cohesion (c) and friction angle (ϕ) parameters. The contact shear strength is then computed according to the following equation: Regarding the joint behavior in the normal direction, a joint tensile strength (ft mo ) was specified, while compressive failure was not considered. As shown in Figure 2, a new contact model recently proposed by Pulatsu et al. 26 was employed in this research. This model accounts for softening regimes in tension and shear through a fracture energy approach. Shear (G s ) and tensile (G tj ) fracture energy values were specified to control the contact post-peak behavior.
The employed shear contact model also allows to account for dilatancy. Despite the fact that the post-peak shearcompression behavior of masonry has been demonstrated to be significantly affected by dilatancy at least at a microscale level, 28,29 dilation angle usually decreases as compressive stresses and shear displacements increase, resulting in a smaller effect at the macroscale level. Moreover, the employed shear contact model allows only to define a constant dilation angle, which has been demonstrated to lead to unconservative overestimations of URM wall shear resistance. 29 In light of these observations, in this research DEM models were developed with a zero dilation as also adopted in similar numerical works available in literature. 22,24,26 However, the authors are aware that this assumption might lead to some conservative results especially when simulating the in-plane lateral response of strengthened masonry elements.
DEM has been extensively applied to simulate the response of blocky structures, especially when the damage, and the associated failure mode, is concentrated at the joint level. However, the proper modeling of shear-compressive and flexural-splitting behavior of masonry units, which are not usually accounted for in the standard formulation of DEMbased codes, remains a topic that has received less attention. Indeed, to accurately reproduce with DEM the complex nonlinear stress-deformation behavior exhibited by each masonry unit subjected to shear-compression, a fine FD mesh is usually required. Meanwhile, to account for flexural-splitting failure in units, a widespread approach consists in the introduction of a physical interface at the unit mid-length, which divides each block into two separate parts. 17 The highlyrefined mesh and the necessity of dividing the blocks in multiple regions (essentially doubling the number of blocks of the model), inevitably increase the computational effort and analysis time, especially for structures consisting of a large number of blocks. These aspects usually lead researchers to perform only monotonic pushover analysis, which have been demonstrated to be capable of correctly capturing the stiffness, strength, and failure mechanisms. However, available literature results for the in-plane response of URM structures usually come from experimental quasi-static cyclic tests, which provide additional information compared to monotonic tests, such as the progressive strength and stiffness degradation, and the energy dissipation capacity. Therefore, although very few works 22,24 have used DEM to simulate the cyclic response of URM structures, the direct comparison with experimental cyclic tests represents an important step to assess the reliability of the numerical models, validating the adopted modeling strategy.
In this contribution, the simplified numerical strategy proposed by Malomo et al. 22 was employed. Units were modeled as continuum deformable blocks discretized into multiple FD hexahedral regions, each constituted by two sets of five overlapping tetrahedral uniform-strain zones ( Figure 2). This discretization strategy is referred as mixed discretization and greatly improves the solution when plastic deformation occurs. 27,30 A Mohr-Coulomb plasticity model (MPM) was then assigned to the blocks to account for both masonry crushing and unit flexural-splitting failure. Therefore, the proposed strategy combines a discrete modeling approach, with masonry units modeled independently and constitutive laws describing their mutual interaction, and a smeared modeling approach to account specifically for masonry crushing and unit failure.
A linearized version of the Feenstra and De Borst 31 compression model originally developed for concrete, was implemented in 3DEC to simulate masonry crushing by assigning a strain-softening version of the MPM to the blocks. As discussed in Malomo et al., 24 by rearranging the failure envelope of a standard MPM in terms of principal stresses and assuming a zero block internal friction (ϕ b ), a simple correlation can be obtained between the cohesion assigned to the MPM blocks (c b ) and the masonry compressive strength (fc m ), resulting in c b = fc m / 2. A bi-linear stress-strain relationship is assumed up to the peak deformation ε p = (4 fc m ) /(3 E b ), where E b is the masonry unit elastic modulus. The first elastic branch is followed by a reduced one when the stress exceeds fc m / 3, as depicted in Figure 2. After attaining the peak value, the strength linearly reduces to a small residual value (herein assumed as 0.10 fc m to ensure numerical stability), in correspondence of the ultimate compressive strain defined as: where G c is the fracture energy of masonry in compression, d c is the ductility index parameter assumed as 1.6 mm consistently with a masonry compressive strength lower than 12 MPa, [32][33][34] while h is the crack bandwidth that is function of the adopted mesh size.
In the same plasticity model, a strain-softening tension model was also implemented in 3DEC to take into account the unit tensile failure. A linear stress-strain relationship is assumed up to the block/unit tensile strength (ft b ); after reaching ft b , the tensile strength linearly decreases to zero in correspondence of the ultimate tensile strain defined as function of the tensile fracture energy of masonry units (G tb ) 35 : The implementation of a fracture energy-based criterion including the crack bandwidth h, significantly decreased the mesh influence in simulating masonry compressive and unit tensile behavior. 24 However, the mesh size, especially for relatively coarse meshes, might affect the numerical flexural-splitting strength of units depending on the specific problem. This last aspect will be addressed in detail in the following sections.

Quasi-static modeling
The topic of this research is the numerical simulation of the in-plane quasi-static response of masonry structures. In 3DEC, quasi-static phenomena are modeled by adopting an approach conceptually similar to the dynamic relaxation. 36 By defining a suitable mechanical damping, equations of motion are damped to reach a force equilibrium state as quickly as possible under the applied initial and boundary conditions. In this work, adaptive global damping 37 was employed, which is a numerical servo-mechanism used to adjust the damping constant automatically. Viscous damping forces are used, and the viscosity constant is continuously adjusted to keep constant the ratio (R) between the power absorbed by damping (P) and the change of system kinetic energy (E k ). The adjustment to the viscosity constant is made to keep the ratio R equal to a target value herein assumed equal to 1.0. The solution scheme employed by 3DEC is conditionally stable. 27 A limiting timestep (∆t) that satisfies the stability criterion for both the calculation of internal block deformation and inter-block relative displacement is determined based on the specific numerical model. Since a higher value of numerical timestep corresponds to a smaller computational cost, specific assumptions were made in this research in order to reduce the analysis time, without losing accuracy of numerical outcomes.
Specifically, size, density, and time-scaling techniques were employed simultaneously in combination with the adaptive global damping, to obtain an acceptable compromise between the accuracy of results and the computational cost. 22,23 Size scaling was accounted for by adopting a reasonably coarse mesh for FD blocks, as function of the modeled masonry typology, positively affecting the numerical timestep. Time scaling was implemented by introducing a multi-stepped loading history: the horizontal velocity was increased gradually, leading to a displacement of 1 mm in 0.135 s; after that, velocities were set to zero, and simulation was continued until the average unbalanced nodal force resultant divided by the average load acting on the nodes became smaller than 10 −5 . Finally, density scaling was automatically performed by selecting adaptive global damping: an increase of density (ρ) corresponds to an increase of the critical timestep, proportional to √ρ.

Derivation of numerical material properties and mesh dependency considerations
This work refers to an extensive experimental campaign which included characterization tests performed not only on small-scale masonry samples (e.g., wallettes, triplets, etc), but also on units and mortars separately aiming at providing more specific ranges for mechanical properties, allowing an easier calibration of DEM models. All the details related to the performed characterization tests, carried out in 2021 and 2022 at the laboratories of University of Pavia (Italy), can be found in the research report by EUCENTRE. 6 However, a reliable characterization of some parameters, such as the elastic modulus (E b ) and the direct tensile strength (ft b ) of units, is challenging since it requires a specific experimental setup, which was not available. Therefore, some empirical correlations available in literature were employed for the derivation of the above-mentioned properties. Two different masonry typologies were considered in this work: the first one consists of solid clay bricks and lime mortar, arranged in a header bond pattern (identified in the following as S-HEA), while the second one of hollow clay "Doppio UNI" units and cement-lime mortar, assembled in a Flemish bond pattern (named as H-FLE).
The elastic modulus of units was evaluated using empirical formulations available in literature and already employed by several researchers, as a function of the unit compressive strength (fc b ). Different formulations 38,39 were employed accounting for the two different unit typologies considered, as summarized in Figure 3. The unit shear modulus (G b ) was then derived assuming material isotropy, with a Poisson coefficient (ν) of 0.25.
Then, starting from the assumed E b and the masonry elastic modulus (E m ) obtained through experimental compression tests, common homogenization formulae 22 were used to estimate the elastic and shear moduli of the mortar, required to compute the contact stiffnesses. Specifically, in this contribution the normal (k n ) and shear (k s ) stiffness assigned to the zero-thickness springs were evaluated using Equation (5), which should be applied only if unit deformability is adequately accounted for by the employed mesh 22 : where E mo and G mo are the mortar Young's and shear modulus respectively, while t mo is the thickness of mortar joints. For S-HEA masonry, made of solid isotropic bricks, the mortar shear modulus (G mo ) was calculated as for units, assuming material isotropy. Differently for H-FLE masonry, the presence of hollows in the units makes them behave as an anisotropic material. 5,6 Since the 3DEC strain-softening version of the employed MPM can only be employed for isotropic materials, to use the same modeling strategy also for H-FLE masonry, the anisotropy was accounted for by a proper reduction of the shear modulus of the mortar (G mo ), and therefore of the contact shear stiffness. The following equation, which considers brick and mortar as springs arranged in series, was employed: where G m is the masonry shear modulus and h b is the height of masonry units. In this research, G m was evaluated through experimental diagonal shear compression tests.

F I G U R E 3 Derivation of masonry unit properties and mesh influence on flexural behavior.
As anticipated, in the reference experimental campaign only indirect tensile splitting tests were performed on masonry bricks. Since the implemented MPM requires definition of block direct tensile strength, the empirical correlation 40 ft b = 0.85 ft b,sp , was employed to derive ft b as function of the experimental splitting strength (ft b,sp ).
Since one of the aims of this work is to reduce the computational expense usually associated with large DEM models, the FD discretization of blocks played a fundamental role. In this work, the mesh size is identified as A × B × C(R), where A, B, C, are respectively the number of hexahedral regions along the length, thickness, and height of masonry units as shown in Figure 3, where a 4 × 2 × 2(R) block is displayed. In the developed models, the mesh size was selected based on the specific problem (i.e., masonry typology, unit dimensions, bond pattern), and the smallest number of FD regions has been used to ensure the accuracy of the achieved results.
Equations (3) and (4) (see Section 2.1) include fracture energy-based criteria leading the simulation of the block softening regimes to be independent of the adopted mesh size. Moreover, fracture energies employed in the joint constitutive models have been demonstrated by Pulatsu et al. 26 to be essentially mesh-independent.
However, even if block compressive and pure-tensile behavior, as well as joint shear and tensile behavior, are not affected by the adopted mesh, block numerical flexural strength (ft b,flex ) might increase with the reduction of the number of FD regions. To investigate this aspect, S-HEA and H-FLE units were oriented according to the related masonry bond pattern and numerically tested in bending using three different levels of discretization, 2 × 1 × 1(R), 4 × 1 × 2(R), and 8 × 1 × 4(R) corresponding respectively to 2, 8, and 32 FD hexahedral regions. As can be observed in Figure 3, the use of a coarse mesh led to the overestimation of the target experimental ft b,flex , which in absence of available tests was evaluated as function of the available brick splitting strength. 41 To overcome this modeling issue, an equivalent tensile strength (ft b,eq ) should be iteratively selected to match the target value of the brick flexural strength.
It is important to point out that at least a 2 × 1 × 1(R) mesh is required to capture unit bending failure. Moreover, since the in-plane lateral response of piers is essentially a two-dimensional problem, in this research no unit discretization was applied along the masonry thickness. 22 However, it is worth to notice that by adopting a proper unit discretization along the wall thickness, the described numerical approach can be employed as well for the simulation of the out-of-plane response of masonry members. TA B L E 1 Experimental and inferred material properties for S-HEA and H-FLE masonry.

Validation on small masonry specimens
In this section, a preliminary verification of the proposed DEM modeling strategy against experimental results is presented, to verify the proper selection of material properties for the employed constitutive models. Cyclic compression tests on masonry wallettes, and shear-compression tests on masonry triplets were replicated numerically for both the considered masonry typologies S-HEA and H-FLE. Specifically, S-HEA solid clay bricks presented average dimensions of 120 × 250 × 55 mm and were assembled in a single-wythe header bond pattern. H-FLE hollow clay units with average dimensions of 250 × 120 × 120 mm were instead arranged in double-wythe Flemish-bond pattern. For both the masonry typologies and for all the samples, mortar joints had an average thickness (t mo ) of 10 mm. Exhaustive details about the experimental layouts, instrumentation setups and loading protocols can be found in the dedicated EUCENTRE research report. 6 All the properties assigned in the numerical models, obtained experimentally and inferred as discussed in the previous section, are listed in Table 1 (where ρ m is the masonry density). In S-HEA masonry, the unit dimensions and the header-bond pattern led to a significant number of blocks in the related DEM models. To reduce computational expense, a relatively coarse 2 × 1 × 1(R) mesh was adopted 22 and, according to results of numerical bending tests showed in Figure 3, an equivalent tensile strength ft b,eq ≈ 0.50 ft b, was assigned to the blocks to correctly simulate their bending response. Differently, for H-FLE masonry, the larger dimensions of units resulted in a significantly smaller number of blocks compared to S-HEA masonry, allowing the use of a more refined 4 × 1 × 2(R) mesh, which did not require the definition of a reduced block tensile strength ( Figure 3).
To shorten computation time, compression tests were simulated using a monotonic quasi-static loading as opposed to the cyclic one employed in the experiments. As shown in Figure 4, numerical models accurately represented the stiffness, strength, and softening regimes of masonry wallettes in compression, validating the adequacy of the block strain-softening compressive model and the employed empirical correlations. Note that the masonry axial stiffness was captured correctly, confirming that block deformability was adequately accounted for, even using a relatively coarse mesh.
Shear-compression tests were useful to assess the correct behavior of the spring layers, and especially to calibrate the shear fracture energy (G s ) to define the post-peak softening branch. For S-HEA masonry, a value of 50 N/m was selected, which represents a suitable value if compared to the ranges of G s commonly adopted in literature. 26,42 On the other hand, H-FLE masonry required different considerations because of the presence of vertical hollows in the units, which influenced the shear behavior of the brick-to-mortar bond as noticeable in the results of shear tests on H-FLE masonry triplets. As shown in Figure 4, the mortar that penetrated in the hollows restrained the formation of a smooth sliding surface, resulting in a large softening behavior. Accordingly, a value of fracture energy equal to 750 N/m was selected to match the mean slope of the softening branch. Since hollows were located only along the unit vertical plane, 5,6 consequently affecting only the shear behavior of bed joints, a lower value of 100 N/m was assigned to head and internal joints as a function of the experimental cohesion. 26,42 Finally, regarding the joint tensile fracture energy, the values of G tj reported in Table 1 were estimated referring to typical ranges commonly employed in research works on similar masonry typologies available in literature, 26,32,39,42 in the absence of dedicated experimental tests.

SUMMARY OF EXPERIMENTAL IN-PLANE CYCLIC TESTS
An experimental campaign aimed at studying the effects of the newly proposed retrofit solution was carried out in 2020−2021 at the EUCENTRE laboratory, Italy. A series of full-scale unreinforced and retrofitted masonry piers, belonging to different masonry typologies (i.e., S-HEA and F-HEA), were subjected to quasi-static shear-compression loading with different combinations of aspect ratios, overburden stresses (σ v ) and boundary conditions. Three loading cycles were performed for each target displacement (δ). The installed retrofit system, called "Resisto 5.9," consisted of modular steel frames connected to each other by means of 12-mm-diameter class 8.8 bolts and to the masonry through chemical anchors realized with 14-mm-diameter class 8.8 threaded bars. Each frame was made up of L-shaped border members and of horizontal and diagonal bracing members realized with plate sections, with a nominal yield (f y ) and tensile (f u ) strength of respectively 250 and 330 MPa. Figure 5A displays the retrofit layout with the typical dimensions of the modular frame. Further details on specimens, test setup and loading protocol can be found in Albanesi et al. 5 Among the several tested piers, two masonry specimens were considered in the present contribution, along with their strengthened counterpart. Note that, no connection was provided between the retrofit frames and the reinforced concrete (RC) top and bottom elements in the retrofit layout addressed in this paper ( Figure 5B,C). The considered specimen details are listed in Table 2, while the associated mechanical properties were already reported in Table 1.  Regarding the S-HEA masonry, unreinforced (S-HEA-U) and retrofitted (S-HEA-R) specimens exhibited a crack pattern typical of shear failures, characterized by diagonal cracks spreading from corner to corner of the piers. Cracks were mainly located in the mortar bed and head joints with a limited crushing of clay bricks, especially at the pier corners ( Figure 5B). The shear damage started developing from a drift ratio (θ = δ/h) of 0.10% for S-HEA-U and at θ = 0.15% for S-HEA-R. In both tests, the shear cracks grew in number and width with the increase of the lateral displacements, leading to a progressive reduction of lateral strength and stiffness up to specimen ultimate conditions. Tests were stopped to prevent the collapse at θ = 0.25% and θ = 1.00% for the S-HEA-U and S-HEA-R, respectively. As it is possible to observe from Figure 5B, the retrofit system allowed the specimen to increase significantly its ultimate displacement capacity, without affecting the URM pier initial stiffness and maximum base shear.
Similarly, both the H-FLE masonry specimens, H-FLE-U and H-FLE-R, displayed a shear failure with the development at θ = 0.10% of diagonal cracks initially located in the mortar joints. As the horizontal displacements increased, the damage in the units became more severe, with a progressive growth of the cracks in number and width, followed by the crushing of some units at θ = 0.30%, which corresponded to H-FLE U ultimate conditions. On the other hand, the strengthened specimen H-FLE-R reached its ultimate conditions at θ = 0.50%. As shown in Figure 5C, the retrofit did not affect the URM pier initial stiffness but allowed an appreciable increase in both maximum base shear and ultimate displacement capacity.
Overall, the retrofit guaranteed a significant improvement of the in-plane lateral response of the tested piers, especially considering the ultimate displacement capacity. Although H-FLE piers denoted an appreciable increment of the lateral strength, a smaller benefit was observed in terms of ultimate deformation in comparison to S-HEA specimens. This may be attributed to the poor tensile strength and brittle behavior of H-FLE units (Table 1), which led to higher damage within the masonry units ( Figure 5B). 5,6

SIMULATION OF THE IN-PLANE QUASI-STATIC RESPONSE OF URM PIERS
DEM models of the considered masonry piers were developed by directly reproducing the actual experimental setup. 6 The RC foundation and top beam were explicitly represented as an assembly of linear elastic hexahedral regions with an elastic modulus of 40 GPa. A relatively refined mesh was employed for these elements to increase the number of contact points along the masonry-to-RC interfaces, avoiding both stress-localization and interpenetration phenomena. Consistently with the experiments, a cohesion and a tensile strength of respectively 0.45 and 0.15 MPa were assigned to the joints between masonry and RC elements. The level of vertical overburden (σ v ) applied at the top of the piers was simulated in the numerical models by assigning an equivalent density to the RC top beam. Lateral in-plane displacements were imposed to the specimens by incrementally applying multi-stepped loading histories (see Section 2.2) to the RC top beam. Moreover, to ensure double-curvature boundary conditions, the rotation of the RC top beam was restrained, while its vertical displacement was left free to guarantee a constant level of applied vertical load throughout the simulation, consistently with the test setup.

Mesh refinement strategy
Before modeling the masonry specimens described above, some fundamental considerations about the impact of unit discretization (i.e., mesh size) on the numerical outcomes needed to be addressed. Despite all the adopted measures discussed in Section 2.2, DEM analysis still requires a significant computational effort. Because of the different studied masonry typologies, the two specimens presented a significant difference in terms of number of masonry units, which is proportional to the required analysis time. Specifically, the S-HEA pier was made of 684 bricks while the H-FLE pier was only made of 370 bricks. In this work, a significant effort was made to define a proper mesh size for each masonry typology, aiming at reducing the analysis time without losing accuracy on the final outcomes. Numerical evidence showed that the use of a coarse mesh, 22 even with a proper definition of ft b,eq , makes it difficult to correctly capture the behavior of masonry piers in the corners, especially when in these regions the failure mechanism is mainly influenced by the splitting of masonry units (see Figure 5B). This aspect might require increasing the number of block FD regions that, especially for S-HEA pier, prohibitively increases the computational expense.
To overcome this issue, a new mesh-refinement criterion based on the macro-distinct element model (M-DEM) discretization proposed by Malomo et al. 24 was employed. As shown in Figure 6, by defining the average slope (φ) of the lines connecting consecutive head joints as function of the masonry texture, it is possible to divide a masonry element into four different regions. Then, by subdividing each determined diagonal line (with length L d ) into three parts, a rectangular region can be defined in correspondence of each corner of the pier. The four identified "corner" regions represent the pier portions in which a refined FD mesh is required to correctly capture the unit flexural-splitting behavior. On the other hand, for the units belonging to the "central" region of the pier, a coarse mesh (with a proper ft b,eq , if necessary) can be employed.

F I G U R E 6
Mesh refinement strategy employed for S-HEA and H-FLE masonry.

F I G U R E 7
Influence of mesh refinement on pushover curve and damage pattern.
Accordingly, for the S-HEA pier an 8 × 1 × 4(R) mesh was adopted for the "corner" regions, assigning to these FD blocks the full tensile strength ft b , while a 2 × 1 × 1(R) discretization with ft b,eq ≈ 0.50 ft b was assigned to the units in the "central" region ( Figure 6). The values of tensile strength assigned to the blocks in the different regions of the pier were inferred from the results of numerical bending tests performed on masonry units with different mesh sizes (Figure 3). For H-FLE masonry, instead, considering the lower amount of employed bricks compared to S-HEA, a satisfactory compromise between computational effort and accuracy of the numerical outcomes resulted by employing a 4 × 1 × 2(R) mesh for the whole pier (Figure 6), without any distinction between "central" and "corner" regions. By using such a discretization, no ft b,eq was defined, since numerical flexural strength of H-FLE units was slightly overestimated if compared the inferred experimental one. However, as shown in Figure 3, no appreciable difference can be observed between a 4 × 1 × 2(R) and 8 × 1 × 4(R) mesh. Figure 7 shows the influence of the mesh refinement on the numerical monotonic "pushover" curve and the related damage pattern of pier S-HEA-U. In particular, the proposed mixed-discretization (8 × 1 × 4(R) + 2 × 1 × 1(R), see Figure 6) and a coarse mesh (2 × 1 × 1(R) for the whole specimen) are compared. As it can be noticed, no appreciable difference in the initial stiffness, peak strength and shear diagonal cracking mechanism can be observed at the varying of mesh refinement. However, if a coarse mesh is adopted, when a tensile failure occurs in a unit, the damage is rapidly spread among adjacent blocks resulting in an overestimated damage in the units at the pier corners also considering the necessary use of ft b,eq . Because of the spread and not localized damage in the corners, the 2 × 1 × 1(R) mesh essentially delayed the opening of the diagonal crack, underestimating the severe strength degradation exhibited by the specimen in the experiments. Differently, the proposed mesh refinement was able to satisfactorily reproduce the damage at the pier corners, capturing the progressive strength degradation with the increasing of lateral displacements. On the other hand, the benefits of adopting a refined 8 × 1 × 4(R) mesh in the corners of pier H-FLE-U was less obvious.

Numerical results
This section provides a comparison between the experimental and numerical response of the investigated URM piers. The capability of numerical models to reproduce the experimental response was assessed by performing both quasi-static cyclic and monotonic (i.e., pushover) in-plane analyses. Note that cyclic tests were simulated by performing only a single cycle for each target drift ratio, instead of the three performed in the experiments, to reduce the computational expense. Specifically, the S-HEA-U pier was tested by applying the numerical loading to progressively attain a drift ratio of 0.05%, 0.10%, 0.15%, 0.20%, and 0.25%. As in the experimental testing protocols, an additional cycle at θ = 0.30% was performed for specimen H-FLE-U. In the following, the obtained numerical curves are presented and compared against the experimental results. A comparison between the numerical failure mechanisms and their experimental counterparts is also presented, diversifying failures occurred in the joints, or in the units. Figure 8 compares the numerical and experimental in-plane responses of pier S-HEA-U. Although the numerical model underestimated the initial peak strength in both cyclic and monotonic analyses, satisfactory results were obtained in terms of initial stiffness, post-peak response and progressive strength and stiffness degradation ( Figure 8A,B). The observed differences between the numerical and experimental peak strength can be likely attributed to possible discrepancies in the numerical simulation of the vertical displacement imposed to the RC top beam in the experiments, 5,6 affecting the axial load and the lateral strength of the pier, or even to the intrinsic variability of the masonry mechanical properties (which were inferred experimentally and assigned in the models without any iterative adjustment). For this pier, the monotonic pushover curve slightly underestimated the magnitude of strength degradation observed in the experiments at large displacements (after 5 mm), while this was satisfactorily captured by numerical cyclic results. Moreover, a good agreement was found between the predicted in-plane failure mechanism and the experimental crack pattern. As shown in Figure 8C, the bi-diagonal stair-stepped cracks were adequately simulated by the DEM model. The refined 8 × 1 × 4(R) mesh employed in the "corner" regions enable effective simulation of the splitting of masonry units located at pier corners observed in the experiment. Note that in the presented numerical crack patterns (and in the other that follow), the joint opening is displayed as a function of the magnitude of the joint relative displacements (shear and normal). Consequently, hairline numerical cracks, as those formed in the initial stages of the pier response, cannot be totally appreciated.
The numerical and experimental results related to pier H-FLE-U are compared in Figure 9. The pier initial stiffness and strength, as well as the experimental diagonal shear failure were adequately reproduced by the numerical model. Although pushover curve matched closely with the experimental envelope ( Figure 9B), some discrepancies can be observed between numerical and experimental cyclic results ( Figure 9A). In particular, although the experimental shear-failure mechanism was correctly simulated, the numerical cyclic response presented wider hysteretic loops and a more pronounced strength degradation, if compared to the experimental one. Similar to specimen S-HEA-U, the numerical hysteresis denoted a more severe strength degradation compared to the pushover curve. In this case, the numerical pushover closely captured the magnitude of experimental strength degradation. A good match can also be observed between the predicted in-plane damage and the experimental crack pattern. As shown in Figure 9C, the DEM simulation resulted in a bi-diagonal crack pattern with damage spread in both mortar joints and masonry units. The severe damage of masonry units observed experimentally was satisfactorily captured.   Table 3 summarizes some of the most relevant parameters associated with the cyclic response of the investigated URM piers. Specifically, the initial stiffness (k el ), the maximum base shear (V max ) and the base shear in correspondence of the last testing cycle (V θmax ) are compared. k el was evaluated as a secant value at 70% of the pier maximum base shear, while values of V max and V θmax were reported for both the positive and negative loading directions. For each parameter, the percentage error between numerical and experimental values is also reported.

SIMULATION OF THE IN-PLANE QUASI-STATIC RESPONSE OF RETROFITTED PIERS
After validating the proposed strategy to model the in-plane lateral response of URM specimens, the numerical modeling of the retrofit system in DEM framework is addressed in this section. In the numerical models of retrofitted piers, the masonry contribution was modeled exactly as discussed in the previous sections. Any variation in the mechanical properties was not necessary since URM and strengthened piers were built using materials (i.e., units and mortars) coming from the same batches. Retrofit frames were then modeled by using structural-element and structural-link tools already available in 3DEC. 27

Retrofit system modeling
The modular steel frames were numerically reproduced by explicitly modeling each component as a one-dimensional beam FE with an isotropic linear elastic behavior, with no failure limit. Although this represents a simplified modeling assumption, it was assumed reasonable since experimental evidences showed that the retrofit behavior was mainly controlled by the failure of the retrofit-to-masonry anchors and by the elastic buckling of retrofit braces rather than the attainment of the axial or flexural strength of steel components. 5,6 FE beams were defined according to the retrofit scheme, connecting the actual positions of the retrofit-to-masonry anchors ( Figure 10A). Cross-sectional areas and moments of inertia were properly assigned to reproduce the actual axial and flexural stiffness of the steel members. Steel bolts employed to connect adjacent modular frames, were modeled as rigid links to reproduce the experimental retrofit layout shown in Figure 10A. An elastic modulus of 210 GPa, a Poisson coefficient of 0.30, and a density of 7850 kg/m 3 were assigned to steel retrofit components. It is important to note that in addition to border elements with a L-shaped cross-section of 88 × 80 × 3 mm, retrofit frames also consist of thin plate elements employed for the diagonal and horizontal braces (i.e., 40 × 3 and 50 × 3 mm sections, respectively). These members, as evidenced in the experiments, can be susceptible to buckling phenomena that can significantly limit their effective compressive strength. Thus, a sinusoidal shape initial deformation (i.e., imperfection) was assigned in the numerical models to each diagonal and horizontal brace to enable simulation of buckling effects 27 ; specifically, a small displacement of 1 mm was pre-imposed along the y-direction (i.e. perpendicular to masonry wall plane, see Figure 10), in correspondence of the crossing points of diagonal braces and of the mid-length of the horizontal ones.
As anticipated above, the interaction of the retrofit system with the URM panel was mainly affected by the behavior of retrofit-to-masonry anchors. Therefore, the following strategy was employed to characterize the connections between the FE beam frames and the assembly of FD blocks. As shown in Figure 10A,B, structural links 27 were defined at the actual position of chemical anchors. Specifically, a shear-yield deformable model was assigned to each link for the F I G U R E 1 0 Retrofit system modeling strategy: (A) retrofit numerical layout; (B) retrofit-to-masonry anchors model; (C,D) experimental and numerical anchors shear behavior along x and z directions. three translational degrees of freedom in directions x, y, and z, while rotational degrees of freedom were considered as free. Since the aim of this research is to model the retrofit contribution to the in-plane response of URM piers, particular attention was given in the definition of the structural link behavior along the x and z directions. 3DEC allows the definition of nonlinear behavior in the link by providing an input series of F i -δ i pairs, as shown in Figure 10B, where F i is the link force, while δ i is the link relative displacement. The structural link behavior along x and z directions was calibrated against the results of shear tests performed on retrofit-to-masonry anchors for both the masonry typologies during the same experimental campaign. 5,6 A comparison between the experimental shear force-displacement behavior of anchors and the numerical rule assigned to the structural links is reported in Figure 10C,D, respectively for S-HEA and H-FLE masonry: specifically, F i -δ i pairs were defined to match the median curve (50th percentile) obtained by the fitting of experimental results. Similarly, referring to the experimental results of extraction tests on retrofitto-masonry anchors, 5,6 a simplified bilinear force-displacement rule, with an elastic stiffness k y = 25 kN/mm and a yielding force F y = 11 kN was assigned to the link translation degree of freedom along the y-direction (i.e., perpendicular to the wall).

Numerical results
In this section the numerical responses of the retrofitted piers, namely S-HEA-R and H-FLE-R, are compared against the experimental results and against their numerical unreinforced counterparts. The addition of the retrofit frames in the models, combined with the larger target drift ratios reached in the experimental cyclic loading protocols for retrofitted piers, significantly increased the computational effort resulting in a prohibitive analysis time required for the simulation of the entire experimental cyclic tests. The numerical results presented in Figures 8 and 9 demonstrated that monotonic pushover analysis satisfactorily captured the fundamental parameters of in-plane behavior, such as stiffness, strength, failure mechanism, and progressive strength degradation with the increasing of lateral displacements. Therefore, in this section the numerical performance of the investigated retrofit solution was investigated through monotonic pushover analysis. However, for the sake of completeness, experimental cyclic tests were also numerically simulated, but only up to the drift ratios corresponding to the ultimate conditions of the unreinforced specimens (i.e., 0.25% and 0.30% for S-HEA and H-FLE masonry, respectively). In Figure 11A,B, the numerical pushover curves for the retrofitted piers are compared with experimental cyclic results and related backbone curves. As in the experimental results, the pier initial stiffness was not affected by the addition of the retrofit frames in the numerical simulation results. Regarding the strength, the S-HEA-R numerical model was not able to fully capture the specimen maximum strength of 238 kN, instead predicting a 15%-lower maximum base shear of around 200 kN, as consistently observed for the unreinforced numerical model ( Figure 8B). Also, for pier H-FLE-R, an underestimation of around 10% of the experimental peak strength of 254 kN can be observed. However, since in this case the H-FLE-U model correctly captured the specimen maximum base shear ( Figure 9B), the underestimation can be potentially attributed to the force-displacement behavior assigned to structural links. Indeed, although in the numerical models a unique rule was equally assigned to all the retrofit-to-masonry connections, experimental tests revealed (especially for H-FLE masonry, Figure 10D) a significant variability in the shear behavior of anchors, depending on their location on the wall surface, that is, in the masonry units or in the mortar joints. Moreover, especially for the strengthened piers, the numerical underestimation of the experimental peak strength might be attributed to the conservative zero-dilation assumption. Indeed, since the retrofit frames are applying a sort of confinement effect to the masonry, the dilatancy may influence the pier shear response resulting in a higher peak strength. 29 However, while the strength of retrofitted piers was underestimated, an important overall behavior change caused by the retrofit solution was still captured. Indeed, by comparing the numerical pushover curves of unreinforced and retrofitted specimens showed in Figure 11A,B, the numerical models predict that the retrofit system was able to increase the pier ultimate displacement capacity, delaying the strength degradation observed in the unreinforced models without significantly affecting the pier initial stiffness and strength, as observed experimentally. The ultimate drift ratio of piers subjected to increasing monotonic in-plane loads was identified on pushover curves by a strength reduction equal to the 20% of the maximum attained base shear; this resulted in a drift ratio of 0.41% for S-HEA-U and of 1.40% for S-HEA-R, and a drift ratio of 0.36% and 0.64% for H-FLE-U and H-FLE-R, respectively.
The benefit of the tested solution can also be appreciated considering the numerical failure mechanisms reported in Figure 11C,D. Herein, the damage pattern at the ultimate conditions of the unreinforced piers is compared, at the same drift ratio level, with the one of retrofitted piers. Moreover, the damage pattern of the retrofitted piers at ultimate conditions is also reported. Note that to clearly appreciate retrofit benefits, joint relative displacements are displayed with the same order of magnitude for both unreinforced and retrofitted specimens. As can be observed, the retrofit solution did not alter the initial diagonal shear failure mechanism exhibited by the unreinforced piers. However, once the shear mechanism formed, the steel frames restrained the relative displacement of masonry portions identified by the diagonal crack, reducing the joint opening and spreading the damage over the masonry pier surface, for the same imposed drift ratio. Figure 12A,B compare the numerical and experimental in-plane cyclic responses of the retrofitted piers, which were simulated only up to the ultimate conditions of the unreinforced specimens, as mentioned above. The numerical backbones of the unreinforced specimens are also reported in the same figures. Despite a lower maximum strength if compared to the experiment, S-HEA-R pier presented a numerical response that closely matched the experimental one. Differently, H-FLE-R numerical response showed similar discrepancies already observed for the unreinforced counterpart, resulting in a more pronounced strength degradation and in a wider hysteresis if compared to experimental loops. Consistently with pushover analysis, the pier initial stiffness and failure mechanism were not affected by the retrofit frames.
The retrofit performance can be clearly appreciated in Figure 12C,D, which compare the numerical and experimental damage patterns of unreinforced and retrofitted specimens at the same level of drift corresponding to the ultimate conditions of URM piers. If compared to unreinforced models, the retrofitted pier damage patterns presented only hairline shear-diagonal cracks with a decreased number of failures in the masonry units, as consistently observed in pushover analysis. In the same figure, the numerical crack patterns are compared also with the corresponding experimental ones: the DEM prediction of retrofit effects on the specimen lateral behavior closely simulated the damage-reduction and spreading effects observed in the experimental tests. 5,6 Although the promising results obtained by the partial simulation of retrofitted specimen cyclic responses, the required computational cost did not allow to perform the cyclic analysis until the ultimate conditions of retrofitted piers. Therefore, it was decided to move the attention on monotonic analysis, which, within a reasonable computational time window, allowed to study the retrofit performance for increasing levels of imposed drift ratios with an acceptable degree of accuracy.
Since the primary benefit of the retrofit was to delay and counteract the opening of shear cracks enhancing pier lateral displacement capacity, the ultimate drift ratio, usually estimated for masonry panels as the drift ratio associated with a  strength reduction equal to the 20% of the maximum attained base shear, can be employed as a metric to assess the retrofit performance benefits. In Table 4, some of the most relevant parameters of the in-plane lateral response 43 of the investigated specimens, that is, the maximum base shear (V max ), the associated drift ratio (θ Vmax ) and the ultimate drift ratio (θ 20% ), are summarized, providing a comparison between experimental and numerical results. The non-dimensional ratio between retrofitted and unreinforced parameters is also provided to highlight the retrofit performance improvement. As can be noticed, although the numerical ultimate conditions of retrofitted piers were studied using monotonic loading rather TA B L E 5 Retrofit system performance benefits. than the experimental cyclic one, the experimental response improvement was satisfactorily simulated by the numerical models, especially considering the displacement capacities. Moreover, reflecting experimental results, a more significant enhancement was observed for S-HEA masonry. However, the 20%-strength-drop criterion cannot be generally applied, especially when the retrofit solution induces mechanisms, like sliding-shear modes, which are not commonly associated with a severe strength degradation. In order to overcome this issue, three alternative metrics are proposed herein to quantitatively assess the retrofit system performance benefits: the maximum normal or shear joint relative displacement (δ j,max ), selected as the parameter representative of the damage in mortar joints, and the maximum plastic tensile (ε t,max ) and shear (ε s,max ) strain of blocks, selected to give an effective indication of the damage experienced by the masonry units. Table 5 summarizes and compares the above mentioned metrics, evaluated through pushover analysis at certain specific drift ratios. In particular, the damage metrics of unreinforced and retrofitted specimens were evaluated first at the ultimate drift ratio of URM piers observed in the experimental tests (Section 3), and then at the drift ratios corresponding to numerical ultimate conditions of unreinforced piers inferred through pushover analysis as discussed above. Results at the ultimate drift ratio (θ 20% ) of retrofitted piers are also reported.
The obtained results revealed a significant reduction of joint and block damage metrics provided by the application of the retrofit solution. For example, considering S-HEA masonry at θ = 0.41%, the retrofit was able to reduce the joint opening (joint relative displacement) by 38%, the block maximum tensile strain by 62%, and the block maximum shear strain by 57%. For H-FLE masonry at a drift ratio of 0.36%, the retrofit reduced the joint opening by 44%, the block maximum tensile strain by 61%, and the block maximum shear strain by 37%. Moreover, looking at the ultimate drift ratios of retrofitted piers, the retrofit allowed the piers to experience higher values of δ j,max , ε t,max , and ε s,max before reaching the 20% strength reduction. The proposed set of damage metrics can also be employed to assess the ultimate conditions of retrofitted panels without relying on the strength drop of 20%, by evaluating the drift ratio at which the metrics of retrofitted piers become equal to the ones associated with URM piers ultimate conditions. By using this approach, a good prediction of ultimate drift ratio was obtained, resulting in θ = 1.10% and θ = 0.52% for S-HEA-R and H-FLE-R, respectively (see Table 5). Figure 13 shows the performance of the retrofit solution as a function of the lateral displacement (or drift ratio) imposed to S-HEA and H-FLE specimens. Specifically, the retrofit performance benefit was evaluated at consecutive values of θ and up to the ultimate conditions of unreinforced piers as the ratio between the proposed damage metrics for unreinforced piers versus retrofitted piers. As it can be noticed, because of the lower stiffness of the retrofit system compared to the masonry pier, no appreciable benefits can be observed for θ < 0.10%, denoting that the retrofit was not able to modify the onset of masonry damage and was effectively engaged only once some damage had already occurred, reflecting experimental observations. Indeed, just like numerical models, also in the experimental tests the diagonal shear failure formed at a drift ratio of around 0.10% for both the masonry typologies (see Section 3). Then, after masonry shear cracking (i.e., θ > 0.10%), the retrofit contribution became appreciable, as denoted by a ratio higher than 1.0.

CONCLUSIONS
This paper presented a computational modeling investigation, based on the DEM, of the in-plane lateral response of unreinforced and retrofitted masonry structures. The studied innovative retrofit solution consists in modular steel frames connected to the masonry external surface by means of chemical anchors and among them through steel bolts. The effectiveness of the retrofit solution was investigated through in-plane quasi-static experimental tests performed on identical masonry piers in both unreinforced and strengthened conditions at the EUCENTRE laboratory. Within this framework, a strategy to model URM structures was developed and validated against experimental tests. The proposed model has the ability of simulate all the possible failure mechanisms of masonry components, that is, mortar joint tensile and shear failure, unit flexural and splitting failure due to tensile stresses, and unit crushing due to highcompressive stresses. Moreover, particular attention was given to the simulation of the post-peak softening behavior and to the selection of the discretization of masonry units and its influence on numerical accuracy. To find a compromise between the analysis time and the numerical accuracy, a new mesh-refinement strategy was proposed and was found to be essential in some cases to effectively capture the actual masonry failure mechanism.
Two different masonry typologies, tested within the same experimental campaign, were studied in this paper: the first consisted of solid clay bricks and lime mortar, arranged in a header bond pattern, while the second one was made up of typical Italian "Doppio UNI" hollow clay units and cement-lime mortar, assembled according to a Flemish bond pattern. The simulation of the lateral response of URM piers with DEM revealed a good numerical prediction of the experimental damage pattern, hysteretic response, and progressive lateral strength and stiffness degradation.
An innovative modeling strategy to include the contribution of the investigated steel-framed retrofit solution in DEM framework was then presented. Despite some differences in terms of maximum lateral strength, the employed model for the strengthening system was able to adequately capture the experimentally observed performance improvement for the lateral response of investigated piers. Numerical outcomes demonstrated the capability of the implemented retrofit in spreading the damage without modifying the URM initial failure mechanism, consequently delaying the strength degradation and enhancing the pier displacement capacities, in accordance with the experimental observations. A novel set of metrics was proposed to assess the retrofit performance improvement, accounting for the level of damage in both mortar joints and masonry units.
In order to generalize the experimental results and to further investigate the effectiveness of the retrofit system, the herein proposed modeling strategy and damage metrics will be employed to perform a series of parametric analyses aimed at assessing the retrofit performance for varying pier aspect ratios, imposed vertical loads, boundary conditions, bond patterns, and masonry typologies, as well as retrofit system details. Moreover, the effectiveness of the retrofit solution in the more realistic scenario of multi-storey masonry façades will be addressed, with a particular focus on the connections between retrofit modules and floor diaphragms, as well as to the interaction between retrofit frames belonging to consecutive storeys.

A C K N O W L E D G M E N T S
The first author gratefully acknowledges the financial support provided by the University School for Advanced Studies IUSS Pavia. The experimental and numerical research activities described in this paper were performed at the EUCENTRE Foundation in Pavia, Italy, and were funded by Progetto Sisma s.r.l. The study of the modeling approach proposed for unreinforced masonry specimens was partially supported and funded by DPC-ReLUIS Work Package 10 Subtask 10.1 (2022−2024) "Modelli di capacità locali e globali per la definizione degli stati limite, definiti in funzione del metodo di analisi". The received financial and technical support is gratefully acknowledged. Special thanks go also to F. Graziotti, G. Magenes, and D. Malomo whose valuable advice was essential to the study.