# An Adaptive Model Order Reduction Method Based on the Damage Evolution for Nonlinear Seismic Analysis

Authors: Jian Wang; Ming Fang; Hui Li
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Nonlinear seismic analysis, an approach to evaluate the seismic performance of a structure, is facing the challenge of computational efficiency for large-scale and high-fidelity simulation. This paper proposes an adaptive model order reduction (MOR) method based on the damage evolution among the overall structure to alleviate the computational burden. The damage state of each component during seismic loadings is distinguished as the initial-elastic phase, the plastic-damage phase, and the residual-elastic phase. In order to exploit the potential of model order reduction based on the damage evolution, a duration spectrum analysis is utilized to evaluate the characteristics of the residual-elastic phase for SDOF systems with bilinear hysteretic behaviour. Thus, an adaptive MOR method has been proposed to handle the nonlinear dynamic analysis of structures during different damage evolution phases. The overall structure is adaptively partitioned into linear substructures and nonlinear substructures on the basis of the time-varying damage distribution. The model order of linear substructures is reduced using the initial stiffness-based vibration modes, while nonlinear substructures that keep in the residual-elastic phase are reduced using the tangent-stiffness-based vibration modes. The residual displacements of nonlinear substructures are treated as the initial deformation during the residual-elastic phase. Compared with the traditional time step integration method, the proposed adaptive MOR method is able to increase the computational efficiency as yielding comparative results.


---

## Body

## 1. Introduction

Taking advantage of the enormous growth of computational capacities, finite element methods have become the most popular simulation techniques for large and complex structures with high fidelity. However, the steep increase in the demand for complex nonlinearity and high resolution of large problems yields undesired facts of expensive time consuming. The model order reduction (MOR) methods attract significant attention on the capacity of balancing the accuracy and efficiency.The principal of MOR methods is to projecting the original model onto a reduced dimensional space. On the basis of the definition of the reduced dimensional space, the MOR methods mainly include the proper orthogonal decomposition (POD) [1–3], proper generalized decomposition (PGD) [4–6], reduced basis (RB) method [7–10], component mode synthesis (CMS) [11, 12], and machine learning approaches [13, 14]. Such projection-based methods have been widely applied in linear systems and abbreviate an enormous computational burden. Recently, these projection-based methods are encouraged to cope with large and complex structures in nonlinear analysis.The POD constructs the reduced order basis in a posteriori manner by solving the original problem or its modified version. Based on the POD technique, some modified methods have been proposed to lower the dimension of nonlinear models, such as the GNAT [15, 16], Hyper-reduction [17, 18], or DEIM [19, 20]. Recently, the POD integrates domain decomposition to solve nonlinear structural problems subjected to dynamic loads [21]. The entire mechanical system is decoupled into a group of subdomains and POD is applied in the elastic subdomains. Nonetheless, once a subdomain enters the plastic stage, its nonlinear behaviour is simulated by the original nonlinear model and the computational saving vanishes. The PGD, which is more recently developed in past decade, generates the reduced order basis in an a priori manner by a powerful iteration algorithm during the resolution of the problem. When the reduced order basis changes often in nonlinear problems, aforementioned techniques such as DEIM become less suitable in the context of PGD than POD. To tackle this issue, PGD has integrated the Reference Point Method to make model order reduction [22]. Contrary to an interpolation technique, this strategy utilizes an approximation technique based on the integrals for the construction of the reduced order model. And the cost of this stage is decreased by one order of magnitude to make it represent less than 10% of the total computational time. The RB method utilizes an offline-online framework to make model order reduction. The reduced basis functions are selected in the offline stage, and the reduced basis coefficients are obtained in the online stage. The RB method has been developed to handle nonlinear issues, such as Static Condensation Reduced Basis Element method [23–25] and the Reduced Basis Element Method [26, 27]. The RB method recently integrated the Gaussian processes’ regression method (GPR) in the online stage to solve nonlinear problems [28]. The GPR method builds a bridge between the projection coefficients and the system parameters and utilizes the reduced basis functions to construct the solution. Nonetheless, this strategy is only effective when applying to a large-scale structure with a predefined linear-nonlinear interface. Driven by the rise of machine learning methods, regression-based MOR methods have been developed to solve nonlinear problems, such as the GPR [29] and artificial neural networks [14]. This approach takes advantages of the data-driven features and utilizes the input-output regression process to decouple the ofﬂine-online stage. Thus, the solution of nonlinear problems is constructed in the online stage from system parameters and reduced basis functions obtained in the offline stage.In nonlinear seismic analysis, the distribution of damaged structural components is usually spatially sparse. Thus, the substructure-based MOR method has a potential to alleviate the computational burden because it is flexible to simulate sparsely distributed damage through a substructure modelling strategy. The CMS method is a very popular substructure-based projection method and the strategy is to retain a set of substructure modes to simulate the dynamic response of the entire structure. Recently, the substructure-based MOR method has been applied to improve the nonlinear structural analysis, such as the landing simulation of the nonlinear aircraft system [30], the soil-structure interaction simulation [31], the reinforced concrete (RC) frame response simulation [32], and the crack and damage simulation [33]. When applying the CMS method to nonlinear seismic response analysis, two issues arise from the process of model order reduction. Firstly, the distribution of damaged structural components cannot be a priori known due to the randomness of seismic excitations. Thus, the substructure interface cannot be previously determined and should change automatically according to the time-varying nonlinearities. Secondly, the increasing ground motion intensity results in the spatially expanding distribution of damaged structural components, which leads to low level of the model order reduction for linear substructures. A nonlinear structural simulation strategy that integrates the adaptive substructure modelling with the nonlinear model order reduction is needed to tackle these issues. However, the abovementioned investigations neglect either adaptive substructure modelling or nonlinear model order reduction, or both.An adaptive modified Craig–Bampton (AMCB) method [34] has been previously proposed, which is inspired by the ideas of adaptive tracking control methods for nonlinear systems [35–38]. The AMCE method incorporated adaptive substructure partitioning technique and the CMS method together for the nonlinear seismic analysis of tall buildings. During a nonlinear seismic analysis, the AMCB method is able to perform model order reduction onto the time-varying linear substructures for tall buildings in which the damage distribution is a priori unknown. However, the increasing earthquake intensity declines the number of linear substructures, which declines computational efficiency. Actually, structural components may undergo various damage phases and the damage state of each component is time-varying during seismic loadings. There is a great potential to perform deeper model order reduction according to the damage evolution among the overall structure.The objective of this paper is to propose and validate an efficient nonlinear seismic analysis method with the nonlinear model order reduction capacity based on the damage evolution. The overall structure is adaptively partitioned into linear substructures and nonlinear substructures on the basis of the time-varying damage distribution. The damage state of each component during seismic loadings is classified as the initial-elastic phase, the plastic-damage phase, and the residual-elastic phase. The simulation speed has been accelerated during the residual-elastic phase through the model order reduction for nonlinear substructures. From the duration spectrum analysis, the duration of the residual-elastic phase occupies a major proportion of the duration of ground motions for common civil structures. Therefore, the proposed method has a great advantage for the efficient nonlinear seismic analysis of civil structures.Section2 introduces the damage evolution of structural components during seismic excitations and the characteristics of the residual-elastic phase are studied through the duration spectrum analysis. In Section 3, the adaptive MOR method is extended to consider the damage evolution phases. The governing equations are formulated in different damage states for the nonlinear structure. In Section 4, a RC frame structure is utilized to test the numerical performance of the proposed adaptive MOR method. The conclusions are summarized in Section 5.


