Study on structural response under vibration load with multiple degrees of freedom

Deyu Kong

Home /

Archives /

Civil and Energy Research, Volume 1, Issue 1, 2025 /

Original Research

Study on structural response under vibration load with multiple degrees of freedom

Deyu Kong ¹ ^*

1. University of Leeds, School of Civil Engineering, LS2 9JT, UK

Correspondence: alexkong0825@outlook.com

Received

26 December 2024
Revised

28 February 2025
Accepted

18 April 2025
Published

23 November 2025

Multi-degree-of-freedom (MDOF) structures Stochastic Incremental Dynamic Analysis (IDA) Seismic excitation Monte Carlo simulation Peak Ground Acceleration (PGA) Maximum interstory drift ratio Nonlinear time-history analysis

Abstract

In recent years, the frequency of seismic activity has increased, highlighting the importance of evaluating the response of structures to varying seismic intensities to ensure their safety. This study explores the stochastic incremental dynamic analysis (IDA) of multi-degree-of-freedom (MDOF) structures under unsteady seismic excitation. To address the complexity of unsteady seismic events, stochastic precesses are integrated into the IDA to generate stochastic seismic records that adhere to the elastic response spectra specified in Eurocode 8. A Monte Carlo simulation approach was employed to generate seismic waves using MATLAB, followed by a nonlinear time-history analysis in ABAQUS to evaluate the structural response, particularly focusing on the maximum interstory drift ratio. The findings indicate that (1) at low Peak Ground Acceleration (PGA) levels, the structure exhibits a minimal risk of failure. However, the risk of structural failure escalates significantly as the PGA increases, particularly beyond 0.4g. (2) The study also identifies gaps in current seismic analysis practices, especially the need for more robust stochastic IDA applications and the consideration of non-smooth excitations. This research offers a more comprehensive understanding of the seismic performance of MDOF structures and provides valuable insights for enhancing seismic design and risk assessment. Nevertheless, the study acknowledges certain limitations, such as the use of simplified structural models and the constraints imposed by computational time and suggests that future research should focus on more sophisticated modelling and simulations with larger sample sizes.

<h5>1. Introduction</h5>
On February 6, 2023, at 04:17 local time (04:17 UTC+3), a Mw 7.8 earthquake struck southern Turkey, characterized by an 80-second duration of strong ground motion. Nine hours later, a second Mw 7.8 earthquake occurred approximately 96 km north of the initial epicentre, lasting 30 seconds with distinct rupture dynamics. Within 24 hours, the region experienced 150 aftershocks of M≥4.0, including a Mw 7.5 event as the largest aftershock. Casualties and economic impacts included:
>50,000 fatalities across Turkey and Syria
100,000 injuries
2,800+ building collapses
USD 100 billion in economic losses
This seismic sequence represents Turkey's most devastating disaster since the 1939 Erzincan earthquake (Mw 7.8). Notably, newly constructed buildings exhibited unexpected severe damage patterns. Gurbuz attributed these failures to underestimated non-stationary ground motion characteristics in structural design codes, particularly the time-varying spectral acceleration and cumulative energy input during the dual mainshocks <a href="#reference_list_0">[1]</a>.
In the field of engineering, structures are often classified into single-degree-of-freedom (SDOF) and multi-degree-of-freedom (MDOF) systems when conducting dynamic studies. SDOF systems are typically used for simplified analysis and allow for easier evaluation of a building's response to seismic forces. In contrast, MDOF systems provide a more accurate representation of the dynamic behaviour of complex structures. As noted by Makarios, the study of MDOF systems is crucial in structural engineering, particularly when assessing the impact of extreme loads <a href="#reference_list_1">[2]</a>. In recent years, increased tectonic activity has led to an increase in the frequency of earthquakes worldwide. It is often reported in the news that some buildings in earthquake-stricken areas are destroyed, while others remain largely intact, with only specific beams or columns suffering significant deformation or damage. These disparate results not only affect the cost of post-earthquake repairs but also have serious implications for the safety of people in the affected areas. This concern underscores the importance of research on the "dynamic analysis of multi-degree-of-freedom structures."
<img src="https://file.luminescience.cn/6923b759cdf3c1763948377" alt="" width="500" height="333" />
                                                                      Figure 1. Completely failure (left) and partial failure (right) <a href="#reference_list_2">[3]</a>
