Nonlinear Dynamic Analysis of a Masonry Arch Bridge Accounting for Damage Evolution

Abstract

This study investigates the nonlinear dynamic response of the masonry bridge ‘Ponte delle Torri’ in Spoleto, aiming at assessing the seismic performance of the structure and evaluating the occurring damaging mechanisms. A 3D Finite Element (FE) macromechanical procedure implemented in the FE program FEAP is adopted to model the bridge. To reproduce the typical nonlinear microcracking process evolving in masonry material when subjected to external loads, an isotropic damage model is used. This is based on a scalar damage variable introduced in the stress-strain constitutive law and equally degrading all the components of the elastic constitutive operator. A nonlocal integral definition of the damage associated variable, that is the equivalent strain measure governing its evolution, is adopted to overcome the mesh dependency problems of the FE solution typically occurring in the presence of strain softening behavior. Based on the results of a recent study by some of the authors, a single equivalent pier is analyzed, whose geometry and boundary conditions are selected so that its response can provide useful information on the out-of-plane dynamic behavior of the overall bridge. To perform the seismic assessment, a set of recorded accelerograms is properly selected to simulate the seismic history of the Spoleto site. The nonlinear dynamic response of the structure is evaluated and monitored in terms of top displacement time history, evolution of the global damage index, and distribution of the damage variable. First, a set of analyses is performed by imposing the selected ground motions one by one on the initial undamaged configuration for the structure with the aim of emphasizing the damaging effects on its dynamic response. Then, the accelerograms are arranged in sequence to reproduce the seismic history of the site and analyze the influence of accumulated damage on the dynamic amplification of the response. A critical comparison of the bridge response to the sequence of accelerograms and the single records is made, and the interaction between the damaged structure dynamic response and the signal characteristic is highlighted, as well.

1. Introduction

Masonry arch bridges are widespread in European roadway and railway networks. Most of them represent ancient constructions and constitute an invaluable part of the world architectural and historical heritage. Roman bridges are still popular in many European countries, as well as medieval bridges, distinguished by the different geometric shapes of the arch and other architectural details [1]. Due to the old building period, in many cases, masonry material constituting the arch bridges experienced severe deterioration processes, which increased the structure vulnerability to seismic events. Great efforts have been devoted in the last few decades to develop numerical procedures to assess the structural performances of these constructions and evaluate their load-carrying capacities and actual safety level [2,3,4], accounting for the change in the life service conditions and the material deterioration that have occurred since they were built.

Various methodologies were adopted in this context based on nonlinear static analyses (pushover) or nonlinear dynamic analyses, differently describing masonry constitutive response. Although pushover analyses are widely used at present to evaluate the structural global response and nonlinear collapse mechanisms, these were satisfactorily employed for the seismic assessment of framed structures [5], while being much less explored for masonry arch bridges [3]. The main issues include difficulty in defining an equivalent single degree of freedom system for the masonry bridge, that instead is a widely employed procedure for other classes of buildings [6], the selection of the most appropriate lateral load distribution, and the choice of the control point. For instance, in References [7,8], the seismic capacity of masonry bridges was evaluated by selecting both the top and center of mass as control nodes for the pushover analysis. Nonlinear dynamic analyses are certainly preferred when the aim is to accurately evaluate the seismic performances of the structures accounting for the interaction between the nonlinear damage and plastic mechanisms evolving during the earthquake and their dynamic response [9]. This method, although computationally expensive, has recently taken hold to study the complex dynamic behavior of masonry arch bridges, as this represents the most advanced deterministic analysis tool available. Usually, artificial or natural earthquake time histories are selected on the basis of the seismic characteristics of the specific site and its design response spectrum [3,7,10,11]. Moreover, first attempts to apply the Incremental Dynamic Analysis (IDA) to masonry bridges should be mentioned [7]. This methodology provides for the application to the structure of multiple ground motions, each scaled in order to obtain various levels of intensity, and then the creation of parametric curves, which relate the chosen damage and intensity measures. In this way, a comprehensive forecast of the structural response produced by the increasing ground motion severity is obtained.

A number of different Finite Element (FE) procedures were developed to evaluate masonry structural response, differing for the level of accuracy in modeling the geometric and mechanical properties of the real heterogeneous material and the related computational burden [12]. As concerns the micromechanical approaches, although these permit to accurately describe the masonry response modeling in detail the constituents geometry and mechanical response [13,14], they are computationally very cumbersome and are rarely adopted for real large structures. Focusing specifically on masonry bridges, macromechanical models are the most widespread [3,15,16]. These substitute the real heterogeneous masonry material with a fictitious homogenized medium, adopting a phenomenological constitutive law. The latter is formulated in order to be capable of approximately reproducing the main masonry nonlinear mechanisms related to the formation, propagation, and coalescence of brittle microcracks and to friction slidings occurring in mortar joints. Usually, damage-based [9,17,18] and plasticity-based [19] formulations are used. Multiscale procedures were also proposed representing a modern and attractive approach to face the mechanical modeling of microstructured materials [20,21,22], allowing to obtain a cost-effective tool when parallel computing sources are employed. Lastly, discrete models were successfully developed and used to evaluate masonry bridge structural response [23,24].

