Seismic vulnerability of structures: application to the Civil Protection building in Andorra

. This works deals with the seismic vulnerability of buildings in the Pyrenees mountains region where almost a thousand earthquakes are recorded each year in the border area. The challenge is twofold: first to detect the damage due to seismic events and then to localize it inside studied buildings. Operational Modal Analysis (OMA) coupled with numerical modelling by Finite Element (FE) constitutes an interesting approach to address these issues. Here we intend to apply such methodology on a strategic building located in Andorre-la-Vieille whose structure is complex, irregular and heterogeneous. The structural behaviour of the building is studied through frequency computation method in order to identify its undamaged behaviour. A seismic event is next simulated by a non-linear dynamic computation method which creates damage within the structure. Numerical results (natural frequencies, modal shapes and damage location) allow highlighting damaged zones induced by the earthquake and quantify degradation level in these areas. Accordingly, some guidelines may be given in view of the future instrumentation of the building (accelerometers and RAR).


Introduction
Every year, a thousand earthquakes are recorded in Pyrenees, 90% of them in a strip about 50 km around the border [1].Since 2013, there have been about fifteen seismic events of magnitude greater than three in the Pyrenees each year [2].To protect populations that may be heavily affected in this area, it is important to be prepared for future seismic crises and especially to identify buildings with high vulnerability.This work specially deals with the case of strategic buildings (housing authorities, emergency services, hospitals, etc.) whose safety is even paramount for emergency management.
Assessing buildings structural integrity includes the two first Structural Health Monitoring (SHM) levels, namely damage detection and damage localization.It can be achieved by several kinds of non-destructive techniques such as visual inspection, stress-wave methods (acoustic emission, impact echo testing...), electro-magnetic methods (ground penetration radar, electrical resistivity measurements...), thermal methods (infrared thermography) or Vibration Based Damage Detection Methods (VBDDM, dynamic/vibration testing).Among them, the VBDDM are the most commonly used for measuring the structural integrity of structure [3].Determining if a structure is damaged or not can be achieved experimentally * e-mail: olivier.dalverny@enit.frusing VBDDM in which Operational Modal Analysis (OMA) helps for studying the changes of its dynamical characteristics [4,5].Indeed, the damage of a structure induced by an earthquake leads to the modification of its mechanical behaviour.In particular, a damaged structure shows a least rigidity and thus a higher damping.This results in a modification of the dynamic characteristics such as the reduction of natural frequencies and modification of modal shapes [6,7].Studying these changes and determining the state of a structure at a given time can be achieved through the processing of experimental data captured by a few of sensors.Defining the accurate location of a damage is however quite less trivial.Several OMA methods allow localizing a structural damage with various success.As example, methods linked to modal shape variations (mode shape curvature and curvature damage factor methods) are sensitive to measurement directions while flexibility variations based methods are sensitive to sensors location [8].Applying these methods in a specific order allow overcoming aforementioned drawbacks [9], but requires a huge amount of data in order to accurately determine damage location [10].Moreover, if the structure is not equipped with a continuous monitoring system, using such a methodology needs to make post earthquake measurements that can be dangerous to perform.
To deal with this issue, a relevant solution consists in coupling numerical modelling by Finite Element (FE) and OMA.Few studies are dedicated to the detailed representation and investigation of existing complex buildings.Generally, FE models are carried out to quickly identify natural frequencies and modal shapes.In order to account for seismic damage, several methodologies are employed, either considering artificial reduction of materials properties [11], applying horizontal forces considered as limits for the building's capacity (static pushover analysis [12,13]) or using accelerations recorded during seismic events as boundary conditions [14][15][16].The last option makes it possible to follow the progressive character of damage evolution, more representative of the dynamic load consequences.Most of these works deal with historical buildings that exhibit a quite homogeneous architectural design, namely of masonry type.Moreover, calibration of FE models is often supported by several experimental measurements or in situ pathologies observations [17][18][19].
This study intends to assess the structural integrity of the building housing the Civil Protection of the Principality of Andorra (Edifici administratiu del Prat del Rull).This building stands out by its structural heterogeneity regarding both the geometry of structural elements (floors, beams, columns) and the materials involved (mixed concrete-metal structure).Moreover, this study intends to provide decision help arguments before an experimental campaign.In this way, the building modelling requires accurate and realistic representation of the mechanical behaviour of different parts so as to guide future building instrumentation.After the description of the geometry and complex structure of this existing structure (Sect.2), the numerical assumptions used to model the building and simulate its dynamic behaviour are described (Sect.3).The two computations methods (frequency and dynamic) are firstly introduced following by the material description and the generation of artificial earthquake representative of the Pyrenees area.Results obtained are finally discussed and guidelines are given in view of the future instrumentation of the building (Sect.4).

