Seismic Vulnerability Evaluation of a Historical Masonry Tower: Comparison between Different Approaches

Abstract

Throughout the last few decades, the scientific community has paid great attention to the structural safety of historical masonry constructions, which have high vulnerability with respect to seismic activities. Masonry towers are very widespread in Italy and represent an important part of the built heritage to be preserved. Different numerical methods with different levels of refinement were developed in the literature to evaluate their seismic performance. The present study shows a practical application of the seismic vulnerability evaluation of a masonry tower using different approaches. The aim is to provide practical suggestions to engineers for the successful evaluation of the performance of masonry towers under seismic loads. An in situ survey was performed to characterize the geometry of the structure and its constitutive material. All the collected information was introduced in a building information model, later used to generate different finite element models for the structural analyses. The global capacity of the structure was evaluated using three different models with different levels of complexity: the first simplified model is made of beam elements with cross-sections discretized in fibers; the second model is made of shell elements and uses a concrete damage plasticity model to describe the nonlinear masonry behavior; the third model adopts solid elements with a concrete smeared crack constitutive law. A preliminary eigen-frequency analysis is performed on the shell model to obtain some basic information about the structural behavior. Nonlinear static analyses were carried out for each model to understand the response of the tower under seismic loads, highlighting the main differences between the approaches. The behavior factor was evaluated on the basis of the analyses results and compared with the ones suggested by the Italian building code. The results showed that the towers do not satisfy the seismic demand required by the standards for all the considered models. Furthermore, the behavior factor calculated according to the Italian design code is overestimated, while the one evaluated by the simplified model is underestimated due to the neglection of the shear behavior. From all the analyzed configurations, the shell model resulted as a good compromise between reliable results and computation efficiency.

1. Introduction

Historical masonry buildings represent the significant part of the Italian historical heritage. Since these buildings were originally designed to resist gravity loads only, they often suffered serious damage during earthquakes [1,2,3,4]. The importance of protecting their historical and cultural value and the need to guarantee their structural safety to avoid loss of life during earthquakes has increased the interest of researchers and engineers in the study and assessment of their seismic behavior [5,6,7,8,9,10,11] in order to propose the most appropriate retrofitting interventions [12,13,14,15,16,17].

Despite the great advancements in research made in the last years and the introduction in Italian design codes of specific guidelines for built heritage [18,19,20], the evaluation of the seismic vulnerability of masonry buildings is still a challenging issue. The main difficulties arise from the heterogeneity of the material and its low tensile strength, when compared with its compressive strength, that results in a complex constitutive behavior.

Historical masonry towers are widespread in Italy and represent an important part of the built heritage to be preserved. The experience of past seismic events has proven how these structures are particularly vulnerable to earthquakes [21,22,23]. Indeed, their unique peculiar morphologic and typological features (e.g., high slenderness, presence of irregularities, leaning phenomena) might significantly affect their behavior under horizontal loads [24]. A realistic estimation of their structural behavior is fundamental to evaluate both the safety level and the most appropriate retrofitting interventions, but a compromise should always be found between required computational resources and accuracy of results.

Different approaches to predict the masonry towers’ seismic response were proposed in the scientific literature. The different methods can be grouped into two main categories: the continuum method and the discrete approach [25].

Continuum methods assume that stresses and strains over the structure are approximated by continuous functions. Among the continuum approaches, the finite element (FE) method proved to be a suitable tool to study the seismic behavior of towers [26,27,28,29,30], and the implementation of full nonlinear solid FE models is highly recommended to obtain accurate results. A macro-modeling strategy is usually adopted [31] together with sophisticated constitutive models (i.e., elasto-plastic with softening and damaging models) suitable for an insight into the nonlinear structural response. The main disadvantage of this approach is that it requires skilled users with strong mechanical backgrounds and powerful and expensive FE codes, which are different from the commercial ones commonly used in engineering practice. As alternative to 3D solid models, simplified approaches were formulated, as suggested in [18], in order to reduce the computational effort and to propose practice-oriented methods. The structure of the masonry tower, thanks to its simple geometry, can be represented by a continuous one-dimensional model (i.e., cantilever beam), and the main mechanical characteristics of the material in all the cross-sections along the height evaluated by means of a nonlinear elastic constitutive law formulated in terms of generalized stress and strain [32,33]. This procedure gives more reliable results for slender towers with a predominant bending behavior. When the tower is not particularly slender, flexural damage is often associated to shear cracks, which cannot be taken into account by the cantilever model, leading to wrong evaluations and inadequate retrofitting interventions.