## 2. Damage Evolution and Duration Characteristic

At the beginning of a seismic loading period, all structural components keep elastic as a result of the low acceleration amplitudes. This phase is considered as the initial-elastic phase. Afterwards, the seismic loading intensity increases and the structure deformation becomes larger. At severe ground shakings, the structural components may undergo an inelastic phase through yielding plastic strain. The plastic deformation causes structural damage. This phase is considered as the plastic-damage phase. The plastic-damage phase continues until the seismic loading intensity drops back to a low level. After severe shakings, structural components enter a damage phase in which the residual deformation remains without incremental plastic strain. Structural components behave in an elastic manner with a residual deformation. This phase is considered as the residual-elastic phase. The plastic-damage phase and residual-elastic phase alternate until the end of seismic loading. The plastic deformation has been utilized to define damage model for structures [39–41]. Herein, the aforementioned damage phases are represented by the plastic strain εp and the plastic strain rate ε˙p, respectively. Three damage phases are shown in Figure 1.Figure 1
Schematic diagram of the residual-elastic phase.During the residual-elastic phase, structural components retain residual deformation and behave in an elastic manner. Thus, it is feasible to make model order reduction using the tangent-stiffness-based vibration modes to avoid the iteration integral. If the duration of the residual-elastic phase lasts for a sufficiently long period, the model order reduction for these components is worthwhile and beneficial to the nonlinear seismic analysis.


### 2.1. Duration Definition

To evaluate the characteristics of residual-elastic phase, this investigation puts forward a definition of the duration of the residual-elastic phase as follows:(1)Der=1Ds∫t0t1stdt,st=0,for:εp≠0,ε˙p≠0,st=1,for:εp≠0,ε˙p=0,where Der is the duration of residual-elastic phase. The subscript “er” denotes the “residual-elastic phase.” Ds is the significant duration of ground motion [42] and herein is defined as the time interval between the 5% and 95% of the root-mean-square acceleration arms, which is shown in equation (2). t0 and t1 are the moments of reaching 5% and 95% of arms, respectively. st represents the residual-elastic phase.(2)arms=1tE∫0tEa2tdt,where tE is the total duration of the seismic event and at is the acceleration time history.


### 2.2. Duration Spectrum Analysis