This paper focuses on the seismic assessment of ‘Ponte delle Torri’ in Spoleto, a Roman bridge whose current configuration results from several constructive phases over the centuries and dates back to the early 13th century. The dynamic response of ‘Ponte delle Torri’ has already been investigated in a recent work by some of the authors [10], considering a set of natural earthquake records properly selected in order to match the defined design spectra. The main aim of the present study is to analyze the structural dynamic response of the masonry bridge under a sequence of recorded accelerograms simulating the seismic history of Spoleto. This is an intriguing and not explored task that allows for highlighting of some interesting issues concerning the bridge structural behavior in presence of damaging mechanisms, arising from the critical comparison of the response to single excitations and the sequence of these. The main interest is to give insight into the most severe damage state attained by the structure during the sequence and investigate how the previously accumulated damage influences the maximum level reached. A 3D macromechanical model implemented in the FE code FEAP [25] is adopted. To describe the nonlinear damaging mechanisms occurring in masonry material due to the onset and propagation of microcracks, an isotropic damage constitutive law is adopted [13]. A scalar damage variable is introduced in the stress-strain relationship, equally affecting all the components of the elastic constitutive operator. This evolves according to an exponential law as a function of an equivalent strain measure accounting for the principal strains, and satisfies thermodynamics irreversibility constraint. A nonlocal integral definition is introduced for the latter to overcome the pathological mesh dependency of the FE solution due to the strain-softening nature of the constitutive law. The described numerical procedure is used to investigate the nonlinear dynamic response of a single pier, whose geometry and boundary conditions are selected so that its response can provide useful information on the out-of-plane dynamic behavior of the overall bridge [10]. This is subjected to the single selected records and to the sequence composed of these to evaluate the effect of the accumulated damage on the overall dynamic response. The definition of a reduced simplified model to investigate the structural response of the bridge takes its inspiration from the procedures widely adopted for the seismic vulnerability evaluation of other classes of buildings and masonry structures. Among others, this is a well-established procedure to evaluate the response of unreinforced masonry walls subjected to out-of-plane seismic excitation aiming to assess their structural stability [6,26]. In this context, commonly a single-degree-of-freedom idealization of the masonry walls subjected to out-of-plane action is made. By the way, the study here presented, which is focused on an equivalent single pier, allows for obtaining information only on the out-of-plane response of the overall bridge and is a suitable choice if the attention is mainly focused on this behavior. As shown in Reference [10], a very good match is obtained when the response of the equivalent pier is compared to that of the bridge, both in the linear elastic and nonlinear field, depending on the frequency content of the applied seismic signal. When the first out-of-plane mode of the bridge is mostly excited, the analysis performed on the single pier gives suitable information in terms of out-of-plane displacement history, acceleration, and damaging paths. At the same time, a considerable computational saving is obtained. If the interest is moved to other aspects of the bridge structural response, different reduced or more sophisticated models are needed. As for example, if other modal shapes and frequencies are investigated, as those involving torsion of the pier, a different reduced model should be conceived. Eventually, for seismic excitations characterized by a rich and large frequency content, simplified models are hard to be identified and adopted, and modeling of the overall bridge has to be resorted to.

Section 2 and Section 3 contain a brief description of the ‘Ponte delle Torri’ geometry, its historical background, and the seismic history of the site. The adopted damage constitutive law and the FE procedure are illustrated in Section 4, while Section 5 contains the results of the nonlinear dynamic analyses performed on the equivalent single pier in terms of displacement time history and damage index evolution. Finally, in Section 6, some concluding remarks are given.

2. Structural Description and Historical Background

‘Ponte delle Torri’ is one of the most prominent monuments of the Municipality of Spoleto, connecting Sant’Elia hill, where the city lies, with Monteluco (Figure 1).

It is a masonry bridge consisting of the abutments and nine hollow piers (‘Torri’) founded on emerging rock and linked at their top by pointed arches. The geometric characteristics of the bridge are shown in Figure 2 [27].

The present configuration of the bridge, dating back to the early 13th century, is the result of several constructive phases over the centuries, as witnessed by the different geometric proportions of piers and arches [28,29,30]. Romans built a smaller bridge to carry the water from the springs of Monteluco to the thermae of Spoleto. The base and the lower part of the piers on Monteluco side, in stonework with units coarsely dressed, is what remains of the Roman bridge. The upper part of these piers, similar to most piers on the side of Sant’Elia hill, are in stonework with dressed units and sharp corners [31]. It can be presumed that the arches at intermediate height between piers 2, 3, and 4, as well as the flying buttress between piers 4 and 5, were built during the same constructive phase. In the subsequent centuries, the bridge underwent many interventions of strengthening and restoration [30]. At present, the bridge shows a diffuse state of cracking, localized cracks in some arches, and a sub-vertical crack in the West elevation of pier 4 [31].

3. Seismic History of the Site

Italy is one of the countries with the richest macroseismic observations catalogues. Using the most up-to-date [32], starting from 1000 AD, for the town of Spoleto, it is possible to find a list of 36 observations whose locally felt intensity is above the damage threshold ($I_s\ge$ V Mercalli-Cancani-Sieberg, MCS). The maximum intensities at the site are $I_s=$ IX-X, recorded in 1328, and $I_s=$ IX, recorded in 1730. A set of 20 earthquakes was selected as representative of the seismic history of the site and reported in

Table 1. List of earthquakes representative of the seismic history of the site.

Date	Time UTC	Earthquake	${\mathit{M}}_{\mathit{w}}$	d [km]
1246		Spoleto	$5.33$	$0.7$
1277		Spoleto	$5.56$	$0.7$
30 April 1279	18:	Appennino umbro-marchigiano	$6.20$	$41.4$
1 December 1298		Monti Reatini	$6.26$	$21.8$
1 December 1328		Valnerina	$6.49$	$26.3$
6 November 1599	01:25	Valnerina	$6.07$	$22.7$
28 April 1646	07:	Monti della Laga	$5.90$	$43.2$
14 January 1703	18:	Valnerina	$6.92$	$26.9$
27 June 1719	06:30	Valnerina	$5.59$	$29.4$
12 May 1730	05:	Valnerina	$6.04$	$30.8$
5 June 1767	01:30	Valle Umbra	$5.45$	$9.8$
3 September 1815	23:	Valnerina	$5.58$	$24.7$
13 January 1832	13:	Valle Umbra	$6.43$	$29.7$
22 August 1859		Valnerina	$5.73$	$30.6$
23 February 1879	18:30	Valnerina	$5.59$	$24.7$
12 September 1898	14:14:10	Valnerina	$5.48$	$27.1$
31 December 1948	03:32	Monti Reatini	$5.42$	$26.7$
26 September 1997	09:40:02	Appennino umbro-marchigiano	$5.97$	$32.5$
24 August 2016	01:36:32	Monti della Laga	$6.18$	$39.6$
30 October 2016	06:40:18	Valnerina	$6.61$	$31.8$