In the discrete approach, the structure is considered as an assembly of distinct bodies, which interact along their boundaries [34]. Between the discrete approaches, the Rigid Body and Spring Model (RBSM) [35], which considers the masonry discretized as 2D rigid elements interconnected by nonlinear axial and shear springs, was adopted to perform the nonlinear dynamic analyses of towers [36] with a good approximation of the results and lower computational costs compared to nonlinear FE models.

More recently, the distinction between continuum and discrete methods has become somewhat blurred by their continuous evolution as they borrowed features from each other. This is the case of the combined finite-discrete element method [37], in which the structure is modelled with deformable solid finite elements with zero-thickness contact elements located at their boundaries to simulate the cracking phenomena. The application of this method to masonry, both in general and specifically for towers, can be found in [38,39].

The choice of method depends on several factors (i.e., complexity of the structure, available data, available computational resources, etc.) but, in general, the use of multiple approaches with different levels of complexity and the combination of the obtained results [40,41] may be useful to achieve a better and more exhaustive understanding of the tower’s seismic behavior.

In this paper, the seismic behavior of a historical masonry tower located in Northern Italy is analyzed utilizing different FE approaches with different levels of refinement and required computational resources. After a survey campaign aimed at characterizing the structure, the constitutive material, and the geometry of the tower, a 3D Building Information Model (BIM) of the building was implemented and used to define three different FE models: (i) a cantilever beam model, (ii) a shell model, and (iii) a solid elements model. Preliminary, an eigen-frequency analysis was performed on the shell model to obtain some basic information about the structural behavior. Then, the nonlinear static analyses of each model were performed, with the aim of highlighting the differences obtained with the different approaches. The main goal is to provide useful suggestions to practitioners that are called to perform the evaluations of the seismic vulnerability of masonry towers.

2. Description of the Tower

The masonry tower (Figure 1) is located in Travo (Northern Italy) at 465 m above mean sea level. It was built in the XI century and was part of an architectural complex also including a medieval castle, destroyed in 1255, and the nearby San Michele Arcangelo church.

The tower has a truncated pyramidal base on which four perimetral load-bearing walls develop with a total height of about 21.96 m. The walls are characterized by a variable thickness ranging from 3.15 m (base) to 1.20 m (top), with the following plan dimensions (Figure 2):

• floor −1 (base): 11.47 m × 12.45 m;
• ground floor: 10.35 m × 10.83 m;
• floors 1, 2, 3: 9.24 m × 9.72 m.

Floor −1 is made from a wooden slab, and the other floors are characterized by the presence of masonry vaults, oriented in an east–west direction, except for the vault of floor 2, oriented in a north–south direction.

A geometric and diagnostic survey campaign was carried out to characterize the structure. Ground penetrating radar investigation [42,43] and video-endoscopes surveys [44] were performed to determine the stratigraphy and the thickness of both masonry walls and vaults, while static penetration tests [45] were conducted to evaluate the quality of the mortar. In addition to on-site tests, compression tests [46] on the masonry blocks and mineralogical and petrographic analyses were performed. The walls are made of rubble stone masonry and vaults are characterized by rough-hewn masonry. The performed diagnostic tests did not give any useful information about the mechanical properties of the material to be used in the analyses, so their evaluation was carried out according to the Italian Building Standard (NTC18) [19,20], as reported in

Table 1. Mechanical properties of masonry.

	E	ν	G	γ	fm	τ0
	[N/mm2]	[-]	[N/mm2]	[kN/m3]	[MPa]	[MPa]
Walls	870	0.2	362.5	19	1.25	0.021
Vaults	1230	0.2	512.5	20	1.67	0.036