The duration of the residual-elastic phase is defined according to the summation of time intervals in the residual-elastic phase compared to the ground motion duration. The ground motion duration may be a key influencing factor since the duration of the residual-elastic phase is also affected by the damage extent of structural components. The damage extent depends on the seismic intensity, period of vibration, damping ratio, and postyield stiffness. All these mentioned influencing factors should be considered. To assess the duration of the residual-elastic phase, this investigation proposes the duration spectrum of the residual-elastic phase, which is expressed byDer versus the natural period of structures. This duration spectrum can directly account for the aforementioned influence factors.The duration spectrum of the residual-elastic phase is computed for SDOF systems with bilinear hysteretic behaviour. The strength reduction factor [43], which is shown in equation (3), is usually used to measure the seismic intensity relative to the lateral strength of the system:(3)R=mSaFy,where m is the mass of the system, Sa is the seismic spectral acceleration, and Fy is the lateral yield strength of the system.To compute the duration spectrum of the residual-elastic phase, a group of earthquake records for nonlinear time-history analyses are selected to be representative of different ground motion duration according to [44]. The ground motion duration index ID [45] is employed to classify the ground motions into three sets of 20 records, i.e., the set of small ground motion duration (ID≈5), the set of moderate ground motion duration (ID≈14), and the set of large ground motion duration (ID≈22). The reason for using ID is that it has been proven to be an acceptable ground motion duration index for predicting the plastic damage demand of structural components.Figure2 shows the mean duration spectrum of the residual-elastic phase corresponding to three sets of earthquakes with different ground motion durations. Herein, the vibration periods of the SDOF system range from 0.1 s to 6.0 s, and four strength reduction factors R = 2, 4, 6, and 8 are considered for possible structure strengths. The damping ratio ξ=0.05 and postyield stiffness ratio α=0.01. It can be seen that, in general, the strength reduction factor R has a dominant influence on the duration spectrum. When R is smaller, the duration spectrum value is larger. This is because a small R corresponds to the high strength of the SDOF system and yields less plastic damage to the structure. Therefore, the duration of the residual-elastic phase will be long. On the contrary, for a constant R value, the duration spectrum generally shows the same trend regardless of the ground motion duration. This phenomenon indicates that the ground motion duration cannot account for the characteristics of the residual-elastic phase of structures, and the duration of the residual-elastic phase is consistent for seismic excitations with different ground motion durations. To give a specific description, the duration spectrum of the residual-elastic phase is strongly dependent on the period of vibration. The duration spectrum of the residual-elastic phase increases sharply in the short period region and declines mildly in the medium-long period region. Throughout the entire period region, there is a peak in the duration spectrum that appears at approximately 0.3 s, while the duration spectrum is larger than 0.6 for most of the period region.Figure 2
Duration spectrum for different ground motion durations. (a)ID = 22. (b) ID = 14. (c) ID = 5. (d) Ratio.
(a)(b)(c)(d)Figure3 shows the mean duration spectrum of the residual-elastic phase corresponding to three damping ratios ξ = 0.02, 0.05, and 0.10 for a total of 60 earthquake records. It can be seen that the duration spectrum decreases as the damping ratio increases. This is because the large damping ratio will suppress the structure deformation, so less plastic damage is generated throughout the seismic excitation. However, the damping ratios have a mild effect on the duration spectrum of a residual-elastic phase throughout the entire period region.Figure 3
Duration spectrum for different damping ratios. (a)ξ = 0.02. (b) ξ = 0.05. (c) ξ = 0.10. (d) Ratio.
(a)(b)(c)(d)Figure4 shows the mean duration spectrum of the residual-elastic phase corresponding to three postyield stiffness ratios α = 0.01, 0.05, and 0.10 for a total of 60 earthquake records. It can be seen that the duration spectrum decreases as the postyield stiffness ratio increases. This is because the large postyield stiffness ratio corresponds to less plastic damage of the SDOF system. However, the postyield stiffness ratios have little effect on the duration spectrum of a residual-elastic phase throughout the entire period region, which is assumed to be neglected.Figure 4
Duration spectrum for different postyield stiffness ratios. (a)α = 0.01. (b) α = 0.05. (c) α = 0.10. (d) Ratio.
(a)(b)(c)(d)Above all, the duration spectrum of a residual-elastic phase is dominated by the strength reduction factorR, while the ground motion duration ID has a minor effect on the duration spectrum. The effect of the damping ratio ξ is to lower the duration spectrum, and the postyield stiffness α has a negligible influence. The duration spectrum of the residual-elastic phase is validated to be larger than 0.6 when the natural period of the structures is greater than 0.3 s. This fact facilitates the model order reduction of the structural components using the tangent-stiffness-based vibration modes during the residual-elastic period.


## 2.1. Duration Definition

To evaluate the characteristics of residual-elastic phase, this investigation puts forward a definition of the duration of the residual-elastic phase as follows:(1)Der=1Ds∫t0t1stdt,st=0,for:εp≠0,ε˙p≠0,st=1,for:εp≠0,ε˙p=0,where Der is the duration of residual-elastic phase. The subscript “er” denotes the “residual-elastic phase.” Ds is the significant duration of ground motion [42] and herein is defined as the time interval between the 5% and 95% of the root-mean-square acceleration arms, which is shown in equation (2). t0 and t1 are the moments of reaching 5% and 95% of arms, respectively. st represents the residual-elastic phase.(2)arms=1tE∫0tEa2tdt,where tE is the total duration of the seismic event and at is the acceleration time history.


## 2.2. Duration Spectrum Analysis

The duration of the residual-elastic phase is defined according to the summation of time intervals in the residual-elastic phase compared to the ground motion duration. The ground motion duration may be a key influencing factor since the duration of the residual-elastic phase is also affected by the damage extent of structural components. The damage extent depends on the seismic intensity, period of vibration, damping ratio, and postyield stiffness. All these mentioned influencing factors should be considered. To assess the duration of the residual-elastic phase, this investigation proposes the duration spectrum of the residual-elastic phase, which is expressed byDer versus the natural period of structures. This duration spectrum can directly account for the aforementioned influence factors.The duration spectrum of the residual-elastic phase is computed for SDOF systems with bilinear hysteretic behaviour. The strength reduction factor [43], which is shown in equation (3), is usually used to measure the seismic intensity relative to the lateral strength of the system:(3)R=mSaFy,where m is the mass of the system, Sa is the seismic spectral acceleration, and Fy is the lateral yield strength of the system.To compute the duration spectrum of the residual-elastic phase, a group of earthquake records for nonlinear time-history analyses are selected to be representative of different ground motion duration according to [44]. The ground motion duration index ID [45] is employed to classify the ground motions into three sets of 20 records, i.e., the set of small ground motion duration (ID≈5), the set of moderate ground motion duration (ID≈14), and the set of large ground motion duration (ID≈22). The reason for using ID is that it has been proven to be an acceptable ground motion duration index for predicting the plastic damage demand of structural components.Figure2 shows the mean duration spectrum of the residual-elastic phase corresponding to three sets of earthquakes with different ground motion durations. Herein, the vibration periods of the SDOF system range from 0.1 s to 6.0 s, and four strength reduction factors R = 2, 4, 6, and 8 are considered for possible structure strengths. The damping ratio ξ=0.05 and postyield stiffness ratio α=0.01. It can be seen that, in general, the strength reduction factor R has a dominant influence on the duration spectrum. When R is smaller, the duration spectrum value is larger. This is because a small R corresponds to the high strength of the SDOF system and yields less plastic damage to the structure. Therefore, the duration of the residual-elastic phase will be long. On the contrary, for a constant R value, the duration spectrum generally shows the same trend regardless of the ground motion duration. This phenomenon indicates that the ground motion duration cannot account for the characteristics of the residual-elastic phase of structures, and the duration of the residual-elastic phase is consistent for seismic excitations with different ground motion durations. To give a specific description, the duration spectrum of the residual-elastic phase is strongly dependent on the period of vibration. The duration spectrum of the residual-elastic phase increases sharply in the short period region and declines mildly in the medium-long period region. Throughout the entire period region, there is a peak in the duration spectrum that appears at approximately 0.3 s, while the duration spectrum is larger than 0.6 for most of the period region.Figure 2
Duration spectrum for different ground motion durations. (a)ID = 22. (b) ID = 14. (c) ID = 5. (d) Ratio.
(a)(b)(c)(d)Figure3 shows the mean duration spectrum of the residual-elastic phase corresponding to three damping ratios ξ = 0.02, 0.05, and 0.10 for a total of 60 earthquake records. It can be seen that the duration spectrum decreases as the damping ratio increases. This is because the large damping ratio will suppress the structure deformation, so less plastic damage is generated throughout the seismic excitation. However, the damping ratios have a mild effect on the duration spectrum of a residual-elastic phase throughout the entire period region.Figure 3
Duration spectrum for different damping ratios. (a)ξ = 0.02. (b) ξ = 0.05. (c) ξ = 0.10. (d) Ratio.
(a)(b)(c)(d)Figure4 shows the mean duration spectrum of the residual-elastic phase corresponding to three postyield stiffness ratios α = 0.01, 0.05, and 0.10 for a total of 60 earthquake records. It can be seen that the duration spectrum decreases as the postyield stiffness ratio increases. This is because the large postyield stiffness ratio corresponds to less plastic damage of the SDOF system. However, the postyield stiffness ratios have little effect on the duration spectrum of a residual-elastic phase throughout the entire period region, which is assumed to be neglected.Figure 4
Duration spectrum for different postyield stiffness ratios. (a)α = 0.01. (b) α = 0.05. (c) α = 0.10. (d) Ratio.
(a)(b)(c)(d)Above all, the duration spectrum of a residual-elastic phase is dominated by the strength reduction factorR, while the ground motion duration ID has a minor effect on the duration spectrum. The effect of the damping ratio ξ is to lower the duration spectrum, and the postyield stiffness α has a negligible influence. The duration spectrum of the residual-elastic phase is validated to be larger than 0.6 when the natural period of the structures is greater than 0.3 s. This fact facilitates the model order reduction of the structural components using the tangent-stiffness-based vibration modes during the residual-elastic period.