Research on earthquakes predates the 21st century, and over time, our understanding of how structures respond to seismic events has evolved—from the equivalent static method to pushover analysis. To further investigate the relationship between earthquake intensity and structural response metrics, the incremental dynamic analysis (IDA) method was introduced in 2001, which establishes a significant link between these two factors <a href="#reference_list_3">[4]</a>. However, this method typically relies on one or more specific seismic inputs and cannot fully capture the uncertainties associated with earthquake events. As computing power has increased and the demand for structural safety has grown, existing earthquake records have proven insufficient for large-scale simulations. For this reason, stochastic analysis has been incorporated into earthquake engineering, in which earthquake records are generated through computer simulations. These simulated records, combined with the incremental loading approach of incremental dynamics, mitigate the limitations imposed by the scarcity of earthquake records while introducing additional uncertainties into the analysis, providing a more comprehensive framework for the seismic assessment of buildings.
Building on these advancements, seismic analysis has gained momentum. Researchers, after examining the relationship between structural response and seismic intensity, have begun to explore the application of statistical methods in incremental analysis. To derive higher-order statistics of engineering demand parameters (EDP), such as probability density functions, researchers have extended the outcomes of incremental dynamic analysis and introduced the concept of the IDA surface <a href="#reference_list_4">[5]</a>. This development not only encapsulates the relationship between intensity measures (IM) and EDP but also visualizes the probability of various EDPs occurring under different IMs, further refining seismic analysis methodologies.
Compared to conventional methods, the equivalent static method simplifies seismic loads as static forces and neglects time-dependent effects and nonlinear material responses <a href="#reference_list_5">[6]</a>. In contrast, Stochastic Incremental Dynamic Analysis (SIDA) explicitly captures the degradation of structural performance under progressively amplified seismic actions by incorporating non-stationary ground motion inputs <a href="#reference_list_6">[7]</a>. For instance, static methods exhibit prediction errors of 20–40% in the evaluation of collapse mechanism assessments, whereas SIDA significantly enhances accuracy through dynamic analysis.
Pushover analysis, on the other hand, assumes a lateral load distribution dominated by a single mode of vibration (e.g., the first vibration mode) and does not take into account the contributions of higher modes of vibration or load redistribution after plastic hinge formation <a href="#reference_list_7">[8]</a>. SIDA addresses this limitation by employing scaled accelerograms to excite multi-mode responses, thereby revealing critical weak points in irregular structures, such as the development of plastic hinges at beam-column joints in steel frames.
Recent advancements in the integration of stochastic incremental dynamic analysis (IDA) and meta-modelling techniques have significantly enhanced the efficiency of seismic response analysis for complex structures. For instance, Chen demonstrated that multi-input multi-output nonlinear autoregressive (MIMO-NARX) models enable efficient response prediction for stochastic nonlinear multi-degree-of-freedom (MDOF) systems, reducing computational costs while maintaining high accuracy <a href="#reference_list_8">[9]</a>. Mitseas and Beer introduced a first-excursion stochastic incremental dynamics methodology, which incorporates Monte Carlo simulations and hysteretic models to quantify time-dependent failure probabilities of MDOF systems under non-stationary seismic excitations <a href="#reference_list_9">[10]</a>. Additionally, Xu proposed a deep reinforcement learning framework for generating non-stationary ground motions that align with target response spectra, offering a novel approach to simulate realistic seismic inputs for infrastructure resilience assessment <a href="#reference_list_10">[11]</a>. These methodologies collectively address critical gaps in conventional seismic analysis, particularly in capturing uncertainties and time-varying characteristics of earthquake events.
Despite these advancements, most existing studies focus on analysing the structural response in stationary seismic excitations and neglect the impact of non-stationary excitations. Additionally, although IDA has been widely used in seismic analysis, research on stochastic incremental dynamic analysis (Stochastic IDA) is still insufficient, particularly when applied to MDOF systems. Given the non-stationarity of actual seismic excitations and the complexity of MDOF structural systems, there is an urgent need to perform stochastic incremental dynamic analysis. With this method, the behavior of structures under complex seismic conditions can be simulated more accurately, providing more reliable data for seismic design and structural safety assessment.
<h5> 2. Methodology</h5>
This study aims to perform a stochastic incremental dynamic analysis of multi-degree-of-freedom (MDOF) structures under non-stationary seismic excitations using Monte Carlo simulation method, with the goal of establishing the relationship between Intensity Measures (IM), Engineering Demand Parameters (EDP), and probability density. To achieve this, the research starts with the identification of the relevant IM and EDP. A combined approach utilising programming software and finite element analysis is then applied. A structural model is developed in ABAQUS, and, following the guidelines of EUROCODE 8, random seismic waves are generated in MATLAB using Monte Carlo simulation. These waves are then input into ABAQUS for nonlinear time-history analysis to determine the maximum inter-story drift ratio of the structure.
<img src="https://file.luminescience.cn/6922a3714eb591763877745" alt="" width="374" height="711" />
                                                                             Figure 2. Overview flow chart of the methodology
2.1 Parameter selection for incremental dynamic analysis
Incremental Dynamic Analysis (IDA) establishes a connection between IM and EDP. In this study, Peak Ground Acceleration (PGA) is selected as the Intensity Measure (IM), and the maximum inter-story drift ratio is chosen as the Engineering Demand Parameter (EDP).
2.2 Model
2.2.1 Introduction to ABAQUS
ABAQUS is an important tool for finite element analysis in earthquake engineering, which is used for static and quasi-static problems, as well as complex nonlinear analysis. Static and quasi-static analyses are employed to simulate the deformation and internal force distribution of structures before and after an earthquake. These analyses help engineers understand the residual deformation and potential damage patterns of structures after an earthquake. On the other hand, dynamic analysis in ABAQUS can deal with complex material nonlinearity and geometric nonlinearity (such as plasticity, fracture, buckling, etc.), which makes it particularly effective for simulating large deformations and complex contact problems.
In dynamic analysis, ABAQUS's implicit solver method demonstrates extremely high stability when dealing with nonlinear large deformation problems caused by earthquakes, especially when analysed over long-duration and when integrating large time steps, making it an ideal choice for detailed structural analysis in earthquake engineering.
2.2.2 Model development
<img src="https://file.luminescience.cn/6922a736a1ed01763878710" alt="" width="690" height="57" />
                                                      Figure 3. Flow chart of modelling in ABAQUS3.3. Generate Seismic Wave
