Seismic Assessment of Existing Masonry Buildings Using Damage Mechanics

Abstract

This paper presents research concerning the numerical simulation of existing masonry buildings when subjected to pushover analysis. A nonlinear static analysis is undertaken using the commercial software ABAQUS standard, in which masonry structures are modelled using damage mechanics. To validate the chosen input parameters, this study compares two different approaches for static nonlinear modelling, the Finite Element Method (FEM) and the Equivalent Frame Method (EFM), for a simple masonry building. The two methods are compared using the guidelines from Part 3 of Eurocode 8. This study identifies the advantages and disadvantages of various modelling approaches based on the results obtained. The results are also compared in terms of capacity curves and damage distributions for the simple case study of a masonry building created to compare numerical methods. Subsequently, nonlinear pushover analyses with ABAQUS (FEM) were performed on the North Tower of Monserrate Palace, Portugal, in which the material parameters were calibrated by considering the results of dynamic characterisation tests conducted in-situ. Regarding the circular body of Monserrate Palace, the damage distribution of the structure is analysed in detail, aiming to contribute to the modelling of such structural configurations through the Equivalent Frame Method.

1. Introduction

Not only in underdeveloped countries but also in industrialised nations, the masonry construction typology is extensively employed in construction practices and makes up a sizeable share of the stock of residential and commercial buildings. According to [1], the percentage of masonry buildings in some Asian countries is around 90%, in South America it is around 70%, and in some countries in Europe it passes the mark of 60%. However, masonry has been losing market share globally despite being a sustainable construction choice due to its thermal and acoustic efficiency, fire resistance, durability, and simple construction technology. The primary cause of this is the emergence of new, less vulnerable options for low- to medium-rise structures, including steel and reinforced concrete, which are less susceptible to earthquakes than masonry structures in seismic zones. Even so, masonry construction is still widely used in areas susceptible to earthquakes [2].

Portugal, a country with a medium to high seismic hazard, is now facing a housing crisis due to the significant influx of people into large cities [3]. With the rising cost and reduced construction of new properties, there has been a growth in the rehabilitation of existing buildings, typically masonry structures. Standardisation and the growing economic interest in this activity resulted in the publication of Decree-Law No. 95/2019, which introduced significant changes, including repealing national structural regulations in favour of the Structural Eurocodes and the imposition of seismic vulnerability reports for rehabilitation interventions.

A large part of the built masonry stock is composed of old structures with enormous walls of weak bricks or stones, weak mortar, and weak connections between orthogonal walls and floors. Because of these characteristics and the common presence of flexible diaphragms that are only sporadically attached to walls, old masonry structures typically exhibit localised out-of-plane damage mechanisms [4]. Unlike ancient masonry buildings, modern structures have limited wall thicknesses and a regular brick-masonry layout. Independently of the type of masonry, it is necessary to guarantee strong floor-to-wall connections and rigid floors in order to ensure that box behaviour is used to accomplish the global seismic response [5].

Existing structures in seismic regions, especially those with a high level of anticipated seismic risk, are increasingly required to undergo seismic assessment. Due to their low tensile strength, low ductility, and limited energy dissipation capacity, unreinforced masonry structures have a significant seismic susceptibility, as demonstrated by previous seismic occurrences [6]. The building structures have a comparatively poor level of seismic safety compared to the requirements of the current standard earthquake engineering practice.

The significant seismic vulnerability of masonry buildings is also a result of the absence of seismic design guidelines, as these structures were frequently created primarily for vertical loads. Research studies have been conducted to enhance masonry structural systems under seismic actions and develop guidelines and tools for their seismic design [7,8]. However, experimental tests are costly and time-consuming. Therefore, different numerical approaches have been developed to represent the complex behaviour of masonry structures and perform advanced computational nonlinear analyses.

It is well-recognised that nonlinear dynamic analysis is the most accurate method for simulating and evaluating a structure’s seismic reaction. However, because of the intricacy of its use in engineering practice, it necessitates a high level of expertise, a significant investment in computational time, and a high calibration cost for the cyclic constitutive laws. The seismic input employed in the dynamic analysis significantly impacts the response of structures. Moreover, the absence of uniform verification protocols makes it challenging to assess a building’s seismic reaction based solely on dynamic analysis results. However, because the material response is significantly nonlinear, regardless of the low amount of loading, linear elastic analysis does not accurately describe the behaviour of a masonry construction. Consequently, nonlinear static (pushover) analysis has been frequently chosen for the seismic design and assessment of structures.

In earthquake engineering, pushover analysis is a widely used method and one of the most adequate instruments for determining the seismic safety of any structure, new or old. It must be implemented correctly and adhere to several requirements. Numerous methods exist for implementing and applying pushover analysis, considering factors like force versus displacement control, different load patterns, higher mode inclusion, and adaptive load patterns. As a result, the parameters selected directly impact the outcomes. The type of lateral load applied, whether a force or a displacement, is crucial to the accuracy of pushover techniques [9]. The goal is to determine a capacity curve that most closely resembles the seismic behaviour of the structure.