($M_w$ = moment magnitude, d = epicentral distance) [33]. The expected maximum horizontal ground acceleration $a_g$ specified by the Italian Building Code [34] for the site and type A soil is reported in Figure 3 as function of the return period $T_R$. In particular, $a_g=0.217g$ for $T_R=475$ years, corresponding to medium seismicity.

The seismic history was simulated by means of recorded accelerograms, selected from the Engineering Strong Motion Database [35], which have occurred in Central and Southern Italy since 1980, with magnitude and epicentral distance close to those of the historical earthquakes, normal style of faulting, which is typical of the earthquakes of Central Italy Apennines, and recorded on soil type A. The list of the recorded accelerograms reproducing the seismic history is reported in

Table 2. List of recorded accelerograms used to reproduce the seismic history.

Date	Time UTC	Earthquake	${\mathit{M}}_{\mathit{w}}$	Station	d [km]	Component	$\mathit{PGA}$ [cm/s${}^2$]	$\mathit{PGV}$ [cm/s]	${\mathit{I}}_{\mathit{H}}$ [cm]
12 October 1997	11:08:35	Valnerina	$5.19$	FORC	4.1	EW	82.37	3.81	7.27
7 April 2009	17:47:37	Aquilano	$5.54$	RM03	3.4	EW	90.93	4.65	12.51
24 August 2016	01:36:32	Monti della Laga	$6.18$	SPM	39.6	EW	65.98	2.38	8.93
26 October 2016	19:18:06	Valnerina	$6.07$	T1212	18.8	NS	86.64	6.37	20.13
30 October 2016	06:40:18	Valnerina	$6.61$	MNF	26.0	EW	130.74	6.43	19.35
26 October 2016	19:18:06	Valnerina	$6.07$	ACC	25.4	EW	88.36	4.35	14.61
24 August 2016	01:36:32	Monti della Laga	$6.18$	SPM	39.6	EW	65.98	2.38	8.93
23 November 1980	18:34:53	Irpinia-Basilicata	$6.81$	ALT	23.4	EW	56.34	6.28	22.05
14 October 1997	15:23:09	Valnerina	$5.62$	FHC	29.5	NS	73.68	2.47	7.41
26 October 2016	19:18:06	Valnerina	$6.07$	MZ19	30.4	EW	93.91	3.86	11.58
26 October 2016	17:10:36	Valnerina	$5.51$	CLO	8.4	NS	125.23	4.79	16.42
14 October 1997	15:23:09	Valnerina	$5.62$	FHC	29.5	NS	73.68	2.47	7.41
30 October 2016	06:40:18	Valnerina	$6.61$	SPM (2)	31.8	EW	92.71	6.80	16.70
14 October 1997	15:23:09	Valnerina	$5.62$	FHC	29.5	NS	73.68	2.47	7.41
14 October 1997	15:23:09	Valnerina	$5.62$	FHC	29.5	NS	73.68	2.47	7.41
14 October 1997	15:23:09	Valnerina	$5.62$	FHC	29.5	NS	73.68	2.47	7.41
26 October 2016	17:10:36	Valnerina	$5.51$	MMO	16.8	NS	71.24	2.48	7.51
26 October 2016	19:18:06	Valnerina	$6.07$	MZ19	30.4	EW	93.91	3.86	11.58
24 August 2016	01:36:32	Monti della Laga	$6.18$	SPM	39.6	EW	65.98	2.38	8.93
30 October 2016	06:40:18	Valnerina	$6.61$	SPM (2)	31.8	EW	92.71	6.80	16.70

, along with their component (NS or EW), peak ground acceleration $PGA$, peak ground velocity $PGV$, and Housner intensity $I_H$. For each record, the component associated with the maximum $PGA$ was considered. A total of 12 different records was used to reproduce the 20 historical events, with some records reproducing more than one historical event with similar characteristics. The accelerograms recorded on top of Monteluco were used for the two last events that occurred in 2016. The magnitude of the earthquakes generating the records was compared with that of the historical events (Figure 4a), resulting in very high correlation ($R=$ correlation coefficient). In a similar manner, the epicentral distance of the station was compared with the distance between the site and the epicenter of the historical event (Figure 4b), again, resulting in very good correlation.

The spectral compatibility of the selected accelerograms was investigated by comparing the average response spectrum of the 12 selected records with the spectrum provided by the Italian Building Code [34] for the site. The 20 seismic events range from since 1246 up to 2016, for a return period equal to 40.5 years. In Figure 5, the acceleration spectra of the 12 records, together with their average spectrum, are compared with the acceleration response spectrum of the Italian Building Code for a return period of 35 years and soil type A. In the range between $T=1.6\mathrm{s}$ and $T=0.6\mathrm{s}$ (significant for the analyzed structure), the shape of the average spectrum is similar to that of the Code spectrum. Although the ordinates of the average spectrum are 20–30% lower, the set of the 12 selected records was considered sufficiently accurate to reproduce the seismic history of the site on the basis of the correlations on magnitude and epicentral distance.

4. Macromechanical Damage Model

Masonry material response is strongly dominated by the onset and evolution of brittle damaging mechanisms which cause strain-softening response and the formation of localized damage regions in masonry structures. Among the different approaches proposed in literature, here, a macromechanical model is adopted considering the real heterogeneous material as an equivalent homogenized medium [9,17,18]. The phenomenological isotropic constitutive law, originally proposed in References [13,36], is used. This is based on the introduction of a scalar damage variable in the stress-strain law relating vectors $\sigma ={[{\sigma }_x,{\sigma }_y,{\sigma }_z,{\tau }_{xy},{\tau }_{yz},{\tau }_{zx}]}^T$ and $\epsilon ={[{\epsilon }_x,{\epsilon }_y,{\epsilon }_z,{\gamma }_{xy},{\gamma }_{yz},{\gamma }_{zx}]}^T$, as follows:

$$\begin{array}{cc}\hfill \sigma & = \left(1-D\right)\mathbf{E}\epsilon ,\hfill \end{array}$$

(1)