Monte Carlo Simulation is a numerical method that utilizes random sampling to estimate the outcomes of complex systems or processes. By repeatedly performing random sampling and iterating numerous times, Monte Carlo Simulation can approximate solutions to problems that are typically challenging to solve directly through analytical methods. The core concept of Monte Carlo Simulation is to model the behaviour of a system or process through extensive random sampling, thereby estimating the probability distribution of its outcomes. In simple terms, if researchers wish to understand the output results of a complex system under varying input conditions, they can randomly generate these input conditions, simulate the system, and then observe the distribution of the resulting outputs. This method is valuable for estimating statistical quantities such as the expected value and variance of the system and is widely used in risk assessment and decision analysis.
Mathematically, Monte Carlo Simulation relies on the Law of Large Numbers and the Central Limit Theorem. Suppose we have a random variable  with a probability density function f(x), and we want to estimate the expected value IE(x):
<img src="https://file.luminescience.cn/6922a89bcbccc1763879067" alt="" width="187" height="68" />
Since directly calculating the above integral can be very challenging, the Monte Carlo method estimates it through the following steps:
(1) Generate Random Samples: Randomly generate  N samples x1, x2, ..., xn from the distribution f(x).
(2) Calculate Sample Mean: Calculate the mean of these samples as an estimate of the expected value:
<img src="https://file.luminescience.cn/6922ab9897a131763879832" alt="" width="144" height="77" />
According to the Law of Large Numbers, when the sample size  is sufficiently large, the sample mean IE(x)  will converge to the true expected value.
Additionally, by calculating the sample variance and standard error, the uncertainty of the simulation results can also be estimated.
In earthquake engineering, the response of multi-degree-of-freedom (MDOF) systems is often complex and highly nonlinear. Monte Carlo simulation effectively addresses this complexity, particularly under non-stationary seismic excitations. By randomly generating various seismic inputs, it captures the dynamic behaviour of the system across different potential earthquake scenarios <a href="#reference_list_11">[12]</a>. This approach offers distinct advantages for analysing the dynamic response of MDOF systems, especially when seismic motion exhibits significant time-varying characteristics.
Monte Carlo simulation quantifies the response uncertainty of MDOF systems under different seismic excitations by generating a large number of random samples. This capability is crucial in earthquake engineering, as it enables engineers to assess potential risks to structures during future seismic events, providing a quantitative basis for design and retrofit decisions. The method is highly flexible and can adapt to various types of seismic excitations and structural characteristics. Whether dealing with short-period or long-period non-stationary seismic waves or considering complex factors such as material and geometric nonlinearities in MDOF systems, Monte Carlo simulation offers accurate response predictions.
In this study, Monte Carlo simulation is primarily employed to simulate seismic uncertainty by generating random peak ground accelerations (PGA)  ag following a lognormal distribution and random phase angles ∅ with uniform distribution. Other parameters, such as the natural vibration period T, are deterministic due to structural configuration and height variations, while soil parameters lack randomness as they are governed by region-specific geotechnical conditions. By integrating these parameters (excluding the random phase angle ∅ ) with Equations (6)-(8) in Section 3.2.1, the elastic response spectrum of the structure is derived. The randomly generated phase angles ∅ are then utilized to synthesize frequency-domain seismic records in MATLAB. Through N independent repeated experiments (where N is researcher-defined based on computational resources and precision requirements), the methodology outputs N sets of ground motion time histories compliant with Eurocode 8 specifications <a href="#reference_list_12">[13]</a>. The complete MATLAB implementation is provided in Appendix A.
The inverse Fourier transform (IFT) is implemented in MATLAB to synthesize time-domain ground motion accelerograms from stochastic frequency-domain signals. The procedure comprises the following steps:
A.Frequency-domain signal construction:
The target elastic response spectrum  Se(f) is derived from Eurocode 8, incorporating site coefficients (S) and damping correction factors (η) :
<img src="https://file.luminescience.cn/6922ae2b5b1bd1763880491" alt="" width="273" height="110" />
Where ag follows a lognormal distribution to represent seismic intensity variability.
B. Phase randomization:
A stochastic phase vector ∅∈ [0,2π) is generated for each frequency component to emulate natural earthquake randomness
C. Complex spectrum synthesis:
The frequency-domain signal combines spectral amplitudes Se(f) and randomized phases:
<img src="https://file.luminescience.cn/6922af368a6bc1763880758" alt="" width="181" height="42" />
ensuring conjugate symmetry for real-valued time-domain outputs.
D. Time-domain conversion via IFT:
The discrete inverse Fourier transform is executed with the 'symmetric' flag to enforce Hermitian symmetry
E. Non-stationary modulation:
A trapezoidal envelope ω(t) modulates the synthesized signal to introduce time-dependent energy characteristics:
<img src="https://file.luminescience.cn/6922afbeebaa71763880894" alt="" width="251" height="46" />
Additional MATLAB code details for implementing the inverse Fourier transform are provided in the Appendix.
<img src="https://file.luminescience.cn/6922afee7107e1763880942" alt="" width="677" height="72" />
                                                                                  Figure 4. Flow chart of generating a random seismic wave
2.3 Non-linear Time History Analysis (NTHA)
Nonlinear Time History Analysis (NTHA) is a detailed and complex seismic engineering analysis method used to evaluate the dynamic response of structures under seismic excitation. Unlike linear analysis methods, NTHA considers material nonlinearity, geometric nonlinearity, and the interactions between components, allowing for a more accurate simulation of the real behaviour of structures under extreme loads. The general analysis steps are as follows:
Seismic Input: In NTHA, the structural model is subjected to actual or artificially generated seismic inputs. The seismic input is provided in the form of acceleration time histories, and multiple seismic records are usually employed to capture the uncertainty of the earthquake.
Modelling of nonlinear behaviour: The nonlinear factors considered in the analysis include: (i) Material nonlinearity: Such as the yielding of steel structures or cracking of concrete. (ii) Geometric nonlinearity: Such as large deformation effects and P-Δ effects.
Numerical solution: The structural equations of motion are solved using numerical integration methods (such as the Newmark method or the Hilber-Hughes-Taylor method) to obtain the structural response at each time step.
Through NTHA, detailed information about the structural response under seismic action can be obtained. This study mainly focuses on:
(1) Displacement and deformation: The time history of displacement and deformation at various points of the structure, especially the displacement curves of key nodes. These data help understand the deformation pattern and identify the maximum deformation locations of the structure during an earthquake.
(2) Interstory drift ratio: The inter-story drift ratio of each floor is calculated to evaluate whether the structure meets the seismic code requirements and to identify potential mechanisms of story drift or collapse.
In certain cases, by analysing the energy distribution (such as kinetic energy, potential energy, and dissipated energy), the structure's ability to absorb and dissipate energy under seismic action can be assessed, providing further insights into the seismic resilience of the structure.
<h5>3. Simulation of multi-degree-of-freedom structures under non-stationary seismic excitations</h5>
3.1 Overview of the structure
The structural model in this study is based on the Sentinel Tower dormitory at the University of Leeds. It is simplified in ABAQUS as a five-story steel frame with five spans in the x-direction and three spans in the y-direction. The model is shown in the figure below.
<img src="https://file.luminescience.cn/6922b11a284a11763881242" alt="" width="666" height="275" />
                                                            （a）Unrendered model                                          (b) rendered model
                                                                                               Figure 5. Schematic of the model
