Analysis of Seismic Action on the Tie Rod System in Historic Buildings Using Finite Element Model Updating

: Historic buildings have a high architectural value and their maintenance, repair and rehabilitation require a special approach. This approach is mainly based on the buildings’ performance under non-destructive tests such as operational modal analysis (OMA). Under extreme loads, such as earthquakes, these buildings require representative numerical models to simulate their expected response. In historic buildings, tie rods transfer axial loads and are typically used to balance horizontal trust due to static and dynamic loads associated with seismic actions. It is very important to determine the possibility of exceeding their load-bearing capacity under extreme loads, such as an earthquake. In this context, this paper presents an approach for the analysis of seismic action on the tie rod system in a historic building. The analysis was performed by combining the on-site experimental testing and the finite element model updating (FEMU) of the local models of tie rods and the global model of the structure. It was shown that the combination of analyzing local and global structural models, experimental on-site testing and FEMU is a viable solution for assessment of historic buildings’ load bearing capacity.


Introduction
The use of stone or masonry elements in buildings is one of the oldest methods still in use today. The durability of historical buildings is proven by the number of such structures which are still in use. Because of their high mass, the problem of high compressive forces occurs in the elements of structures built from stone. This problem is particularly interesting in buildings whose load-bearing systems are arches, buttresses, vaults, walls, and domes. The load-bearing system is mainly supported by masonry walls. The walls must support this load in such a way that only compressive stresses occur in them. This is achieved by making the load-bearing system much heavier or thicker. In order to keep these dimensions within certain limits characteristic of stone and masonry, the horizontal components of the compressive forces are balanced using tie rods [1]. These structural elements control the horizontal forces and displacements caused by static and dynamic loads associated with seismic actions. In extreme cases, they can reach the maximum loadbearing capacity by being subjected to high stresses and pulled out of their anchorage ( Figure 1). Therefore, the value of the axial force in such systems is often a topic of discussion. There are several uncertainties in determining the axial force value in this type of structural element, which include the complex boundary conditions [2][3][4][5] and geometrical and material properties [6]. Several approaches are used to estimate the force values in tie rods-static [7], dynamic [8], and mixed approaches [9]. In this paper, the previous approaches are only mentioned, while their detailed description can be found in [10]. In addition to the above-mentioned approaches, a combination of experimental onsite testing and numerical analysis based on the model updating approach can also be used. Model updating was developed in the 1990s as a topic of great importance to the design, construction and maintenance of mechanical systems and structures [11]. It refers to updating the finite element (FE) numerical model to adopt more precise and accurate structural dynamics parameters based on the known dynamics or statical parameters obtained through experimental investigations of a structure [12][13][14][15][16] (Figure 2). During the model updating procedure, the most commonly used model updating parameters are the material [16][17][18][19][20][21][22][23] and geometrical characteristics [23][24][25][26] and boundary conditions [27][28][29], but their selection is generally based on performing the sensitivity analysis [30]. The selection of materials and the definition of their properties cause great difficulties in the development of a numerical model due to their specific mechanical properties. For historical buildings, it is difficult to obtain documentation that provides exact, reliable information about the used material and its properties. Except for the problem of the availability of documentation, an experimental approach to determine the properties is challenging because it requires quite extensive non-destructive testing [31][32][33][34]. Regarding the geometric issues, historic buildings are also known for their complex geometry, resulting from the use of different types of structural elements such as vaults, arches, columns and domes. Moreover, the integrations of very slender and massive elements, for example, columns and walls, are common. The aforementioned causes of complex geometry require special attention when modelling these structural elements and their connections. The complexity of finite elements and the ways in which they can be connected, as well as the modelling of connections, have been addressed in several studies [35][36][37][38]. Although for some buildings there is no project documentation, nowadays the geometry is not a major problem, as there are several modern methods that can be used to determine the geometry very easily [39][40][41][42]. Besides the correct definition of the geometrical and material properties of the numerical model, the boundary conditions are also very important and have an impact on the dynamic behavior of the structure. When developing a numerical model, their definition can be very demanding and challenging, especially since the scatter of results for different boundary conditions can be very high [43][44][45]. Experimental investigations are the best control and help to define the boundary conditions of the model. Using experimentally determined structural dynamic parameters (natural frequencies, mode shapes, damping ratios) together with the model updating methods, it is possible to predict the unknown system properties of a real structure (geometry, material properties, boundary conditions) with quite high accuracy and reliability [46][47][48]. This paper is an extended version of the paper published in [49]. We extend our previous work by performing the operational modal analysis of the cathedral (Figure 3) to determine the mode shapes of the entire structure, that we used to validate the FE model and to update the finite element model. Based on the experimental testing and the application of operational modal analysis (OMA), the structural dynamics parameter, natural frequencies and mode shapes of the whole structure and tie rods were determined. The initial values of the forces in the tie rods were taken from authors' previous paper [10]. The initial numerical model of the cathedral was iteratively updated based on the experimentally determined mode shapes and natural frequency values of the cathedral. In order to select the model updating parameters, a sensitivity analysis was performed. Based on this analysis, the stone masonry elasticity modulus and the boundary conditions were selected as the model updating parameters. Since there were few selected updating parameters, the FEMU was performed with manual updating. After the initial state of the tie rods was obtained, and the global numerical model of the cathedral was updated, an analysis of the internal forces and stresses levels in the tie rods caused by the characteristics of the seismic action for the Šibenik area was performed. On this basis, it was determined whether the tie rods will reach their maximum load-bearing capacity during an earthquake.
In the following text, Section 2 describes the procedure of the experimental determination of the structural dynamic parameters of tie rods and the cathedral; Section 3 describes the numerical analysis of the initial numerical model of the cathedral; Section 4 describes the procedure of model updating; Section 5 presents the results of the analysis of internal forces and stresses levels in the tie rods due to seismic action; and Section 6 presents the resulting discussion and conclusions.

