1 Introduction

Before the personal computer became widespread, structural analysis was limited to simplified models requiring reasonable computational effort and equivalent graphical methods [1]. Thanks to the development of powerful computers and refined software, which are able to quickly return a solution for different engineering problems, the attention in structural analysis has shifted to the calibration of the model and the validation of its results.

As a consequence of these technological and scientific advancements, some methods have found limited applications and development (i.e. graphical statics, the force method approach [2], kinematic analysis), while finite element analyses based on the displacement method have found great success and employment in several fields.

However, the representation of complex physical phenomena with elaborated finite element models requires the adoption of numerous input parameters that introduce uncertainties in the solution. To improve the accuracy of the results, different strategies are used to reduce disagreement with the experimental behavior. Model updating is a widely used technique in the calibration of finite element models (FEM), in particular considering global properties such as the dynamic behavior of the structure.

In general, the methods to update the FEM fall into two main categories: direct methods and iterative methods. Direct methods try to reproduce the data obtained from the structure by making small changes to the stiffness and mass matrix, without considering the variation of physical parameters; there are different examples in the literature where this method is applied to small structures in laboratory [3, 4] and slender steel buildings [5]. Iterative methods update the physical parameters of the model until it reproduces the measured data to a sufficient degree of precision, and the error between the experimental data and the numerical results is reduced to a tolerated value [6, 7]; this method is more adaptable, because allows the updating of the physical properties of the model. The development of the update relies on experimental data from dynamic testing, the earliest of these tests on masonry buildings were forced vibration tests [8] or ambient vibration tests on slender structure such as bell towers [9] and bridges [10, 11]. Experimental data acquired with the ambient vibration test (AVT) are usually processed according to the operational modal analysis (OMA) and successfully implemented in different model updating procedures [12]. The selection of the variables used to update the model are often identified by global and local sensitivity analysis [13] that drive the significant changes in the model. The calibration of large and complex models can follow a manual iterative procedure [14, 15] or approximated semiautomated methods to reduce the trial analysis. The Douglas–Reid method [16] is based on the approximation of the natural frequencies with quadratic relationships of the unknown structural parameters; this method has been successfully implemented on civil structures [17] and cultural heritage [18,19,20]. To properly apply this technique, it is necessary to develop an accurate numerical model of the structure, consider the proper sensitive parameters of the problem, and set reliable ranges of variation of these parameters.

In the present paper, different options for the optimizations of three finite element models (FEM1, FEM2, FEM3) of San Giovanni in Macerata were carried out through the Douglas–Reid method. In addition to the dynamic identification and the local tests on the masonry (double flat jack test, sonic tomography and video endoscopy), the information derived from the geotechnical campaign were considered to define the proper range of variation of the material parameters and the boundary conditions. This research aims to select the best set of parameters of the problem after evaluating different solutions for the model updating.

2 San Giovanni in Macerata

2.1 The history of the building

The church of San Giovanni, together with the bell tower and the collegiate is part of an aggregate of buildings located in the south-west side of the well-preserved historic center of Macerata (central Italy), between “Piazza Vittorio Veneto” and the main street “Corso della Repubblica”, as shown in Fig. 1.

Fig. 1
figure 1

Urban plan of Macerata according to G. Lauro in Heroico Splendore delle Città del Mondo (1642), highlighted in red the Church of San Giovanni [17]

Full size image

San Giovanni was one of the last projects of the renowned local Architect Rosato Rosati, who inspired its design based on his earlier work of San Carlo ai Catinari in Rome [21] and died before the completion of the church in 1622. The architecture respects the typical style of the Jesuit’s order: a Latin cross-plan, single wide nave, three deep side chapels on each side, imposing dome with a tall lantern, and a lateral bell tower [22].

The construction of the current framework of San Giovanni commenced in 1610 on top of the previous building attempts of the Jesuit order to erect the church. The nave represents the first core of the building constructed in 1620 [23], and it is covered by a large barrel vault and buttressed on each side by two walls in the transverse direction which shape the six lateral chapels. The sides of the transept and the chancel in addition to the classical hemispherical semi-dome above the apse all maintain the same vault typology, from a second construction stage starting in 1640 when the surrounding properties were acquired, and the 40-m bell tower was erected. In 1712 with the last construction stage, the façade and the collegiate were added (Fig. 2). The construction of the dome suffered numerous delays not only due to technical and managerial issues, but also due to the earthquake of VII grade (Mercalli–Cancani–Sieberg scale) in 1741 [24]. Finally in 1769, the dome, 12 m in diameter with the typical top lantern, was completed [25]. Repairs to the dome were recorded in 1784 and in 1835, while in 1854, the exterior was restored [26]. In Table 1 is recorded the historical evolution of the church considering its three main construction phases.

Fig. 2
figure 2

View of San Giovanni from Corso della Repubblica before the 2016 Central Italy earthquake [18]

Full size image
Table 1 Historical evolution of the church
Full size table

Two significant seismic events of VII grade (MCS) affected Macerata on the 25 August 1809 and on the 1 September 1951 [24, 27]. The most recent major restoration work started in 1995, but it was interrupted two years later due to the earthquake of Umbria and Marche [24]. In 2002 and 2006, the works were restarted and stopped again by the 2006 and 2009 seismic events. Finally, in 2016, the Central Italy Earthquake caused the church to be permanently condemned. Table 2 and Fig. 3 summarize the seismic history of San Giovanni.

Table 2 Relevant seismic events: I0 (epicentral seismic intensity), dist. (distance from Macerata to the epicenter of the earthquake), Mw (epicentral moment magnitude)
Full size table
Fig. 3
figure 3