The structure uses S355 grade steel, with a floor height of 2.7 meters. Total length of structure 17.5 metres, width 10.5 metres, height 13.5 metres. The columns have a box section, and the beams have an I-shaped section. The specific dimensions are shown in the table.
                                                                                   Table 1. Member size of the ABAQUS model
<img src="https://file.luminescience.cn/6922b2229bb431763881506" alt="" width="549" height="80" />
For S355 steel, considering its density, elastic stress, and plastic stress/strain, the specific data are provided in the table below.
                                                                                         Table 2. Elastic data of ABAQUS model
<img src="https://file.luminescience.cn/6923b8bb3da9d1763948731" alt="" width="542" height="58" />
                                                                                         Table 3. Plastic data of ABAQUS model
<img src="https://file.luminescience.cn/6923b84be86301763948619" alt="" width="618" height="114" />
3.2 Seismic wave
The process of generating random seismic waves using Monte Carlo simulation through MATLAB has already been presented in the previous chapter with a flowchart. Here, only the explanation and selection of the key parameters in the code are provided; the full code can be found in the appendix.
3.2.1 Parameters
According to EC8, Clause 3.2.2.2 (2), assuming the soil type is Class B with a viscous damping ratio of 5%, a Type 1 response spectrum is used. The parameters describing the response spectrum are provided in the table below.
                                                                                Table 4. Spectrum parameters according to EUROCODE 8
<img src="https://file.luminescience.cn/6922b3ade25f41763881901" alt="" width="551" height="57" />
The horizontal elastic response spectrum of the structure is as follows:
<img src="https://file.luminescience.cn/6922b4a5b11481763882149" alt="" width="386" height="201" />
3.2.2 Parameters for Monte Carlo Simulation
The structure has a floor height of 2.7 meters and a total height of 13.5 meters. The period of the first mode is calculated using the empirical formula as follows:
<img src="https://file.luminescience.cn/6922b5111af4d1763882257" alt="" width="103" height="29" />
Where Ct and x are coefficients; for steel structures, they are selected as 0.085 and 0.75, respectively.
<img src="https://file.luminescience.cn/6922b55e6f3a71763882334" alt="" width="258" height="37" />
To capture the high-frequency characteristics during analysis, the lower limit of the period is set to T = 0.1  seconds, and to ensure the accuracy of low-frequency information, the upper-frequency limit is set to 4T = 2 seconds. Considering computational efficiency and convergence, 1000 sampling points are selected for the period.
During an earthquake, PGA follows a lognormal distribution, but the elastic response spectrum in EC8 does not directly provide PGA. Therefore, ag is assumed to follow a lognormal distribution, and by changing the logarithmic mean of ag, PGA can be controlled. Changing the standard deviation of ag allows adjustment of the dispersion of acceleration at different moments in the seismic record and the probability of generating extreme values.
3.2.3 Time modification function
Amin and Jennings proposed using a time modulation function to simulate the non-stationary characteristics of an earthquake, specifically to represent the initial and decay stages of an earthquake through the following equation <a href="#reference_list_13">[14]</a>:
<img src="https://file.luminescience.cn/6922b67c839741763882620" alt="" width="208" height="113" />
Where ta and tb are the start and end times of the stationary phase of the seismic motion, respectively, and α is the control parameter for the decay phase. For an earthquake record with a total duration of 30 seconds, ta is set to 8 seconds, tb is set to 18 seconds, and α.
3.3 Mesh control
Mesh generation is a critical step in ABAQUS analysis, where a continuous component is divided into a finite number of elements, allowing for the individual solution of these elements. This discretization process is particularly effective for analysing geometrically complex shapes. However, a denser mesh requires solving more finite elements, which can significantly increase computation time. For the model used in this study, which features a regular shape and simple components, it is especially important to balance computational accuracy and efficiency.
In this section, to minimise the influence of seismic records on the results, the north-south component of the El Centro earthquake is used for the simulation. All other settings remain unchanged except for the mesh size. The mesh sizes and corresponding simulation times are presented in the table below, and the displacement diagrams are illustrated in the accompanying figure.
<img src="https://file.luminescience.cn/6922b73d541e21763882813" alt="" width="734" height="281" />
                                                                                     Figure 6. Structural displacements for mesh size=100
<img src="https://file.luminescience.cn/6922b7933e79a1763882899" alt="" width="758" height="296" />
                                                                                Figure 7. Structural displacements for mesh size=250
<img src="https://file.luminescience.cn/6922b7de55a2c1763882974" alt="" width="756" height="305" />                                                                                 Figure 8. Structural displacements for mesh size=500
                                                                                      Table 6. Analysis time for different mesh sizes
<img src="https://file.luminescience.cn/6922b850140c41763883088" alt="" width="495" height="102" />
Typically, the smaller the mesh size, the more accurate the results. However, in Group C, with a mesh size of 500, the displacement of the structure at the same measurement points became highly distorted. Moreover, during the 30-second earthquake simulation, displacement was only recorded in the first 4 seconds, indicating that the analysis results are divergent. In contrast, the displacement diagrams for Groups A and B are nearly identical, but the computation time for Group B was 30% longer than for Group A. Therefore, a mesh size of 250 will be uniformly applied for all subsequent analyses in this study.
3.4 Selection of monitor point
When seismic motion occurs, the maximum inter-story drift in a regular structure generally occurs at the corners of the structure. Therefore, the measurement points are selected as shown in the figure below:
<img src="https://file.luminescience.cn/6922b89f0556b1763883167" alt="" width="371" height="281" />
                                                                             Figure 9. Monitoring points (MP) for structural modelling
<h5>4. Data analysis</h5>
<img src="https://file.luminescience.cn/6922b904ebba21763883268" alt="" width="517" height="480" />
                                                                              Figure 10. Overview flow chart of data analysis
