Research Article  Open Access
Francisco L. SilvaGonzález, Sonia E. Ruiz, Alejandro RodríguezCastellanos, "NonGaussian Stochastic Equivalent Linearization Method for Inelastic Nonlinear Systems with Softening Behaviour, under Seismic Ground Motions", Mathematical Problems in Engineering, vol. 2014, Article ID 539738, 16 pages, 2014. https://doi.org/10.1155/2014/539738
NonGaussian Stochastic Equivalent Linearization Method for Inelastic Nonlinear Systems with Softening Behaviour, under Seismic Ground Motions
Abstract
A nonGaussian stochastic equivalent linearization (NSEL) method for estimating the nonGaussian response of inelastic nonlinear structural systems subjected to seismic ground motions represented as nonstationary random processes is presented. Based on a model that represents the time evolution of the joint probability density function (PDF) of the structural response, mathematical expressions of equivalent linearization coefficients are derived. The displacement and velocity are assumed jointly Gaussian and the marginal PDF of the hysteretic component of the displacement is modeled by a mixed PDF which is Gaussian when the structural behavior is linear and turns into a bimodal PDF when the structural behavior is hysteretic. The proposed NSEL method is applied to calculate the response of hysteretic singledegreeoffreedom systems with different vibration periods and different design displacement ductility values. The results corresponding to the proposed method are compared with those calculated by means of Monte Carlo simulation, as well as by a Gaussian equivalent linearization method. It is verified that the NSEL approach proposed herein leads to maximum structural response standard deviations similar to those obtained with Monte Carlo technique. In addition, a brief discussion about the extension of the method to mutidegreeoffreedom systems is presented.
1. Introduction
The method of stochastic equivalent linearization (SEL) is one of the most common methods within the approximate approaches for stochastic dynamic analysis of nonlinear systems. The SEL method has proven to be a versatile and efficient technique from the computation point of view. A great impulse to the method was given by Atalik and Utku [1], who demonstrated that for Gaussian responses the calculation of linearization coefficients can be performed in a very simple way. However, an important deficiency of the method has been found when the response of systems with nonlinear inelastic behavior subjected to intense excitations has been considered as Gaussian [2–6]. The authors [7] have found that when the hypothesis related to the response of singledegreeoffreedom (SDOF) systems with hysteretic behavior has Gaussian distribution, the standard deviation of the displacement response can be underestimated in the order of 45%; consequently, the probability of failure widely deviates from the Monte Carlo simulation results, especially for systems with high displacement ductility demands. This is mainly due to the fact that a Gaussian probability distribution is assumed for all variables; however, the restoring force should lie on a finite region, implying that its probability density function (PDF) is nonGaussian. Several attempts have been made in the last years to overcome the problem mentioned above [8–10]. For example, Asano and Iwan [11] obtained alternate linearization results for a hysteretic system by using equations that were identical to the usual form on all subsets having a nonzero probability density in the nonlinear system but differed in the zero probability subsets. An ample review related to applications of statistical and equivalent linearization in the analysis of structure and mechanical nonlinear stochastic dynamic systems is found in [12].
A nonGaussian equivalent linearization criterion (NSEL) that uses mathematical expressions for the linearization coefficients is presented herein. In order to derive the nonGaussian equivalent linearization coefficients, a model appropriately representing the time evolution of the joint probability density function of the structural response of softening systems subjected to seismic ground motions is used. The parameters of the joint PDF depend on a function which measures the nonGaussianity of the response process.
The structural softening models are of particular interest in earthquake engineering because they are commonly used to analyze buildings made with softening materials such as steel, reinforced concrete or masonry, which present saturation of the restoring force when deformation increases [13]. The nonlinearities considered in the present study are only due to material behavior.
It is noticed that the NSEL methods mentioned in the first paragraph of this section do not take into account the time evolution character of the joint PDF of the structural response in the way that is formulated here.
The mathematical expressions corresponding to the nonGaussian equivalent linearization coefficients derived here are applied to calculate the response of hysteretic SDOF systems with different vibration periods and design displacement ductility demands, subjected to seismic ground motions that are representative of a nonstationary ergodic random process. The comparison of the results with those obtained by means of Monte Carlo simulation shows that the proposed approach predicts with acceptable accuracy the maximum standard deviation of the structural response of hysteretic SDOF systems.
An extension of the proposed method to multidegreeoffreedom (MDOF) systems is outlined at the end of the paper.
2. Equation of Motion
The equation of motion of hysteretic SDOF system studied here is where , , and represent the relative acceleration, velocity, and displacement, respectively. is the fraction of viscous critical damping; is the circular frequency of vibration of the system; is the mass; is the coefficient of viscous damping; is the initial stiffness. is the ratio between the postyielding stiffness and the initial stiffness; represents the ground acceleration; is the time and is the hysteretic component of displacement represented here by the BoucWen model [14, 15] which is defined by the following nonlinear differential equation: where , , , and are parameters that define the size and the shape of the hysteretic cycles, controls the smoothness of the transition from elastic to plastic structural behavior in such a way that the smaller , the smoother the transition, and corresponds to a bilinear model. and define the softening or the hardening of the system; the first case corresponds to while the second corresponds to . Casciati and Faravelli [16], recommend considering for steel and for reinforced concrete. In this study softening systems with smooth transition represented by are considered.
3. Stochastic Model of the Seismic Excitation
The seismic motion is modeled as a nonstationary Gaussian random process generated from a filtered white noise modulated in amplitude. The transfer function proposed by Clough and Penzien [17] is used here; therefore, the power spectral density of the seismic excitation is given by The square of the amplitude modulating function is [18] In practical applications, the amplitude of the power spectral density of the white noise and the filter parameters (, , , and ) are determined with a nonlinear fitting of the power spectral density estimated by means of the Fourier spectrum of the basic seismic ground motion representative of the excitation process. The parameters , , , , and of the modulating function are determined with a nonlinear fitting of the Arias intensity curve of the basic accelerogram [19]. In this study, the excitation process is modulated only in amplitude; to study the effect of the frequency modulation on the nonlinear systems response the reader is referred to [20–23].
In the analysis presented below, the seismic record obtained in Mexico City at the Ministry of Communications and Transportation (SCT) during the September 19, 1985 earthquake is used to obtain the parameters of the stochastic model of the seismic excitation. The amplitude of the spectral density of the white noise is cm^{2}/s^{3}, the filter parameters are , , , and , and those corresponding to the modulating function are cm^{2}/s^{4}, , , , and . The bimodal power spectral densities of the seismic motion and its fitted model are shown in Figure 1. The cumulative Arias intensity curve and its fitted model are in Figure 2.
4. Modeling the Joint Probability Density Function of the Structural Response
The NSEL approach presented here is based on an appropriate representation of the time evolution of the joint PDF of the structural seismic response of softening systems with smooth transition. Results from Monte Carlo simulation analyses have shown that the PDF of the displacement, , and of the velocity, , of the structural system is Gaussian, and their evolution in time is as shown in Figures 3 and 4, respectively. However, the probability density of the restoring force and of the displacement hysteretic component, , changes with time because it depends on the nonlinearity of the structural behavior. At the beginning and at the end of the excitation, when the structure is under low seismic intensities, tends to be Gaussian because the structure presents linear behavior; however, as the seismic intensity increases and the structural behavior becomes nonlinear (hysteretic), then becomes nonGaussian, as shown in Figure 5. It can be seen that when the system is vibrating within the nonlinear range equals its maximum value and its PDF presents concentrations at + and −. Figure 5 also shows that as the excitation intensity grows (up to s) and consequently the displacement ductility demand increases, turns into a bimodal PDF. It is observed that the support of is − , where the maximum value of is given by [13, 19] Considering the concentration of values of in the vicinity of its maximum value and the shape that such PDF adopts (see Figure 5 for s), the following mathematical model is used: where is a Gaussian density function with zero mean and standard deviation , is a Gaussian density function with mean (or ) and standard deviation , and is a weighting factor. The idea is that when the nonlinear behavior becomes important, the weighting of , functions should increase while that corresponding to should decrease. It is noticed that (6) uses the Gaussian density functions and and not Dirac’s pulses as has been proposed by other authors [9]. The advantage of using Gaussian density functions, instead of Dirac’s pulses, is that it is possible to model more appropriately the time evolution of the joint PDF and, as a consequence, the NSEL approach leads to results that better approximate the Monte Carlo results. Furthermore, it is possible to get mathematical expressions of the equivalent linearization coefficients.
Using the marginal PDFs , , and the information about the correlation structure of the random variables, it is possible to compute the joint PDF [24, 25]. Here, the joint PDF is obtained using the conditional PDF of and given , , as follows: Based on the simulation results, it is reasonable to assume that is jointly Gaussian. The parameters of , that is, conditional means and variances, are estimated assuming that the former are linear functions of , and the latter do not depend on , as it happens in a mean square estimation of normal variables [26]. Additionally, it is considered that the response process has zero mean.
Based on (6) and (7), the following expression for the joint PDF is obtained: where , , and are the standard deviations of , , and , respectively; , , and are the corresponding correlation coefficients; and .
4.1. Evaluating the Parameters , , , and
The four parameters , , , and of the joint PDF given by (8) should satisfy that the area under the PDF given by (6) is close to unity: where is the error function given by The above mentioned condition is necessary because (6) is not a bounded function. If (9) is satisfied, the second moment (in this case, equal to the variance ) can be approximated as The four parameters , , , and could be computed by means of their first four moments; however, here they are determined by means of nondimensional coefficients that depend on the system response nonGaussianity level which is represented by [13] From (12) it is observed that is close to zero when the response is Gaussian, that is, when the system strength is high (where is high) and the excitation or response is low (and as a consequence is low). On the other hand, if the response is nonGaussian, that is, if the system strength is low (where is low) and the excitation or response is high (and consequently is high), then is high.
Here, the following mathematical expressions for the normalized parameters , , , and are proposed: where , , , and , and 4, are constants to be determined. The general behavior of (13) is presented in Figures 6(a)–6(d), where it can be seen that and increase with (Figures 6(a) and 6(b)), meaning that when the nonlinearity of the response increases, the function moves far from and its weight increases; hence, (given by (6)) becomes a bimodal PDF. On the other hand, when the structural behavior is linear is close to zero and the weight of vanishes; therefore becomes Gaussian. It can also be seen (Figures 6(c) and 6(d)) that and decrease with , but it should be noticed that as increases, grows as well.
(a) versus
(b) versus
(c) versus
(d) versus
Equations (13) show that the parameters , , , and are time dependent; thus, the proposed approach considers the time evolutionary character of the joint PDF of the response, while other methods proposed in the literature do not take it into account. Figure 7 shows an example of the time evolution of given by (6).
4.1.1. Case of NarrowBand Seismic Inputs
Parameters , , , and , and 4 of (13) corresponding to the response of singledegreeoffreedom systems (SDOF) systems subjected to the action of narrowband seismic inputs, were calculated by means of Monte Carlo simulation analysis. The critical case, where the natural vibration period of SDOF system is equal to the dominant period of seismic input (in this case s, see Figure 1) and the SDOF system has a high design ductility demand , was considered. The design ductility factor is defined here as the ratio between the expected maximum and the yield displacement of the system. In this study, the first one is found by means of a stationary Gaussian equivalent linearization analysis, and the latter is estimated by means of the initial stiffness and the yield force of the system.
The SDOF system was subjected to the action of 50,000 artificial accelerograms based on the narrowband motion recorded in SCT on September 19, 1985 (see Section 3). For this case, it was found that lies within the interval and that leads to good results. The parameter values obtained from the analysis are shown in Table 1.

