Inelastic Analysis of Mdof Systems Damaged by Earthquakes, Posteriorly Subjected to Wind Load

This paper deals with the analysis of the inelastic response of buildings originally damaged by earthquakes and subjected to earthquake aftershock and wind loading. The overall aim is to establish the effect of wind actions on structural stability. To that end, one four-story bare frame benchmarked by the European Laboratory for Structural Assessment, is subject to various levels of winds and earthquake joint load while monitoring changes on the ductility demand. In this paper is shown that the combined action of strong winds and earthquakes, however its low probability of occurrence, would cause a decrease of strength reduction factors and considerably increase the ductility demand of damaged infrastructure hence inducing additional risks that would otherwise remain unquantified. The paper examines the nonlinear performance of Multi-degree of freedom systems subject to various levels of winds and earthquake load and deals with the estimation of strength reduction factors. This is a relatively unexplored area of research which builds on past developments whereby inelastic performance of buildings has been discussed. It also links to various other paths of development such as structural reliability, forensic and control systems engineering.


Introduction
Most structures designed according to current code provisions will sustain damage in the event of a design-level earthquake occurrence, which to a good extent, is expected. It is well known that structural damage is directly related to ductility demands [1] and therefore evaluation of their relationship is very important for accurately assessing structural performance. Past research shows that current design standards for seismic design and seismic assessment, can improve significantly through the explicit account of ductility demands [2].
The evaluation of ductility demands for structures is abundantly reported in past research, see for example Yi et al. (2007) study [2]. According to these, the extreme design values of ductility relate to rare events. Earthquake design of structures generally allows for inelastic deformation depending on the target performance level. Proper structural performance thus relies on inelastic deformations that correspond to expected ductility or drift levels.
The structural demand due to an earthquake ground motion can be approximated by performing nonlinear time history analyses with computational models subject to a given ground motion. However, it seems appropriate not only consider isolated 'design earthquakes' but take one step forward to account for the inelastic response of buildings originally damaged by earthquakes that are subject to moderate or strong aftershocks. In such scenarios, if not during the original design, once can assume that wind is flowing. This depict multi-hazard scenarios not been studied in the past. To the authors' knowledge, only a few studies have examined the aftershocks effects on buildings while only a few revised combined earthquake and wind joint load. Moreover, to date some studies still report on simplified modelling techniques that omit integrating hysteretic beam or column behavior. Hatzigeorgiou and Beskos (2009) and Hatzigeorgiou (2010) [3,4] focused on the response of Single-of-Degree-Of-Freedom (SDOF) systems assuming bilinear elastoplastic hysteresis with no stiffness and strength degradation. MDOF systems can more efficiently integrate beams, columns, and beam-column joints that can be deform in different manners due to cyclic loading reversals propagated by earthquake ground motion. It is to note that, global ductility depends on inter-story drift at higher rates that it does on roof drift [5]. Hysteretic models with stiffness degradation better simulate peak displacement of short period of structures than those assuming non-degrading hysteresis behavior [6]. Chopra and Chintanapakdee (2001) [7] studied far-fault and near-fault ground motions to determine qualitative differences on structural performance. That study showed that near-fault motions tend to induce greater displacement demand than far-fault motions in the constant acceleration region. This is apparently due to the fact that for near-fault ground motions, the average corner period (the period at which the constant velocity region begins) of the response spectra shifts to longer periods. The constant velocity region for near-fault motions is also generally narrower than that of far-fault motions.
All these studies assume, however, that earthquake and wind events are uncorrelated. This assumption seems reasonable when accepting that the structure's lifetime assumed for either design condition does limit their joint probability of occurrence. However little to no research has been produced for cases in which major earthquake events clash with low to moderate winds whose probability of occurrence is significant or when notwithstanding its low probability of occurrence, major seismic and wind events coincide. The latter does not seem unlikely after all if we consider the number of aftershocks that usually follow major earthquake events. For example, the earthquake that hit Nepal in 2015 (Ms =7.8) killing more than 8,000 [8] was followed by 30 aftershocks of Ms < 5 occurring within three weeks and killing 200 more. The earthquake that hit the Sichuan Province in China in 2008 (Ms =7.9) killing over 87,000 [9] was followed by 12 weeks with 42 aftershocks ranging in magnitude between 5 < Ms < 6.4. If we add this to the fact that wind is constantly flowing with a minimum speed that is equal to the average local wind velocity, then the current design assumptions concerning the unrelated action of wind and earthquakes do not seem conservative [10]. It thus appears that new research avenues generate through consideration of near-or far-fault earthquake events followed by moderate to strong aftershocks that occur simultaneously to wind actions .This paper attempts to address this knowledge gap through the estimation of ductility demands (µ) of structures subject to the combined effect of earthquake and wind load. The latter, flowing at mean speeds ranging between 5 ms -1 and 50 ms -1 . The subject of the investigation is the four-story bare frame standardized by the European Laboratory for Structural Assessment (ELSA) [11,12]. This building is thus subject to multiple near-and far-fault seismic ground motions. The structural response is assessed by nonlinear dynamic (time history) analysis, to estimate the overall response. This paper is organized mainly in different sections. The introduction portion is followed by an overview of the ductility demand, the structure configuration and material properties. Then the seismic input, wind filed simulation, wind and earthquake loading, which is followed by the results of the combined effect of earthquake and wind load. The latter, flowing at mean speeds ranging between 5 ms -1 and 50 ms -1 is conducted to obtain expressions for the ductility demands, in terms of the force reduction factor, and finally, conclusions and recommendations are provided.