, where E is Young’s modulus, ν is Poisson’s modulus, G is the shear modulus, γ is the unit weight, f_m is the masonry compressive strength, and τ₀ is the masonry shear strength.

3. Numerical Analyses

3.1. Introduction

All information collected during the survey campaign was reported in a BIM model of the tower (Figure 3), realized with Revit Structure software [47], including the georeferencing of the tests and the related results. Starting from the BIM model and taking advantage of the interoperability between the BIM and the FE software, three different FE models with different levels of complexity were implemented and developed to evaluate the seismic behavior of the masonry tower: a fiber element model, a shell element model, and a solid element model (Figure 4). Isotropic materials were considered for all models because a clear assemblage of blocks and layers of mortar has not been identified for the studied masonry, so orthotropic behavior cannot be expected. The masonry is instead a rubble material made with irregular stones, embedded in thick lime mortar. For each model, nonlinear static analyses were carried out considering the appropriate constitutive laws.

3.2. Shell FE Model

The shell FE model (Figure 4a) was implemented in Midas GEN [48]. Before performing the nonlinear analyses, an eigenvalue analysis was conducted to evaluate the fundamental vibration modes of the tower. Figure 5 presents the first three vibration mode shapes. A rather standard dynamic behavior can be observed: the first two vibration modes involve a high value of participant mass and are mainly characterized by a bending behavior with negligible torsional effects, while the third vibration mode shows a predominant torsional behavior. The value of the participant mass of the fundamental vibration mode is affected by the presence of the basement, which is characterized by large wall thickness when compared with the upper walls.

The nonlinear behavior of masonry is simulated using the Concrete Damage Plasticity (CDP) model. The adopted compression and tensile stress–strain inelastic curves are illustrated in Figure 6a,b, while

Table 2. Values of the parameters adopted to define the CDP model.

Ψ	ε	fb0/fc0	K	µ
[°]	[-]	[-]	[-]	[-]
10	0.1	1.16	0.667	0.002

reports the adopted CDP parameters, evaluated according to [49,50], where ψ is the dilation angle, ε the eccentricity, f_b0/f_c0 is the ratio of initial equibiaxial compressive yield stress to initial uniaxial compressive yield stress, K is the ratio of the second stress invariant on the tensile meridian, and µ the viscosity parameter.

The seismic action is represented by means of an equivalent static loads system evaluated as:

$$F_i=F_h·z_i·\frac{W_i}{{{\displaystyle \sum }}_j{z}_j{W}_j}$$

(1)

$$F_h=S_d\left(T_1\right) ·W· \alpha /g$$

(2)

where F_i is the force applied to the i-th mass, W_i and W_j are the weight of the i-th mass and j-th mass, z_i and z_j are the height of the i-th mass and j-th mass measured from the tower base, S_d(T₁) is the value of the spectral acceleration from the design response spectrum for the first natural period T₁, W is the total weight of the tower, α is a coefficient equal to 0.85 if T₁ < 2T_c or equal to 1 in all the other cases, and g is the gravity acceleration. The load profile and the involved parameters are reported in

Table 3. Equivalent static loads distribution used to represent the seismic action.

h [m]	zi [m]	Wi [kN]	Wi·zi [kNm]	Fi [kN]
−2.24	0	0	0	0
0	2.24	3729	8354	99
3.62	5.86	3990	23,382	278
8.92	11.16	5394	60,199	716
13.89	16.06	4759	76,431	909
19.72	21.96	5479	120,326	1431

.

The response of the tower was evaluated in all the principal directions because of the non-symmetric distribution of the openings in the perimeter walls. Figure 7 shows, for each principal direction of the XY plane, the areas of the tower where the tensile stress exceeds the tensile strength, together with the respective value of the scale factor λ of the equivalent static lateral loads. The damaged areas are mainly located near the openings or on the second and third vaults.

3.3. Fiber Model