with $\mathbf{E}$ denoting the $6\times 6$ isotropic elastic constitutive matrix of the virgin material depending on the material elastic parameters, i.e., Young’s modulus E and Poisson ratio $\nu$. This can be also written as:

$$\begin{array}{cc}\hfill \sigma & =\tilde{\mathbf{E}}\epsilon ,\hfill \end{array}$$

(2)

where, relying on the principle of the equivalent strain [37], the effective constitutive matrix $\tilde{\mathbf{E}}=(1-D)\mathbf{E}$ is introduced. The damage variable D, according to its physical definition as the ratio between the damaged and total area of the material representative volume element, is bounded in the range $[0,1]$, and it satisfies the thermodynamics irreversibility condition $\dot{D}\ge 0$. The evolution process of the damage variable during the loading history is governed by the following law, accounting for the above physical and thermodynamics constraints:

$$D=\underset{history}{max}\left\{0,min\left\{\tilde{D},1\right\}\right\},$$

(3)

with the variable $\tilde{D}$ defined as:

$$\tilde{D}=1+\frac{1}{{\overline{\epsilon }}_{eq}{\left({\epsilon }_t-{\epsilon }_u\right)}^3}{e}^{-\beta \left({\overline{\epsilon }}_{eq}-{\epsilon }_t\right)}{\left({\overline{\epsilon }}_{eq}-{\epsilon }_u\right)}^2\left({\overline{\epsilon }}_{eq}{\epsilon }_u+{\overline{\epsilon }}_{eq}{\epsilon }_t-2{\epsilon }_t^2\right).$$

(4)

In Equation (4), ${\epsilon }_t$ denotes the tensile strain threshold, ${\epsilon }_u$ the ultimate strain value corresponding to full damage, and $\beta$ is a parameter governing the shape and slope of the softening branch in the constitutive law, defined as:

$$\beta ={\beta }_t+\frac{1}{1+exp\left(\alpha \left(I_1-I_{1m}\right)\right)}\left({\beta }_c-{\beta }_t\right),$$

(5)

where ${\beta }_c$ and ${\beta }_t$ are the values of $\beta$ for a mostly contracted or elongated strain state, respectively, $\alpha$ governs the rate of the variation of $\beta$ from ${\beta }_c$ to ${\beta }_t$, and vice-versa, $I_1$ is the strain first invariant, and $I_{1m}$ its value corresponding to $\beta =\left({\beta }_c+{\beta }_t\right)/2$. According to Equation (4), a scalar strain measure controls the damage evolution. This is defined on the basis of the principal strain values ${\epsilon }_i$, as:

$${\epsilon }_{eq}=\sqrt{{〈\sum _{i=1}^3{〈{\epsilon }_i+{\epsilon }_o〉}_{+}^2-\kappa \sum _{i=1}^3\sum _{j=1}^3\frac{\left(1-{\delta }_{ij}\right)}{2}{〈{\epsilon }_i〉}_{-}{〈{\epsilon }_j〉}_{-}〉}_{+}}-{\epsilon }_o.$$

(6)

The Macaulay’s brackets ${〈•〉}_{+/-}$ evaluate the positive/negative part of the strains, and ${\delta }_{ij}$ denotes the Kronecker’s symbol. To be noted is that the second term under the square root, depending on the negative part of the principal strains, allows for reduction of the equivalent tensile strain due to the presence of compressive strains. The effect of this reduction is ruled by the parameter $\kappa$, while the coefficient ${\epsilon }_o$ regularizes the shape of the limit domain.

The damage constitutive law described was implemented in a displacement-based FE procedure in the FEAP program [25]. An 8-node hexahedral solid element was developed, where the three unknown displacement fields are interpolated by tri-linear shape functions on the basis of 24 kinematic degrees of freedom, 3 at each node. The Gauss integration technique is used to perform the integration of the element stiffness and mass matrices, as well as the nodal residual vector, by adopting a $2\times 2\times 2$ quadrature rule. To avoid mesh localization problems related to the presence of strain-softening in the constitutive law evaluated at each Gauss integration point, the nonlocal integral procedure is adopted [38,39], and the nonlocal definition of the strain measure in Equation (6) is introduced as:

$${\overline{\epsilon }}_{eq}\left(\mathbf{x}\right)=\frac{1}{{\int }_{V_r}\psi \left(\mathbf{x}-\mathbf{y}\right)dV}{\int }_{V_r}\psi \left(\mathbf{x}-\mathbf{y}\right){\epsilon }_{eq}\left(\mathbf{y}\right)dV,$$

(7)

where $V_r$ is the reference volume. The weighting function $\psi$ ruling the influence of the generic point $\mathbf{y}$, located in the neighborhood $V_r$ on the analyzed point $\mathbf{x}$, is given by:

$$\psi \left(\mathbf{x}-\mathbf{y}\right)=e^{-{\left(\frac{\parallel \mathbf{x}-\mathbf{y}\parallel }{l_c}\right)}^2},$$

(8)

where $l_c$ is the material characteristic length, usually set on the basis of the microstructural size, that is the block size in masonry material.

Finally, a global damage variable is introduced to monitor the damage state at structural level. Adopting the definition in Reference [10], the Global Damage Index ($GDI$) is evaluated on the basis of the local values D, given by Equation (3), at each integration point of the FE model, as follows:

$$GDI=\frac{1}{{\int }_{V}dV}{\int }_{V}DdV,$$

(9)

only accounting in the numerator of the integral for the region where $D\ge 0.1$, with V being the structural volume. To be noted is that the introduced global index $GDI$ does not provide information on the spatial damage distribution but is a scalar measure of the damage severity able to follow the progression of the degrading process in the structure. Other definitions could be adopted as an alternative to Equation (9), which provides an average measure of the overall damage and does not detect the maximum levels reached in localized zones. An evaluation based on the 95% fractile seems to be an effective option because it is closer to maximum local damage. Its evaluation is schematically depicted in Figure 6. Here, the values of damage $D_k$ ($0\le D_k\le 1$) evaluated at all integration points of the FE mesh are reported in abscissa and related to the corresponding ratio between the accumulated volume, where that damage level is reached, and the total structural volume V. Hence, the global damage index $D_{0.95}$ is determined as the value that is not exceeded in the 95% of the entire volume.