Many practices follow the general method of considering a pushover analysis finished if a capacity curve is established to consider the capacity up until a minimum 20% reduction in the maximum obtained base shear. An alternative strategy is to consider the point at which a single pier or a subset of piers reach certain drift thresholds and then activate a soft-storey mechanism in a wall/level. When performing pushover analysis, nonlinear static procedures are the most practical way to determine the target displacement of structures (i.e., the maximum inelastic displacement requirement corresponding to a particular level of seismic activity).

Nevertheless, the methods are mostly suitable for standard buildings, and their response does not account for higher mode impacts. Included in the Eurocode 8 [10], the N2 technique [11] is a nonlinear static procedure that produces findings acceptable for steel-framed and reinforced concrete structures, often categorised as medium- to high-period structures. Concerns have been raised about the method’s capacity to estimate findings for short-period masonry constructions accurately [12]. The fundamental nonlinear static procedure assumes that a building experiences gradual, monotonic lateral loading, which changes the building’s displacement response. Moreover, a nonlinear static procedure transforms a multi-degree-of-freedom structure into a single-degree-of-freedom system, in which the primary translational mode shape alone determines the dynamic behaviour.

Using sophisticated and precise methods [13] is essential when considering that many built-stock structures are part of the historical and architectural legacy [14,15]. This essay’s subjects are masonry towers, regular bell towers, fortifications, and chimneys. Nonetheless, by considering the commonalities and the universal method of pushover analysis implementation, practitioners can apply many of the treated elements of numerical techniques to a wide variety of structural typologies.

Masonry towers are extremely unusual structural typologies, and like all historical masonry projects, they are usually designed to support only vertical loads [16,17]. Therefore, they are anticipated to be highly vulnerable to seismic events, and special consideration should be given to their seismic assessment in the context of a preservation strategy. The two primary factors contributing to the anticipated high susceptibility are the poor mechanical qualities of the masonry material under strain and the geometric characteristics, such as apertures, imperfections, and slenderness.

Estimating the uncertain factors (material properties and boundary conditions) influencing the structural behaviour is crucial, regardless of the analytical methodology employed [18]. To identify the seismic vulnerabilities in a historical structure, a multidisciplinary approach using cutting-edge research and simulations is often advised [19,20,21]. A thorough approach to seismic safety would involve progressing from sophisticated numerical models, including in-plane and out-of-plane local studies, to diagnostic examinations.

As previously discussed, at the moment, the use of nonlinear analysis to perform structural design is gaining popularity among structural designers. The possibility of using several different structural software packages (ABAQUS [22], ADINA [23], ANSYS [24], ATENA [25], DIANA [26], 3DEC [27], 3Muri [28], HiSTrA [29]) and a set of structural design guidelines [30,31] supports extra confidence in their use by structural engineers. The use of nonlinear analysis is not intended to solve simply supported beams but complex structures with unknown stress distributions, which are found in the case of masonry buildings [32].

The abovementioned software systems use three methodologies to perform nonlinear analyses, classified as macro-, simplified, and detailed micro-modelling, as discussed by several authors in the literature [33,34,35]. In macro-modelling, masonry is represented as a homogeneous continuous medium without differentiation between the units, making it highly attractive due to its lower computational requirements. Simplified micro-modelling defines masonry units but excludes mortar joints by lumping them in discontinuous elements. This method leads to a good compromise between computational effort and accuracy. On the other hand, detailed micro-modelling considers the definition of all masonry elements (units, mortar, and unit–mortar interfaces), which leads to a high computational load for large-scale structures. A very widely used macro-modelling approach in masonry structures is the Equivalent Frame Method (EFM) because of its simplicity in terms of input requirements and very low computational consumption. The method was developed specifically to represent masonry behaviour; however, due to its simplicity, it is limited in terms of representing irregular geometries, such as curved surfaces, and only analyses in-plane behaviour. Even so, this method has been applied to large-scale historical masonry structures with overall good results [36,37,38]. Within the macro-modelling approach, it is important to highlight the Discrete Macro-Element Model (DMEM) introduced by Caliò et al. (2012) [39]. Originally, this element was designed to simulate the nonlinear in-plane behaviour of unreinforced masonry walls. It was subsequently upgraded for application to monumental masonry buildings, enabling the simulation of both in-plane and out-of-plane behaviour [40]. The DMEM strategy was recently implemented in commercial software that is able to automatically perform safety verification processes [29].