4.1 Interstory drift ratio
4.1.1 Limit value
According to EC8 Clause 4.4.3.2, consider the worst-case:
<img src="https://file.luminescience.cn/6922b989561401763883401" alt="" width="123" height="32" />
Where dr is the inter-story drift, and v is a reduction factor related to damage limitation requirements. However, due to varying practical requirements, the story heights of a structure often differ across floors. To eliminate the influence of varying story heights and more intuitively quantify structural deformation per unit height under seismic loads, it is more rational to perform failure probability analysis using the ratio of the maximum interstory displacement to the corresponding story height for each floor.
Therefore, the inter-story drift ratio:
<img src="https://file.luminescience.cn/6922b9d4cb39b1763883476" alt="" width="346" height="55" />
Which hi is the height of each floor, considering the low recurrence period of seismic activity. It is generally taken as 0.4 or 0.5, depending on the importance level of the structure. In this study, considered the worst case,  is used.
4.1.2 Analysis of incremental dynamic analysis surface
After performing kernel density estimation on the simulation results, the resulting IDA surface is shown in the figure below:
<img src="https://file.luminescience.cn/6922ba3f877e91763883583" alt="" width="682" height="382" />
                                                                    Figure 11. IDA surface (IM: PGA, EDP: Max inter-story drift ratio)
The red plane represents the maximum elastic inter-story drift ratio limit. The probability density peak in the figure is concentrated between an inter-story drift ratio of 0.003-0.004, which is more likely to occur during lower-intensity earthquakes (with PGA ranging from 0.1g to 0.3g).
The IDA surface is presented in a three-dimensional format, illustrating the relationship between PGA (one horizontal axis) and inter-story drift ratio (the other horizontal axis), with probability density as the third dimension (vertical axis). This 3D surface reflects the inter-story drift response of the structure under different earthquake intensities (PGA) and represents the probability density through colour. The yellow areas in the figure indicate that, at specific combinations of PGA and inter-story drift ratio, the structure has a higher probability of reaching that drift ratio. This implies that as earthquake intensity increases, the structure is more likely to exhibit a particular range of inter-story drift ratios.
The red plane represents a limit value of 0.005 for the drift ratio between the projecticles. This is the limit value for the elastic inter-story drift ratio for steel frame structures as specified by Eurocode 8. This limit value is used to assess whether the structure remains in a safe state under seismic action. Exceeding this value generally indicates structural failure or the potential for significant deformation. The area on the right side of the IDA surface, beyond the red plane, indicates that the structure is very likely to exveed the Interstory Drift Ratio limit at these specific PGA values. As the PGA value increases, the IDA surface demonstrates a corresponding increase in the inter-story drift ratio. This trend aligns with physical intuition: as earthquake intensity increases, structural deformation also increases.
The figure illustrates that at lower PGA values (e.g., 0.1g to 0.2g), the inter-story drift ratio remains relatively small, and the surface appears fairly flat. This indicates that in these cases, most structural deformation is concentrated within a narrower range, with a lower risk of exceeding the specified limit. Conversely, at higher PGA values (e.g., 0.4g to 0.5g), the inter-story drift ratio increases significantly, and the surface becomes steeper, suggesting that structural deformation is more widespread, and the risk of exceeding the limit is greater.
On the right side of the red plane, the probability density is at its lowest when the PGA reaches 0.2g, followed by a general upward trend as the PGA reaches or exceeds 0.35g. For structural engineers, this indicates that once the PGA surpasses this threshold, the risk of structural failure becomes significant. The shape of the surface also reflects the deformation distribution under different combinations of PGA and inter-story drift ratio. For instance, the probability density associated with higher PGAs and larger inter-story drift ratios is lower, suggesting that while such scenarios are less likely, they could result in severe structural failure if they do occur.
Engineers might use an IDA surface like this to identify critical PGA values for design purposes. For example, if a building is situated in a high seismic risk area where the expected PGA could reach 0.4g or higher, it is essential to ensure that the structure maintains sufficient elasticity at these levels to avoid exceeding the limit. Based on this figure, engineers could also calculate the cumulative failure probability of the structure at different PGA levels by integrating the surface above the red plane, which would help quantify the overall risk of failure at a specific PGA.
4.1.3 Estimation of failure probability
The failure probability of the structure at different PGA levels was calculated using the cumulative distribution function (CDF). The results are shown in the table below.
                Table 8. Means, standard deviations, and failure probabilities of maximum inter-story displacement ratios of structures under different PGAs
<img src="https://file.luminescience.cn/6922bacb6761d1763883723" alt="" width="497" height="154" />
The data includes the mean (μ), standard deviation (σ), and failure probability of the structure's inter-story drift ratio at different PGA levels.
The mean value of the inter-story drift ratio gradually increases with increasing PGA. This is due to the fact that the higher the seismic intensity, the greater the deformation of the structure. At PGA=0.1g, the mean inter-story drift ratio is 0.00063, while it increases to 0.00421 at PGA = 0.5g. The standard deviation represents the degree of dispersion in the data. As the PGA increases, the standard deviation of the inter-story drift ratio also increases, indicating that the uncertainty (or variability) of the structural deformation increases with stronger seismic actions. For example, the standard deviation at PGA=0.1g is 0.00020, while it increases to 0.00114 at PGA=0.5g.
At PGA=0.1g, 0.2g, and 0.3g, the probability of failure is 0. This means that at these earthquake intensities, the inter-story drift ratio of the structure is not to exceed the set failure limit (i.e., 1/200 = 0.005). When the PGA reaches 0.4g, the failure probability rises to 0.03254 (approximately 3.3%). This indicates that at a seismic intensity of PGA=0.4g, there is a 3.3% probability that the structure will fail. At PGA=0.5g, the failure probability significantly increases to 0.24281 (approximately 24.3%). This means that at this level of seismic intensity, there is a high risk of structural failure, with nearly a one-in-four chance that the inter-story drift ratio will exceed the set limit.
As the PGA increases, the failure probability shows a nonlinear growth trend. Particularly beyond PGA=0.4g, the failure probability begins to rise sharply, indicating a higher risk of failure. This nonlinear growth may be related to factors such as the nonlinear behaviour of structural materials and geometric nonlinear deformation. The data suggests that PGA=0.4g is a critical point. Beyond this point, the failure probability begins to increase significantly, indicating that structures approaching PGA=0.4g may require higher safety design or reinforcement measures.
Therefore, for buildings located in high seismic risk areas, particularly those that may experience seismic intensities of PGA=0.4g or higher, special consideration should be given to enhancing the seismic capacity of the structure in the design to reduce the failure probability.
4.1.4 Data distribution fitting
To validate the statistical characteristics of the inter-story drift ratio, a normal distribution fitting was performed. The probability density function (PDF) and cumulative distribution function (CDF) graphs are shown below:
<img src="https://file.luminescience.cn/6922bb5edeafb1763883870" alt="" width="633" height="402" />
                                                                                               Figure 12. PDF for PGA=0.5g