Significant earthquake events recorded in the area surrounding Macerata

Full size image

It is not possible to properly date and connect the present damage of the building complex to the past seismic events. Due to the frequent interruptions, the extension of the site, and the unclear interventions the main restoration works that were spaced out over two decades starting from 1995, were not able to overcome the safety problems of the construction. In its present state, the church reports two main crack patterns across the nave and by the apse (as shown in Fig. 4). In particular, deeper cracks and concentrated damage are localized on the nave following the keystone of the lateral arches and towards the façade, in proximity of the abutments of the vaults of the lateral chapels, on the North and West side of the drum of the central dome, at the base and top of the lantern, and on the top exterior portion of the façade as reported in Fig. 5.

Fig. 4
figure 4

Crack survey of the church in red the crack extension, in blue the detached material

Full size image
Fig. 5
figure 5

Details of the relevant damages of San Giovanni

Full size image

3 On-site testing campaign

From November 2018 to July 2019, the PRiSMa Lab. (Laboratory of Testing and Research on Structures and Materials) from the Architecture Department of Roma Tre University, carried out extensive onsite testing measurements on San Giovanni. During the onsite campaign different tests were performed, including 4 double flat jack tests, 5 ambient vibration tests (AVT), 47 video endoscopy inspections and 2 sonic tomographies. In the present work, the results from the onsite campaign are implemented in the updating of the numerical model.

3.1 Ambient vibration test

Ambient vibration measurements and output-only dynamic testing are well known techniques for dynamic characterization of historical masonry constructions [7, 28, 29]. The tests were conducted on the 2 and 3 July 2019. The response accelerations of the construction were recorded in 5 different runs, and every run employed 12 uniaxial accelerometers and lasted for about 45 min as reported in Table 3. A total of 16 different points located at four different heights were measured: nave level (14 m), drum level (22 m), dome level (29 m) and lantern level (36 m) as shown in Fig. 6. The sensors were connected to a 24-bit acquisition system equipped with 2 signal conditioners, each acquisition was recorded with 441 Hz sample frequency with dedicated anti-aliasing filter using 12 uniaxial ICP piezoelectric accelerometers (no. 10 PCB 393B12, 10 V/g sensitivity, range ± 0.5 g, frequency range 0, 15–1000 Hz, 8 μg rms resolution and no. 2 PCB 393B31, 10 V/g sensitivity, range ± 0.5 g, frequency range 0, 10–100 Hz, 1 μg rms resolution).

Table 3 Recording information
Full size table
Fig. 6
figure 6

Location of the accelerometers in the AVT

Full size image

OMA is a technique used to process the results from output-only measurements (AVT) based on the hypothesis that white noise conditions allows for not recording the input and in this circumstance, all the natural frequencies of the structure are assumed to be excited [30,31,32].

In this work, the extraction of the modal parameters was carried out using the polyreference least-squares complex frequency-domain algorithm implemented in the software Simcenter Testlab [33]. To optimize the number and the position of the accelerometers, the multi-run operational modal analysis method was adopted, this technique allows to scale the modes based on the most appropriate and complete run, finally, the results are merged to return the global behavior of the construction [34].

A kinematic geometric model was implemented to represent the modal deformations: the nodes with the letters from A to T correspond to the points where the accelerations were measured (master nodes), while the nodes with the numbers are points dependent from the measurements (slave nodes) to improve the representation of the construction. In this work, all five acquisitions have in common at least two horizontal measurements in the reference node B. Being that the church is part of a large irregular building aggregate that requires a significant amount of energy to be excited, only the first three modes of vibration were considered and reported in Fig. 7.

Fig. 7
figure 7

Experimental results from multi-run operational modal analysis

Full size image

3.2 Double flat jack

The double flat jack test is a consolidated semi-destructive technique frequently used in private practice to characterize the mechanical properties of the masonry onsite.

In this work, the tests were performed according to the ASTM Standard C 1197-91, the four locations were decided based on the accessibility and minor impact to the site. The results are reported in Fig. 8.

Fig. 8
figure 8

Double-flat jack test plan and results

Full size image

3.3 Sonic tomography

Sonic tests were performed in two locations in the church (Fig. 9). The speed of the sonic wave was measured on a mesh of points along the two opposite sides of the pillar. A piezoelectric instrumented hammer was used to generate a low-frequency sonic impulse, while the arrival of the propagating wave was captured by a wide spectrum accelerometer and recorded with a CMS HLF-P series SG02_0128 acquisition system equipped with a signal conditioner which allows an optimum amplification of the signal. The propagation of the wave is influenced by the physical–mechanical characteristics of the material, in the case of homogeneous and isotropic material, the dynamic modulus of elasticity (Ed) is evaluated according to the following equation:

$${E}_{{\rm d}}={V}^{2}\rho\frac{\left(1+\nu \right)\left(1-2\nu \right)}{\left(1-\nu \right)},$$
(1)

where \(V\) is the speed of the sonic impulse, \(\rho\) is the density of the material (in this case, estimated as 18 kN/m3 according to the Italian Code [35]) and \(\nu\) the Poisson’s ratio (previously calculated in the same location with the double flat jack test).

Fig. 9
figure 9

a Sonic tomography test plan. b Results of TG02. c Results of TG01

Full size image