Another widespread approach for the modelling of masonry structures is the Finite Element Method (FEM), which offers the possibility of performing either macro- or micro-modelling. Macro-modelling is, again, the most commonly used approach due to its already mentioned benefits; plus, with the FEM, there is the possibility of considering the out-of-plane behaviour of walls, the damage propagation is much more detailed, and it is possible to model highly irregular structures. However, it is also worth noting that when compared to the EFM, the FEM, even when using the macro-modelling approach, requires the introduction of material parameters that are usually unknown to the engineers of masonry structures, and the time-consuming and computational efforts increase significantly. Therefore, for the modelling of masonry structures, the FEM is mainly applied to highly irregular or complex structures where the EFM cannot be used [14,41]. In addition, lately, simplified micro-modelling approaches have started to be adopted as an alternative to detailed micro-modelling, such as Discrete Element Modelling (DEM). Different authors have demonstrated it to be a suitable tool for the assessment of masonry structures, presenting studies with regular and irregular structures, including curved surfaces, and even the application of retrofitting techniques [42,43,44]. Moreover, some authors have started to compare the described approaches applied with different software to irregular masonry structures [38,45], which is essential for practitioners when deciding which methodology to apply.

The work presented here emphasises methods that any Finite Element Method (FEM) software could use to conceptualise a nonlinear static analysis of rubble stone masonry. It examines how shape, material characteristics, and load patterns affect results. Given the high likelihood of brittle collapse in regular masonry constructions, one might anticipate that realistic simulation behaviour will show some discernible softening. In this context, the use of numerical simulations becomes extremely important. The total modelling and analytical procedure is complex due to the unique properties of the masonry material, which include a brittle failure in shear and a relatively low tensile strength [46]. Full 3D solid elements have been used to perform static nonlinear analyses of existing masonry walls [47], in which brittle material behaviour was modelled using damage mechanics [36].

The main purpose of this study was to carry out a seismic safety assessment of a simple masonry building using the Equivalent Frame Method (EFM) with 3Muri [28] and the Finite Element Method (FEM) with ABAQUS standard commercial software [22]. The aim was to explore and compare modelling approaches, identifying the advantages and disadvantages of the two software packages. The seismic safety analyses were performed according to the earthquake assessment requirements of Part 3 of Eurocode 8 [48], applying the pushover methodology.

In addition, pushover analyses were carried out on the North Tower of Monserrate Palace, a structure of notable architectural value with a circular profile, using the Finite Element Method. This model was calibrated using the results of dynamic characterisation in-situ tests of the structure. A detailed analysis of damage distribution in the circular tower was also performed to investigate thoroughly how damage propagates throughout the circular structure under various loading scenarios. This study aimed to advance the modelling of similar structural configurations through the EFM. By scrutinising the damage distribution, insights can be gained into the behaviour of the tower’s components and their interactions, thus enhancing the accuracy and reliability of structural predictions using the EFM.

Even though some works have already been carried out, applying nonlinear static analysis for existing heritage masonry buildings and considering out-of-plane mechanisms and circular façades still requires more development towards defining simplified methods that are attractive to practitioners. This work planned to fill this knowledge gap by providing outputs with the ABAQUS standard, which offers a wide range of modelling tools for simulating the nonlinearity of masonry buildings. The numerical simulations obtained some important information on the seismic assessment of masonry towers. Therefore, structural designers can readily apply the discoveries presented here for more realistic and accurate outcomes.

2. Model Validation for a Simple Masonry Building

The study starts with the modelling of an existing residential building using different approaches and two software systems. The case study is a hypothesised masonry building with a regular geometry encompassing three floors.

2.1. Modelling with the EFM

This method proposes dividing the vertical resisting elements into equivalent frames, based on years of studies on damage to mixed and masonry structures caused by seismic activity. This division applies to deformable elements, which are divided into piers and spandrels (where the nonlinear response of the structure is concentrated), and nondeformable elements, designated as rigid nodes. The piers are responsible for bearing most of the loads, while the spandrels limit the rotation of the piers. Both the piers and the rigid nodes are formed due to the presence of openings, with the rigid nodes appearing between the ends of the openings and the spandrels occupying the space between the rigid nodes. The software used in this study is called 3Muri [28].

2.1.1. Description of the Building and Seismic Characterisation

The building presents front and side façades with 0.45 m and 0.40 m thicknesses, respectively, and side walls 0.25 m thick, while the interior walls present a thickness of 0.20 m. The interior doors throughout the building are consistent on all floors, with a height of 2 m and a width of 1 m. However, the exterior openings vary by floor. On the ground floor, the openings correspond to doors that are 2.5 m high and 1.35 m wide. For the other floors, the external masonry features openings with thresholds 0.7 m high and 1.35 m wide and lintels of 0.9 m. A 3D representation of the structure can be found in Figure 1. This building is situated in Lisbon, in Seismic Zone 1.3/2.3, on a type C site, as defined in the Portuguese National Annex to Part 1 of Eurocode 8 [49]. Earthquakes of type 1 and type 2 were simulated separately.