Points in Figures 6(a)–6(d) represent the results of Monte Carlo simulation analyses.
It is noticed that parameters in Table 1 are proposed for narrowband seismic inputs and are independent of any other parameter, such as structural vibration period and design ductility factor.
5. NonGaussian Equivalent Stochastic Linearization
In the equivalent linearization method, (2) is replaced by a linear differential equation as follows [5, 27]: where , , and are linearization coefficients. In order to minimize the expected value of the square of the difference of (2) and (14), the linearization coefficients must satisfy [5, 28, 29] where is the response vector , is the vector of equivalent linearization coefficients given by , and is given by (2).
When the joint PDF proposed in this paper (see (8)) is substituted in (15), mathematical expressions for the linearization coefficients are derived (see Appendix A). In particular, for , which corresponds to SDOF systems with softening behavior, the equivalent linearization coefficients are as follows: where and .
6. Covariance of the Response of the Hysteretic System
The covariance matrix of the response of the hysteretic SDOF system, where , is calculated by solving the following equation [28, 29]: where is a matrix depending on the mechanical and geometrical properties of the system, the linearization coefficients, the modulating function coefficients, and the filter parameters. In this study, is given by
is a matrix of zeros, except the element which is equal to the amplitude of the white noise spectral density .
In order to compute , it is essential to start the integration process of (17) when the CloughPenzien filter response has reached the stationary stage. Thus, the initial condition for should correspond to the covariance of the stationary response of the CloughPenzien filter as described in Appendix B [30].
7. Verification of the Accuracy of the Proposed NSEL Approach
In this section the standard deviation of the response of several hysteretic SDOF systems is calculated by means of the proposed NSEL criterion using the mathematical expressions obtained above ((13) and (16) using values of Table 1). The SDOF systems have vibration periods s and s and design ductility factors and , and they are subjected to a nonstationary random process with the statistical properties of the seismic motion recorded in Mexico City at the SCT during the September 19, 1985 earthquake (see Section 3).
The response of the SDOF hysteretic systems is calculated by means of SEL and alternatively by Monte Carlo simulation using 50,000 simulated accelerograms. Two SEL criteria are used: (a) assuming Gaussian response [1] and (b) assuming nonGaussian response using the criterion and the mathematical expressions proposed herein.
The statistical response of the SDOF system with s and is shown in Figures 8, 9, and 10, and that corresponding to the system with s and is in Figures 11, 12, and 13. Figures 8 to 13 show that the proposed approach predicts with reasonable accuracy the maximum response; however, as well as the Gaussian method, it does not predict the permanent drift of the hysteretic system especially for high values of the design ductility (see Figures 8 and 11). The statistical response of the SDOF system with s and is shown in Figures 14, 15, and 16, and that corresponding to the system with s and is in Figures 17, 18, and 19. As in the case of SDOF system with s, the proposed approach predicts with reasonable accuracy the maximum response of the SDOF system with s; however it does not predict with good accuracy the permanent drift.
The relative error () in the calculation of the maximum standard deviation of the response (, and ) is obtained as follows:
The relative errors of the maximum response, corresponding to all the cases analyzed, are shown in Tables 2 and 3 corresponding to structural vibration periods of 2.1 and 4.0 s, respectively. A negative value of indicates that the equivalent linearization underestimates the response. Tables 2 and 3 show that the proposed NSEL method leads to results close to those obtained by means of Monte Carlo simulation.