In this particular case study, the solid brick masonry with lime mortar was well manufactured and preserved, thus its characteristics can be compared to a low-strength concrete. According to Mehta et al. [36], the dynamic elastic modulus of low-strength concrete is generally 40% higher than the static elastic modulus. Lydon et al. [37] suggests a ratio of 0.83 between the static (\({E}_{{\rm s}}\)) and dynamic (\({E}_{{\rm d}}\)) modulus of elasticity. The value of the static modulus was therefore approximated as an average of these two methods.

The results of TG2 (Fig. 9) showed a fairly homogeneous section with an average velocity of 1450 m/s. The results of TG1 (Fig. 9) instead, highlighted two areas of the section with different velocities, the west part with a lower average velocity of 700 m/s and the east part with a higher average velocity of 1200 m/s. The results and output of the sonic tests are summarized in Table 4.

Table 4 Modulus of elasticity from sonic tomography
Full size table

3.4 Endoscopy

To validate the results from the sonic tomography and reduce the uncertainties in the wall morphology, three video-endoscopy inspections were performed and shown in Fig. 10.

Fig. 10 
figure 10

Video-endoscopy inspection plan and results

Full size image

Inspection E01 shows two slightly different states of conservation of the material: a partially disconnected material at the beginning from 15 to 40 cm and a compact masonry from 45 to the end of the bore hole as confirmed by the velocities of the sonic tomography TG02 in Fig. 9 (from 1200 to 1400 m/s in the first section and from 1400 to 1600 m/s in the second part).

Inspection E02 shows a substantial homogeneous deteriorated masonry: the first section of render, brick and lime mortar of 5 cm corresponds to a sonic velocity of about 800 m/s followed by layers of discontinuous masonry and cavities with sonic velocity between 500 and 700 m/s (Fig. 9).

Inspection E03 shows again an essentially fragmented masonry, the sonic tomography TG1 in Fig. 9 shows a constant velocity of 500–600 m/s for most of the length of the bore hole (which is confirmed by the last section of 95 cm), the first segment of about 15 cm has slightly more compact material (presence of a solid brick) and velocity between 600 and 700 m/s.

The video endoscopy inspections returned an overall similar morphology of the elements composed by solid organized masonry of excellent bricks with congruous amount of lime mortar, the presence of small cavities and discontinuities are in compliant with the construction technology adopted at the time and the state of conservation of the elements.

3.5 Geotechnical survey

With the support of the “Confraternita delle Sacre Stimmate di San Francesco” in Macerata, the owner of the Church of San Giovanni, it was possible to carry out a geotechnical survey of the subsoil inside and outside the church [38].

The analysis showed that the ground conditions do not present any particular geotechnical problem. The foundation soil is mainly made of sands and silty-sands, with silt–clay intercalations. From a global point of view, it is, therefore, possible to assign the soil a mainly granular behavior, only partially cohesive.

The foundations of the church consist of cemented masonry, comparable to the masonry of the elevated structure. These were made at different depths due to the natural slope of the soil towards the south-west of the original surface trend. Based on the carried out geotechnical investigations, the North and South longitudinal sections of the church foundations are schematized in Fig. 11.

Fig. 11
figure 11

Schematic of the geotechnical survey results

Full size image

Three different groups of stiffnesses have been identified based on the geometry of the foundations and their soil interaction by assigning higher stiffness to the elements with greater depth and embedment in the silty sand. The stiffness of the elastic springs supports was estimated according to the Winkler method and the soil parameters suggested by Bowles [39]. The imprint of the church of San Giovanni is shown in Fig. 12 with three different areas representing the different stiffnesses due to the soil structure interaction.

Fig. 12
figure 12

Representation of the three groups of stiffnesses

Full size image

4 Model updating

The model updating is an iterative process made of different phases, which aims at reducing the error between the onsite outcome and the numerical results.

At first, it is suggested to improve the models through the manual tuning of the parameters to control and understand the changes in behavior of the model. Considering that the stiffness of the elastic spring support is the most influential parameter in the dynamics of the model, this parameter was initially varied and analyzed.

Subsequently, the influence of each parameter of the model was quantified with sensitivity analysis.

Once these preliminary analyses were completed, the semi-automatized model updating was carried out through the Douglas–Reid method. The process is schematized by the workflow in Fig. 13.

Fig. 13
figure 13

Flowchart of the model updating process

Full size image

Three different finite element models were implemented with the software Midas Gen [40]. The basic model (FEM1) includes the church of San Giovanni and the convent on the South side of the building complex, it is formed by 18,596 nodes and 9165 three-dimensional brick finite elements with fixed external constraints at the base. The advanced model (FEM2) has the same characteristics previously described for the basic model, but it is externally constrained with three-dimensional elastic spring supports considering three different stiffnesses according to the scheme defined in Fig. 12.

In Fig. 14 are shown the two preliminary models FEM1 and FEM2.

Fig. 14
figure 14

On the left the basic model (FEM1) with all the degrees of freedom of the base nodes constrained, on the right the advanced model (FEM2) constrained with three-dimensional elastic spring supports

Full size image

At the end, a conclusive optimization was performed on the final finite element model (FEM3), derived from the advanced model but with intermediate boundary conditions and proper weight coefficients for the objective function of the optimization.

The model updating of two models with different constraints was performed to relate the evolution of the optimization to the variation of the boundary conditions and tune up the proper spring support [41,42,43].

Both models were divided into macro-volumes considering the structural elements and the construction phases which may influenced their characteristics. Different material parameters were assigned to the macro-volumes according to the results from the double flat jack tests, the sonic topographies and the video endoscopies.

The values of the unknown Young’s modules were estimated through appropriate averages of the onsite testing outputs. In Fig. 15 are reported the material properties assigned to the models.

Fig. 15
figure 15