2.1.2. Material Properties and Dead Loads

The exterior walls are constructed from rubble stone masonry, while the interior walls are made of brick masonry. The mechanical properties of the masonry were determined using the Italian NTC-2018 code [50] and are presented in

Table 1. Masonry mechanical properties.

	E [GPa]	G [MPa]	w [kN/m3]	fm [MPa]	τ [MPa]	γm [Safety Factor]
Stone masonry	1	0.33	19	2.0	0.017	1.5
Brick masonry	1.2	0.40	18	2.6	0.138	1.5

. $E$, $G$, $w$, $f_m$, $\mathsf{\tau }$, and ${\gamma }_w$ are, respectively, the Young’s and shear moduli, the specific weight, the average compressive strength, the average shear strength without axial action, and the safety factor. Following Part 3 of Eurocode 8, a knowledge level factor of 1.2 was used. During nonlinear static analysis, materials undergo cracking and stiffness loss. Therefore, the Young’s modulus was reduced by 50% in the EFM model. The roof was also modelled as a horizontal timber floor, but its weight was applied with linear loads at the top-floor external walls supporting the roof timber structure. The material properties for the timber were defined as pine, according to LNEC’s N2 Datasheet [51], with an average density of 580 kg/m³, a mean Young’s modulus ($E)=\mathrm{12,000}$ MPa, a mean shear modulus ($G)=750$ MPa, and a mean tensile strength parallel to the fibres ($f_t)=18$ MPa. The floor-structure equivalent shear modulus was calculated according to New Zealand norms (NZSEE-2006) and the experimental tests performed by Giongo et al. (2014) [52], resulting in a $G_{eq}=11$ MPa.

2.1.3. Equivalent Frame Modelling

The final model is illustrated in Figure 1c, in which the green elements are the spandrels, the blue elements are the rigid nodes, and the orange elements are the piers. The floor elements are modelled as orthotropic membranes.

2.2. Modelling with FEM

This section describes how the numerical model was assembled with the ABAQUS standard for the same case study.

2.2.1. Geometry Assembly

In ABAQUS, the modelling process was organised into parts based on materials. Initially, the exterior walls (Figure 2c) were modelled, and, subsequently, the openings were assembled using the Cut/Extrude feature. The interior walls (Figure 2b) were modelled floor by floor, with each part later copied to a different floor. Finally, the timber components (Figure 2a) were designed to fit the gaps on top of the interior walls. These timber elements were then inserted into the Assembly environment, and after performing an Instance/Union operation, the complete geometry of the building was obtained. The results of this process can be observed in Figure 2d, in which all elements are continuously connected in the end, being allowed to generate a conforming mesh without displacement discontinuities between nodes. In Figure 2d,e the different colours represent the different materials considered for the walls.

2.2.2. Type of Analysis

The model used a physically and geometrically nonlinear analysis with force control applied at all nodes that present mass due to self-weight. The geometrically nonlinear analysis was necessary since the high vertical displacements could potentially increase the stress in the foundations, thereby increasing the level of nonlinearity. It was expected that, due to the brittle behaviour of the masonry structure, softening may have occurred in the static nonlinear analysis; therefore, to control the decrease in force in the final structural response, a variation of the arc length method known as the modified Riks algorithm [53] was applied.

2.2.3. Constitutive Relations

The main material properties, including the density and elasticity, are the same as the ones admitted in 3Muri. Due to the existence of high spurious modes during the modal analysis due to the vibration of the wood floors, the mass of these was artificially reduced to a value of 10⁻⁵ kg/m³ to prevent their existence, and the missing corresponding load was applied to the adjacent masonry walls.

The mechanical properties of the masonry materials were assigned using the simplified Concrete Damaged Plasticity (CDP) hypothesis. For tensile and compressive behaviour, the constitutive relations adopted by Malcata et al. [36] were adapted and applied to both brick and stone masonry (Figure 3). It is important to point out that since these are written in terms of stress vs. strain, some mesh dependency may occur, and these were only calibrated for an average equivalent finite element size of around 20 cm. Although methods like the Turnsek–Cacovic method [54] may provide more conservative values, a separate tensile and compressive constitutive relation was adopted for the modelling due to the choice limitations in the ABAQUS materials library (

Table 2. Mechanical properties for the stone and brick masonry used in ABAQUS.

Material	ε’eloc	ε’cm	fc [MPa]	ft [Mpa]
Stone masonry	0.00125	0.003	1.67	0.025
Brick masonry	0.00135	0.0037	2.17	0.173

). This study examined the influence of variable d, which defines the damage variable for masonry materials. The damage was applied for one model with both compressive and tensile damage set to a maximum of 0.7.