The fiber FE model of the tower was implemented in Midas GEN [48] using only beam elements (Figure 4b). Each level is characterized by an appropriate hollow squared cross-section, according to the drawings reported in Figure 2. The masses of the vaults, slab and roof were taken into account as nodal loads applied at different levels. Dead and live loads—evaluated according to the Italian Standard (NTC18) [19]—acting on the vaults, slab and roof were also considered as nodal loads. The openings were neglected. The fiber model was calibrated to be equivalent to the shell FE model in terms of total weight, fundamental natural periods, and vibration mode shapes. The site response spectrum was evaluated according to the NTC18 (Figure 8). The constitutive law proposed by Kent and Park [51] was chosen to describe the masonry compressive behavior with the following parameters: compressive strength f’_c = 1.4 MPa, strain at maximum stress ε_c0 = 0.004, confinement coefficient K = 1 (no lateral confinement considered), ultimate strain ε_cu = 0.014, and slope of the strain-softening branch Z = 533. The tensile strength of the material was ignored.

Two kind of load profiles were defined for the pushover analyses: one profile was proportional to the first (or second) vibration mode shape, and the other one was a uniform acceleration distribution. Figure 9 shows the obtained capacity curves in terms of base shear and top displacement for both X and Y directions.

The results show that the uniform acceleration profile leads to a greater value of collapse base shear with respect to the modal load case. In the X direction, the collapse base shear is 5247 kN for the modal distribution and 9532 kN for the uniform acceleration, with an increase of 82%, while in the Y direction, the collapse base shear is 5265 kN for the modal distribution and 9720 kN for the uniform acceleration with an increase of 85%. Therefore, as expected, the tower behavior is more similar to that of a cantilever than that of a frame structure, which usually collapses for soft story mechanisms associated to a uniform acceleration distribution.

Starting from the obtained results, it was possible to evaluate the values of the behavior factor q, according to the N2 method proposed by [52], by transforming the multi degrees of freedom (MDOF) capacity curves into the equivalent single degree of freedom (SDOF) capacity curves. Nonlinear force–displacement relations were simplified into elastic-plastic relations, whose main parameters are reported in

Table 4. Parameters obtained from the bilinearization of the capacity curves with modal distribution.

Direction	d*u	d*y	F*y	T*
[-]	[mm]	[mm]	[kN]	[s]
X	34.27	24.10	3340	0.59
Y	34.48	22.59	3472	0.57

in which T* is the period of the equivalent bilinear system, F*_y is the plastic limit force with respective displacement d*_y, and d*_u is the ultimate displacement of the bilinear system. Only the modal distribution was considered because, in the case of the uniform acceleration profile, the intersection between the capacity curve and the capacity spectrum occurs in the elastic branch for both the X and Y directions. The obtained values of the behavior factor for the modal distribution are q = 1.37 (X direction) and q = 1.39 (Y direction).

As alternative to the N2 method, the Capacity Spectrum Method (CSM) proposed in [53] was also used to evaluate the behavior factor q, adopting the procedure A of [54]. A comparison of the results is reported in

Table 5. Comparison of results obtained with N2 and CSM methods.

	Modal Distribution X Direction		Uniform Acceleration X Direction		Modal Distribution Y Direction		Uniform Acceleration Y Direction
	CSM	N2	CSM	N2	CSM	N2	CSM	N2
Sa [g]	0.305	0.276	0.411	0.412	0.307	0.280	0.411	0.412
Sd [mm]	35.14	33.30	17.23	18.77	35.78	31.80	16.49	18.01
V [kN]	4901	4495	6597	6700	5030	4645	6740	6836
d [mm]	52.23	44.82	25.61	25.26	52.55	42.54	24.22	24.09
q	1.08	1.37	1	1	1.07	1.39	1	1

, where S_a is the spectral acceleration evaluated in correspondence to the performance point, S_d is the related spectral displacement, V is the base shear of the inelastic structure evaluated considering the performance point, and d is the related maximum displacement of the inelastic structure. It has to be noted that the values of V and d obtained with the CSM method are related to the MDOF structure system, while those obtained by the N2 method consider a SDOF system. For this reason, to make a comparison possible, it was necessary to convert the results of the N2 method from SDOF to MDOF systems.