Ductility Demand: An Overview
An adequate design is accomplished when a structure is dimensioned and detailed in such a way that the local (story and member) ductility demands are smaller than their corresponding capacities. Thus, during the preliminary design of a structure, there is a need to estimate the lateral strength (lateral load capacity) of the structure that is required in order to limit the global (structure displacement) ductility demand to a certain pre-determined value, which results in adequate prediction of localized damage.
The inelastic deformation experienced by a system under a given ground motion can be determined through the displacement ductility ratio, µ, which is defined as the ratio of maximum absolute relative displacement to its yield displacement this definition is given by Miranda (1993) [13].
The governing equation for SDOF inelastic system subjected to horizontal ground acceleration ü ( ) is For a given excitation ̈( ), the deformation ( ) depends on three systems parameters: , ζ, and , in addition to the form of the force-deformation relation. This becomes evident if Equation 2 is divided by m to obtain: Where is the natural vibration frequency, = 2 ⁄ is the natural vibration period, and ζ is the damping ratio of the system vibrating within its linear elastic range (i.e., ≤ ). Equation 3 is rewritten in terms of ( ) = ( )⁄ in order to identify the parameters that influence the ductility factor µ, Equation 1, the peak value of ( ). Substituting ( ) = ( ), ̇( ) = ( ), ̈( ) = ( ) in Equation 3 and dividing by gives: Where = ⁄ may be interpreted as the acceleration of the mass required to produce the yield force .
For given ̈( ), µ depends on three system parameters: , , [Chopra, 2001]. The maximum force response of a linear elastic system can be denoted by , while the yield strength of a nonlinear elasto-plastic system can be denoted by . Thus, the force reduction R factor can be defined as  This method has been widely used to determine inelastic response spectra considering at least three non-linear load-deformation models: elastoplastic, bilinear, and stiffness degradingsee for example Miranda and Bertero (1994), Chopra and Goel (1999), and Riddell et al. (2002) [14][15][16]. In the present research the elastoplastic model is adopted.

Structural Configuration and Material Characteristics
The present study adopts the four-story reinforced concrete prototype shown in Figure 2. The concrete frame has been tested at the European Laboratory for Structural assessment (ELSA) [11,12] under two subsequent pseudodynamic loadings. The two reinforced concrete frames tested at (ELSA) laboratory can be considered representative of the design and construction common practice until the late 1970's in southern European countries such as Italy, Portugal and Greece. They were designed to primarily withstand vertical loads such that its lateral resistance is of 5% of the frame weight [12]. The reinforcement details were specified in accordance to the normative available and to the construction practice at that time. Thus, no specific seismic detailing provisions were considered, preferential inelastic dissipation mechanisms were not assumed and no specific ductility or strength provisions were provided.