Subdivision and initial values of the mechanical properties of San Giovanni

Full size image

4.1 Manual tuning

The manual tuning is a well-known initial model optimization [44]; during this procedure, the main parameters are changed manually to understand the behavior of the model and estimate the approximated values of the parameters that reduce the gap between the results.

The constraints are the structural parameters with greatest influence on the dynamic behavior of the structure [42]. Consequently, it was decided to create a preliminary model with uniform elastic spring supports (FEM 0) and to focus on the tuning of foundation stiffnesses. The modulus of elasticity of the horizontal and vertical spring supports have been homogeneously varied considering a fixed ratio of 1–10 [15]. Eigenvalue numerical analyses were performed for each variation, Fig. 16 shows the frequency error (2) between the AVT output and the numerical results. The average minimum error corresponds to the spring stiffness of 480 kN/mm, which was then adopted for the advanced model (FEM2).

Fig. 16
figure 16

Frequency error with the variation of the stiffness of the elastic spring supports

Full size image
$${\rm Error}=\frac{\left|{f}_{{\rm exp}}-{f}_{{\rm FEM}0}\right|}{{f}_{{\rm exp}}}\times 100.$$
(2)

4.2 Sensitivity analysis

The sensitivity analysis refers to the identification of the most influential variables in the response of the model. The dependence of the input parameters on the overall behavior of the FEM is related to the heterogeneity of the masonry, the different levels of damage, and the interaction between the structural elements.

The sensitivity coefficient (3) quantifies the influence of each parameter on the dynamic behavior of the structure, in this analysis, it was considered a variation of 5% from the initial values of Fig. 15 [44, 45].

$$S_{i,j} = 100 \cdot \frac{{X_{j} }}{{R_{i}^{{{\text{FEM}}}} }} \cdot \frac{{\Delta R_{i}^{{{\text{FEM}}}} }}{{\Delta X_{j} }}\begin{array}{*{20}c} {i = 1, \ldots ,M} & {M \, = {\text{ Frequencies}}} \\ {j = 1, \ldots ,N} & {N \, = {\text{ Parameters}}} \\ \end{array}$$
(3)

  • \({X}_{j}\) represents the j-th model parameter;

  • \({R}_{i}^{{\rm FEM}}\) represents the i-th output of the analysis (in this case the natural frequencies);

  • \(\Delta {X}_{j}\) represents the range of variation of the model parameter;

  • \(\Delta {R}_{i}^{{\rm FEM}}\) represents the variation of the output of the analysis.

The parameters with the highest sensitivity coefficient are those that most influence the structure considering the first three modes of vibration. The density, being the parameter that influences all the elements of the church, revealed as the most significant on the modal behavior of the finite element model. The sensitivity coefficients were averaged over the first three modes and proportionated to the global, since the sensitivity coefficient of the density influences 69% of the dynamic response, it was not considered in the graphical layout. Column 4 resulted with the lowest sensitivity coefficients (Fig. 17), while the rest of the parameters reported comparable influences. To reduce the variables of the model updating, the modulus of elasticity from the local testing was assigned to column 3 and column 4 being the less influential and the most certain. Thus, the parameters selected for the model updating are E1, E2, E6, E7, E8, E9, in addition to the more obvious K1, K2, K3 and ρ.

Fig. 17
figure 17

Sensitivity analysis results

Full size image

4.3 The Douglas–Reid method

The Douglas–Reid method is a simplified procedure introduced in 1982 to minimize the error between two sets of data considering different combinations of parameters [46]. In this work, the Douglas–Reid method was applied to minimize the error between both the numerical natural frequencies and modal deformation of the FEM and the estimated frequencies and modal deformation of the OMA.

The modal parameters derived from the OMA, indicated as \({R}_{i}^{{\rm EXP}}\) with 1 ≤ i ≤ M (M represents the number of extracted modes of vibration), are compared with the related \({R}_{i}^{{\rm FEM}}\) values obtained from the analysis of the FEM (function of N significant structural parameter indicated with \({X}_{k}\)) as described by the relation (4).

$$R_{i}^{{{\text{FEM}}}} = R_{i}^{{{\text{FEM}}}} \left( {X_{1} ,X_{2} ,X_{3} , \ldots ,X_{k} , \ldots ,X_{N} } \right)\begin{array}{*{20}c} {i = 1, \ldots ,M} \\ {K = 1, \ldots ,N} \\ \end{array}$$
(4)

To reduce the error between \({R}_{i}^{{\rm FEM}}\) and \({R}_{i}^{{\rm EXP}}\), the proper structural parameters \({X}_{k}\) are selected from the set range of values using a batch analysis procedure.

The Douglas–Reid method uses an approximated model to interpolate with a quadratic function (5) the proper structural parameters within a range of coordinate points \(({X}_{1}, {X}_{2},{X}_{3},\dots ,{X}_{k},\dots ,{X}_{N})\) [41, 47, 48].

$$R_{i}^{{{\text{DR}}}} = C_{i} + \mathop \sum \limits_{K = 1}^{N} \left( {A_{iK} X_{K} + B_{iK} X_{K}^{2} } \right)\begin{array}{*{20}l} {i = 1, \ldots ,M} & \quad {M \, = {\text{ Frequencies or modal deformations}}} \\ {K = 1, \ldots ,N} & \quad {N \, = {\text{ Significant structural parameters}}} \\ \end{array}$$
(5)

  • i represents the i-th target parameter;

  • K represents the k-th significant variable in the FEM;

  • \({R}_{i}^{{\rm DR}}\) represents the structural response of the approximated model;

  • \({C}_{i}\), \({A}_{Ki}\), \({B}_{Ki}\) are 2N + 1 coefficient of the interpolating quadratic function.