Experimental Analysis
Within the experimental investigation of St. James Cathedral in Šibenik (Croatia), operational modal analysis (OMA) was used to determine the structural dynamic parameters of the whole structure and the tie rods [10]. The description of the structure and the experimental determination of the dynamic parameters of the tie rods can be found in [10], while the experimental investigation of the cathedral will be described in the following text.

Experimental Investigation of the Structure
The experimental investigation of the structure was performed by using piezoelectric accelerometers (PCB Piezotronics, type 393B31 with nominal sensitivity of 10 V/g, PCB Piezotronics, Depew, NY, USA), an analyzer (Bruel and Kjaer, type 3560c, Bruel and Kjaer, Naerum, Denmark) and associated software. The measurement points were defined in 9 axes in the transverse direction (axes 1-9), 4 axes in the longitudinal direction (axes A-D) (Figure 4a), and in 5 axes along the height of the structure (axes R1-R5) (Figure 4b) according to accessibility within the structure. The measurement was performed at each of the measuring points in two perpendicular directions, x and y ( Figure 4a). In this way, 216 measuring points were obtained, the sum of which corresponds to two horizontal directions. The reference measuring point was set at point 7 C at level R3.  The mode shapes ( Figure 5) and natural frequency values (Table 1) were obtained during ambient excitation. Frequency domain decomposition (FDD) and enhanced frequency domain decomposition (EFDD) were used to estimate the modal parameters. The values of the experimentally obtained natural frequencies (Table 1)

Experimental Investigation of Tie Rods
In addition to the natural frequencies and mode shapes of the whole structure, dynamic parameters of the tie rods were also obtained. Since the tie rods system consists of aluminum and steel tie rods, for simplicity one aluminum and one cast steel tie rod were taken as a reference. A more detailed analysis was performed on these tie rods, as described in [10]. The analysis is based on the use of the experimentally determined natural frequency values and mode shapes to update the finite element model of the tie rods. In addition, the existing equations for determining the axial force values in tie rods were corrected according to the updated finite element model of the tie rods. As the results are shown in [10], for simplicity they will not be described in this paper.

Numerical Analysis
The experimentally obtained results were used to update the numerical model of the tie rods and the cathedral. Both the local model of tie rods and the global model of the cathedral were developed in SAP 2000 (Computers and structures, Walnut Creek, CA, USA). The numerical model of tie rods was modeled as the Euler Bernoulli beam element with the released bending and torsional moment. In the global model of cathedral columns, arches and tie rods were developed as a beam element, while walls, roof, and floors were developed as a tick shell element. A detailed description of the numerical model of tie rods can be found in [10], while the global model of the cathedral is described in the next subsection.