## 3. Adaptive MOR Method Based on Damage Evolution

It has been demonstrated that the damaged structural components will undergo a residual-elastic phase after the severe ground shakings. The duration spectrum value of the residual-elastic phase is larger than 0.6 when natural period is greater than 0.3, which means the duration of the residual-elastic phase is sufficiently long for model order reduction. Therefore, the previously proposed AMCB method is extended by conducting a further model order reduction in the residual-elastic phase. In the AMCB method, linear and nonlinear substructures are automatically modelled on the basis of the time-varying damage distribution. The DOFs of the linear substructures are condensed by a modified CB method, while the nonlinear substructures are maintained in the original form. A small set of vibrational modes of linear substructures are retained to simulate the dynamic response of the structure. Herein, this newly developed method is named as an adaptive model order reduction method based on the residual-elastic phase (AMCB-E) method. Figure5 illustrates the flowchart of the AMCB-E method. The corresponding MOR methods are applied for each different damage state.Figure 5
Flowchart of the AMCB-E method.The nonlinear governing equation of a lumped-mass structure is formulated as follows:(4)Mgu¨g+Cgu˙g+Kgug+Rug,u˙g=fg,where Mg, Cg, and Kg represent the mass matrix, the initial damping matrix, and the initial stiffness matrix, respectively. ug and fg are the displacement vector and external force vector, respectively. Rug,u˙g is the nonlinear restoring force vector. The subscript “g” denotes the degrees of freedom (DOFs) of the global structure.It is noted that the stiffness matrix of a nonlinear structure is composed of elastic stiffness elements and plastic stiffness elements. The whole structure is initially elastic, and the initial stiffness matrixKg is composed of elastic stiffness elements. While the nonlinear restoring force vector Rug,u˙g is composed of plastic stiffness elements. No conflicts exist between the definition of Kg and Rug,u˙g, and such a stiffness matrix partition facilitates the formulation of nonlinear governing equations based on substructures in the following context.


### 3.1. Initial-Elastic Phase

The governing equation in the initial-elastic phase is formulated by modal coordinates as follows:(5)Mdeq¨de+Cdeq˙de+Kdeqde=fde,where Mde=ΦdTMgΦd, Cde=ΦdTCgΦd, Kde=ΦdTKgΦd, and fde=ΦdTfg are the mass matrix, the damping matrix, the stiffness matrix, and the external load vector, respectively. Φd is the reduced basis selected from vibration modes. qde is the corresponding modal coordinates. The superscript “e” represents the initial-elastic phase, and the subscript “d” represents the dominant vibration modes.Normally, a small set of vibration modes are sufficient to simulate the seismic response of the structure. The model order can be substantially reduced for large-scale structures. By solving such a reduced governing equation (5), the nonlinear seismic analysis is remarkably accelerated during the initial-elastic phase.


### 3.2. Plastic-Damage Phase