<img src="https://file.luminescience.cn/6922bbab1e3b51763883947" alt="" width="652" height="376" />                                                                                                Figure 13. CDF for PGA=0.5g
<img src="https://file.luminescience.cn/6922bbfa5e7c91763884026" alt="" width="653" height="388" />
                                                                                                   Figure 14. PDF for PGA=0.4g
<img src="https://file.luminescience.cn/6922bc35b8db41763884085" alt="" width="661" height="373" />
                                                                                                    Figure 15. CDF for PGA=0.4g
<img src="https://file.luminescience.cn/6922c0c22447f1763885250" alt="" width="696" height="385" />                                                                                                  Figure 16. PDF for PGA=0.3g
<img src="https://file.luminescience.cn/6922c143be6981763885379" alt="" width="658" height="379" />
                                                                                                      Figure 17. CDF for PGA=0.3g
<img src="https://file.luminescience.cn/6922c192c53a21763885458" alt="" width="657" height="377" />
                                                                                                     Figure 18. PDF for PGA =0.2g
<img src="https://file.luminescience.cn/6922c1f04d27e1763885552" alt="" width="651" height="412" />                                                                                                  Figure 19. CDF for PGA=0.2g
<img src="https://file.luminescience.cn/6922c22a8fb501763885610" alt="" width="662" height="377" />
                                                                                                   Figure 20. PDF for PGA=0.1g
<img src="https://file.luminescience.cn/6922c29491fa81763885716" alt="" width="605" height="362" />
                                                                                                   Figure 21. CDF for PGA=0.1g
Observing the probability density function graph, as the PGA increases, the probability density of generating a larger interstory drift ratio also increases, which is a natural response of the structure under stronger seismic excitation. At lower PGA levels (e.g., 0.1g and 0.2g), the normal distribution can fit the inter-story drift ratio data well. However, as the PGA increases, at PGA=0.3g, the fitting curve shows deviations from the actual data in the extreme value regions, i.e., the tails of the probability density distribution. This deviation becomes more pronounced at PGA=0.4g and PGA=0.5g, where the normal distribution starts to struggle in accurately describing the data distribution, especially in the high inter-story drift ratio region. Similarly, the cumulative distribution function graph exhibits similar characteristics, with its variability becoming more pronounced at PGA=0.5g.
4.2 Peak ground acceleration
As mentioned in Chapter 4, PGA follows a lognormal distribution <a href="#reference_list_14">[15]</a>. To verify the fit of the simulated data to a lognormal distribution and to confirm the data variability resulting from increased PGA as discussed in Section 5.1, the logarithm of PGA was taken in MATLAB, and a normal distribution fitting was performed. The Q-Q plot is shown below:
<img src="https://file.luminescience.cn/6922c30ec306e1763885838" alt="" width="720" height="410" />
                                                                                             Figure 22. Q-Q plot for PGA=0.1g
<img src="https://file.luminescience.cn/6922c37c8f0251763885948" alt="" width="680" height="376" />
                                                                                                  Figure 23. Q-Q plot for PGA=0.4g
<img src="https://file.luminescience.cn/6922c413107491763886099" alt="" width="667" height="369" />
                                                                                                  Figure 24. Q-Q plot for PGA=0.5g
Observing the Q-Q plot, whether for small PGA (0.1g) or large PGA (0.4-0.5g), the middle-range data fits the curve quite well, with 0.1g being the closest and 0.4g showing slightly more fluctuation. However, the data set for PGA=0.5g shows significant deviations at the ends and tails, which could be due to the smaller sample size, making it more susceptible to random fluctuations. Therefore, a second verification was conducted using the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). The results are shown in the table below.
                                                                                          Table 9. AIC/BIC values at different PGAs<img src="https://file.luminescience.cn/6922c46e0371d1763886190" alt="" width="390" height="109" />