Initial Numerical Model of the Cathedral
Based on the existing drawings and available 3D scans, a numerical model of St. James Cathedral in Šibenik was developed. From 3D scans and the obtained cross sections of the cathedral, the geometry of the initial numerical model was first developed in Auto-CAD (Autodesk AutoCAD 2019, educational, Portland Oregon), after which it was imported into SAP. After importing the model, the elements are associated with the corresponding geometric and material characteristics. As part of the FEMU, static and dynamic analysis was performed by associating the mechanical characteristics of the material and determining the natural frequencies and mode shapes of the whole structure. The values of elastic modulus and material density of cast steel are adopted from the previous paper [10] while the material properties of the aluminum were determined using the same procedure as for the cast steel tie rods. The obtained properties of the aluminum tie rods were: 70 GPa for elastic modulus and 2660 kg/m 3 for material density. For the stone masonry, material density was 2490 kg/m 3 and the modulus of the elasticity was assumed as 20 GPa. The material behavior of stone was described by linear theory and homogeneous behavior was assumed. The sum of the material characteristics assigned to different cross-section properties of the global model of St. James Cathedral is shown in Table 2. The numerical model was described using the finite element method (FEM). The size of the walls' finite elements was chosen according to the approximate dimensions of the stone blocks with which the cathedral was built, and it was 0.5 m × 0.5 m. Thickness values of stone blocks were read from the 3D scans of the cathedral and are different for different parts of the cathedral (0.45-1.00 m). In the area around the opening, portal, stained glass, rosette, the size of the finite elements was adjusted so that their shape can better describe the real situation and not affect possible errors in calculation. Since the various roof slabs and a half-domes are connected to the walls of the east facade, the finite elements were adjusted. The finite elements in this part of the structure were meshed more densely and mostly described by triangular finite elements. The complete numerical model is shown in Figure 7. The tie rod system is put into the three levels, and it was modeled as a series of beam elements clamped into walls and columns with the cross Section 55 mm × 55 mm to 60 mm × 60 mm for cast steel tie rods and 30 mm × 50 mm for aluminum tie rods. Their initial stress was obtained using the previously described method presented in [10]. The results of the initial stress calculation of the tie rod for each level of the structure can be seen in Figure 8. According to the obtained values of the stress level, in the global numerical model of the cathedral, the tie rods are subjected to axial tension. Since the foundation conditions are not known, for the initial numerical model, instead of the foundation strips under the walls and columns, hinge boundary conditions were defined.

Results of the Initial Numerical Model of the Cathedral
On the initial numerical model for which the hinge boundary condition was proposed and elasticity modulus of stone masonry as 20 GPa, the numerical modal analysis was performed through which the natural frequency values and mode shapes were obtained. To assess how well the initial numerical model describes the real state of the structure, the numerically obtained natural frequency values (Table 3) were compared with the experimentally obtained one (Table 1). Based on the obtained natural frequency error (Table 3), it can be concluded that the initial numerical model deviates significantly from the real state of the structure, where the deviation in the natural frequencies values is up to 16.5%, but the mode shapes are visually well correlated. Bending around x axis, antimetrical 5.14 5.55 6.5 3 Second bending mode shape around x axis 6.24 5.85 6.3 Given the natural frequency error (Table 3) in order to approach the determination of stress levels in tie rods due to the action of seismic load, the global numerical model of the cathedral needs to be updated. The next section describes the selection of the parameters used for updating the model, and the results of the model updating performed on the global numerical model of the cathedral.

Finite Element Model Updating
First, the sensitivity analysis was performed to determine the influence of different structural parameters on structural dynamic behavior. The geometrical characteristics of the structure are well known from the 3D scan. As for the material characteristics, the only unknown parameter is the elasticity modulus of the stone masonry. Besides, the principle of the foundation is also unknown, but geotechnical soil testing data are available [50]. Therefore, the soil stiffness obtained by geotechnical testing varies significantly for certain parts of the area where the cathedral is located, so the boundary conditions are taken as additional model updating parameters. As there were only two updating parameters selected, the manual model update was performed in such a way that the value of the elasticity modulus of stone masonry E n num was iterated for two boundary conditionshinge and clamped ( Table 4). The interaction between the structure and the ground is described by the same boundary conditions. The modulus of the elasticity of stone masonry E n num was iterated with the step of approximately 0.5 GPa to 2 GPa within the lower bound of 20 GPa and upper bound of the 35 GPa.