The AMCB method can be used to make model order reduction in the plastic-damage phase. The displacement vectorug is divided and formulated as ug=ulubunT. The subscripts “l,” “b,” and “n” represent the linear substructures, the linear-nonlinear interface, and the nonlinear substructures, respectively. Equation (4) can be written as follows:(6)Mll000Mbb000Mnnu¨lu¨bu¨n+CllClb0CblCbbCbn0CnbCnnu˙lu˙bu˙n+KllKlb0KblKbbKbn0KnbKnnulubun+00Rnun,u˙n=flfbfn.The model order reduction is made through coordinate space transformation as follows:(7)ulubun≈Φ¯dΨ¯b00Ib000Inqdpubpunp,where Φ¯d represents the dominant vibration modes and Ψ¯b represents constraint modes. qdp is the corresponding coordinates for linear substructures in modal space. The superscript “p” represents the plastic-damage phase. The governing equation is written in the hybrid coordinate space, which composes of physical and modal coordinates:(8)Mhpu¨hp+Chpu˙hp+Khpuhp+Rhp=fhp,where(9)Mhp=Φ¯dTMllΦ¯dΦ¯dTMllΨ¯b0Ψ¯bTMllΦ¯dΨ¯bTMllΨ¯b+Mbb000Mnn,Khp=Φ¯dTKllΦ¯d000Kbb−KblKll−1KlbKbn0KnbKnn,Rhp=00Rnun,u˙n,uhp=qdpubpunp,Chp=Φ¯dTCllΦ¯dΦ¯dTCllΨ¯b+Clb0Ψ¯bTCll+CblΦ¯dΨ¯bTCllΨ¯b+CblΨ¯b+Ψ¯bTClb+CbbCbn0CnbCnn,fhp=Φ¯dTflΨ¯bTfl+fbfn,where uhp is the hybrid displacement vector. The subscript “h” represents the hybrid coordinates.The model order reduction is only made for linear substructures in the governing equation (8). When the seismic loading intensity is low, such a MOR method will achieve enormous computational savings through the projection of a large number of linear DOFs onto a reduced basis. However, when the seismic intensity increases, the expanding damage distribution results in a decreasing number of linear DOFs, which weakens the MOR effect for linear substructures. Therefore, the MOR effect of the entire structure will steeply decline. One way to break through this bottleneck is to make model order reduction in the residual-elastic phase.


### 3.3. Residual-Elastic Phase

In the residual-elastic phase, the subvectorun is calculated by summarizing the residual displacement and an incremental displacement Δun. To make model order reduction, a coordinate space transformation is conducted for the incremental displacement vector of the structure. The displacement vectors of both the linear substructures and the nonlinear substructures are transformed from physical space to modal space, while only the interface remains in the physical space:(10)ΔulΔubΔun≈Φ¯dΨ¯b00Ibb00Ψ˜bΦ˜dΔqlerΔuberΔqner,where Φ˜d is the dominant tangent-stiffness-based vibration modes and Ψ˜b is the constraint modes. Δqler and Δqner are the incremental coordinates for linear and nonlinear substructures, respectively. The governing equation is formulated in incremental form as follows:(11)MherΔu¨her+CherΔu˙her+KherΔuher=Δfher,where(12)Mher=Φ¯dTMllΦ¯dΦ¯dTMllΨ¯b0Ψ¯bTMllΦ¯dΨ¯bTMllΨ¯b+Mbb+Ψ˜bTMnnΨ˜bΨ˜bTMnnΦ˜d0Φ˜dTMnnΨ˜bΦ˜dTMnnΦ˜d,Δfher=Φ¯dTΔflΨ¯bTΔfl+Δfb+Ψ˜bTΔfnΦ˜dTΔfn,Kher=Φ¯dTKllΦ¯d000Kbb−KblKll−1Klb−KbnKnn−1Knb000Φ˜dTKnnΦ˜d,Δuher=ΔqderΔuberΔuner,Cher=Φ¯dTCllΦ¯dΦ¯dTCllΨ¯b+Clb0Ψ¯bTCll+CblΦ¯dΨ¯bTCllΨ¯b+CblΨ¯b+Ψ¯bTClb+Cbb+Ψ˜bTCnnΨ˜b+CbnΨ˜b+Ψ˜bTCnbΨ˜bTCnn+CbnΦ˜d0Φ˜dTCnnΨ˜b+CnbΦ˜dTCnnΦ˜d,where Δuher is the incremental displacement vector in hybrid coordinate space. Equation (11) is solved for the structure during the residual-elastic phase until the structure turn backs to the plastic-damage phase.In the governing equation (11), the model order reduction has been made for both the linear substructures and nonlinear substructures. When the seismic loading intensity drops down to a low level after severe ground shakings, such an AMCB-E method will achieve a further DOF reduction than the AMCB method. The longer the duration of the residual-elastic phase is, the more efficient the AMCB-E method can be.


## 3.1. Initial-Elastic Phase

The governing equation in the initial-elastic phase is formulated by modal coordinates as follows:(5)Mdeq¨de+Cdeq˙de+Kdeqde=fde,where Mde=ΦdTMgΦd, Cde=ΦdTCgΦd, Kde=ΦdTKgΦd, and fde=ΦdTfg are the mass matrix, the damping matrix, the stiffness matrix, and the external load vector, respectively. Φd is the reduced basis selected from vibration modes. qde is the corresponding modal coordinates. The superscript “e” represents the initial-elastic phase, and the subscript “d” represents the dominant vibration modes.Normally, a small set of vibration modes are sufficient to simulate the seismic response of the structure. The model order can be substantially reduced for large-scale structures. By solving such a reduced governing equation (5), the nonlinear seismic analysis is remarkably accelerated during the initial-elastic phase.


## 3.2. Plastic-Damage Phase