5. Masonry Bridge Dynamic Response

In this section, the dynamic response of ‘Ponte delle Torri’ is studied, adopting the FE procedure described above. First, the structural modal characteristics investigated in Reference [10] are briefly recalled. Then, the response to the records listed in Table 2 is analyzed, following two different ways. As a first step, a set of analyses is performed by applying, one by one, the selected natural earthquakes on the initial undamaged configuration of the structure. After that, the structural behavior is studied arranging in sequence all the input motions with the aim of evaluating the effect of the seismic history on the most relevant aspects of the response. The Newmark implicit time integration scheme, with coefficients $\gamma =1/2$ and $\beta =1/4$, is used, together with the Newton-Raphson iterative algorithm, to solve the dynamic equilibrium equations governing the nonlinear FE procedure. The mass and damping matrices are defined according to the consistent formulation and the Rayleigh method, respectively. A damping factor equal to 5% is set on the basis of the two first vibration modes.

5.1. Modal Analysis

Figure 7a–d show the first four modal shapes of the bridge, numerically evaluated with the 3D FE model described in detail in Reference [10], assuming the following material parameters: Young’s modulus $E=8000\mathrm{MPa}$, Poisson ratio $\nu =0.2$, and mass density $\rho =2067\mathrm{kg}/{\mathrm{m}}^3$. These parameters are set on the basis of the results reported in Reference [31], where an experimental investigation was conducted to identify the masonry mechanical characteristics. The frequencies $f_i$, corresponding to the modal shapes in Figure 7, are reported in

Table 3. ‘Ponte delle Torri’: numerical and experimental first four natural frequencies $f_i$.

	${\mathit{f}}_1$ [Hz]	${\mathit{f}}_2$ [Hz]	${\mathit{f}}_3$ [Hz]	${\mathit{f}}_4$ [Hz]
numerical	0.605	1.039	1.489	1.819
experimental	0.634	1.021	1.508	2.003

and compared with the experimental values derived by ambient test vibration [29,40], showing a very good match. The results indicate that the first modes involve the out-of-plane bridge behavior. Hence, focusing on the out-of-plane response and relying on the numerical evidences reported in Reference [10], the modal characteristics of an equivalent single pier are also analyzed. It was proved that the analysis of the dynamic response of the equivalent pier can provide useful information on the bridge overall behavior with significant computational saving. Obviously, the match between bridge and pier responses depends on the frequency content of the applied records. Indeed, setting proper boundary conditions (fully restrained base and symmetry conditions on the left and right vertical transverse planes) and geometry, the first pier natural frequency, $0.603\mathrm{Hz}$, is very close to that of the overall bridge, $0.605\mathrm{Hz}$. Figure 8a shows a schematic of the analyzed pier, whose height is the average height of all the piers, whereas Figure 8b,c contain its first two modal shapes. It is interesting to note that the first four modes of the complete model are mainly governed by the horizontal bending of the deck: the first mode is approximately symmetric with respect to the centerline of the bridge, with the piers following a “first mode” deformation along the height; the second mode is approximately skew-symmetric with respect to the centerline, with the piers that still follow a first mode deformation; the third mode is approximately symmetric with respect to the centerline, with two central piers following a second mode deformation in phase with each other and the remaining ones deformed according to their first mode; and the fourth mode is approximately skew-symmetric with respect to the centerline with two central piers that follow a second mode deformation in phase opposition with each other, while the remaining ones are deformed according to their first mode.

5.2. Response to Natural Earthquakes Simulating the Seismic History

5.2.1. Individual Events

The dynamic response of ‘Ponte delle Torri’ to the accelerograms listed in Table 2 is simulated using the equivalent pier model introduced in the previous section. The masonry compressive strength, $f_c$, is set equal to $6\mathrm{MPa}$ according to the experimental results collected in Reference [31], whereas the tensile strength, $f_t$, is assumed equal to $0.35\mathrm{MPa}$. These values lead to a typical ratio $f_t/f_c$ for masonry material, usually ranging between $\frac{1}{20}$ and $\frac{1}{10}$. The characteristic length $l_c$ is set equal to $4\mathrm{m}$. The material mechanical parameters are reported in

Table 4. Material parameters adopted for the numerical model.

E [MPa]	$\mathit{\nu }$	${\mathit{\epsilon }}_{\mathit{t}}$	$\mathit{\kappa }$	${\mathit{\epsilon }}_{\mathit{u}}$	${\mathit{\epsilon }}_0$	${\mathit{\beta }}_{\mathit{t}}$	${\mathit{\beta }}_{\mathit{c}}$	$\mathit{\alpha }$
8000	$0.2$	$0.25\times {10}^{-4}$	$0.03$	$15\times {10}^{-4}$	$1\times {10}^{-5}$	14,000	200	8000

. Two-stage analyses are performed: first, the structure is subjected to self-weight, and then the ground acceleration is applied along the bridge out-of-plane direction. The structural response is monitored in terms of out-of-plane displacement of point A, $v_A$, located at the top side of the pier, as indicated in Figure 8a. To be noted is that the displacement induced by the application of self-weight is not represented in the time histories of the displacement $v_A$ reported in the following. Figure 9a, Figure 10a, Figure 11a, Figure 12a, Figure 13a, Figure 14a, Figure 15a, Figure 16a, Figure 17a, Figure 18a, Figure 19a and Figure 20a show the significant part of the response obtained by applying each accelerogram, starting from the pier undamaged configuration, with the aim of highlighting the effect of the onset of degrading mechanisms and better understanding the influence of accumulated damage during the seismic sequence. To this end, the responses evaluated assuming a linear elastic constitutive behavior for the material (blue lines in Figure 9a, Figure 10a, Figure 11a, Figure 12a, Figure 13a, Figure 14a, Figure 15a, Figure 16a, Figure 17a, Figure 18a, Figure 19a and Figure 20a) are compared with those obtained for the damaged structure considering the damage model described in Section 4 (red lines in Figure 9a, Figure 10a, Figure 11a, Figure 12a, Figure 13a, Figure 14a, Figure 15a, Figure 16a, Figure 17a, Figure 18a, Figure 19a and Figure 20a).