The smaller the values of AIC and BIC, the more precise the fit. Using this method, it was found that the data for PGA=0.4g conforms more closely to a lognormal distribution, which differs from the results obtained from the Q-Q plot. This discrepancy may be due to the following reasons:
(1) Difference Between Visual Judgment from Q-Q Plot and Numerical Differences in AIC/BIC: The Q-Q plot is an intuitive tool used to compare the fit between the actual data distribution and the theoretical distribution. By observing whether the data points align along the theoretical straight line, one can judge the similarity between the data and the theoretical distribution. The conclusion from a Q-Q plot relies on the comparison of data quantiles, particularly the performance of the data's tail and middle parts.
AIC and BIC are numerical indicators calculated based on the log-likelihood function. They evaluate models by combining goodness of fit and complexity. AIC/BIC considers the fit of the model to the entire dataset and includes a penalty term to prevent overfitting. Therefore, the results from AIC/BIC might differ from the visual judgment of the Q-Q plot.
(2) Different Weights on Data Characteristics: When judging the goodness of fit, the Q-Q plot might focus more on the match of the distribution shape, especially the match of quantiles with the theoretical distribution. It visually demonstrates the deviation of the data from the theoretical distribution, particularly in extreme values (tails).
AIC and BIC evaluate the model's log-likelihood as a whole and pay less attention to local distribution shapes. Thus, even if the Q-Q plot shows that some local fits (such as the middle part) are good, AIC/BIC might give a poorer score due to a lower overall log-likelihood.
(3) Impact of Model Complexity: AIC and BIC include a complexity penalty term. If a model has a good fit to the data but is complex, the penalty term increases the AIC and BIC values. This could result in models with good fit performing poorly in AIC/BIC scores. Conversely, simpler models have smaller penalty terms, potentially leading to better AIC/BIC scores.
(4) Sample Size and Data Characteristics: With a smaller sample size, the Q-Q plot may be more susceptible to random fluctuations, potentially leading to misleading results. AIC/BIC considers sample size and adjusts through complexity penalties, which might explain the inconsistency between the two results. If there are significant differences in distribution characteristics between datasets (e.g., different frequencies of extreme values), the Q-Q plot may be more sensitive in such cases, while AIC/BIC may focus more on overall fit.
(5) Bias in Model Assumptions: AIC/BIC relies on the assumptions of the chosen model. If the assumed model does not fully suit the data distribution, it may lead to higher AIC/BIC values. On the other hand, the Q-Q plot is a visual comparison against a specific model distribution, so its sensitivity to model assumption bias might differ.
<h5>5. Conclusion</h5>
This study first describes a method for generating stochastic non-stationary seismic records in MATLAB that align with the elastic response spectrum outlined in EUROCODE 8. Monte Carlo simulations were employed to generate ground accelerations (ag) following a lognormal distribution. An inverse Fourier transform was then applied to produce a smooth seismic excitation, which was subsequently modified with a time modulation function to simulate its non-smooth characteristics. Statistical estimation and distribution fitting of the PGA from the generated seismic recordings demonstrated that:
(1) Its statistical characteristics are largely consistent with a lognormal distribution. Subsequently, the generated series of seismic records were input into ABAQUS to perform a nonlinear time-dependent analysis of a five-degree-of-freedom steel frame model to obtain the probability density of the maximum inter-story drift ratio of the structure, i.e., the IDA surface, at different PGAs. It is found that for a limit of 1/200 of the interstory drift ratio.
(2) At low Peak Ground Accelerations (PGA) of 0.1g, 0.2g, and 0.3g, the interstory drift ratio of the structure is minimal, indicating that the relative movement between floors is limited. This suggests that the building undergoes very little deformation small deformations under these seismic loads, which minimizes stress and strain on the structural components. As a result, the risk of structural failure is almost negligible, so the integrity of the building is maintained even under such low-level seismic events. Consequently, the structure is considered safe and stable under these conditions, providing confidence in its ability to withstand minor ground motions without any significant damage.
(3) Under medium Peak Ground Acceleration (PGA) of 0.4g, the probability of structural failure begins to increase, signaling a shift in the performance level of the building. At this level of seismic intensity, the interstory drift ratio becomes more pronounced, indicating greater relative movement between floors. This increased movement can result in heightened stress and strain on the structural components, such as beams, columns, and joints, which may compromise their ability to perform optimally. Consequently, there is a growing need for careful monitoring and assessment of the structure's response to ensure that any signs of potential damage or failure are promptly identified and addressed. Special attention should be given to the design, detailing, and material properties of the structural elements to ensure they can adequately resist the forces and deformations imposed by such seismic events. Reinforcement strategies, retrofitting measures, or other forms of intervention may be necessary to maintain the desired level of safety and performance.
(4) The failure probability increases significantly at high PGA (0.5g), indicating that the structure has a high risk of failure at this seismic intensity. This set of data can help engineers in seismic risk assessment to better understand the performance of structures under different seismic intensities, and in particular to identify which PGA levels have a significantly higher risk of structural failure so that targeted protective measures can be developed.
It is important to note that this study has certain limitations, primarily manifested in limited computational resources and the use of simplified structural models. The computational power constraints of a single workstation restrict the number of simulations that can be performed within a reasonable timeframe, potentially resulting in a limited number of samples that may not fully capture the statistical characteristics of seismic responses. Furthermore, modern architectural designs increasingly adopt irregular geometries to enhance aesthetic and functional performance. However, such irregularities can amplify or attenuate seismic waves, significantly increasing uncertainty in structural responses. Consequently, the simplified models employed in this study lack generalizability for practical applications involving complex architectural configurations.
To address these limitations, a potential solution is proposed for future research. First, artificial intelligence (AI) could be leveraged to predict and model irregular building geometries, enabling realistic representation of modern architectural forms. Subsequently, cloud computing platforms (e.g., AWS, Azure) could be utilized to perform large-scale Monte Carlo simulations. This approach not only mitigates the oversimplification of building geometries but also overcomes computational bottlenecks, thereby reducing costs while improving simulation accuracy and efficiency. Such advancements would enhance the applicability of the methodology to real-world seismic risk assessments.
Overall, this study addresses the limitations of conventional stationary ground motion models in capturing critical features of extreme seismic events, such as long-period pulses, while revealing the nonlinear relationship between the probability of structural failure and seismic intensity. Specifically, the results demonstrate that under low Peak Ground Acceleration (PGA) conditions (e.g., 0.1–0.3g), the structure remains predominantly elastic, aligning with current code-based safety requirements for minor earthquakes. However, a sharp nonlinear increase in failure probability is observed at high PGA levels (e.g., ≥0.5g), highlighting the necessity for optimized ductility design strategies, including joint reinforcement and energy dissipation mechanisms, to enhance structural resilience under extreme loading. By integrating Monte Carlo simulations for probabilistic ground motion generation with ABAQUS-based nonlinear time-history analysis, this research establishes a probabilistic framework for seismic risk assessment. This framework enables the identification of high-risk zones (e.g., 0.5g PGA thresholds) to prioritize retrofitting efforts and quantifies expected economic losses under varying seismic intensities, thereby informing data-driven resource allocation and infrastructure resilience planning. The findings collectively advance performance-based seismic design by bridging non-stationary excitation modeling, nonlinear structural response prediction, and probabilistic risk management.
<h5>Conflicts of interest</h5>
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the analysis or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
<h5>Funding </h5>
The authors declare that no specific funding or external grants were received from any public, commercial, or not-for-profit sectors for the research, authorship, and/or publication of this article.
<h5>Nomenclatures</h5>
Roman symbols
c: Damping coefficient, use N·s/m.
k: Structural stiffness, use N/m.
m: Mass, use kg.
Greek symbols
ζ: Damping ratio.
ω: Angular frequency, use rad/s.
Abbreviations
IDA: Incremental Dynamic Analysis. A method for assessing structural performance by scaling ground motions to increasing intensity levels.
MDOF: Multi-Degree-of-Freedom. A structural system requiring multiple independent coordinates to describe its motion.
MIMO-NARX: Multi-Input Multi-Output Nonlinear Autoregressive with Exogenous Inputs. A surrogate modelling technique that predicts nonlinear dynamic responses of multi-parameter systems using historical outputs and external inputs, enabling efficient seismic analysis of MDOF structures.
NTHA: Nonlinear Time History Analysis. A computational approach simulating structural responses using time-varying earthquake loads and nonlinear material behaviour.
PGA: Peak Ground Acceleration. The maximum horizontal ground acceleration observed during an earthquake, measured in g (gravity units).
SDOF: Single-Degree-of-Freedom. A simplified structural model where motion is described by a single coordinate.
SIDA: Stochastic Incremental Dynamic Analysis. An extension of IDA that incorporates stochastic ground motion simulations to quantify uncertainties in structural fragility estimates.
General Terms
Non-stationary excitation: Ground motion with time-varying frequency/amplitude.