The AMCB method can be used to make model order reduction in the plastic-damage phase. The displacement vectorug is divided and formulated as ug=ulubunT. The subscripts “l,” “b,” and “n” represent the linear substructures, the linear-nonlinear interface, and the nonlinear substructures, respectively. Equation (4) can be written as follows:(6)Mll000Mbb000Mnnu¨lu¨bu¨n+CllClb0CblCbbCbn0CnbCnnu˙lu˙bu˙n+KllKlb0KblKbbKbn0KnbKnnulubun+00Rnun,u˙n=flfbfn.The model order reduction is made through coordinate space transformation as follows:(7)ulubun≈Φ¯dΨ¯b00Ib000Inqdpubpunp,where Φ¯d represents the dominant vibration modes and Ψ¯b represents constraint modes. qdp is the corresponding coordinates for linear substructures in modal space. The superscript “p” represents the plastic-damage phase. The governing equation is written in the hybrid coordinate space, which composes of physical and modal coordinates:(8)Mhpu¨hp+Chpu˙hp+Khpuhp+Rhp=fhp,where(9)Mhp=Φ¯dTMllΦ¯dΦ¯dTMllΨ¯b0Ψ¯bTMllΦ¯dΨ¯bTMllΨ¯b+Mbb000Mnn,Khp=Φ¯dTKllΦ¯d000Kbb−KblKll−1KlbKbn0KnbKnn,Rhp=00Rnun,u˙n,uhp=qdpubpunp,Chp=Φ¯dTCllΦ¯dΦ¯dTCllΨ¯b+Clb0Ψ¯bTCll+CblΦ¯dΨ¯bTCllΨ¯b+CblΨ¯b+Ψ¯bTClb+CbbCbn0CnbCnn,fhp=Φ¯dTflΨ¯bTfl+fbfn,where uhp is the hybrid displacement vector. The subscript “h” represents the hybrid coordinates.The model order reduction is only made for linear substructures in the governing equation (8). When the seismic loading intensity is low, such a MOR method will achieve enormous computational savings through the projection of a large number of linear DOFs onto a reduced basis. However, when the seismic intensity increases, the expanding damage distribution results in a decreasing number of linear DOFs, which weakens the MOR effect for linear substructures. Therefore, the MOR effect of the entire structure will steeply decline. One way to break through this bottleneck is to make model order reduction in the residual-elastic phase.


## 3.3. Residual-Elastic Phase

In the residual-elastic phase, the subvectorun is calculated by summarizing the residual displacement and an incremental displacement Δun. To make model order reduction, a coordinate space transformation is conducted for the incremental displacement vector of the structure. The displacement vectors of both the linear substructures and the nonlinear substructures are transformed from physical space to modal space, while only the interface remains in the physical space:(10)ΔulΔubΔun≈Φ¯dΨ¯b00Ibb00Ψ˜bΦ˜dΔqlerΔuberΔqner,where Φ˜d is the dominant tangent-stiffness-based vibration modes and Ψ˜b is the constraint modes. Δqler and Δqner are the incremental coordinates for linear and nonlinear substructures, respectively. The governing equation is formulated in incremental form as follows:(11)MherΔu¨her+CherΔu˙her+KherΔuher=Δfher,where(12)Mher=Φ¯dTMllΦ¯dΦ¯dTMllΨ¯b0Ψ¯bTMllΦ¯dΨ¯bTMllΨ¯b+Mbb+Ψ˜bTMnnΨ˜bΨ˜bTMnnΦ˜d0Φ˜dTMnnΨ˜bΦ˜dTMnnΦ˜d,Δfher=Φ¯dTΔflΨ¯bTΔfl+Δfb+Ψ˜bTΔfnΦ˜dTΔfn,Kher=Φ¯dTKllΦ¯d000Kbb−KblKll−1Klb−KbnKnn−1Knb000Φ˜dTKnnΦ˜d,Δuher=ΔqderΔuberΔuner,Cher=Φ¯dTCllΦ¯dΦ¯dTCllΨ¯b+Clb0Ψ¯bTCll+CblΦ¯dΨ¯bTCllΨ¯b+CblΨ¯b+Ψ¯bTClb+Cbb+Ψ˜bTCnnΨ˜b+CbnΨ˜b+Ψ˜bTCnbΨ˜bTCnn+CbnΦ˜d0Φ˜dTCnnΨ˜b+CnbΦ˜dTCnnΦ˜d,where Δuher is the incremental displacement vector in hybrid coordinate space. Equation (11) is solved for the structure during the residual-elastic phase until the structure turn backs to the plastic-damage phase.In the governing equation (11), the model order reduction has been made for both the linear substructures and nonlinear substructures. When the seismic loading intensity drops down to a low level after severe ground shakings, such an AMCB-E method will achieve a further DOF reduction than the AMCB method. The longer the duration of the residual-elastic phase is, the more efficient the AMCB-E method can be.


## 4. Numerical Example

The AMCB-E method has been implemented in the finite element program OpenSees [46, 47].The implementation of the AMCB-E method requires program modifications over class inheritance at both the low level and the high level of OpenSees framework. At the low level of OpenSees framework, a new abstract layer is created to represent the time-varying substructures and implement the modelling process of substructures. At the high level of OpenSees framework, a new transient analysis class is created to implement the solving process of substructure-based governing equations. The implementation of the AMCB-E method has been presented with details in a recent publication [48].


### 4.1. Description of the Test Structure

A 12-story RC frame structure is employed to test the performance of the AMCB-E method through nonlinear seismic analysis, as shown in Figure6.Figure 6
Elevation view and section information.The 2D model has been developed in the OpenSees program. The displacement-based beam-column elements are selected to simulate the beams and columns and each component is discretized with 5 elements. The component cross-sections are divided into 50 fibers. The uniaxial Kent–Scott–Park concrete model is employed to model the concrete fiber for the degraded linear uploading/reloading stiffness and an assigned tensile strength. The compressive strength and strain of a cover concrete fiber is 30 MPa and 0.002. For a core concrete fiber, the compressive strength and strain are calculated according to the Kent–Part confined concrete model. On the other side, the uniaxial bilinear model is employed to simulate the steel rebar fibers for the kinematic hardening. Young’s modulus is 200 GPa, the strain hardening ratio is 0.01, and the yield strength is 400 MPa. The top three vibration periods of the 2D model are 1.83 s, 0.63 s, and 0.35 s. The Rayleigh damping method has been used in the nonlinear dynamic analysis.The El Centro, Kobe, Northridge, and Parkfield earthquakes are employed as seismic excitation inputs. The peak ground accelerations (PGA) range from 100 gal to 500 gal with a spacing interval of 100 gal. The normalized time histories of the seismic input accelerations are illustrated in Figure7.Figure 7
The normalized seismic input time histories.The nonlinear seismic analysis results are used to assess the performance of the proposed adaptive MOR method. All numerical simulations are run on a workstation equipped with a CPU of Intel Xeon E5-2640 @ 2.4 GHz and RAM of 64 GB.