Manual FEMU
Considering the different types of boundary conditions-hinge and clamped-the elasticity modulus value was updated. The obtained changes of the first three natural frequencies for different values of the dynamic elasticity modulus and for a different type of boundary condition are summed in Table 4. The experimentally obtained natural frequency values are indicated in the graph by the dashed line, where the first experimentally obtained natural frequency is shown in black, the second in red and the third in green, respectively. Following the known experimental values of natural frequency from the graph (Table 4) For the assumed value of the elasticity modulus, an analysis of the changing natural frequency values depending on the change of the spring stiffness was performed. Soil stiffness was defined by line springs where the spring stiffness was iterated between 5,000-20,000 kN/m, which corresponds to the soil characteristics determined by geotechnical testing. Due to the simplicity of the FE model of the cathedral, the same spring stiffness values are defined for each degree of freedom. Based on the performed iteration procedure, it can be concluded that the second and third modes are insensitive to spring stiffness.
The sum of the performed FEMU iteration steps and results can be found in the following table (Table 4).

Evaluation of the Updated Numerical Model Using Modal Assurance Factor Criteria (MAC)
In this section previously updated numerical models based on the variation of elasticity modulus and boundary conditions will be compared. The first numerical model (1.HBC) describes hinge boundary conditions with an adopted elasticity modulus E = 26.3 GPa. The second numerical model (2.CBC) describes the clamped boundary conditions with an adopted elasticity modulus E = 23.2 GPa. The third numerical model (3.ELS) describes boundary conditions approximated by elastic springs with adopted stiffness k = 15,000 kN/m and elasticity modulus E = 25 GPa. The comparison will be based on the MAC factor criteria Equation (1) of experimentally obtained mode shapes of the structure with the mode shapes of the corresponding updated numerical models. The characteristics of the estimated numerical models are summarized in Table 5. For simplicity, the comparison of the experimentally, ϕ n,i exp , and numerically, ϕ n,j num , obtained mode shapes were performed at seven measurement points in the middle of the arches which separate the roof slabs of the middle nave ( Figure 9). In this way, the MAC Equation (1) factor values of the observed mode shapes were obtained (Table 5). Based on the MAC factor, the model that the beast corresponds to the real state of the structure (the model whose MAC factor assumes the highest values for all three mode shapes) will be selected. A graphical representation of the experimentally obtained normalized mode shape vector, ϕ n,i exp of the seven selected points is given in Figure 10. Since in each of the mode shapes the dominant horizontal displacement is in the y direction, the displacement in the x and z directions is negligible. The comparison of mode shapes will be performed based on the dominant direction y.

Results of the Evaluation of the Numerical Model
The numerical mode shape vector, ϕ n,j num , for each developed model (Table 5) for three mode shapes were read from the numerical model and afterward normalized to the unit (Table 6a). Using MAC factor the numerically obtained mode shape vector ϕ n,j num was compared to the experimentally obtained one ϕ n,j exp (Table 6b). The graphical representation of the obtained MAC factor for each of the developed numerical models is shown in Table 6. From the obtained MAC matrix (Table 6b) an excellent correlation of numerical model 1 (HBC) with the real state of structure can be noticed. The MAC factors for all three mode shapes take extremely large values. From the obtained MAC matrix for the numerical model 2, CBC, it can be seen that the MAC factors for first and third mode shapes have very high values, while for the second mode shapes this value is much lower. As for the third numerical model (3. ELS), it can be seen that the first numerically obtained mode shape vector, of all developed numerical models best correlate with the experimentally obtained one. However, considering the second and third numerically obtained mode shape vectors there is a significant deviation. This is supported by the fact that only the first numerically obtained mode shape vector was achieved for the selected spring stiffness. For the second mode shape vector, the natural frequency value deviated significantly from the experimentally determined value. In addition to the deviation of the natural frequency, the MAC factor is very low and is only 0.17. In the third mode shapes, for assumed soil stiffness, a frequency value slightly closer to the experimentally determined value was achieved, but still insufficient to make the MAC factor value satisfactory [51]   (c) Graphical representation of MAC matrices for comparison of the numerically and experimentally obtained normalized mode shape vectors for the selected measurement points (1-7) As the MAC factor values for numerical model 1.HBC achieve the highest values for all three mode shapes, it was selected as the model on which an analysis of force in tie rods for seismic action will be performed.