2.2.4. Loads

Two types of load distribution were applied (see Figure 4): a uniform load proportional to the mass and a pseudo-triangular load, usually associated with equivalent horizontal static forces, where it is assumed that, during an earthquake, higher floors present higher acceleration.

2.2.5. Finite Element Modelling

Due to the case study’s apparent regularity, hexahedral elements (C3D8) were initially adopted. However, the analysis would not converge, and the building’s resistance would decrease suddenly at a precocious stage due to mesh dependency. Therefore, a 10-node modified tetrahedron with hourglass control (C3D10M) was used, as presented in Figure 2e. This choice was related to the structure’s having better mesh generation and fewer elements than linear solid brick elements like C3D8. Also, using C3D10M allowed a better convergence rate than C3D8, since the latter, in some tested cases, did not converge to a final solution. The maximum size of the mesh side was defined as the thickness of the thinnest walls of the building, corresponding to the 20 cm thick brick masonry partitional walls. In some areas, such as window thresholds, the number of elements was increased by reducing their size to enable linear approximations that would be unfeasible in sections defined by just one element. The mesh size was determined considering the need to balance the accuracy of the results and reasonable time-consuming analyses in such large models.

The finite element model used in this study did not include contact elements.

2.3. Comparing the EFM with the FEM

The nonlinear behaviour of the structure was characterised by its capacity curve in terms of base shear versus top displacement until the ultimate displacement.

The nonlinear analysis in 3Muri was fast and efficient. The displacement was defined as the weighted average of the top displacements proportional to the masses of the walls. In this software, a total of eight analyses were conducted: four uniform and four pseudo-triangular lateral force distributions. The analyses were carried out in the X (longitudinal) direction, along the main façade, and the Y (transverse) direction, along the side façades. The positive directions of the building in the plan were considered in the X direction from left to right and in the Y direction from bottom to top. In 3Muri, the ultimate displacement was defined after the occurrence of a mechanism leading to partial or complete collapse.

The nonlinear analysis in ABAQUS required considerably more time and computational effort than in 3Muri, with each analysis taking between 3 and 5 days to complete. Consequently, the decision was made to focus only on the conditioning cases, specifically those with lower safety values (−X Uniform and +Y Triangular). In ABAQUS, a total of four analyses were conducted, with two of them accounting for the damage variable (d). To define the ultimate displacement in this software, the Italian code NTC 2018 [50] provides two options: either a displacement corresponding to 80% of the maximum resistance of the structure or a displacement that creates an equivalent system with acceptable ductility.

2.3.1. Calibration with Modal Analysis

Taking into account the fact that the EFM and the FEM use different formulations to simulate masonry buildings, it was necessary to see if they presented the same vibration modes to compare the structural response later. The modal analysis results showed that the models were calibrated in the Y direction. However, in the X direction, even though the difference in frequency values was within an admissible range (17.7% for the FEM model), the vibration mode did not appear in the same order for both models. The FEM model presented lower stiffness in the Y direction, while the opposite occurred for the EFM model. The difference in results should be due to the fact that 3Muri, the EFM strategy adopted in this work, does not consider the out-of-plane stiffness of the walls. As previously mentioned, if the DMEM had been used, it would have been possible to account for both the in-plane and out-of-plane behaviour of the masonry walls, even though it is framed within the simplified approaches.

Thus, the values from ABAQUS tended to describe the actual behaviour of the structures more accurately in comparison with the values obtained by 3Muri. However, since the results of the nonlinear analyses were to be compared in this theoretical case study, it was decided to calibrate the ABAQUS model to match the modes obtained in 3Muri.

This calibration was accomplished by initially attributing the floor’s elasticity to rigid properties and then through various reductions in the masonry’s Young’s modulus in ABAQUS. The reduction applied was approximately 35% of the original stiffness. The results of the calibration are shown in

Table 3. Frequencies for the first mode after the calibration.

1° Mode	3Muri (Hz)	ABAQUS (Hz)	Error (%)
X	3.616	4.254	17.7%
Y	3.731	3.624	−3%

and Figure 5.

2.3.2. Capacity Curves

As can be observed in Figure 6, even with a reduction in the Young’s masonry modulus, the ABAQUS analyses exhibited higher stiffness values than the 3Muri analyses in the elastic phase, with a significant difference in the X direction. The obtained results suggest that the values of the masonry’s Young’s modulus in ABAQUS should have been further reduced, as recently suggested by Parisse (2024) [55]. These curves illustrate that 3Muri also presents more conservative values in terms of resistance for both directions.