8. Extension to MultipleDegreeofFreedom Systems
In this section an extension of the proposed method to multipledegreeoffreedom (MDOF) systems is presented.
Consider a nonlinear dynamic MDOF structural system subjected to a vector input with an dimensional response satisfying the nonlinear differential equation: Here bold letters are related to MDOF systems. The equivalent linear system of (20) is defined as where , , and are defined as the equivalent mass, viscous damping, and stiffness matrices of . They are determined minimizing the mean square error , where is the dimensional vector: Denoting by and the dimensional response vector and the equivalent structural matrix given, respectively, by it can be shown that can be obtained from [5, 28, 29] Herein, the random excitation is modelled as a linear combination of the responses of linear filters to white or shot noise. With the aim of calculating the statistical secondorder response of the MDOF structure, it is convenient to express (21) in state space form: where is the state vector and the timedependent system matrix given by and the matrix contains the coefficients of the linear combination of the response of the filters. The dynamics of the latter is governed by where the matrix contains the coefficients of the filters and the vector the white or shot noise. Equations (25) and (27) can be rewritten as where Matrix has the same structure as but with higher dimension. Thus, (28) controls the response of an augmented system consisting of the linear filters and the structural system in series [5]. Finally, the covariance matrix of the response of the MDOF system is calculated by solving where is a matrix of zeros, except the last element which is equal to the amplitude of the white noise spectral density .
The mathematical expressions given in this section correspond to a MDOF structural system with dimensional response, while those given in the previous sections correspond to the special case of a SDOF with 2dimensional response . For example (20) took the particular form of (1) and (2); (21) took the particular form of (14); (24) was reduced to (15) while (30) took the particular form of (18).
9. Conclusions
A model appropriately representing the time evolution of the joint PDF of the structural response of softening systems subjected to seismic ground motions was proposed. The displacement and velocity was assumed jointly Gaussian and the marginal PDF of the hysteretic component of the displacement was modeled by a mixed PDF which is Gaussian when the structural behavior is linear and turns into a bimodal PDF when the structural behavior is hysteretic.
The NSEL method proposed in this study uses mathematical expressions for the equivalent linearization coefficients (16) corresponding to the nonGaussian response of inelastic nonlinear systems. The mathematical expressions depend on the parameters of the joint PDF of the structural response, and they take into account the time evolutionary character of the PDF. It is noticed that other NSEL methods found in the literature have proposed mathematical expressions for the equivalent linearization coefficients; however, they do not consider the time evolution of the joint PDF.
The mathematical expressions derived here (16) can be applied to obtain the nonGaussian response of any kind of SDOF softening system with smooth transition from elastic to plastic behavior (). It is possible to derive expressions of the linearization coefficients for other values of .
The NSEL method proposed was applied to calculate the statistical response of SDOF systems with different vibration periods and different design ductility factors. The results show that, for the SDOF systems analyzed, the proposed criterion estimates with good accuracy the time evolution of the structural response; however, it does not predict the permanent drift exhibited by hysteretic systems with high design ductility values.
It is concluded that the proposed NSEL approach is useful to estimate the maximum standard deviation of the structural response of softening systems subjected to seismic ground motions. The relative errors in the calculation of the maximum standard deviation of the displacement, velocity, and hysteretic component of displacement as compared with those obtained by means of Monte Carlo technique are low and in general they are smaller than those obtained by the Gaussian equivalent linearization method.
The parameters , , , and , , in (13) can be found for other types of seismic excitation (e.g., wideband ground motions).
Appendices
A. Linearization Coefficients
This appendix presents the derivation of the equivalent linearization coefficients given by (16).
The linearization coefficients are obtained from where is the response vector and is the vector of linearization coefficients, given by and .
The inverse of the covariance matrix iswhere . Substituting (A.2) in the second factor of (A.1), the following is obtained: where Considering that is modeled by (8), , , and can be easily obtained, thus making it possible to obtain mathematical expressions for the linearization coefficients. For the case of , (A.4) are as follows: where where , , is the exponential function, and is the error function given by Substituting ((A.5)(A.6)) in (A.3), the value of for is obtained. This result together with the inverse of the covariance matrix (A.2) is substituted in (A.1) in order to obtain the linearization coefficients corresponding to . Consider
B. Initial Condition of the Covariance Matrix of the Response
Following the state space approach, the state equation of a linear system subjected to a nonwhite excitation modeled as a filtered white noise is [29] where is the state vector, the matrix is the system matrix, and the matrix contains the coefficients of the linear combination of the response of the filter. The response of the filters is controlled by where the matrix contains the coefficients of the filters and is a vector with zeros except the last element which correspond to white noise. (B.1) and (B.2) can be rewritten as where Matrix has the same structure as but with a higher dimension. To avoid introducing the effect arising from the transient response of the filter, the output of the filter must be allowed to reach the stationary phase before it is multiplied by the modulating function [5]. Analytically, it is reached by selecting the suitable initial conditions for the covariance matrix of the overall state variable vector , which is governed by the following differential equation: where is a matrix of zeros except the last element which is equal to the amplitude of the white noise spectral density .
It is convenient to introduce a partition in in agreement with (B.4). Thus, where , , and . If the original system is at rest, at , then the response will be zero and all the elements of and will be zero; however the elements of , at , will correspond to the stationary response of the filter which is governed by where has the same form as but with a lower dimension. For the case when the CloughPenzien filter is used, the solution of (B.8) is where