Focusing on the elastic responses, various levels of dynamic amplification emerge (note that different ranges for y-axis are used in Figure 9a, Figure 10a, Figure 11a, Figure 12a, Figure 13a, Figure 14a, Figure 15a, Figure 16a, Figure 17a, Figure 18a, Figure 19a and Figure 20a). These can be interpreted on the basis of both the effects of $PGA$ (see Table 2) and displacement response spectrum (see Figure 9c, Figure 10c, Figure 11c, Figure 12c, Figure 13c, Figure 14c, Figure 15c, Figure 16c, Figure 17c, Figure 18c, Figure 19c and Figure 20c) of each earthquake. For instance, although the ground acceleration recorded at ALT station is characterized by low $PGA$ intensity ($56.34\mathrm{cm}/{\mathrm{s}}^2$), this induces the largest response amplification, as the displacement spectrum shows high values in correspondence of the pier first natural period (indicated with dashed vertical line in Figure 15c).

The lowest displacement amplifications are detected in cases of the FORC, FHC, and MMO records. In fact, although characterized by high values of the $PGA$ ($82.37\mathrm{cm}/{\mathrm{s}}^2$, $73.68\mathrm{cm}/{\mathrm{s}}^2$, and $71.24\mathrm{cm}/{\mathrm{s}}^2$, respectively), their response spectra show low values (lower than $1\mathrm{cm}$) in correspondence of the structural periods (see Figure 9c, Figure 16c, Figure 20c). In general, it emerges that earthquakes with comparable response spectrum ordinates (consider also RM03-SPM and T1212-MNF-ACC-SPM${}_2$) lead to similar maximum displacements, regardless of their $PGA$ values. This confirms that the ratio between the forcing and natural frequencies plays a significant role on the dynamic response amplification, and the overall behavior is influenced by several factors highly correlated with each other. Even more complex is the interpretation of the structural response when the nonlinear phenomena are accounted for: displacements can be reduced or amplified with respect to the elastic case, depending on the combination of various phenomena [10,41]. Indeed, the onset and evolution of damage modify the structure mechanical properties, making it more flexible; thus, the structural natural periods increase, as a consequence. Correspondingly, the ordinates of the displacement spectrum may either increase or decrease, depending on the frequency content of the accelerogram.

Schematic representation of the degrading mechanisms evolution during each record is shown in Figure 9b, Figure 10b, Figure 11b, Figure 12b, Figure 13b, Figure 14b, Figure 15b, Figure 16b, Figure 17b, Figure 18b, Figure 19b and Figure 20b by the time histories of the $GDI$ defined in Equation (9), whose initial value represents the damage state induced in the structure by self-weight.

The comparison of the maximum attained displacements in the elastic and nonlinear cases can, again, be better understood referring to the displacement spectra in Figure 9c, Figure 10c, Figure 11c, Figure 12c, Figure 13c, Figure 14c, Figure 15c, Figure 16c, Figure 17c, Figure 18c, Figure 19c and Figure 20c. Let us focus on the responses to some representative earthquakes taken as example. Under the ALT accelerogram, the pier response significantly involves the first mode of vibration, whose period ($T_1=1.658\mathrm{s}$) increases at the onset of damage, likely moving the structure towards the decreasing branch of the displacement spectrum. This causes the reduction of the dynamic amplification with respect to the elastic case. Conversely, increment of the maximum response emerges in case of T1212 record when damage is activated. As a matter of fact, the corresponding displacement response spectrum (Figure 12c) shows increasing branches near the main periods of the structure. To better clarify, in Figure 21a,b, the variation of the first two periods of the damaging structure ($T_d^i$), normalized with respect to the elastic ones ($T_e^i$), are plotted for ALT and T1212 earthquakes, respectively. These are derived from a step-by-step modal eigenvalues analysis, performed on the basis of the damaged stiffness matrix. Progressive increments of the natural periods emerge in the first part of the loading histories, and then these stabilize at constant values when the evolution of damage stopped. The final values of the damaged periods are also indicated with dashed blue vertical lines in the most significant part of the displacement response spectra shown in Figure 21, further validating the reported considerations.

When the earthquake intensity induces moderate structural damage and/or the branches of the elastic response spectra do not show high gradients, no significant differences emerge between elastic and nonlinear response, as evident, for instance, in cases of SPM, ACC, FHC, and MMO. Obviously, the interpretation of the phenomenon is not straightforward, considering the continuous evolution of the degrading processes during the loading histories. Indeed, a succession of time intervals in which the response of the damaged structure is alternatively amplified or reduced with respect to the elastic one could occur.

Finally, Figure 22 gives a direct comparison between the maximum values attained by the displacement for all analyzed events, where the blue and red bars refer to the elastic and damaged cases, respectively.

To further investigation, Figure 23 shows the correlations, in terms of regression lines depicted in red, between the maximum Global Damage Index $GD{I}^{max}$ evaluated for each event with (a) maximum response displacement $|v_A{|}^{max}$, (b) $PGA$, (c) $PGV$, and (d) Housner intensity $I_H$. It appears that the $GD{I}^{max}$ quite well correlates with the Housner intensity, both indexes being evaluated as average quantities. Even more satisfying is the agreement obtained with $|v_A{|}^{max}$ and $PGV$, as opposed to $PGA$, which shows the lowest level of correlation.