The analysis was carried out according to the following steps:

Definition of the lower and upper bound of variation for each parameter (6), in this case, a variation of ± 50% was adopted on all the selected parameters with the exception of the density which was set within a range of ± 10%.

$${X}_{K}^{\rm L}\le {X}_{K}^{\rm B}\le {X}_{K}^{\rm U}.$$
(6)

(2) Each variable was normalized over its nominal initial value (7):

$$\frac{{X}_{K}^{\rm L}}{{X}_{K}^{\rm B}}\le \frac{{X}_{K}^{\rm B}}{{X}_{K}^{\rm B}}\le \frac{{X}_{K}^{\rm U}}{{X}_{K}^{\rm B}}.$$
(7)

(3) The identification of the 2N + 1 coefficients \({A}_{K,i}\), \({B}_{K,i}{ , C}_{i}\) of the approximated model was carried out by equaling the modal parameters of the Douglas–Reid approach (\({R}_{i}^{{\rm DR}}\)) and those obtained from the FEM (\({R}_{i}^{{\rm FEM}}\)).

The Douglas–Reid method was applied both for the modal frequencies and the modal deformations and it was finally combined in the minimization of the error of a single objective function.

4.3.1 The Douglas–Reid for frequency

The system of equations is generated considering the output frequencies from the FEM, where the \({X}_{K}\) variables assume the basic values \({X}_{K}^{\rm B}\), and the output referrers to the 2N models considering once the lower bound value \({X}_{K}^{\rm L}\) and then the upper bound value \({X}_{K}^{\rm U}\). The procedure is represented by the following equation:

$$\begin{array}{l} {f}_{i}^{{\rm DR}}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)={f}_{i}^{{\rm FEM}}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right) \\ {f}_{i}^{{\rm DR}}\left({X}_{1}^{\rm U},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)={f}_{i}^{{\rm FEM}}\left({X}_{1}^{\rm U},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right) \\ {f}_{i}^{{\rm DR}}\left({X}_{1}^{\rm L},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)={f}_{i}^{{\rm FEM}}\left({X}_{1}^{\rm L},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right) \\ \vdots \\ {f}_{i}^{{\rm DR}}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm U}\right)={f}_{i}^{{\rm FEM}}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm U}\right) \\ {f}_{i}^{{\rm DR}}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm L}\right)={f}_{i}^{{\rm FEM}}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm L}\right) \end{array}$$
(8)

(4) The system of Eq. (8) and the relation in (4) and (5) are expressed in matrix form in (9) according to the matrices defined in (10):

$$\left\{{f}_{i}^{{\rm FEM}}\right\}=\left[C\right]\cdot \{{K}_{i}\}$$
(9)

  • \(\left\{{K}_{i}\right\}\) is the vector containing the unknown coefficients \({C}_{i}, { A}_{K,i}\), \({B}_{K,i}\);

  • \(\left[C\right]\) is the matrix containing the combinations of the variables;

    $$\left[C\right]=\left[\begin{array}{ccc}\begin{array}{ccc}1& {X}_{1}^{\rm B}& {X}_{1}^{{B}^{2}}\\ 1& {X}_{1}^{\rm U}& {X}_{1}^{{U}^{2}}\\ 1& {X}_{1}^{\rm L}& {X}_{1}^{{L}^{2}}\end{array}& \cdots & \begin{array}{cc}{X}_{N}^{\rm B}& {X}_{N}^{{B}^{2}}\\ {X}_{N}^{\rm B}& {X}_{1}^{{B}^{2}}\\ {X}_{N}^{\rm B}& {X}_{1}^{{B}^{2}}\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{ccc}1& {X}_{1}^{\rm B}& {X}_{1}^{{B}^{2}}\\ 1& {X}_{1}^{\rm B}& {X}_{1}^{{B}^{2}}\end{array}& \cdots & \begin{array}{cc}{X}_{N}^{\rm U}& {X}_{1}^{{U}^{2}}\\ {X}_{N}^{\rm L}& {X}_{1}^{{L}^{2}}\end{array}\end{array}\right] \quad \left\{{K}_{i}\right\}=\left\{\begin{array}{c}{C}_{i,1}\\ {A}_{i,1}\\ \begin{array}{c}{B}_{i,1}\\ \vdots \\ \begin{array}{c}{A}_{i,N}\\ {B}_{i,N}\end{array}\end{array}\end{array}\right\}$$
    (10)

(5) The system of equations is solved according to (11). Therefore, a matrix of the constants \(\left[K\right]\) containing all the vectors \(\left\{{K}_{i}\right\}\) \({\rm is defined}\) in (12):

$$\left\{{K}_{i}\right\}={\left[C\right]}^{-1} \{{f}_{i}^{{\rm FEM}}\}$$
(11)
$$\left[K\right]=\left[\begin{array}{c}\left\{{K}_{1}\right\}\\ \begin{array}{c}\left\{{K}_{2}\right\}\\ \vdots \\ \left\{{K}_{M}\right\}\end{array}\end{array}\right]$$
(12)

Finally, the approximated frequency \({f}^{{\rm DR}}\) are defined in (13) according to the matrices defined in (14).