Analysis of Seismic Action of the Tie Rod Systems
After the numerical model and its properties were obtained, the seismic performance assessment of the global numerical model was performed using the modal response spectrum analysis according to the EN 1998-1:2004 [52]. For the defined the load combination, which includes the dominant action of the earthquake in the x-direction (ERTHSTX) and y-direction (ERTHSTY), the values of the force in the tie rod were determined. The load combinations for dominant seismic actions in the x (Equation (2)) and y directions (Equation (3)) include the following load patterns: self-weight (dead load) (SW), additional constant load (G), service load (Q), seismic action in the x-direction (ERTHX), and seismic action in the y-direction (ERTHY). The load combination for the dominant seismic action in the x-direction (ERTHSTX) contains all load patterns multiplied by the scale factor 1.0, except for the load pattern for the seismic action in the y-direction and the service load. The load pattern for the seismic action in the y-direction and the service load were multiplied by the scale factor 0.3. The load combination for the dominant seismic action in the y-direction (ERTHSTY) contains the same load pattern, except for the service load, and the load pattern for seismic action in the x-direction was also multiplied by the scale factor 0.3 instead of the y-direction. ERTHSTX= 1.0·SW+0.3·Q+1.0·G+1.0·ERTHX+0.3·ERTHY (2) According to the available map of the seismic area in Croatia and Šibenik, expected horizontal peak soil acceleration was read. For the return period of 95 years, its value is 0.094 g, while for the return period of 475 years, its value is 0.189. According to the defined soil type category C, the values of elastic response spectrum type 1 were determined. Considering the type of construction, a behavior factor q of 1.5 was assumed for the project spectrum. The seismic action was analyzed for two axes-the x axis and the y axis. The deformation of the 3d model of the updated finite element model of the cathedral due to seismic action in the x and y directions is shown in the following (Figure 11). After the definition of the seismic action the stress levels in the tie rods were determined (Tables 7 and 8). The graphical comparison of the stress levels in the aluminum tie rods at the level R2 is shown in Figures 12 and 13, while for the cast steel tie rods graphical comparison is shown in Figure 14. Table 7. Values of forces and stress levels in the aluminum tie rods in the axis at the level R2 for the seismic action in the x and y-direction.

Discussion and Conclusions
This research presents the method for the analysis of the possibility of exceeding the load bearing capacity of tie rods in a historic building (Cathedral of St. James in Šibenik) under seismic action. The performed analysis is based on the combination of structural dynamic parameters obtained by the experimental on-site investigation and the finite element model updating (FEMU) procedure. In order to obtain the initial state of each aluminum and cast steel tie rod, the FEMU of each representative tie rod was performed by using the approach elaborated by the authors in [10]. Then, by using the elastic modulus of the stone masonry and the boundary conditions as the update parameters, the global numerical model of the cathedral was updated in order to obtain a representative global numerical model describing the real state of the structure. By subjecting the updated global numerical model to seismic load according to EN 1998-1:2004, the stress level in the tie rods was determined. Based on the obtained stress levels, it can be concluded that due to the seismic action, the axial forces in the tie-rods increase. For aluminum tie rods this increase was 13.25% with st. deviation 0.85% for seismic action in x direction. For seismic action in y direction, this increase was 10.80% with st. deviation of 0.63%. The increase in the stress levels for cast steel tie rods was 16.15% with st. deviation of 1.15% for seismic action in x direction. For seismic action in the y direction, the increase in the stress level was 18.72% with a standard deviation of 3.99%. For the seismic action in the x direction, the tie rods in x direction were activated more than the tie rods in y direction. An analogous situation was observed in the case of seismic action in the y direction. Compared to their tensile strength, the maximum stress level in the tie rods induced by the seismic action was 36.1% in the case of the cast steel tie rod and 21.31% in the case of the aluminum tie rod. Based on the results of conducted research, it can be concluded that the proposed method is a viable solution for determination of the load bearing capacity of tie rods in historical buildings. Since the updating of the tie rods local model and the global cathedral model was performed by manual tuning, as well as the linking of the experimental and numerical data sets, it is recommended to perform the automatic updating of the finite element model by using some kind of computer intelligence algorithm in future. Additionally, since a linear model was used in this study, it is recommended to perform a nonlinear analysis in the following studies. Due to the possibility that the force values in the tie rods may change due to settlement, the authors recommend that this influence also be considered in the following studies.