References

[1]Gurbuz T, et al. Damages and failures of structures in İzmir (Turkey) during the October 30, 2020 Aegean Sea Earthquake. Journal of Earthquake Engineering. 2022;27(6):1565-1606. doi:10.1080/13632469.2022.2086186.

[2]Makarios KT. The Equivalent Non-Linear Single Degree of Freedom System of Asymmetric Multi-Storey Buildings in Seismic Static Pushover Analysis [Internet]. In: Earthquake-Resistant Structures - Design, Assessment and Rehabilitation. InTech; 2012. Available from: http://dx.doi.org/10.5772/27245.

[3]Carocci CF. Small centres damaged by 2009 L'Aquila earthquake: On site analyses of historical masonry aggregates. Bulletin of Earthquake Engineering. 2012;10(1):45-71.

[4]Vamvatsikos D, Cornell CA. Incremental dynamic analysis. Earthquake Engineering & Structural Dynamics. 2002;31(3):491-514.

[5]Li J, Yan Q, Chen JB. Stochastic modeling of engineering dynamic excitations for stochastic dynamics of structures. Probabilistic Engineering Mechanics. 2012;27(1):19-28.

[6]FEMA. Prestandard and Commentary for the Seismic Rehabilitation of Buildings. FEMA-356. Washington, DC; 2000.

[7]FEMA. Seismic Performance Assessment of Buildings. FEMA P-58. Washington, DC; 2012.

[8]Chopra AK, Goel RK. A modal pushover analysis procedure for estimating seismic demands for buildings. Earthquake Engineering & Structural Dynamics. 2002;31(3):561-582.

[9]Chen M, Gao X, Chen C, Guo T, Xu W. A Comparative Study of Meta-Modeling for Response Estimation of Stochastic Nonlinear MDOF Systems Using MIMO-NARX Models. Applied Sciences. 2022;12(22):11553. doi:10.3390/app122211553.

[10]Mitseas IP, Beer M. First-excursion stochastic incremental dynamics methodology for hysteretic structural systems subject to seismic excitation. Computers & Structures. 2021;242:106389.

[11]Xu K, Zhang L, Wang Y. A deep reinforcement learning approach for non-stationary ground motion synthesis. Earthquake Engineering & Structural Dynamics. 2023;52(6):1789-1810. doi:10.1002/eqe.3845.

[12]Musson RMW. Determination of design earthquakes in seismic hazard analysis through Monte Carlo simulation. Journal of Earthquake Engineering. 1999;3(4):463-474.

[13]CEN. Eurocode 8: Design of structures for earthquake resistance – Part 1: General rules, seismic actions and rules for buildings. BS-EN-1998-1-2004. 2004.

[14]Amin M, Ang AHS. Nonstationary stochastic models of earthquake motions. Journal of the Engineering Mechanics Division. 1968;94(2):559-584.

[15]Li J, Chen J. Stochastic Dynamics of Structures. John Wiley & Sons; 2010.

Download PDF

Download XML

How to Cite

Kong, D. (2025). Study on structural response under vibration load with multiple degrees of freedom. Civil and Energy Research, 1(1), 39–60. Retrieved from https://ojs.luminescience.cn/CER/article/view/344

Download Citation

Issue

Civil and Energy Research, Volume 1, Issue 1, 2025

Section

Original Research

License

This work is licensed under a Creative Commons Attribution 4.0 International License.

Copyright licenses detail the rights for publication, distribution, and use of research. Open Access articles published by Luminescience do not require transfer of copyright, as the copyright remains with the author. In opting for open access, the author(s) should agree to publish the article under the CC BY license (Creative Commons Attribution 4.0 International License). The CC BY license allows for maximum dissemination and re-use of open access materials and is preferred by many research funding bodies. Under this license, users are free to share (copy, distribute and transmit) and remix (adapt) the contribution, including for commercial purposes, providing they attribute the contribution in the manner specified by the author or licensor.

Home /

Archives /

Civil and Energy Research, Volume 1, Issue 1, 2025 /

Original Research

Study on structural response under vibration load with multiple degrees of freedom

Abstract

How to Cite

Issue

Section

License

Follow Luminescience

Further Information

Home / Archives / Civil and Energy Research, Volume 1, Issue 1, 2025 / Original Research

Study on structural response under vibration load with multiple degrees of freedom

Abstract

How to Cite

Issue

Section

License

Follow Luminescience

Further Information

Home /

Archives /

Civil and Energy Research, Volume 1, Issue 1, 2025 /

Original Research