Additionally, there was a notable variation in the capacity for deformation of the structure, with higher values in ABAQUS in the Y direction. This discrepancy occurred primarily because, in 3Muri, the damage was concentrated on a single floor F, whereas in ABAQUS the damage was distributed along the height of the walls in the Y direction, as is discussed in more detail in Section 2.3.4. The differing behaviour in terms of capacity of deformation observed in the X direction was influenced by the higher stiffness values in ABAQUS; this difference would nearly disappear if the Young’s modulus of the masonry walls were reduced in this software. In the ABAQUS curves, both with and without the d variable, only minor differences were observed due to some convergence errors in both cases. However, in the X direction, the case with damage showed a significant decrease, explained by the materials’ different capacities for redistributing loads.

2.3.3. Damage Distribution for Uniform Load in the −X Direction

In ABAQUS, the plastic deformations showed concentrated damage to the vertical elements of the first and second floors, consistent with the 3Muri observations in Figure 7. Only plastic strains (maximum principal) are plotted in the model without damage, but in the model with damage, both tensile and compressive damage fields are plotted. However, no significant damage was found on the top floor in the finite element software. Ground-floor damage includes the propagation of vertical cracks, affecting upper floors, and rocking in the corners. Mixed damage, combining bending and shear, was identified in the elements, like diagonal cracks in the façade elements and the central wall, as presented in Figure 8. The perpendicular walls exhibited an out-of-plane effect on the side walls, resembling damage from the 2009 L’Aquila earthquake [56].

2.3.4. Damage Distribution for Triangular Load in the +Y Direction

In 3Muri (Figure 9), the blank façades collapse at the first-floor mezzanines for the pseudo-triangular distribution. The ABAQUS plastic deformations (Figure 10a) revealed concentrated damage in the blank façades, indicating potential diagonal and vertical cracks. The interior walls exhibited loaded lintels and rocking at the base of some elements. In ABAQUS, there was more concentrated damage caused by diagonal cracks on the first floor of one blank façade, contrary to the expected behaviour observed in the other façade.

In the ABAQUS model using the damage variable, the predominant damage observed was attributed to tension (Figure 10), with localised compression damage (Figure 10b) in areas heavily loaded with tension. Analysis of elements with plastic deformations revealed concentrated damage in the spandrels and mixed damage in the vertical elements of the ground and first floors. The results between cases with and without damage were very close.

2.3.5. N2 Method Safety Verification

N2 method results are obtained with an equivalent SDOF system, and a bilinear capacity curve, that must be established for masonry structures. The evaluation of a transformation factor is necessary to define the equivalent SDOF system. Unlike 3Muri, ABAQUS does not compute the transformation factor automatically. Instead, it must be performed manually for each direction, which is a laborious process for intricate structures like the Monserrate Palace (Section 4).

In fact, while 3Muri automatically handles the transformation into an equivalent system with a single degree of freedom, in ABAQUS, each step is performed by the user, resulting in different coefficients of transformation ($\Gamma$) between the software systems. Applying the transformation factor revealed close approximations between the two numerical models. The displacements of each floor related to the vibration modes were taken into consideration when calculating the mass associated with each degree of freedom. As mentioned in Part 3 of Eurocode 8 [57] for heritage buildings, various limit states were examined; however, only the outcomes for near collapse, the most demanding limit state, are given.

The safety verification followed the proposal from Part 1 of Eurocode 8, beginning with the bilinearisation of capacity curves and then converting the results into an equivalent system with a single degree of freedom. For bilinearisation, Bondarabadi’s proposal [58] was adopted to account for the softening of the masonry.

Afterwards, comparisons were made between the target displacement (d^*_t) and the ultimate displacement (d^*_u) to assess the safety for both type 1 and type 2 earthquakes according to Eurocode 8 Part 1. In the X direction, both software packages confirmed safety, while in the Y direction, only the ABAQUS models with d considered verified safety. It is possible to observe in Figure 11 the verification coefficient for each software system using different load cases, in which the black line shows the point at which the ratio complies with safety. It is possible to observe that the safety ratios for both 3Muri and ABAQUS are similar when the same load case is compared.

3. Seismic Assessment of Monserrate Palace

The Palace of Monserrate, located in Monserrate Park in the district of Sintra, achieved UNESCO World Heritage status in 2000. Known for its architectural fusion of Gothic, Indian, and Moorish influences, the Palace boasts three distinct towers, two of which have a circular profile. This heritage building was tested with pushover analysis using the ABAQUS standard with physically and geometrically nonlinear analysis.

3.1. Numerical Modelling of the Tower

Due to the circular profile of the North Tower, the decision was made to use the FEM in the ABAQUS software instead of modelling through equivalent frames with the 3Muri software. This choice was supported by the requirement for an appropriate mesh definition, which must be imposed by the user for such irregular structures. Furthermore, the EFM was developed to represent regular structures; therefore, such a simplified method may not represent curved surfaces well. The model in ABAQUS was derived from the BIM model developed by Machete et al. (2023) [59] and was subsequently simplified to focus solely on the structural elements (Figure 12).