### 4.2. Seismic Response of the Test Structure

The simulation results of seismic response are compared between the traditional time step integration (TSI) method and AMCB-E method to validate the performance of the AMCB-E method. To quantify the differences between the global seismic response of the AMCB-E and TSI method, the story-drift ratio is firstly compared and the average relative errors of numerical results are listed in Table1. The maximum errors of the story-drift ratio for the four seismic inputs are approximately 2.76%, 5.15%, 3.47%, and 3.96%, respectively.Table 1
Story-drift ratio errors (%).
Ground motion intensity (gal)Seismic input recordsEl CentroKobeNorthridgeParkfield1001.702.733.472.262002.765.151.463.963001.593.231.462.574002.600.881.122.475002.071.131.641.62The time histories of roof displacement are shown in Figure8, in which the PGA is 500 gal. In addition, Figure 9 plots the moment-rotation curves of the end section of a beam element, in which the PGA is also 500 gal. The simulation results calculated by the AMCB-E and TSI method are comparative with each other. The simulation error is mainly derived from two reasons. Firstly, a half of the yield strength for the curved uniaxial concrete model is considered to trigger the switch of the elastic state and plastic state, which may introduce minor simulation errors. Secondly, the response of substructures is described by a small set of low-order modes, and the ignorance of high-order modes may introduce additional simulation errors. However, the simulation errors presented in the table and figures are quite small and acceptable.Figure 8
Time histories of roof displacements. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)Figure 9
Moment-rotation curves. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)


### 4.3. Model Order Reduction

The MOR effect and computational time are compared between the AMCE method and the AMCB-E method to validate the performance of the AMCB-E method. The effectiveness of the AMCB-E method depends on the characteristics of the residual-elastic phase. When the duration of the residual-elastic phase is longer and the switching frequency of damage evolution phases is lower, the AMCB-E method may have a better model order reduction effect. The reason is that the long duration of the residual-elastic phase provides a large number of time intervals, in which the model order reduction for nonlinear substructures is available. The low switching frequency of the damage evolution phases limits the additional computational time required by modelling substructures and updating governing equations.Figure10 illustrates the time history curves of the damage evolution phase of the test structure. The “i,” “p,” and “e” stand for the three damage evolution phases, i.e., the initial-elastic phase, plastic-damage phase, and residual-elastic phase, respectively. The alternation of the plastic-damage phase and residual-elastic phase can be observed in the time history curves. As the seismic loading intensity increases, the duration of the residual-elastic phase decreases and the switching frequency increases. Since the duration of the residual-elastic phase for all seismic excitations are longer than half of the seismic loading period, the AMCB-E has a great potential to make a model order reduction based on the residual-elastic phase.Figure 10
Time history of the damage state. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)The total number of DOFs is an important index to access the model order reduction effect. This study compares the ratio of the reduced model DOFs to the nonreduced model DOFs between AMCB-E and AMCB. Thus, it is able to investigate the contribution of the model order reduction based on the residual-elastic phase. The DOF reduction ratio of the AMCB method and AMCB-E method is illustrated in Figure11. It can be observed that the AMCB-E makes a further model order reduction in the residual-elastic phase and generates a lower DOF reduction ratio. Because of the long duration of the residual-elastic phase and the low switching frequency, the model order reduction for all seismic excitations is substantially effective.Figure 11
Time histories of DOF reduction ratio. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)


### 4.4. Computational Time

The purpose of the model order reduction in the residual-elastic phase is to further minimize the number of DOFs in governing equations and to accelerate the solution of governing equations. Figure12 compares the simulation time between the AMCB-E and the AMCB. The simulation time ratio refers to the ratio of the simulation time of the MOR methods to the TSI method. When the seismic loading intensities are low, no apparent improvement has been gained from the AMCB to the AMCB-E. That is because the test structure remains in the initial-elastic phase. When the seismic loading intensities increases, the runtime of the AMCB-E for all seismic excitations declines. The reason is that the model order reduction is further conducted in the residual-elastic phase.Figure 12
Comparison of the simulation times. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)


## 4.1. Description of the Test Structure

A 12-story RC frame structure is employed to test the performance of the AMCB-E method through nonlinear seismic analysis, as shown in Figure6.Figure 6
Elevation view and section information.The 2D model has been developed in the OpenSees program. The displacement-based beam-column elements are selected to simulate the beams and columns and each component is discretized with 5 elements. The component cross-sections are divided into 50 fibers. The uniaxial Kent–Scott–Park concrete model is employed to model the concrete fiber for the degraded linear uploading/reloading stiffness and an assigned tensile strength. The compressive strength and strain of a cover concrete fiber is 30 MPa and 0.002. For a core concrete fiber, the compressive strength and strain are calculated according to the Kent–Part confined concrete model. On the other side, the uniaxial bilinear model is employed to simulate the steel rebar fibers for the kinematic hardening. Young’s modulus is 200 GPa, the strain hardening ratio is 0.01, and the yield strength is 400 MPa. The top three vibration periods of the 2D model are 1.83 s, 0.63 s, and 0.35 s. The Rayleigh damping method has been used in the nonlinear dynamic analysis.The El Centro, Kobe, Northridge, and Parkfield earthquakes are employed as seismic excitation inputs. The peak ground accelerations (PGA) range from 100 gal to 500 gal with a spacing interval of 100 gal. The normalized time histories of the seismic input accelerations are illustrated in Figure7.Figure 7
The normalized seismic input time histories.The nonlinear seismic analysis results are used to assess the performance of the proposed adaptive MOR method. All numerical simulations are run on a workstation equipped with a CPU of Intel Xeon E5-2640 @ 2.4 GHz and RAM of 64 GB.


## 4.2. Seismic Response of the Test Structure