The values obtained with the two different methods are comparable. As an example, considering the value of the base shear obtained for the modal load distribution, the difference between the CSM and the N2 method is about 8.2% and 7.65% in the X and Y direction, respectively.

The CSM confirms the results predicted by the N2 method: the tower does not collapse due to a soft story mechanism because, under the uniform acceleration profile, the structure remains in the elastic range (q = 1). In fact, the shear capacity associated to the modal load distribution is always lower than the uniform acceleration one.

3.4. Solid FE Model

The last FE model of the tower was implemented in Midas FEA [55] using solid elements (Figure 4c). The nonlinear behavior of masonry was described using the Concrete Smeared Crack (CSC) approach (Figure 10), adopting the parameters reported in

Table 6. Parameters adopted for the definition of the CSC model.

	fC	GC	h	ft	${\mathbf{G}}_{\mathbf{f}}^{\mathbf{t}}$
	[MPa]	[N/mm]	[mm]	[MPa]	[N/mm]
Walls	1.25	0.1	440	0.15	0.1
Vaults	1.67	0.1	440	0.15	0.1

according to [56]. For the shear behavior, a linear model was adopted with a shear reduction factor β = 0.1. As in the shell model, an equivalent static lateral load distribution was taken into account for the seismic actions (Table 3).

Figure 11 shows the results of the nonlinear static analysis in terms of tensile stress and cracking paths at the beginning of the inelastic phase for the different loading directions. Cracks are mainly located around the openings and on the vaults, especially when considering the −X direction. The stresses configuration and the horizontal loads scale factors λ are comparable with the ones obtained from the shell model previously presented.

A further pushover analysis was performed to obtain the capacity curves of the structure, using a horizontal load profile distribution proportional to the first vibration mode in the X direction and to the second vibration mode in the Y direction (Figure 12). The uniform acceleration load profile was neglected since, from the previous analyses, the tower subjected to a uniform acceleration load profile turns out to be within the elastic range.

According to the N2 method, the behavior factor q was obtained by transforming the MDOF capacity curves into the equivalent SDOF capacity curves. Nonlinear force–displacement relations are then simplified into elastic–plastic relations, whose main parameters are reported in

Table 7. Parameters evaluated to apply the N2 method.

Direction	Γ	M*	Fu,MDOF	Fu,SDOF	du,MDOF	du,SDOF
[-]	[-]	[kN]	[kN]	[kN]	[m]	[m]
X	37.51	9508	5636	3882	0.126	0.087
Y	37.86	9798	5295	3690	0.099	0.069

, where M* is the mass of the equivalent SDOF system, Γ is the transformation factor, and F_u,MDOF and d_u,MDOF are the ultimate load and displacement of the MDOF system, while F_u,SDOF and d_u,SDOF are the ultimate load and displacement of the SDOF system. The bilinear capacity curves parameters are reported in

Table 8. Parameters obtained from the bilinearization of the capacity curves.

Direction	d*u	d*y	F*y
[-]	[m]	[m]	[kN]
X	0.071	0.021	3572
Y	0.069	0.018	3538

, where d*_u, d*_y, and F*_y are the ultimate displacement, the yielding displacement, and the yielding force of the idealized bilinear system, respectively. The curves are plotted in Figure 13.

For both the considered directions, the natural period of the equivalent SDOF system T* is higher than the value T_C = 0.545 s, which characterizes the elastic design response spectrum (Figure 8). Therefore, the demand in terms of displacement of the inelastic system is taken equal to the one of the elastic system with period T*.

Figure 14 highlights the collapse configuration of the masonry tower in both the main directions.

The obtained results are shown in

Table 9. Results obtained from the N2 method.

Direction	Se (T*)	d*e,max	q	de,max	PGAc	PGAd
[-]	[g]	[m]	[-]	[m]	[g]	[g]
X	0.0641	0.0035	1.67	0.0051	0.180	0.163
Y	0.0675	0.0033	1.83	0.0048	0.167	0.163

where S_e(T*) is the spectral acceleration evaluated at T*, d*_e,max is the displacement demand of the elastic system with period T*, q is the behavior factor, d_e,max is the displacement demand of the MDOF system, and PGA_c and PGA_d are the capacity and demand peak ground accelerations, respectively.