$$\left\{{f}^{{\rm DR}}\right\}=\left[K\right]\cdot \{x\}$$
(13)
$$\left[K\right]=\left[\begin{array}{ccc}\begin{array}{ccc}{C}_{{\rm 1,1}}& {A}_{{\rm 1,1}}& {B}_{{\rm 1,1}}\\ {C}_{{\rm 2,1}}& {A}_{{\rm 2,1}}& {B}_{{\rm 2,1}}\end{array}& \cdots & \begin{array}{cc}{A}_{1,N}& {B}_{1,N}\\ {A}_{2,N}& {B}_{2,N}\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{ccc}{C}_{M,1}& {A}_{M,1}& {B}_{M,1}\end{array}& \cdots & \begin{array}{cc}{A}_{M,N}& {B}_{M,N}\end{array}\end{array}\right] \quad \left\{x\right\}=\left\{\begin{array}{c}1\\ {X}_{1}\\ \begin{array}{c}{X}_{1}^{2}\\ \vdots \\ \begin{array}{c}{X}_{N}\\ {X}_{N}^{2}\end{array}\end{array}\end{array}\right\}$$
(14)

4.3.2 The Douglas–Reid for vibration mode

The system of equations for the approximation of the vibration modes is generated starting from the output modes of the FEM, where the \({X}_{K}\) variables assume once the basic value \({X}_{K}^{\rm B}\), then the lower bound value \({X}_{K}^{\rm L}\) and finally the upper bound value \({X}_{K}^{\rm U}\). The procedure is represented by the following equation:

$$\begin{array}{l} {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)\cdots {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)={\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)\cdots {\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right) \\ {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm U},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)\cdots {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm U},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)={\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm U},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)\cdots {\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm U},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right) \\ {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm L},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)\cdots {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm L},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)={\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm L},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right)\cdots {\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm L},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm B}\right) \\ \vdots \\ {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm U}\right)\cdots {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm U}\right)={\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm U}\right)\cdots {\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm U}\right) \\ {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm L}\right)\cdots {\varphi }_{i,j}^{\rm DR}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm L}\right)={\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm L}\right)\cdots {\varphi }_{i,j}^{\rm FEM}\left({X}_{1}^{\rm B},{X}_{2}^{\rm B},\dots , {X}_{N}^{\rm L}\right) \end{array}$$
(15)

  • i = 1, …, M, M = modes of vibration,

  • j = 1, …, S, S = modal displacements,

  • \({1<{\rm X}}_{{\rm K}}<{\rm N}\), N = Significant structural parameters.

(4) The system of equations is also written in matrix form (16), according to the matrices defined in (17):

$$\left\{{\varphi }_{i}^{\rm FEM}\right\}=\left[C\right]\cdot \{{KM}_{i}\}$$
(16)

  • \(\left\{{KM}_{i,j}\right\}\) is the vector containing the unknown coefficients \({C}_{i}\) \(, { A}_{i,j,K}\),\({B}_{i,j,K}\);

  • \(\left[C\right]\) is the matrix containing the combinations of the variables;

    $$\left[C\right]=\left[\begin{array}{ccc}\begin{array}{ccc}1& {X}_{1}^{\rm B}& {X}_{1}^{{B}^{2}}\\ 1& {X}_{1}^{\rm U}& {X}_{1}^{{U}^{2}}\\ 1& {X}_{1}^{\rm L}& {X}_{1}^{{L}^{2}}\end{array}& \cdots & \begin{array}{cc}{X}_{N}^{\rm B}& {X}_{N}^{{B}^{2}}\\ {X}_{N}^{\rm B}& {X}_{1}^{{B}^{2}}\\ {X}_{N}^{\rm B}& {X}_{1}^{{B}^{2}}\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{ccc}1& {X}_{1}^{\rm B}& {X}_{1}^{{B}^{2}}\\ 1& {X}_{1}^{\rm B}& {X}_{1}^{{B}^{2}}\end{array}& \cdots & \begin{array}{cc}{X}_{N}^{\rm U}& {X}_{1}^{{U}^{2}}\\ {X}_{N}^{\rm L}& {X}_{1}^{{L}^{2}}\end{array}\end{array}\right] \quad \left\{{KM}_{i}\right\}=\left\{\begin{array}{ccc}\begin{array}{c}{C}_{i,{\rm 1,1}}\\ {A}_{i,{\rm 1,1}}\\ {B}_{i,{\rm 1,1}}\end{array}& \cdots & \begin{array}{c}{C}_{i,S,1}\\ {A}_{i,S,1}\\ {B}_{i,S,1}\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{c}{A}_{i,1,N}\\ {B}_{i,1,N}\end{array}& \cdots & \begin{array}{c}{C}_{i,S,N}\\ \begin{array}{c}{A}_{i,S,N}\\ {B}_{i,S,N}\end{array}\end{array}\end{array}\right\}$$
    (17)

(5) The system of equations is solved according to (18). Therefore, it is defined a matrix of the constants \(\left[KM\right]\) containing all the vectors \(\left\{{KM}_{i}\right\}\) as shown in (19):

$$\left\{{KM}_{i}\right\}={\left[C\right]}^{-1} \{{\varphi }_{i}^{\rm FEM}\}$$
(18)
$$\left[KM\right]=\left[\begin{array}{c}\left\{{KM}_{1}\right\}\\ \begin{array}{c}\left\{{KM}_{2}\right\}\\ \vdots \\ \left\{{KM}_{M}\right\}\end{array}\end{array}\right]$$
(19)

Finally, it is defined \({\varphi }^{\rm DR}\) as shown in (20) according to the matrixes defined in (21):