An interesting feature is highlighted in Figure 23a: although the general trend is that the higher the maximum displacement, the higher the $GDI$, this is not always true, i.e., an increase of the maximum achieved displacement in some cases corresponds to a decrease of damage. For instance, the highest value of $GD{I}^{max}$ is obtained for MNF, in spite of the highest value of $|v_A{|}^{max}$ corresponding to the ALT record. This is due to the deformed configuration assumed by the structure during the loading history, that influences the areas where maximum tensile strains and, consequently, damage occur. Thus, the pier top displacement $|v_A{|}^{max}$ could be not sufficient alone to fully interpret the structural damage state. Focusing on the responses to MNF and ALT accelerograms, Figure 24a,b show the distribution of the damage variable D on the pier amplified deformed configurations referring to points B and C in Figure 13a and Figure 15a, respectively. It appears that ALT mainly activates mode 1 deformed shape (Figure 24b), whereas MNF seems to combine mode 1 and 2, inducing a more spread damage along the pier height (Figure 24a). Moreover, similar intensities of the maximum local damage D emerge, although these correspond to significantly different values of the pier top displacement, which is $|v_A|=23.4\mathrm{mm}$ and $|v_A|=42.23\mathrm{mm}$ for MNF and ALT, respectively.

The correlations shown in Figure 23a–d are evaluated on the basis of the global damage index defined in Equation (9). If other definitions were adopted, such as that referring to 95% fractile (Figure 6), it is reasonable to assume low variations in correlation levels. Indeed, the comparison between the time histories of the $GDI$ and $D_{0.95}$ in Figure 13b and Figure 15b points out different damage intensities but a similar trend of evolution, as well as analogous ratios $GD{I}^{max}$/$D_{0.95}^{max}$.

In Figure 25, other two correlations are shown: the first (Figure 25a) relates $PGV$ with maximum response displacement $|v_A{|}^{max}$, and the second (Figure 25b) analyzes the relation between Housner intensity $I_H$ and $|v_A{|}^{max}$. These are noteworthy in being characterized by high values of the correlation coefficients, unlike $PGA$ and $|v_A{|}^{max}$ quantities, that are almost unrelated.

Finally, some considerations are made on the single-pier model in comparison with the complete bridge model. The single-pier model neglects several effects, among which the horizontal bending of the “deck” and the torsion of the piers. These effects generally lead to a non-conservative solution. A further effect is related to the non-synchronous motion of the ground, which, however, is outside the scope of this study. The horizontal bending of the deck produces a collaboration between the squat, stiff piers, and the tall, flexible piers. This effect was taken into account by defining an equivalent pier, such as to match the frequency of the first mode of the complete and single-pier model. For the second mode, however, the results of the two models can be significantly different. In fact, it has been noted that, in the complete model, the second mode of the central piers appears in the third and fourth modes, whose periods are equal to 0.672 s and 0.550 s, respectively (see Table 3). In the single-pier model, the second mode has a period of 0.285 s (Figure 8), resulting in a spectral displacement less than that of the complete model for almost all the accelerograms analyzed (Figure 9b, Figure 10b, Figure 11b, Figure 12b, Figure 13b, Figure 14b, Figure 15b, Figure 16b, Figure 17b, Figure 18b, Figure 19b and Figure 20b). An evaluation of the torsional effects in the piers can be performed by means of the complete model by making a comparison between the first mode, which induces the maximum normalized out-of-plane displacement at the top of the central pier (point A) equal to 0.305 mm, and the second mode which induces the maximum torsional rotation at the top of the same pier, which corresponds to the maximum normalized displacements ±0.080 mm at the lateral ends of the pier, equal to 26% of the normalized displacement of the first mode.

5.2.2. Seismic Sequence

As mentioned before, the response of the equivalent pier is here analyzed considering the input motions of the previous section arranged in series. In detail, the seismic sequence is composed by the 12 records listed in Table 2 in order to to reproduce the 20 historical events reported in Table 1. Only the significant part of each event, determined on the basis of the results of Section 5.2.1 in terms of onset of the degrading mechanism and maximum reached damage, was included in the sequence with the aim of reducing the overall time of analysis. Tails of zeroes, properly calibrated for every single accelerogram, were added to dampen the oscillations of each single event. The resulting ground acceleration time history is shown in Figure 26a.

Figure 26b illustrates the result of the simulation through the comparison of the responses corresponding to a linear elastic (blue line) and damaged (red line) constitutive behavior. Figure 26c contains the corresponding time history of the $GDI$. The coupled effect of frequency content of each signal and the accumulated damage during the loading history induces amplification or reduction of the displacements of the damaged structure with respect to the elastic case. To emphasize the influence of pre-existing degradation, Figure 22 summarizes the maximum displacements achieved by the damaging structure when subjected to single seismic events (red bars) and those registered for the sequence of the same signals (yellow bars). It can be noted that displacement peaks might significantly vary, and some comments on this phenomenon are reported in the following, focusing on MFN and ALT events. For MFN record, the difference is not relevant since the damage caused by the individual accelerogram (Figure 13b), i.e., $GD{I}^{max}=0.060$, is slightly lower than that reached during the sequence (Figure 26c), i.e., $GD{I}^{max}=0.062$. Indeed, the structure approaches the MNF earthquake after passing through FORC and RM03 records, which induce a low degradation level, comparable to the initial state of the pier when subjected to the isolated MNF event. The same considerations do not hold for the ALT accelerogram, where different amplification responses are detected in case of seismic sequence and individual earthquake. This discrepancy is explained by the fact that, having experienced a degradation of its mechanical characteristics during the loading history, the pier shows altered dynamic properties.

The values of the maximum displacement experienced by the pier, shown in Figure 22, are also reported in

Table 5. Maximum elastic and damage displacement for individual event (ie) and seismic sequence (ss).

	Elastic Disp. [mm]	Damage Disp. (ie) [mm]	Damage Disp. (ss) [mm]
FORC	8.98	8.42	8.42
RM03	17.56	17.48	17.50
SPM	15.83	14.95	14.75
T1212	19.61	23.78	25.21
MNF	23.59	23.80	24.78
ACC	23.74	22.99	20.96
SPM	15.83	14.95	13.66
ALT	50.11	42.28	36.64
FHC	9.02	9.04	10.46
MZ19	13.80	9.97	9.62
CLO	19.74	17.21	16.72
FHC	9.02	9.04	10.46
SPM${}_2$	20.53	17.42	16.89
FHC	9.02	9.04	10.46
FHC	9.02	9.04	10.46
FHC	9.02	9.04	10.46
MMO	6.05	6.43	6.41
MZ19	13.80	9.97	9.62
SPM	15.83	14.95	13.76
SPM${}_2$	20.53	17.42	16.95