Building description
The studied building is located in Andorre-la-Vieille (Edifici administratiu del Prat del Rull del govern d'Andorra, Fig. 1a and b).The structure built in 2013 is mainly constituted of reinforced concrete and exhibits two underground floors and six upper floors leading to a total height of 20.5 m.The basement (approximately 37 m in width and 41 m in length, Fig. 1c) is supported by a 1-m thick concrete screed and the two underground floors are sustained by concrete walls, beams and columns.Regarding the upper floors, slabs are composed of load-bearing Haircol R material and are sustained by metal beams.The link between the first level (over the ground) and the top of the building is mainly done by metal columns.Finally, a lift surrounded by concrete walls and two metal stairways allows to move within the building.In what follows, all these features, including the different elements types (floors, beams, columns), geometry and materials (metal, concrete), have been accounted in the FE model

Computation methods
Two kind of computation methods implemented in the finite element software Abaq-us R are used in the following sections: a frequency computation method and a dynamic computation method.The first method corresponds to a linear perturbation procedure and allows the extraction of natural frequencies and modal shapes of the structure at a given state [9].Based on Lanczos algorithm, such approach includes initial stress and load stiffness effects due to preloads.This algorithm is widely appreciated in the numerical analysis community as a very powerful tool for extracting some of the extreme eigenvalues of a real symmetric matrix, i.e. to find the largest and/smallest eigen values and vectors of the symmetric eigenvalue problem [20][21][22].It is thus necessary to use a static non-linear computation method in order to integrate the proper weight of the structure before running a frequency computation study.The dynamic method provides a non-linear analysis of the structure and aims to simulate seismic events and induced damage [15,16].Since such procedure uses an explicit solver, the duration of the computation step is linked to the size of the smaller mesh element.Special attention must therefore be paid to mesh elaboration so as to minimize computation time.

Building modelling and numerical assumptions
Concrete elements (screed, slabs, walls, beams, columns) and metal stairways are modelled by shell elements while wire elements are used for metal elements (beams, columns and window frames).These modelling choices are motivated by the geometry of the structural elements (that exhibit higher width and length than their respective thickness) and in view of saving computation time.Indeed, due to the large dimension of the building, using 3D solid elements for all parts would involve a very long computation time.Respective thickness and cross section of shell and beam elements are then taken into account by defining section properties.
Three type of materials law can be used to model the complex and non-linear behaviour of concrete: the concrete smeared cracking model [23], the brittle cracking model [24], and the Concrete Damage Plasticity (CDP) model [25].All of these can simulate damage evolution within a concrete structure, but only the concrete smeared craking model allows to account for crack propagation.In addition, only the two last models can be used with a dynamic loading case [26].The CDP model was used in this study as it is well suited for seismic load.The CDP model is an elasto-plastic damageable model which can continuously represent the degradation of material according to the mechanical loading.This model is governed both in tension and compression by the following equation: where σ i , ε i , and ε pl i , denote, respectively, the stress, total strain, and plastic strain; E 0 is the initial elastic stiffness; d is the damage factor, characterizing the degradation of the elastic properties with values in the range from 0 (undamaged) to 1 (fully damaged); i denotes the type of simulated behaviour (c for crushing in compression and t for tension).Haircol R slab is composed of a concrete layer showing different thicknesses which is sustained by a metallic plate (Fig. 2).In view of further FE calculations, we intend to represent such complex structure by means of a simplified equivalent material, defined by its thickness and weight.To achieve this, an homogenization approach has been implemented upon the analysis of an elementary structure composed of two floors sustained by wire elements (simulating metal columns) embedded to the ground at their bases.Each floor is composed of four slabs linked together by tie interactions.Homogenized properties were then identified consequently to fit the natural frequencies and modal shapes of the 3D Haircol R model.The quality of the homogenization was thus validated by comparing the values of natural frequencies and using the Modal Assurance Criterion (MAC).The MAC criterion allows to check the correlation between two modal shapes using the following equation [29]:  and 2. A maximum relative difference of 6.66% and 5.21% are obtained for the first flexion and torsion mode, which seems acceptable.Moreover, the MAC criterion shows a really good correlation for the two first flexion and torsion modes while there is no correlation for the second flexion mode as the modal shapes are not in the same direction.Homogenizing Haircol R slab seems thus a relevant choice to provide a more easy-handling representation of the building.
Shell elements are then merged together and tie constraints (equivalent to the coupling of degrees of freedom) are finally applied between shells and wire elements to structurally link them together.The final geometry of the numerical model is presented in Figure 4.
The structure is finally meshed using 3D linear beams (B31), linear triangular and quadrilateral reduced integration shells, respectively (S3R) and (S4R) elements (Fig. 5).A particular care is taken on mesh.Indeed, for dynamic computation methods, the duration of computation steps is closely linked to the size of the smaller element (L) and the speed of wave propagation (C 0 ) within material (δt = L C0 ).A convergence study has thus been conducted to ensure the results stability according to the mesh size leading to an average element size of 0.5 m.Precisely, the average aspect ratio (δ) is below five for the whole model which allows to minimize mesh distortion issues during the computation process.Details  of the meshing characteristics (number and type of used elements and distortion criteria) are presented in Table 3. Gravity and non-structural mass are then defined in order to integrate the proper weight of the structure and the initial stress of the structure.Gravity is applied on the whole model while non-structural masses are applied on each upper floor (45 kg • m −2 : representing an approximation of the weight office supplies) and the first underground floor (125 kg • m −2 : representing an approximation of the weight of parked vehicles and operational weight).
Regarding boundary conditions, conventional seismic analysis generally does not take into account the flexibility of the foundation and adjacent soil [13].As a consequence, the evaluated seismic performance may be significantly different from that of the actual building [30].In the present case, data on soil-structure interactions were not available.Yet, in view of the geological context (Pyrenees mountain, location altitude: 1023 m), soil type is considered as generic rock and soil-structures interactions can be neglected [31].For frequencies computation, the outer boundary of the bottom part of the building is assumed to be embedded in the ground, i.e. all horizontal or vertical walls in contact with the ground are considered fixed.Note that the inside walls of the underground are not constrained.In the case of dynamic analysis, the embedding of the base is replaced by base motion accelerations along directions X, Y and Z. Precisely, the software Seismo Artif 2016 [32] was used to generate synthetic accelerograms representative of seismic activity in the Poctefa region (French-Spanish-Andorran border).Pyrenees mountains are located within an interplate region [33], and the maximum recorded magnitude was 6.5 [34].In this study, the duration of the earthquake was set to 5 s during maximum  acceleration peaks acceleration can be observed.Figure 6 shows the three accelerograms generated by Seismo Artif 2016.

Dynamic characteristics of the undamaged structure
Frequency computation is applied to the FE model presented below in order to study natural frequencies and their associated modal shapes.The first 200 modal shapes have been evaluated in the aim to ensure that the Cumulative Effective Mass Participation Factor (CEMPF) corresponds to at least 80% of the total mass of the building [13,35].Indeed, several local modes appear due to the presence of numerous metal columns leading to a value of approximately 9.5 Hz for the 200th modal shape.The most important Effective Mass Participation Factors (EMPF) are obtained with the first modes of the structure (mode no. 1, no. 2 and no.4) for which the EMPF are higher than 7% (Fig. 7 and Tab. 4).EMPF factors decrease then drastically with higher modes and numerous local and mixed modal shapes appear.Though the first modes induce a large movement of the upper part of the structure, the CMEPF reaches only a maximum of 23% of the total mass of the building upon X direction (Fig. 8a) and 22% upon Y direction (Fig. 8b).This result differs from classic results found in the literature [13], but can be explained by the singular constitution of the building.Indeed, the upper part of the building is mainly carried by metal beams and columns while the bottom part is entirely constituted of concrete.Moreover, the bottom part that represents 66% of the total weight is underground and thus highly constrained by boundary conditions.It results thus two distinct mechanical behaviours characterized by a higher flexibility of the upper part than those of the bottom part.Consequently, the upper part is subjected to large displacements while the almost totality of the bottom part (underground floors) does not move.If the embedding of the underground storey is taken into account in calculating the CEMPF, this one reaches 80% of the total mass from the mode no. 4 (Fig. 8, Tab.4).It is thus relevant to consider that the dynamical behaviour of the structure is well described by the numerical model.Modal analysis is a powerful tool to identify the natural frequency behaviour of buildings, and specially to understand how the structure may be sensitive to a particular direction of the mechanical excitation.The analysis of modal shapes thus makes it possible to characterize the resonance phenomenon that amplifies the vibratory response of the structure for a given stress.In the case of the studied structure, we can observe that the main modes with respect to the EMPF factor, namely 1, 2 and 4 (Tab.4), correspond to torsion and/or bending displacements (Fig. 7).The easiest mode to describe is mode 2, which corresponds to bending in the X-direction for the whole construction.Mode 1 and 4 are also global modes   (involving the motion of whole structure), but combining bending and slight torsional effects.Indeed, Mode 1 corresponds to a large bending of the northern part of the building (on the right side of Fig. 7), with displacement in the Y-direction.Mode 4 is also related to a large bending motion in the Y-direction, with movement of the southern part of the structure (on the left side of Fig. 7).For both modes, bending movements induce simultaneously torsional effects on the respective remaining part of the building.The other three modal shapes shown in Figure 7, namely for modes 75, 86 and 171 respectively, exhibit a low EMPF factor and correspond to minor modal shapes.Mode 75 shows a second-order bending mode on the northern part of the building, mode 86 is a global second-order bending mode and mode 171 is related to a local bending mode (Z-direction) on the floors of the southern part of the building.

Damage process
The introduction of damage within the structure is achieved by the simulation of a 6.5 magnitude earthquake using a non-linear analysis.During a short period of time (about half a second) corresponding to the establishment of the seismic amplitude, the kinetic energy of whole model increases slightly (Fig. 9a) and stresses are mainly concentrated in metal beams around windows.At this time, all structural elements remain in the elastic domain (as example: σ metal beams ≈ 26 MPa < σ pl ≈ 450 MPa) and no damage is detected (Figs. 10 and 11).Consequently, no dissipative energy is detected within the model, neither damage energy nor plastic dissipation (Fig. 9b).Beyond this point, the kinetic energy still fluctuates.The initiation of plasticity and damage, namely at around 0.55 s, is specially correlated with a first peak in the amplitude of the acceleration given by the Euclidean norm of the acceleration vector (Fig. 9b).From this, one note an important increase in the plastic dissipation induced in both metal and concrete parts.Precisely, the plastic yield strain is reached in several structural metallic parts (beams and columns) around window frames, as well as in the concrete elements close to these areas.Damage in concrete evolves simultaneously, due to the coupling with plasticity considered in the CDP model.Compared to plastic dissipation, damage energy increases slower as damaged areas are less extended than plastic ones (Figs.9b, 10 and 11).Moreover, it should be underlined that the strong variations in both plastic and dissipations evolutions are correlated, not with high amplitudes of acceleration but with peaks of kinetic energy (see for instance at time 0.75 s, Fig. 9c).Just before 1 s, the seismic acceleration exceeds 1.2 m • s −2 (in absolute value) along East-West (Y) direction (Fig. 9a).Considering the average value over the first four modes (Tab.4), this direction corresponds to the highest mean Effective Mass Participation Factor, that is the most sensitive excitation direction of the structure.Such a strong acceleration applied in the weakest direction of the building significantly increases the stresses and plastic strains in metallic columns on the East side of the building and consequently the damage in concrete slabs attached to them.This leads to an important loss of concrete stiffness related to crack propagation in these areas (98% damage level, see in Figs. 10 and 11).From this point, the computation process is stopped due to the large stresses within metal elements that induce element distortion issues.

Optimization of sensors positioning
In view of the computational analysis, it is possible to optimize the location of sensors.On the one hand, the preponderant modal shapes extracted from the frequency analysis show a large sensibility of the building in regards with torsion and flexion modes, but also a weak motion of the underground floors due to the embedding of this part of the structure within the ground (Fig. 7).On the other hand, results obtained from the dynamical analysis show that stresses and strains are concentrated near the zones where windows are located.Setting-up multi-directional sensors (accelerometers) near windows that is to say near the corners of the building seems thus a relevant solution to capture both torsion and flexion modal shapes, but also accelerations induced by an earthquake.Using sensors in the underground levels of the building is however not necessary because of this weak motion.Indeed, accelerations would be hard to detect in these zones.Recommendations for sensors positioning are presented in Figure 12.

Conclusions and perspectives
This paper has assessed the seismic vulnerability of a complex structure: The Edifici administratiu del Prat del Rull del govern d'Andorra.Addressing such kind of strategic buildings is of crucial interest for seismic crisis management.As for most existing building, experimental data were not available to provide any information on the vibratory response of this mixed concrete-metal structure.Numerical modelling by finite elements appears thus a relevant tool to understand and characterize the dynamic behaviour of the structure.This helps to estimate the damage phenomena that may affect the building during a seismic event and also to provide guiding elements for its future instrumentation.
First, dynamic characteristics of the undamaged structure (natural frequencies and modal shapes) were identified by linear perturbation procedure.Then, a non-linear computation procedure was used to simulate an earthquake representative of the local seismic activity in Pyrenees.
This work focused on the structural behaviour of the structure.It has allowed to identify the most stressed zones of the building and understand the damaging process during the seismic event.We have shown in this study that the analysis of this complex structure was not trivial.It requires the simultaneous monitoring of the solicitation acceleration curve and the kinetic energy evolution curve.While the onset of damage is correlated to an acceleration peak, the evolution transitions are rather correlated to kinetic energy peaks.Moreover, due to the building complexity, specially the heterogeneity of the structural elements (geometry and materials), stress localization and dissipative phenomena (plastic and damage) kinetics can hardly be estimated without an accurate modelling of the structure.Using these results, optimization of the location of sensors (accelerometers) within the structure can now be done to achieve its continuous monitoring.
To go further, the numerical model could be calibrated using experimental data acquired by mean of accelerometers and real aperture radar.Indeed, such approach seems necessary to calibrate the soil/structure interactions and better represent the dynamic behaviour of the building.

Fig. 3 .
Fig. 3. Numerical results of the homogenization procedure: (a, b): first flexion mode, (c, d): first torsion mode with (a) and (c) corresponding to 3D model and (b) and (d) to homogenized model.
2) where ψ lb = [ψ lb1 , ψ lb2 . . .ψ lbj ] and ψ hs = [ψ hs1 , ψ hs2 . . .ψ hsk ] denote respectively the modal shape of the load-bearing floor and of the homogenized structure; n is the number of measured nodes.Diagonal values of the MAC j,k matrix represent the correlation between same modes (j = k), a value of 1 indicates that the modal shapes are exactly the same while a value of 0 denotes the total absence of correlation.Modal shapes obtained for 3D and homogenized models (for a thickness of 0.092 m a Young modulus of E = 15 GPa and a density of ρ = 2018 kg • m −3 ) are compared on Figure 3. Natural frequencies and MAC criterion are given in Tables 1

Fig. 8 .
Fig. 8. Modal analysis results: (a) CEMPF upon X direction, (b) CEMPF upon Y direction (note that most of modal mass is related to the very rigid subsurface part of the building).

Fig. 9 .
Fig. 9. Accelerations data along X, Y, and Z directions and evolution of energies within the model in function of time: (a) kinetic energy and acceleration components, (b) acceleration norm and dissipation energies (c) kinetic and dissipation energies.

Fig. 10 .
Fig. 10.Simulation of damage process following the North and East directions.

Fig. 11 .
Fig. 11.Simulation of damage process following the North and West directions.

Table 1 .
Comparison between the natural frequencies of the 3D and homogenized models.

Table 2 .
Modal assurance criterion (MAC j,k ) of the homogenized Haircol R slab.

Table 4 .
Effective contributions upon directions X, Y and Z for the most important modes.