$$\left\{{\varphi }^{\rm DR}\right\}=\left[KM\right]\cdot \left\{x\right\}$$
(20)
$$\left[KM\right]=\left[\begin{array}{ccc}\begin{array}{ccc}{C}_{{\rm 1,1},1}& {A}_{{\rm 1,1},1}& {B}_{{\rm 1,1},1}\\ {C}_{{\rm 2,2},1}& {A}_{{\rm 2,2},1}& {B}_{{\rm 2,2},1}\end{array}& \cdots & \begin{array}{cc}{A}_{{\rm 1,1},N}& {B}_{{\rm 1,1},N}\\ {A}_{{\rm 2,2},N}& {B}_{{\rm 2,2},N}\end{array}\\ \vdots & \ddots & \vdots \\ \begin{array}{ccc}{C}_{M,S,1}& {A}_{M,S,1}& {B}_{M,S,1}\end{array}& \cdots & \begin{array}{cc}{A}_{M,S,N}& {B}_{M,S,N}\end{array}\end{array}\right] \quad \left\{x\right\}=\left\{\begin{array}{c}1\\ {X}_{1}\\ \begin{array}{c}{X}_{1}^{2}\\ \vdots \\ \begin{array}{c}{X}_{N}\\ {X}_{N}^{2}\end{array}\end{array}\end{array}\right\}$$
(21)

4.3.3 The objective function

The choice of the objective function is a fundamental aspect of the optimization process, this represents the ultimate step to be implemented in the algorithm. In this work, the algorithm “fmincon” (find minimum of constrained nonlinear multivariable function) was implemented on MATLAB [49] to minimize the gap between the numerical and experimental results, objective function (22) represents the error considering both modal frequencies and modal shape vectors.

$$f\left( x \right) = \mathop \sum \limits_{i = 1}^{M} \left[ {\alpha_{i} \left( {\frac{{f^{{{\text{DR}}}} - f^{{{\text{EXP}}}} }}{{f^{{{\text{EXP}}}} }}} \right)^{2} + \beta_{i} \left( {\frac{{1 - \sqrt {{\text{MAC}}_{i} } }}{{{\text{MAC}}_{i} }}} \right)} \right],\quad i \, = \, 1, \, 2, \, 3\;{\text{are modes of vibration}}$$
(22)

where the MAC is defined according to [50] as

$${\rm MAC}_{i}\left({\varphi }_{i,j}^{\rm DR};{\varphi }_{i,j}^{\rm EXP}\right)=\frac{{\left({\left({\varphi }_{i,j}^{\rm DR}\right)}^{T}\cdot {\varphi }_{i,j}^{\rm EXP}\right)}^{2}}{\left({\left({\varphi }_{i,j}^{\rm DR}\right)}^{T}\cdot {\varphi }_{i,j}^{\rm DR}\right)\cdot \left({\left({\varphi }_{i,j}^{\rm EXP}\right)}^{T}\cdot {\varphi }_{i,j}^{\rm EXP}\right)}$$

The factors αi and βi are the weight coefficients used to variate the influence of the modal frequencies and the modal shape vectors on the objective function. Five sets of weight coefficients were considered to investigate the optimization algorithms and tune up the final updating of FEM3 the five sets of coefficients are presented in Model Updating Ver. 0 to Ver. 4 of Table 5.

Table 5 Weight factors in the model updating
Full size table

5 Model updating results

The results of the model updating on the basic model (FEM1 with fixed constraints) and the advanced mode (FEM2 with elastic spring supports) are reported in Tables 6 and 7.

Table 6 Values of the structural parameters of FEM1 after the model updating
Full size table
Table 7 Values of the structural parameters of FEM2 after the model updating
Full size table

The diagrams in Figs. 18 and 19 summarize the percentage of variation of the optimized parameters due to the model updating.

Fig. 18
figure 18

Variation of the structural parameters of the model updating on FEM1

Full size image
Fig. 19
figure 19

Variation of the structural parameters of the model updating on FEM2

Full size image

FEM2 with the springs supports has a greater adaptability than the FEM1 with fixed constrains. In fact, in FEM1, the use of fewer parameters for the optimization, forces the parameters E8 and E9 on the extremes of the range of variation. On the other hand, the optimization of FEM2 yields to more rigid spring supports and reduced the errors with the experimental results. By incrementing the stiffness of the support, the behavior of FEM2 would tend to get closer to the one of FEM1, consequently by performing this comparison it is possible to adjust the range of variation and the initial value of the final optimization on FEM3.

5.1 Frequencies

The frequencies variations of the model updating of FEM1 (fix constrains model) and FEM2 (elastic spring supports model) are reported in Tables 8 and 9.

Table 8 Frequencies variation of the model updating on FEM1
Full size table
Table 9 Frequencies variation of the model updating on FEM2
Full size table

The results are summarized in Fig. 20. After the optimization process, the modal frequencies of FEM1 resulted higher than the experimental frequencies, while the modal frequencies of FEM2 resulted lower (FEM1 positioned above the target line while FEM 2 below), this behavior suggests the study of a final model with average characteristics of the foundations. Moreover, by increasing β (weight factor of the modal deformations), the distance of the numerical frequencies from the target line becomes greater especially for the updating versions 2, 3 and 4.

Fig. 20
figure 20

Frequencies results of the model updating

Full size image

5.2 Vibration modes

Similarly, the same process was performed to control the optimization of the vibration modes. MAC allowed the comparison between the i-th numerical mode with the i-th experimental mode, values close to the unit represent an optimal correlation between the numerical and the experimental modal shape vectors.

In Fig. 21, it is reported the MAC of the first three modes of vibration. Overall, FEM2 returned higher MAC than FEM1 which corresponds to a more accurate approximation of the experimental modal shape vectors in the model with the elastic spring supports. It is relevant that in FEM1 by raising the value of β, the MAC relative to mode 2 and 3 consequently increase while in FEM2, this trend is not evident.