, where the second column refers to the elastic response, while the third and fourth correspond to the damage case for individual event and seismic sequence. Moreover, in

Table 6. Percentage difference between elastic response (er) and damage individual event (die) and seismic sequence (dss) displacement.

	Diff. [%] (er vs. die)	Diff. [%] (er vs. dss)	Diff. [%] (die vs. dss)
FORC	−6.26	−6.26	0
RM03	−0.47	−0.35	−0.12
SPM	−5.60	−6.85	1.34
T1212	21.24	28.54	−5.67
MNF	0.86	5.04	−3.97
ACC	−3.15	−11.71	9.69
SPM	−5.60	−13.72	9.41
ALT	−15.63	−26.88	15.38
FHC	0.26	15.98	−13.56
MZ19	−27.76	−30.26	3.57
CLO	−12.83	−15.31	2.92
FHC	0.26	15.98	−13.56
SPM${}_2$	−15.12	−17.70	3.13
FHC	0.26	16.03	−13.59
FHC	0.26	16.03	−13.59
FHC	0.26	16.03	−13.59
MMO	6.33	5.90	0.40
MZ19	−27.76	−30.26	3.57
SPM	−5.60	−13.65	9.32
SPM${}_2$	−15.12	−17.42	2.79

, the percentage difference for each record is reported, comparing elastic versus damage displacement for individual event (second column) and the seismic sequence (third column). Lastly, in the fourth column, the percentage difference between damage displacement for the individual event versus the value obtained with the seismic sequence is shown. In a few cases, the elastic and damage responses for both cases (single event and sequence) are almost indistinguishable. However, considering the damage responses, displacements are lower or higher in the case of the seismic sequence if compared to the individual event. As described above for some relevant events, this is due to the variation of the mechanical properties caused by the damaging process. The pier subjected to the single events starts from an initial condition characterized by a low damage state induced by self-weight; for the sequence, its properties are affected by the further significant or irrelevant deterioration caused by the previous events. Thus, the variation of the structural stiffness causes the reduction or amplification of the seismic response.

Looking at the evolution of $GDI$ in Figure 26c, it appears that damage mainly grows during FORC, RM03, and MNF records and, then, assumes its constant maximum value $GD{I}^{max}$. The most relevant distributions of the damage variable D on the pier undeformed configurations are shown in Figure 27a–c, referring to the points indicated in Figure 26c with red circle markers. Finally, an interesting feature should be noted: the $GDI$ attains its maximum value, and a similar intensity, under the same seismic event (MNF), both in the case of sequence and individual events.

6. Conclusions

The nonlinear dynamic response of the prominent masonry Italian bridge ‘Ponte delle Torri’, located at Spoleto, was investigated under a set of accelerograms representative of the history of the site. The main 20 seismic events were simulated by means of 12 recorded accelerograms selected from the Engineering Strong Motion Database [35] and characterized by very high correlation level with the historical events in terms of magnitude and epicentral distance of the station. Relying on the results of a previous study [10], the response of the bridge to the chosen records was simulated by adopting a 3D macromechanical finite element model of a properly sized single pier. This allowed for significant reduction of the computational cost, while providing useful information on the bridge’s overall behavior. To account for the nonlinear mechanisms typically characterizing masonry material, the constitutive law proposed in References [13,36] was adopted. This describes the degradation phenomena in a phenomenological fashion by introducing a scalar damage variable in the stress-strain relationship, whose evolution is driven by an equivalent strain measure. A Global Damage Index ($GDI$) was defined to monitor the onset and progression of the damaging process in the overall structure, as well as to provide a measure of the structural integrity.

The performed analyses confirmed several recognized effects caused by degradation on the dynamic response but, at the same time, revealed some new interesting features. The first set of analyses, which involved the pier response to the individual accelerograms starting from the initial undamaged configuration, emphasized that the onset of damage modifies the structural dynamic properties and, in turn, the dynamic amplification of the response [10,41]. The phenomenon is quite complex, as this is influenced by the frequency content of each earthquake. Indeed, the ordinates of the displacement response spectrum may either increase or decrease in correspondence of the periods of the damaged structure, thus inducing amplification or reduction of the displacement with respect to the undamaged case. Correlations between several significant input-output quantities ($PGV$-$GD{I}^{max}$, $PGA$-$GD{I}^{max}$, etc.) were investigated. Among these, the relationship between the maximum Global Damage Index evaluated for each event and the maximum response in terms of displacement of the point A located at the pier top side turned out to be noteworthy. Although the general trend is easily understood (the higher the maximum displacement, the higher the $GD{I}^{max}$), in a few cases, further investigation was required, as an increase in top displacement did not correspond to an increase in damage. This was caused by the deformed configuration assumed by the structure during the loading histories, which can involve higher mode shapes affecting the displacement field and, in turn, intensity and location of damage.

The investigation on the response to the individual earthquakes allowed for the identification of the significant part of each event to be included in the analysis of the seismic sequence. The latter pointed out the effects of the pre-existing damage caused by the deformation history on the structure dynamic response. Hence, different dynamic amplifications emerged with respect to those evaluated in cases of individual events as a consequence of the accumulated damage. Notwithstanding, it appeared that the Global Damage Index attained its maximum value under the same seismic event (i.e., MNF earthquake), with similar intensity in both cases. Thus, the presented study showed that, although the detail of the structural response to the sequence of the seismic records can be different from the response obtained by separately applying each signal, no significant differences emerge as concerns the estimation of the most severe global damage state of the structure. In other words, the analysis performed with the most damaging seismic event gives suitable information also concerning the most severe response to the sequence, at least regarding the Global Damage Index. This interesting feature deserves to be further investigated considering different seismic sequences and constitutive models involving other nonlinear mechanisms, for instance, plasticity and combination of plasticity and damage. Future investigations will also address the combination rule to compute the damage under a sequence starting from the damage under individual events, as well as the effect of modifying the order of the single events in the sequence.