4. Discussion

To make an exhaustive comparison in terms of behavior factors, an evaluation of q was carried out according to the Italian design code NTC18 as:

$$q=q_0·K_R=2.975 ·0.8=2.38$$

(3)

where q₀ is the initial value of the behavior factor, which depends on the ductility class and the type of the structure, and K_R is a coefficient which considers the regularity of the structure. The obtained value is valid for both the X and Y directions.

Table 10. Comparison of the behavior factors evaluated with the different methods.

Approach	q
	X	Y
NTC18	2.38	2.38
Fiber FEM	1.37	1.39
Solid FEM	1.67	1.83

reports all values of q obtained with the different approaches, while

Table 11. Comparison of the horizontal load scale factors λ obtained with the shell and solid models.

Approach	λ
	+X	−X	+Y	−Y
Shell FEM	0.5	0.4	0.4	0.24
Solid FEM	0.6	0.4	0.4	0.3

shows a comparison in terms of collapse scale factor. Figure 15 shows a comparison in terms of the capacity curves for the modal load profile.

The FE model that better approximates the structure of the tower is clearly the solid model, which is associated to the highest computational cost and whose results will be considered as a reference in the comparison, since it represents the closest solution to reality.

The comparison in terms of the behavior factor shows that Italian Technical Standards (NTC18) tends to significantly overestimate the ductility of the tower. The values obtained with the fiber FE model are instead lower than the solid ones in both directions. This is due to the main assumptions underlying the fiber approach, which approximates the geometry (e.g., openings are not modelled) and, more important, neglects the shear behavior. This last hypothesis makes the fiber models more suitable for tower with high slenderness, where the flexural behavior is predominant.

In the solid FE model, a slight difference is observed for q in the X and Y directions due to the asymmetry of thicknesses and openings distribution in the perimetral walls. This difference is less evident in the case of the fiber model, since the approach does not allow the modelling of openings.

The scale factors of the equivalent static lateral loads λ obtained with the shell and solid models are comparable and, in both cases, a damage distributed near the openings or on the vaults was observed.

The comparison in terms of capacity curves shows a maximum value of base shear of about 5000 kN for both X and Y directions and in both the solid and fiber models. However, the solid model shows a higher ductility with respect to the fiber one, which is due to the main simplification hypotheses made by the fiber approach, as previously described.

5. Conclusions

In this paper, the seismic behavior of a historical masonry tower located in Northern Italy is discussed through the implementation of three different FE models, starting from a BIM model realized after the execution of a survey campaign aimed to characterize the structure. Given the impossibility to execute flat jack tests, the masonry mechanical parameters are evaluated according to the Italian Building Code.

The three implemented FE models are a fiber model, in which the tower is represented using beam elements and considering a cantilever behavior; a shell 3D model; and a solid 3D model.

Different analyses were carried out:

(1)

a pushover analysis of the fiber model, performed considering a modal distribution and a uniform acceleration profile, to evaluate an approximate behavior factor q using the N2 and the CSM methods;

(2)

a pushover analysis of the shell model with an equivalent static load distribution to evaluate the collapse scale factor of the horizontal forces;

(3)

a pushover analysis of the solid model with an equivalent static load distribution to evaluate the collapse scale factor of the horizontal forces.

The scale factors, λ, obtained from the different analyses are always less than 1. The tower is not able to resist the seismic forces defined by NTC18 and retrofitting interventions are required. The most critical areas are located on the first and second vaults and near the openings.

The results also show the importance of the correct evaluation of the behavior factor q using appropriate models, since, in this specific case, the evaluation of q made according to NTC18 results in an overestimation of the ductility of the structure, which can lead to inadequate and unsafe retrofitting interventions.

Finally, the use of a simplified cantilever model, even if it implies a significant reduction of the computational effort, it is not adequate for representing the structural behavior of a masonry tower characterized by a low slenderness, since it underestimates its real ductility. The adoption of a shell model can be seen instead as a good compromise between reliable results and computation efficiency.