Fig. 21
figure 21

Modal assurance criterion after the model updating

Full size image

6 Final finite element model

The conclusive optimization was performed on the final finite element model (FEM3). In Fig. 20, the ver. 0 updating of FEM2 reported the overall minimum frequency error, while in Fig. 21, it is reported the irrelevant variation of the MAC. Thus, for the final optimization, the weight coefficients of model updating ver. 0 with α and β equal respectively to 1 and 0.5 were chosen. The initial values of the significant parameters were modified considering the intermediate value of the results relative to the ver. 0 model updating as reported in Fig. 22.

Fig. 22
figure 22

Subdivision and initial values of the structural parameters of FEM3

Full size image

6.1 Final model updating results

The range of variation of the parameters have been carefully modified to assure physical and engineering coherence of the solution, the results are presented in Table 10.

Table 10 Results from the final model updating of FEM3
Full size table

From the trend of variation of the parameters summarized in Fig. 23. An overall reduction of stiffness of the nave, columns 1 and 2, and the convent was observed in respect to the estimated initial values. In the nave, this variation is possibly connected to the underestimation of the initially identified damage or its further extension in inaccessible areas. As highlighted by Table 1, the abutments 1 and 2 and the convent were constructed in different periods, this could have resulted in poor connectivity between these structural parts or localized microcracking of columns 1 and 2 due to their compression behavior. Another explanation lies also in the initial overestimation of the modulus of elasticity of the masonry in the convent, which derives from measurements on the exterior wall of the complex, usually tougher than the interior walls and often consolidated by periodical interventions. No onsite measurements were collected on the masonry of the dome and apse; however, an increment in the stiffness was observed in respect to the approximated initial values. In this case, the estimation was less accurate (overall highest variation of the parameters) and possibly these parts better sustained the seismic actions than initially hypothesized with the low modulus of elasticity. The foundations, which often play a crucial role in the dynamic behavior of a building, were repeatedly refined during the updating process and they were finally modelled with considerably stiff elastic springs that (in agreement with the geotechnical survey [27]) recall the well-preserved conditions and its good interaction with the subsoil.

Fig. 23
figure 23

Variation of the parameter of FEM3 after model updating version 0

Full size image

The final model updating of FEM3 returned in terms of frequencies the best correlation with the experimental results as shown in Fig. 24. The gap of the numerical frequencies from the experimental data is contained in a range of ± 0.13 Hz and an average error of 4.63% (Fig. 24).

Fig. 24
figure 24

Results in terms of modal frequency of the final model updating of FEM3

Full size image

In terms of modal shapes, the final model updating of FEM3 returned overall satisfying results considering the low versatility of the finite element model to approximate the experimental values. A MAC of 0.83 for the first mode and an average of 0.75 for the overall three modes express an improvement of the results compared to the initial results. Moreover, sufficient independency between the modes expressed by the contained MACs between the different modes enhanced the reliability of the results (Fig. 25).

Fig. 25
figure 25

Results in terms of modal shapes of the final model updating of FEM3

Full size image

The final updated finite element model (FEM3) allows a reliable and detailed study of the dynamic behavior of the church of San Giovanni as reported in Fig. 26. The first vibration mode presents a uniform translational along the transversal axis of the church (y-axis), with an evident translational movement of the tower. The second vibration mode presents a uniform translational mode mainly in the longitudinal direction of the church (x-axis), with a small component in the transversal direction which involves both the church and its bell tower. The third vibration mode has a main torsional component (rotation around the z-axis).

Fig. 26
figure 26

Numerical modal analysis of the updated FEM3

Full size image

7 Conclusion

This work aims to create a finite element model of the church of San Giovanni with a dynamic behavior as close as possible to the results from the AVT by integrating the results from the local tests on the masonry and the geotechnical investigations. The ambient vibration tests and operational modal analysis technique returned fundamental global information on the seismic behavior of the church, while appropriate tests such as sonic tomography, video endoscopy, double flat jacks test, and geotechnical survey assessed the relevant characteristics of the structure and supported the estimation of the unknown parameters.

Considering these outputs, two different preliminary models (FEM1, FEM2) were implemented with two different constraint conditions (fixed constraints and elastic springs support), since the type of constraint greatly influences the dynamic behavior of the model.

A preliminary manual tuning approximated the stiffness of the springs at the base by minimizing the average frequency error. Sensitivity analysis evaluated the influence of the model parameters on the dynamic behavior of the church and identified the most important parameters for the model updating.

Subsequently, the model updating was calibrated with the Douglas–Reid method considering five weight coefficients (α, β) for the objective function with the aim of improving both correlation with the experimental modal frequencies and modal shape vectors.

Finally, the last model updating of FEM3 with the appropriate weight coefficients α and β, boundary conditions and input values returned the best approximation of the dynamic identification, the average frequency error has been reduced to 4.63%, while the average MAC increased to 0.75.

The optimized parameters present acceptable values in accordance with the local masonry, the use of the model updating yielded to a small variation of the known modulus of elasticity and supported the estimation of the untested structural elements. The variation of the updated modulus of elasticity can be related to the reduction of stiffness and connectivity due to damage.

The presented study outlines a process to fully exploit the information gained during an exhaustive onsite testing campaign to create a reliable finite element model. This procedure, if properly implemented, can be considered as a tool in the damage assessment of historical large masonry buildings affected by seismic events and it can guide further investigations on the construction.