The simulation results of seismic response are compared between the traditional time step integration (TSI) method and AMCB-E method to validate the performance of the AMCB-E method. To quantify the differences between the global seismic response of the AMCB-E and TSI method, the story-drift ratio is firstly compared and the average relative errors of numerical results are listed in Table1. The maximum errors of the story-drift ratio for the four seismic inputs are approximately 2.76%, 5.15%, 3.47%, and 3.96%, respectively.Table 1
Story-drift ratio errors (%).
Ground motion intensity (gal)Seismic input recordsEl CentroKobeNorthridgeParkfield1001.702.733.472.262002.765.151.463.963001.593.231.462.574002.600.881.122.475002.071.131.641.62The time histories of roof displacement are shown in Figure8, in which the PGA is 500 gal. In addition, Figure 9 plots the moment-rotation curves of the end section of a beam element, in which the PGA is also 500 gal. The simulation results calculated by the AMCB-E and TSI method are comparative with each other. The simulation error is mainly derived from two reasons. Firstly, a half of the yield strength for the curved uniaxial concrete model is considered to trigger the switch of the elastic state and plastic state, which may introduce minor simulation errors. Secondly, the response of substructures is described by a small set of low-order modes, and the ignorance of high-order modes may introduce additional simulation errors. However, the simulation errors presented in the table and figures are quite small and acceptable.Figure 8
Time histories of roof displacements. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)Figure 9
Moment-rotation curves. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)


## 4.3. Model Order Reduction

The MOR effect and computational time are compared between the AMCE method and the AMCB-E method to validate the performance of the AMCB-E method. The effectiveness of the AMCB-E method depends on the characteristics of the residual-elastic phase. When the duration of the residual-elastic phase is longer and the switching frequency of damage evolution phases is lower, the AMCB-E method may have a better model order reduction effect. The reason is that the long duration of the residual-elastic phase provides a large number of time intervals, in which the model order reduction for nonlinear substructures is available. The low switching frequency of the damage evolution phases limits the additional computational time required by modelling substructures and updating governing equations.Figure10 illustrates the time history curves of the damage evolution phase of the test structure. The “i,” “p,” and “e” stand for the three damage evolution phases, i.e., the initial-elastic phase, plastic-damage phase, and residual-elastic phase, respectively. The alternation of the plastic-damage phase and residual-elastic phase can be observed in the time history curves. As the seismic loading intensity increases, the duration of the residual-elastic phase decreases and the switching frequency increases. Since the duration of the residual-elastic phase for all seismic excitations are longer than half of the seismic loading period, the AMCB-E has a great potential to make a model order reduction based on the residual-elastic phase.Figure 10
Time history of the damage state. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)The total number of DOFs is an important index to access the model order reduction effect. This study compares the ratio of the reduced model DOFs to the nonreduced model DOFs between AMCB-E and AMCB. Thus, it is able to investigate the contribution of the model order reduction based on the residual-elastic phase. The DOF reduction ratio of the AMCB method and AMCB-E method is illustrated in Figure11. It can be observed that the AMCB-E makes a further model order reduction in the residual-elastic phase and generates a lower DOF reduction ratio. Because of the long duration of the residual-elastic phase and the low switching frequency, the model order reduction for all seismic excitations is substantially effective.Figure 11
Time histories of DOF reduction ratio. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)


## 4.4. Computational Time

The purpose of the model order reduction in the residual-elastic phase is to further minimize the number of DOFs in governing equations and to accelerate the solution of governing equations. Figure12 compares the simulation time between the AMCB-E and the AMCB. The simulation time ratio refers to the ratio of the simulation time of the MOR methods to the TSI method. When the seismic loading intensities are low, no apparent improvement has been gained from the AMCB to the AMCB-E. That is because the test structure remains in the initial-elastic phase. When the seismic loading intensities increases, the runtime of the AMCB-E for all seismic excitations declines. The reason is that the model order reduction is further conducted in the residual-elastic phase.Figure 12
Comparison of the simulation times. (a) El Centro. (b) Kobe. (c) Northridge. (d) Parkfield.
(a)(b)(c)(d)


## 5. Conclusions

In this paper, an adaptive MOR method based on the residual-elastic phase (AMCB-E) is developed, which can make further model order reduction for nonlinear structures in the residual-elastic phase. The nonlinear substructures in the residual-elastic phase are reduced using the tangent-stiffness-based vibration modes. The duration spectrum of the residual-elastic phase is proposed to study the characteristics of the residual-elastic phase for nonlinear structures under strong ground motions. It is concluded that the duration of the residual-elastic phase generally occupies more than 60% of the duration of ground motion when the structure vibration period is longer than 0.3 s. This fact facilitates further model order reduction for nonlinear substructures under strong seismic excitations. The performance of the AMCB-E method is tested using a 12-story reinforced concrete frame structure. Acceptable agreements of the global and local seismic response are achieved between the traditional time integration method and the AMCB-E method. Both the duration of the residual-elastic phase and switching number dominate the simulation efficiency of the AMCB-E method. The simulation efficiency of the AMCB-E method is demonstrated to be better than the AMCB method.The AMCB-E method has been validated to be an efficient nonlinear dynamic analysis method. The simulation time is saved during the residual-elastic phase through the model order reduction for nonlinear substructures. From the duration spectrum analysis, the duration of the residual-elastic phase occupies a major proportion of the duration of ground motions for common civil structures. Therefore, the AMCB-E method has a great advantage for the efficient nonlinear seismic analysis of civil structures. In the future work, the influence of the switching frequency should be eliminated to further mitigate the additional computational time. Moreover, the proper orthogonal decomposition (POD) method is encouraged to be applied to the model order reduction of linear and nonlinear substructures, instead of the component mode synthesis (CMS) method. The POD method can avoid the simulation errors introduced by the mode truncation of the CMS method. And the POD method is believed to increase the numerical accuracy of the AMCB and AMCB-E methods.



---

*Source: 8865255-2020-09-16.xml*