Figure 2. Four-story, three-bay RC frame geometry (m), (elevation and plan views [12]
The vertical loads, represented in Figure 3, were defined in order to simulate the dead load other than the selfweight of the frame, considering that parallel frames have a distance of 5.0 m Carvalho et al. (1999) [12]. The frame model includes a 4.0 m wide slab, which requires additional vertical load accounting for such a slab portion missing. Vertical distributed loads on beams and concentrated loads on the column were considered in order to simulate the dead load of the frame other than the weight of partitions, finishings and live load The distribution of vertical loads applied to the infilled frame to simulate the dead load other than the self-weight of the frame was identical to the one used for the bare frame, and imposing the same load distribution for all the tests. The scheme of the loads considered can be found in Pinto et al. (1999-c).  Mander et al. (1988) [32] have proposed a unified stress-strain approach for confined concrete applicable to both circular and rectangular shaped transverse reinforcement. The stress-strain model is illustrated in Figure 4 and ′ is based on an equation suggested by Popovic (1973) [33]. For a slow (quasi-static) strain rate and monotonic loading, the longitudinal compressive concrete stress is given by: Where ′ is the compressive strength of confined concrete and x is a ratio of longitudinal compressive concrete strain ( ), r is the ratio of the concrete's initial modulus to the difference of the initial and secant moduli of elasticity. These parameters and their components are mathematically expressed by: As suggested by Richart et al. (1928) [34], where ′ and 0 = the unconfined concrete strength and corresponding strain, respectively (generally 0 = 0.002 can be assumed), and; Is the tangent modulus of elasticity of the concrete (1 MPa = 145 psi), and The characteristic parameters are listed below: = 16300 kPa, = 1900 kPa, = 0.002 m/m, = 18975 MPa

Steel Model
The cyclic behavior of the steel bars is simulated using the classical Menegotto and Pinto (1973) [35] model with kinematic hardening ( Figure 5) and possibility to take into account buckling by introducing a negative modulus slope in compression depending on the transverse steel reinforcement spacing. The monotonic behavior is defined through the initial Young's modulus ( ), the plastic threshold ( , ), the ultimate strength and strain ( , ) and the yielding slope ( ℎ ). The unloading and reloading process, is guided by analytic relations (Equations 12 to 15) corresponding to a set of curves ranging between the elastic and the yielding asymptotes.
Tension: i = 1 and j = 2; Compression: i = 3 and j = 4.( , ) is the studied point; ( 0 , 0 ) is the crossing point of the elastic and yielding slopes; ( , ) are the coordinates of the previous point of load reversion; b is the ℎ / ratio; R is a shape parameter; is the ratio between the maximum reached strain during loading; , 0 , , 0 , , and , are material constants that can be obtained from experimental results.
The Menegotto-Pinto steel model is employed for defining the steel material with the following properties: = 2.00E + 008 kPa, = 343000 kPa, = 0.0024.
The materials considered at the design phase were a low strength concrete class C16/20 (CEN, 1991) and smooth reinforcement steel class Fe B22 k (Italian standards). The latter refers to smooth bars with a yield stress of 235 MPa and ultimate strength of 365 MPa.

Earthquake Record Database
The seismic excitation consists of near-fault and far-fault earthquakes. In this study, the near-fault accelerograms database that it has been used here consists of 80 recorded ground-motion time histories from different fault types and damping ratio of 5% and earthquake magnitudes 5 < Ms < 8 that have occurred in the United States, Canada, Japan, Taiwan and Turkey. These accelerograms present peak ground acceleration (PGA) greater than or equal to 0.10g. These strong ground motions were recorded at stations where the closest to fault rupture is less than or equal to 20km. The complete list of these earthquakes, which were downloaded from the strong motion database of the Pacific Earthquake Engineering Research (PEER) [17] Center, is shown in Table 1 and Figure 6.
Near-fault seismic ground motions are frequently characterized by intense velocity and displacement pulses of relatively long period that clearly distinguish them from typical far-fault ground motions. Generally, forward directivity and permanent translation are the two main causes for the velocity pulses observed in near-fault regions [18].
A total of 80 real far-fault earthquake acceleration time histories from around the world are also used in this study. These accelerograms present maximum ground acceleration greater than or equal to 0.10g and are recorded at sites ranging from hard rock to soft soil conditions. The complete list of these earthquakes, which were also downloaded from the strong motion database of the PEER Center [17] is shown in Table 2 and Figure 7. On the other hand, Gillie et al. [19], Bray and Rodriguez-Marek [20], and Chopra and Chintanapakdee (2001) [7], suggest that Near-fault earthquakes (e.g. ED < 20 km) could have more devastating effects on structures that far-fault ones due to energy contents stored in relatively long-term pulses embedded in the ground motion.

Simulated Wind Field
In the present study, the wind field simulation is based on the algorithm suggested by Vanmarcke et al. (1993) [21] which consists of inferring correlated data series in a number of target points across a region bounded by a number of stations in which the characteristics of the random signal are known. This process is based on the best estimator of Fourier coefficients that can be used to reconstruct power spectra at the target locations which are consistent with realistic correlation laws such as those described by Dyrbye and Hansen (1997) and Simiu and Scanlan (1978) [22,23]. The effectiveness of the method has been reported by a number of authors in the past, see for instance Gurley et al. (1997) [24], Martinez-Vazquez and Rodriguez-Cuevas (2007) [25], Martinez-Vazquez et al. (2010) [26]. Therefore, a brief explanation of the method together with some statistics of simulation results is provided below.
The simulation consisted of generating two uncorrelated time series that represent the wind regime of a suburb at heights of 10 and 250 m above the ground, using Monte Carlo techniques. These time series were taken as recorded data in the algorithm discussed in Vanmarcke et al. (1993) [21] from where partially correlated series at intermediate points could be inferred. The target turbulence intensity at z = 10 m was 0.295 and ~0.2 for the structures, respectively. The former being a standard value for suburban area whilst the latter was the value used in the experimental work reported in Melbourne (1980) and in Thordal et al. (2020) [27]. Table 1 shows the statistics of the simulated series located along the height of 250 for when the wind velocity at 10 m above the ground (U10) is 33 ms -1 .
In Table 3, U and 2 represent mean and variance of the series in ms -1 and m 2 s -2 whereas the subscripts represent target (t) and simulated (s) values, respectively. The mean square error between target and simulated values are of 4.7% for the mean velocity and 0. 9% for the variance. The target and simulated cross-correlation amongst time series is shown in Tables 4 and 5.
The mean square error of the correlation parameter is of 1.135% across all stations. This was considered acceptable and representative of wind events in suburban areas. A comparison between the simulated and theoretical wind power spectrum in the low frequency range is shown in Figure 8 for the position z = 10 m above the ground level.
The simulated power spectrum represents a single realisation and thus the spectral density fluctuates around the Von Karman spectrum. In Martinez-Vazquez and Rodriguez-Cuevas (2007) [25] and Martinez-Vazquez et al. (2010) [26], it is demonstrated that, by increasing the number of realisations, the ensemble of simulated spectral ordinates tends to theoretical values. For the purpose of the present study, the single realisation of the synthetic fields is considered to be appropriate [10].    Figure 8. Simulated and theoretical wind power spectrum at 10 m above the ground

Dynamic Response
The static forces induced by wind were calculated using Equation 16. The reference height (zr) was taken as of 10m.
Where represents the density of air, is a drag coefficient and A is the area exposed to wind. Whilst a value of = 0.22 was used to represent the wind profile in suburban area.

Wind and Earthquake Loading
There is a fundamental difference between wind and earthquake load acting on large areas. The horizontal accelerations induced by earthquakes along the height of structures are assumed to be fully correlated by Chopra whereas the corresponding wind forces are not Dyrbye and Hansen (1997) [22]. The correlation of wind forces depends on the spatial distribution of wind gusts and that is usually taken into account through suitable correlation laws such as those proposed by Vickery (1970) [28] and Tanaka and Lawen (1986) [36] together with admittance functions such as that proposed by Davenport (1967) [29]. For point-like structure the input load can simply be obtained by superimposing both actions thus assuming these acts simultaneously whilst randomly off-phase [10].

Results and Discussion
After designing and detailing the reinforced concrete frame structure, the inelastic response of structure subject to earthquake and wind is performed by the nonlinear analysis software SeismoStruct 2016 [31]. Time history analysis is used to determine the maximum displacement response, and the Pushover analysis to obtain yield displacement of structure.
The pushover analysis consists of the application of gravity loads and a representative lateral load pattern. The lateral loads were applied monotonically in a step-by-step nonlinear static analysis. The applied lateral loads were accelerations in the x direction representing the forces that would be experienced by the structures when subjected to ground shaking. Under incrementally increasing loads some elements may yield sequentially. Consequently, at each event, the structures experience a stiffness change as shown in Figure 9, where IO, LS and CP stand for immediate occupancy, life safety and collapse prevention respectively. An elasto-plastic model for the above curves can be considered, so that the effect of this model would be equal to the effect of the real curve. In this method, drawing the complete elasto-plastic curve is based on the equalling of the structure energy absorption. It means that the area under the curve that indicates the amount of energy exerted to the structure should be equal in the elasto-plastic and the actual case. The yield point of the roof story can be obtained from these curves ( Figure 9).

Displacement, Base Shear and Ductility Demand-µ
In Figure 10a, the time histories (experimental results and numerical calculations) of top-displacement are plotted for Near-Fault input motion performed on the bare frame structure. Figure. 10b shows the top-displacement curves for Far-Fault earthquake, the comparisons between the experimental results and numerical results. The results are in good agreement for the first few second (< 7 s). Figure 11 presents the curves base-shear versus the period of the combined effect of earthquakes occurring simultaneously with the wind flowing.    Figure 12 presents the average of ductility demand-µ for both the Near-fault and Far-fault sets. Looking at the results, it is concluded that the ductility-µ increases with wind velocity and therefore the influence of wind load on the inelastic response of structures increases with respect to that due to the sole action of earthquake load. At low wind velocity values (e.i velocity < 5 m/s) the ductility demand for Near-fault and Far-fault motions are approximate equivalent at zero wind speed. As the velocity of wind increases the ductility demand increases for short and intermediate wind speed. For the average of strength reduction factor show in Figure 13. Observe that a rapid increase in the interval 0 < T < 0.5 followed by a change of direction into a region which tends to be flat until ~ T = 3 s. The strength reduction factor depends not only on the characteristics of system, but also on the ground motion input. For given ground motion, is a function of the period of vibration T of structure, the damping the type of hysteretic behavior and the level of inelastic deformation studies reviewed in the previous section agree that for a given acceleration time history the strength reduction factor is primarily influenced by the period of vibration. Furthermore, Figure 13 show tow linear segments. In the first segment which corresponds to the short-period region, Increases linearly with increasing period from = 1 to value that it is equal to ductility factor. In the second segment, the strength reduction factor maintains a constant value. The de details of the proposed relations depend on the hysteretic behavior and damping of system. This result coincides with comparisons of results in literature [29].

Conclusion
This study presents ductility demands for Multi Degree Of Freedom (MDOF) systems under multiple near-and far-fault seismic ground motions. The main innovation of this work has to do with the influence of near-or far-fault earthquake commonly followed by a number of moderate to strong aftershocks and that wind is constantly flowing on ductility demands, a phenomenon which has not been studied in the past. This paper examines only and deals with the estimation of ductility demand (µ) and the strength reduction factor covering the combined effect of earthquakes occurring simultaneously with the wind flowing at mean speeds ranging between 5 and 50 ms -1 . This investigation shows the impact of considering wind load in the estimation of ductility demands. It is shown that ductility-µ increases with wind velocity and therefore the influence of wind load on the inelastic response of structures increases with respect to that due to the sole action of earthquake load. The estimated strength and ductility demand of wind simulated indicate that under certain conditions these can undergo unforeseen inelastic performance during extreme events. Furthermore, for given target ductility, the strength reduction factor can be exhibit great variations from one ground motion to another. For the design of a structure this means that the lateral strength capacity required to avoid displacement ductility demands larger than a given limit can have important variations from one ground motion to another. Finally, it is found that wind and earthquake combined effects strongly affects the ductility demand. This adds further arguments against current design assumptions which ignore multi-hazard scenarios. It seems therefore necessary to consider the risk to extreme events in the formulation of a more robust design framework for practical use.

Data Availability Statement
The data presented in this study are available in article.

Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.

Acknowledgements
The present work has been developed at the school of Engineering of the University of Birmingham. I gratefully acknowledge the generous help and facilities of Birmingham University and many individuals, for their generous support. The Authors also acknowledges the reviewer suggestions made during the review phase of the publication.

Conflicts of Interest
The authors declare no conflict of interest.