3.1.1. Adopted Geometry and Type of Analysis

To simplify the modelling of the structure, it was divided into three parts, as shown in Figure 13a: one for the entrance, a second for the tower itself, and a third to represent a part connected to the central body of the Palace. The primary challenge encountered in this model was the representation of openings in circular regions, which was addressed by employing temporary walls that were subsequently removed after the openings were created. The divisions and plane parts are presented in Figure 13b and Figure 13c, respectively. All these parts were merged into a final 3D model presented in Figure 13d–f.

A physically and geometrically nonlinear analysis was employed in the model, and force control was performed at each node where self-weight is a source of mass. Given that the high vertical displacements may raise the stress in the foundations and hence increase the degree of nonlinearity, a geometrically nonlinear analysis was required. A modification of the arc length approach known as the modified Riks algorithm [53] was employed in order to regulate the decrease in force in the final structural response, as it was anticipated that softening may have occurred in the static nonlinear analysis due to the brittle behaviour of the masonry structure.

3.1.2. Material Models and Adopted Mesh

The ABAQUS material software’s simplified Concrete Damaged Plasticity (CDP) hypothesis was used to assign the mechanical properties of the masonry materials in

Table 4. Adopted material properties for the Palace.

Material	ε′eloc	ε′cm	fc [MPa]	ft [MPa]	E [GPa]	G [GPa]	ρ [Ton/m3]
Stone masonry	0.00125	0.003	2.0	5.0	1.23	0.41	2.4
Masonry	0.00135	0.0037	7.0	10	2.85	0.95	2.2
Wood in floor	-	-	-	-	1.0	-	1.0 × 10−5
Wood in wood	-	-	-	-	1.3	-	1.0 × 10−5

and Figure 3. It was used in the constitutive relations adopted by Machete et al. (2023) [59] for compressive and tensile behaviour. It is crucial to remember that they were only calibrated for the adopted mesh and that there may have been some mesh dependency, as they were expressed in terms of stress vs. strain. To accommodate the complex geometry of the circular profile of the tower, the mesh elements assigned to this model were 10-node tetrahedral elements (C3D10M), as shown in Figure 13b,c. Since the thickness of the walls was greater than in the hypothetical building, the mesh average size was also greater, around 70 cm.

The finite element model employed in this study did not incorporate contact elements.

3.1.3. Applied Loads and Boundary Conditions

Again, two different distribution lateral load types were used: a uniform load, proportional to the mass and constant throughout the building’s height, which is usually a more demanding load, and a triangular load, typically related to similar horizontal static forces and in which upper floors exhibit greater acceleration during an earthquake. According to the in-situ report, it was acknowledged that there was a near-fixed connection and that the soil was sufficiently rigid to prevent deformation; therefore, fixed boundary conditions were chosen.

3.2. Dynamic Characterisation Tests and Material Calibration

Initially, a dynamic characterisation test of ambient vibration was undertaken to validate the material properties adopted for the Palace. Two high-sensitivity accelerometers (±4 g) and a central data acquisition system were used. While one accelerometer was placed in various locations on the structure to record accelerations, the other remained fixed. ARTeMIS software [60] was used to process the data, enabling the handling of accelerations by applying filters and identifying fundamental frequencies and vibration modes. The test covered the tower under study and the central tower due to its position as a focal point. However, for the purpose of this work, only the case-study tower was modelled in ARTeMIS for the dynamic characterisation.

By observing Figure 14, it is possible to conclude that the deformed shape for the first modes in the X and Y directions is very similar for both the experimental and numerical outputs. As illustrated in

Table 5. Experimental and numerical fundamental frequencies obtained with ARTeMIS and ABAQUS, respectively.

Mode	ARTeMIS (Hz)	ABAQUS (Hz)	Error (%)
1st X	4.941	6.067	22.6
1st Y	4.740	4.972	4.89

, the model had 22.55% and 5% errors in the X and Y directions, respectively, which were considered adequate to assume that the model was calibrated. It was then possible to admit that the initially attributed elasticity and mass density were indeed correct.

3.3. Nonlinear Static Analysis

A total of six analyses were conducted: four for the X direction (+X and −X, with two distributions each) and two for the Y direction (two distributions). This decision was based on the nearly perfect symmetry characterising the structure in the Y direction. Two options were considered to define the ultimate displacement: either the displacement where the structure reaches 80% of its maximum resistance or when the displacement creates an equivalent system with acceptable ductility.

3.3.1. Capacity Curves

As observed in Figure 15, the results obtained show that uniform distributions always lead to higher values of resistant capacity and lower deformation capacities. On the other hand, the triangular force distributions display less resistance than their uniform counterparts in every direction while exhibiting better deformation capacity. For the Y direction, only the positive direction was simulated since the structure presents symmetry.

3.3.2. Final Damage Distribution

The equivalent plastic strain (PEEQ), a scalar variable used to describe a material’s inelastic deformation, indicating where damage occurs, is presented in Figure 16 for the conditioning distribution in each direction. In the case of the Y direction, since the structure is symmetric, the equivalent plastic strain is displayed only in the positive direction. Part of the building adjacent to the tower was hidden in Figure 16a–e for better visualisation of the damage, while in f the plat band entrance was hidden.

Generally, the damage distribution revealed a localised concentration at floor-level 1. Diagonal cracks were observed in the piers on this floor, influenced by shear behaviour. The presence of lintels between the openings on the different floors is not evident due to the arches. Figure 16a,b, for the −X loading direction, indicate a significant concentration of bending damage (rocking) in the interaction between the tower and the section connected to the Palace’s central tower.

In Figure 16c,d, plastic deformations are visible for the +X loading direction, outlining the formation of diagonal cracks propagating towards the corners of the openings, indicating damage caused by shear. There is also a concentration of damage in the spandrels between the ground floor and the first level, at the level of the entrance plat band.

In Figure 16e,f, it is possible to see the plastic deformations that lead to the formation of diagonal cracks that propagate towards the corners of the openings, showing damage caused by shear. It is also possible to see in Figure 16f the concentration of damage caused by the connection to the entrance plat band and extensive diagonal cracking at the two ground-floor piers below the plat band. In all of these studies, it was possible to observe the plastic deformation and collapse mechanism of the arches.

4. Conclusions

The first part of the study compares two modelling approaches in modelling a hypothetical masonry structure: the Equivalent Frame Method (EFM) and the Finite Element Method (FEM). The main conclusions regarding the approaches are as follows:

▪

The EFM demonstrated simplicity and efficiency in the modelling process, with tools facilitating the import of plans from other software. 3Muri stands out for its ability to independently perform modal and nonlinear analyses;

▪

Only 3Muri’s model (EFM) underwent the safety verification process, taking only a few minutes to analyse all directions;

▪

ABAQUS (FEM) showcased its ability to model a complex geometry. However, even with imports from BIM models, it was less efficient than the other two software packages and took around a week to complete a nonlinear analysis;

▪

ABAQUS (FEM) faced challenges in defining loads on flexible floors or masses representing nonstructural elements, which 3Muri could easily handle;

▪

ABAQUS (FEM) requires the definition of many variables to describe masonry behaviour when subjected to compression and tension, which are usually unknown. In 3Muri (EFM), this problem does not arise, thanks to its effective handling of material inputs, as defined in the Italian norm [50];

▪

3Muri (EFM) highlights conditioning damage leading to the collapse of elements. At the same time, ABAQUS (FEM) can display various types of damage to its elements in more detail, but it does not directly indicate which drifts are exceeded.

It is also worth noting that the study was carried out on a hypothetical building that was not calibrated against experimental data. In the future, a simple, real case study, calibrated against experimental data, could be used for the comparison of other methods, including the Discrete Macro-Element Model (DMEM) approach. The DMEM, although framed within simplified approaches, is capable of accounting for both the in-plane and out-of-plane behaviour of masonry walls. It has recently been implemented in the HiSTrA commercial software, which can automatically perform the safety verification process.

In the context of the study on Monserrate Palace’s North Tower, its plastic deformations exhibited anticipated failure, primarily characterised by shear on the piers and spandrels of the arch-like openings. Specifically, plastic hinges developed at these arches’ base, top, and mid-span. Some bending damage was also observed thanks to the resistance of the entrance area and the connection to the middle body of the Palace. The option to model different parts in ABAQUS proved beneficial for segmenting the work. However, this software also revealed some modelling limitations, especially in creating openings in the circular profiled body, which were solved through the required nonconventional means. In addition, the pushover analysis in ABAQUS required approximately one week for each direction, which may not be feasible for practitioners. Therefore, methods that can model the irregular geometry of masonry structures in a practicable time are essential.

It is also worth noting that with the FEM, the adopted constitutive relations for the masonry used in the damage model are limited to the chosen size of the finite element; therefore, with different meshes, it is expected that the capacity curves will present different structural responses that may cause stronger softening, and, hence, some mesh dependency is expected. The mesh size was defined considering the requirement for accurate results and reasonably time-consuming analyses. Thus, an average mesh size equivalent to the smallest wall thickness was deemed sensible for studying the wall damage propagation and the building’s overall capacity curves.

Lastly, to reduce the mesh dependency, an inverse analysis is expected to be performed in the future to estimate the fracture energy to be used when modelling the global behaviour of a masonry wall. This will allow the analysis of future models without any mesh dependency.