On the background and biresonant components of the random response of single degree-of-freedom systems under non-Gaussian random loading

On the background and biresonant components of the random response of single degree-of-freedom systems under non-Gaussian random loading

Engineering Structures 33 (2011) 2271–2283 Contents lists available at ScienceDirect Engineering Structures journal homepage: www.elsevier.com/locat...

2MB Sizes 1 Downloads 6 Views

Engineering Structures 33 (2011) 2271–2283

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

On the background and biresonant components of the random response of single degree-of-freedom systems under non-Gaussian random loading V. Denoël University of Liège, Department of Architecture, Geology, Environment and Constructions, Chemin des Chevreuils, 1, B52/3, 4000 Liège, Belgium



Article history: Received 19 July 2010 Received in revised form 10 March 2011 Accepted 1 April 2011 Available online 11 May 2011 Keywords: Bispectrum B/R decomposition B/bR decomposition Perturbation Integral Non-Gaussian Skewness

abstract The variance of the response of a single degree-of-freedom system subjected to a low frequency excitation is usually decomposed into its background and resonant contributions. In this paper we aim at the formulation of such a decomposition for the third statistical moment, which should in principle sidestep the heavy double integration of the bispectrum of the response. In large finite element models, the estimation of the bispectrum of the loading for a given frequency (ω1 , ω2 ), which is a necessary stage towards estimation of the response, is the most expensive computational task. We therefore formulate the problem with the underlying constraint that the number of estimations of the bispectrum of the force be kept as small as possible. Invoking the perturbation theory in the context of the computation of integrals, we propose to decompose the third moment into background and biresonant components with expressions that are not trivially adapted from the decomposition of the variance. Thanks to the proposed method, the double integration of the bispectrum is avoided and represents accurately the response of lightly to moderately damped structures under a low-frequency loading. The formal expression for the biresonant component still requires an integration that should preferably be avoided. In the second part of the paper, we then investigate the practical implementation of the proposed formula and study the possible application of local, global or hybrid numerical approximations of that remaining integral, so as to further increase the computational efficiency. Finally two numerical experiments illustrate the prospect of the proposed method: the third statistical moment of the response is accurately computed with less than 20 estimates of the bispectrum of the loading, whereas an advanced numerical procedure for the double integration would require a mesh of probably more than (a) thousand(s) points. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Many engineering applications consist in determining the steady state statistics of the output of a linear time invariant system subjected to a Gaussian stationary random input. This well posed problem has been studied meticulously since 1950’s and formal solutions of this problem may nowadays be found in the literature (e.g. [1]). A typical application -commonly called buffeting analysis- in wind engineering deals with the study of gusty wind flows around structures [2,3]. In this field, a series of papers published by Davenport in the early 1960’s [4,5] is usually recognized as the incipient engagement of the beneficial pooling of applied mathematics, statistics, meteorology and structural dynamics. In his work, Davenport formulated the necessary assumptions that allowed his concept of buffeting analysis fall within the typical framework of that class of problems. Even if the framework of this class of problems appears to be restrictive since it requires (i) stationarity, (ii) system linearity and

E-mail address: [email protected] 0141-0296/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2011.04.003

(iii) Gaussianity, the majority of applications still hinge on this concept. The reason is that the application of this theory is rather simple and that the suppression of any of the three limitations brings the problem to a much higher level of complexity, which is typically not affordable to practitioners. Formal solutions to the enlarged class of problems (considering non-stationarity, nonlinearity and/or non-Gaussianity) have however existed since the 1950’s [6]. Ito proposed to tackle this more general problem by means of stochastic differential equations which allow to cope with any or several of the above limitations. Nevertheless, stochastic differential equations are doubtlessly not suitable for everyday applications. The transfer of knowledge from pure mathematics to the engineering community took a couple of decades as it required these solutions to be presented in a simpler and applicable form (e.g. [1]). Actually, several analysis methods emerged concurrently with formulations simpler than Ito’s stochastic differential equations, but encompassing only one of the three refinements at a time. For instance and setting again our sights on wind engineering applications, nonstationary random wind loadings are now used to model downbursts [7,8] and their influence on structures; random nonlinear structural behavior is also sometimes considered [9–12]; and non-Gaussianity of


V. Denoël / Engineering Structures 33 (2011) 2271–2283

the loading and its effect on the structural response is also a topical research subject [13,14,10,15–19]. Other engineering fields have experienced the same developments with a fairly similar timeline. In this paper we consider one of the three semi-extended class of problems where the input is still stationary and applied to a linear system, but features a non-Gaussian character. In this case, an accurate description of the response is given by means of its power spectral density (psd), and similar higher order spectra such as the bispectrum and trispectrum [18]. These generalized spectra are integrated in a multi-dimensional frequency space in order to provide the statistical moments of the response. The numerical integration of these functions is usually not trivial in view of the extreme acute shape of the generalized spectra. In many applications, the characteristic time scale of the input process (the loading) is one or several orders of magnitude above the intrinsic time scale of the system (the natural period of the structure). When this condition is met, a well-known approximation for the computation of the second moment, namely the variance or the integral of the psd, consists in decomposing the total integral into two components, which results in a drastically simplified computational task [20,4]. In this paper we seek to develop the same approximation for the estimation of the third statistical moment, i.e. the integral of the bispectrum. The paper is organized as follows. The problem is formally stated and motivated in Section 2. In Section 3, a particular example with explicit solutions is presented; it is used in a validation procedure in Section 6. In the interval, the proposed approximation is developed in Section 4 along with some implementation issues. Also a degenerated case encountered in many practical applications is exposed in Section 5. 2. Statement of the problem

the instantaneous Gaussian probability density function of the response. It is formally obtained as


Sx (ω) dω.

m2,x =



Notice that we will write the variance as m2,x -unlike the conventional notation σx2 -to better emphasize the parallelism with the third and higher order analyses. When the characteristic time scale of the loading is one or several orders of magnitude above the natural period of the structure, which is usual for wind loaded structures [2,3], the psd of the response presents three distinct peaks, as sketched in Fig. 1(a), corresponding to background and resonant components (B/R). A fair approximation, apparently suggested by Davenport [5] and later on formally rationalized by Ashraf and Gould, see [20], consists in estimating the variance of the response as the sum of two terms corresponding to these components, respectively,

2,x = m2,b + m2,r m



m2,b = |H (0)|2


Sf (ω) dω =



m2,r = Sf (ω0 )


|H (ω)|2 dω = −∞



π ω0 Sf (ω0 ) . 2ξ k2


This indicates that the background component is simply obtained from the variance of the force m2,f and the resonant one is obtained by replacing the actual psd of the force by a constant value equal to that taken at ω = ω0 (which earned this method to be appointed white noise approximation or B/R decomposition). An alternative representation is based on a dynamic amplification factor defined as m2,x

The equation of motion for a single degree-of-freedom system writes

A2 =

mx¨ + c x˙ + kx = f

so that the dynamic response is simply obtained by multiplication of the quasi-static response m2,b by the dynamic amplification factor A2 . The approximate dynamic amplification factor resulting from the B/R decomposition is


where m, c and k represent the mass, viscosity and stiffness; f (t ) and x(t ) are the random processes representing the loading and the structural response. As only stationary responses are considered in this problem, no initial conditions are attributed to (1) and its side-by-side Fourier transform is therefore strictly equivalent. It is written, see e.g. [21], X (ω) = H (ω) F (ω)


where X (ω) and F (ω) are respectively the Fourier transforms of the response and of the loading, whilst the complex transfer function H (ω) given by H (ω) =

1/k 1−

ω2 ω02

+ 2iξ ωω0


is √ here expressed as a function of the natural frequency ω0 = k/m and the damping ratio ξ = c /2mω0 . In a Gaussian context, the force f would be fully described, in a probabilistic manner, by its mean value µf and its power spectral density (psd) Sf (ω), which represents the distribution of the second central moment -the variance- in frequency space. Likewise the response would be fully described by the mean response µx simply obtained as µx = µf /k, and the psd of the response Sx (ω) obtained as, see e.g. [21], Sx (ω) = |H (ω)| Sf (ω) . 2


The second central moment of the response m2,x is of major importance as it allows, together with the mean value, inferring




2 = 1 + A



π ω0 Sf (ω0 ) . 2ξ m2,f


Owing to its high computational effectiveness compared to an expensive numerical integration of (5), this approximation is recommended in many codes and standards and is therefore systematically applied in everyday practice. In a non-Gaussian context, the first two moments are no longer sufficient to fully describe the loading nor the response. There exist quantities similar to the psd that represent the distribution of statistical moments higher than two in the frequency space. In particular, the bispectrum Bf (ω1 , ω2 ) represents the distribution of the third central moment of the force f in a 2-D frequency space [22]. Similarly to (4)–(5) for the second order analysis, the bispectrum of the response is expressed as, see [22], Bx (ω1 , ω2 ) = K2 (ω1 , ω2 ) Bf (ω1 , ω2 ) ,


and the third central moment is obtained as



m3,x =

Bx (ω1 , ω2 ) dω1 dω2 .



Eq. (10) expresses that the dynamic system under investigation, filters out and amplifies (due to the resonance phenomenon) the statistical dissymmetry in the loading process to provide a more or less dissymmetric response. The second Volterra kernel K2 (ω1 , ω2 ) is expressed as, see [22],

V. Denoël / Engineering Structures 33 (2011) 2271–2283


Fig. 1. Typical psd (left) and bispectrum (right) of the response of a single degree-of-freedom system subjected to low frequency turbulence. The psd (left) is commonly decomposed as the sum of background and resonant components. The band widths associated to each component are respectively α and ξ ω0 , see Section 4. This paper aims at the investigation of such a decomposition for the bispectrum.

Fig. 2. Real part of the second Volterra kernel K2 /k3 (for ξ = 0.1).

K2 (ω1 , ω2 ) = H (ω1 ) H (ω2 ) H (ω1 + ω2 )


where the overbar denotes the complex conjugate. As the bispectrum of the force Bf is usually a real function and because of some parity conditions on K2 , the imaginary part of K2 brings no contribution to the third central moment [23]. Only the real part of K2 requires therefore attention. It is represented in Fig. 2 for a rather large value of the damping coefficient ξ = 0.1 in order to better assess the particular shape of that function. It features

• a characteristic basin close to the origin with elliptic level curves, which corresponds to a quasi-static behavior K2 (ω1 , ω2 ) ≃ 1/k3 for (ω1 , ω2 ) = ord(ϵω0 ) with ϵ ≪ 1, • six high peaks centered in (±ω0 , 0) , (0, ±ω0 ) and ± (ω0 , −ω0 ), also referred to as quasi-static biresonance peaks hereafter, as they correspond to resonance in two of the three factors in (12) whilst the third one corresponds to quasi-staticity. The height of these peaks is therefore 1/4ξ 2 k3 , • six low peaks centered in ± (ω0 , ω0 ), ± (2ω0 , −ω0 ) and ± (ω0 , −2ω0 ), also referred to as inertial biresonance peaks, as they correspond to resonance for two of the three factors in (12) whilst the third one corresponds to the inertial range of the response. The height of these peaks is 1/12ξ 2 k3 , • six triangular basins connecting the peaks to one another. They have steep slopes along their perimeter and a rather flat bottom located below the zero level.

Higher order probabilistic analyses are seldom performed today [24]. Moreover, the B/R decomposition (6) offers an affordable access to the second order probabilistic analysis that few would afford if the integration in (5) had to be performed numerically. There is evidence that higher order statistical analyses will not be commonly applied as long as the double integration in (11) has to be performed numerically. In spite of several advanced numerical integration procedures [25–29], authors universally agree that it is difficult to find an acceptable tradeoff between accuracy and computation time. In this connection, there is no doubt that the estimation of the psd, bispectrum, trispectrum, etc. of the loading (typically generalized forces) is the most extensive computational task if the structure is represented by a large finite element model. In practice, it is unlikely that the third order analysis be fluently applied for large structures if the number of integration points in the 2-D frequency space exceeds one hundred. The objective of this paper is to formally derive an approximation of the third central moment of the response, in a fashion similar to the B/R decomposition used for the second central moment, with the firm intention to keep the number of integration points as low as possible. The typical shape of the bispectrum of the response of a single degree-of-freedom system to a low frequency loading is sketched in Fig. 1(b). The topological disjunction of the different peaks suggests the possibility of extending, to the third order,


V. Denoël / Engineering Structures 33 (2011) 2271–2283

the principles of the B/R decomposition. We show next that this operation is not necessarily trivial, but at least that the double integration may be replaced by a single one.

where a is a characteristic frequency, typically one or several orders of magnitude below the natural frequency ω0 . The introduction of (18) into (14) and (16) yields

3. A particular example with explicit solutions

Sf (ω) = m2,f




π + ω2    a 2 a2 + 23 ω12 + ω1 ω2 + ω22       . (20) Bf (ω) = m3,f π a2 + ω12 a2 + ω22 a2 + (ω1 + ω2 )2

3.1. General developments for a quadratic Gaussian process


Let us consider that the loading is expressed in the quadratic form

Then substitution of (19) and (20) into (4), (5), (10) and (11) yields, after some cumbersome calculus and algebra, see [23]

f = γ (1 + u)2 = γ + 2γ u + γ u2

m2,x =


where the low-frequency content process u (regarded here as the turbulence) is modeled as a dimensionless Gaussian process with zero mean and standard deviation σu , and γ is a reference force. For small turbulence intensity (σu ≪ 1) and under some monotonicity conditions on Su , see [23], -which we assume herein to be fulfilled as they do not affect the illustration of the main argument of this paper-, the psd of f is expressed as Sf (ω) = 4γ 2 Su (ω) .


The variance of the force m2,f , obtained by integration of (14) along frequencies is


Sf (ω) dω = 4γ 2 σu2 .

m2,f =



Under similar assumptions, the bispectrum of f is expressed as

+ Su (ω2 ) Su (ω1 + ω2 )

Substitution of (15) and (19) into (6)–(7) yields

m2,r (16)

and the third central moment as



m3,f =

Bf (ω1 , ω2 ) dω1 dω2 = 24γ 3 σu4 .



In the context of a buffeting analysis, the equations developed in this section refer to an ideal point-like structure. Similar expressions exist for more complex structures with dimensions comparable with the turbulence length scales. For such structures, the right-hand side in ((14)) is multiplied by an aerodynamic admittance which represents the imperfect correlation of the wind pressures on the structure. This results in an equivalent psd of force. In principle, the bispectrum of the force, as given by ((16)) should also be modified by a third order admittance function, and an equivalent bispectrum of force should be obtained. However no information is available to date concerning this topic. In the following, Sf (ω) and Bf (ω1 , ω2 ) refer to the actual psd and bispectrum in the case of a point-like structure, but to equivalent quantities, obtained by multiplication by appropriate admittance functions, in all other cases.




4γ 2 σu2 k2

π ω0 Sf (ω0 ) m2,f 1 aω 0 , = = 2 2 2 2ξ k k 2ξ a + ω02

2 = 1 + A


aω 0

2ξ a2 + ω02

σu2 2 π a + ω2


3.4. Rudimentary transposition of background-resonant decomposition to the third order Because it seems so straightforward for the estimation of the second moment, it could be appealing to apply the same method for the third order, i.e. to decompose the third central moment as +∞


m3,f k3

∫ 3,r = Bf ∗ m = Bf ∗


24γ 3 σu4

+∞ ∫ −∞






This function is also represented in Fig. 3 and shows an excellent agreement with the exact result A2 . This is an illustration of the good performance of the B/R decomposition for the estimation of the second order response.

3.2. Application to an Ornstein–Uhlenbeck process Next we further particularize the example by considering that the turbulence u is an Ornstein–Uhlenbeck process [30], i.e. that its psd is given as


which then introduced in (6) provides the usual B/R approximation of the variance of the response. To facilitate comparison with (21), it is written m2,f 2 (ω0 , a, ξ ) 2,x = 2 A m (26) k where the approximate dynamic amplification is readily obtained as

m3,b = K2 (0, 0)

Su (ω) =


3.3. Background-resonant decomposition for the second order

m2,b =

A2 (ω0 , ξ )

k2 m3,f

A3 (ω0 , ξ ) (22) k3 Eqs. (23) and (24) given in Box I where ck (ξ ) , k = 0, . . . , 7, are polynomials in ξ (not reported in this paper for conciseness). Identical results are also obtained by Grigoriu [31] with the method of moments in a time domain analysis. Because m2,f /k2 and m3,f /k3 are precisely the second and third statistical central moments of the response that would be obtained if the response was quasi-static, factors A2 and A3 may be seen as the second and third order dynamic amplification factors. They are represented by solid lines in Fig. 3. A major difference between both factors is that the second order one is unbounded for ξ → 0 whereas the third order one is bounded to a finite value. m3,x =

 Bf (ω1 , ω2 ) = 8γ 3 Su (ω1 ) Su (ω2 ) + Su (ω1 + ω2 ) Su (ω1 )


Bf (ω1 , ω2 ) dω1 dω2 −∞

k3 +∞

Re [K2 (ω1 , ω2 )] dω1 dω2 −∞

8ω02 π 2

3k3 1 + 8ξ 2


V. Denoël / Engineering Structures 33 (2011) 2271–2283

A2 =


ω0 (a + 2ξ ω0 ) 2ξ + 2ξ aω0 + ω02 1



1 1 3 1+8ξ 2

A3 = 

1 + 2ξ

ω0 a


1 + 2ξ

ω0 a


=7  ω0 2 k∑ a

 ω0 2 2  a

ck (ξ )

k =0

4 + 4ξ

ω0 a


 ω0 k a

 ω0 2   a

1 + 4ξ

ω0 a



 ω0 2  a

Box I.

Fig. 3. Second order (left) and third order (right) dynamic amplification factors. Exact results are represented by solid lines; black dots represent the results obtained with the B/R decomposition. The second order dynamic amplification factor is unbounded for ξ → 0 whereas the third order dynamic amplification factor is bounded to 3.

where m3,b is the background component obtained by approximating the Volterra kernel by its value K2 (0, 0) at the origin (a quasi3,r is static behavior). An estimation of the resonant component m obtained, as that resulting from the B/R decomposition, by replacing the actual bispectrum of force by a constant value Bf ∗ , reducing therefore the problem to the estimation of the integral of Re [K2 ]. This computation may be done once for all, with the residue theorem of contour integration [32] for instance and yields the above result. Whatever is chosen for Bf ∗ , the resulting approximate ex3,r is obviously flat for the small values of ξ enpression m3,b + m countered in typical applications. The expected result is however damping dependent as suggested in Fig. 3. This attempt at obtaining an approximate expression of m3,x fails to be efficient because the   double integral of the Volterra kernel behaves as 1/ 1 + 8ξ 2 , i.e. is almost constant in operating values of ξ and is therefore not suitable to model a damping dependent amplification factor. On the contrary the integral of the squared transfer function behaves as 1/ξ which is the expected behavior of the resonant part and makes then possible the replacement of the psd by a constant value. From this preliminary illustration, we may conclude that the development of an accurate estimation of the integral is subject to the consideration of the particular shape of the bispectrum of the loading in the vicinity of the peaks, or even probably along the 2-D frequency space. 4. Approximation of integrals A classical problem of applied mathematics consists in the estimation of the definite integral of a function of a small nonzero parameter ε , while the integral for ε = 0 may eventually be obtained explicitly [33]. The computation of integrals (5) and (11) happens to be such a problem. All the more, both the ratio α/ω0 which quantifies the disjunction of the quasi-static peak from the (bi-)resonant ones and the damping ratio ξ are assumed to be small parameters in our investigation. We formalize straight away the definition of the disjunction of these peaks by considering that the psd of the loading may be written Sf (ω) =




ω α


where S is a density function (with the main property that its integral is equal to one, see [34]) and α is therefore defined as the center of gravity of Sf (ω). For simplicity in the notations, we keep the same symbol and use an alternative definition in the context of the third order analysis, by noticing that the bispectrum of the force may be written Bf (ω1 , ω2 ) =









ω2  α


where B is a 2-D density function and α is now the center of gravity of Bf (ω). The disjunction of the quasi-static peak from the other ones is formalized by stating that α = ζ ω0 with ζ ≪ 1. Hinch presents a series of methods to establish the answers to such perturbed integrals [35]. In particular, one method consists in (i) identifying the zones of the domain of integration that will contribute the most to the integral, (ii) approximating the integrand over these zones, eventually after the introduction of a stretched co-ordinate in order to better highlight the identified region whenever it is small, (iii) subtract the approximate integrand from the original function, (iv) repeat the process iteratively until the integral is represented up to the desired order or the analytical formulation of the response becomes too complex to allow for a supplementary iteration. The methods presented next are inspired from this general procedure, but are adapted because the integrals to be approximated herein are given by means of functions Sf and Bf that are not known analytically and for which we will endeavor to limit the number of numerical estimates. In this section, we first apply the method to the second order analysis. In an innovative fashion, we develop expressions for the background and resonant contributions, which are cheerfully consistent with the usual B/R decomposition (7). Then the method is adapted to the third order analysis and we will derive the approximate expression for the computation of m3,x . 4.1. Estimation of the second moment The second central moment of the response is obtained by (5) where Sx (ω) exhibits three peaks. In this particular case, the domain is thus decomposed into


V. Denoël / Engineering Structures 33 (2011) 2271–2283

Fig. 4. Sketch of the integrand 1Sx (ω) after subtraction of the quasi-static component as a function of (a) the original variable ω and (b) the stretched co-ordinate η introduced to focus on the rightmost peak. The dashed line represents the asymptotic behavior of the integrand for small ξ .

• a background region in the vicinity of the origin. The spread of this region is of order α = ζ ω0 , with ζ ≪ 1. At no extra cost, this region may be extended to almost the resonance peak, at least in a frequency band much wider than α , along which the transfer function is roughly constant. The integral I1 along this extended region is therefore m2,f

ord(1); (31) k2 • two resonant regions on a frequency band of width ord (ξ ω0 ) centered on ±ω0 (resonance peaks). The height of the peaks is Sf (ω0 ) / (2kξ )2 . If we momentarily assume that the psd of the I1 =

force behaves as Sf (ω) ≃ k∞ α2,f ω for ω ≫ α , the integral α I2 along each of these short regions is therefore m

I2 =

ζ q −1 ord 4 ξ 

m2,f k∞ k2

  −q


the psd of the force Sf are smaller than their values in the background region. Depending on the relative smallness of ξ and ζ , the background or the resonant component contributes the most to the integral. So both necessarily have to be computed, keeping in mind though that one of them could eventually be much smaller than the other one. Let us first consider the quasi-static contribution in the origin, where the integrand is approximated as Sx,0 (ω) = |H (0)|2 Sf (ω). The integral of this function is Sx,0 (ω) dω

= |H (0)|



Sf (ω) dω = −∞


1Sx [ω0 (1 + ξ η)] ω0 ξ dη.

R2 =

m2,f k2



  1 1Sx (ω) = Sx (ω) − Sx,0 (ω) = |H (ω)|2 − 2 Sf (ω) k

The width of the peak is now of order 1 and the integrand 1Sx may be approximated by a Laurent series expansion for small ξ . The first factor of 1Sx is

|H (ω0 (1 + ξ η))|2 −

1 k2




4ξ 2 k2 1 + η2

  + O ξ −1


while the second factor is expanded as a Taylor series Sf [ω0 (1 + ξ η)] = Sf (ω0 ) + ξ ηω0 Sf′ (ω0 ) + O ξ 2

 


as Sf is a C1 -function. Only leading order terms in each factor are considered next so that the integrand is approached by Sf (ω0 ) 1 1Sx [ω0 (1 + ξ η)] ≃ 1 Sx [ω0 (1 + ξ η)] = 2 2 4ξ k 1 + η2


which, after substitution into (37) provides an approximation of the resonant component associated to the peak centered on ω0 +∞

1 Sx [ω0 (1 + ξ η)] ω0 ξ dη

ω0 Sf (ω0 ) = 4ξ k2




As sketched in Fig. 4(a) and as expected anyway, 1Sx (ω) presents essentially two peaks in the resonant regions around ω =


π ω0 Sf (ω0 ) . 4ξ k2


Similar developments for the peak around −ω0 yield an identical expression for m2,r − , which is also expected because of symmetry. Finally the resonant component is given as m2,r = m2,r + + m2,r − =


2 −∞ 1 + η


which needs to be integrated along frequencies

1Sx (ω) dω.




After subtraction of this first contribution from the integrand, we are left with


so that (35) becomes

m2,r + =


R2 =

ω = ω0 (1 + ξ η)


m2,b =

small damping coefficient. In order to better focus on the peak centered on ω = ω0 , we introduce a stretched co-ordinate η, see Fig. 4(b), defined as

in which q > 1 for integrability reasons. • two inertial regions ranging above ω0 (1 + ξ ) and below −ω0 (1 + ξ ). In these regions, both the transfer function H and

±ω0 . The spread of those regions is of order ξ ω0 , where ξ is the

π ω0 Sf (ω0 ) 2ξ k2


which is precisely similar to the traditional B/R decomposition, see (7). On top of being a novel formal demonstration of the commonly applied B/R decomposition, this section is essentially aimed at the introduction of the integration method that is next extended to the estimation of the third moment.

V. Denoël / Engineering Structures 33 (2011) 2271–2283

4.2. Estimation of third moment The third central moment of the response is obtained by (11) where Bx (ω) exhibits essentially one peak close to the origin, six plus six biresonance peaks and six distinct triangular basins. The integral is therefore decomposed into

• a background region in the vicinity of the origin. For reasons similar to those mentioned above, this region covers essentially the quasi-static basin in the Volterra kernel, and the first contribution to the integral is J1 =

m3,f k3

a suitable integrand from 1Bx and iteratively repeat the procedure with the processing of all identified components. Next, we decide to first focus on the peak located in (ω1 , ω2 ) = (ω0 , 0), and therefore to introduce the stretched co-ordinates

ω1 = ω0 (1 + ξ η1 ) ω2 = ω0 ξ η2


• six quasi-static regions ranging in frequency areas  biresonant  of order ord ω02 ξ 2 centered on the quasi-static biresonant peaks. The height of the peaks is Bf (ω0 , 0) /4k3 ξ 2 . As a

ω (1). (44) 4k3 This is a major difference with the computation of the second moment.   Indeed, J2 is independent from ξ whereas I2 = ord ξ −1 . This is due to the fact that at most two ‘‘singular lines’’ may intersect in K2 and hence that K2 is at worst ord(ξ −2 ). This observation is the key point that justifies that a basic extension of the B/R decomposition was not successful; • six inertial biresonant regions ranging in frequency areas of  order ord ω02 ξ 2 centered on the inertial biresonant peaks. The height of the peaks is Bf (ω0 , ω0 ) /12k3 ξ 2 . As a first approximation, the integral J3 is J2 =

Bf (ω0 , ω0 )

2 0 ord

ω (1), (45) 12k3 which is probably at least one order of magnitude below J2 according to the (presumably fast) decrease of Bf ; • six inertial bi-quasi-static regions, corresponding to the triangular basins with negative valued bottom. Along these regions, K2 = ord(1)/k3 , and the contribution J4 corresponding to these regions is thus J3 =

J4 =




ω0 , − 23 ω0

 ω02 ord(1)


k3 which is of an order of magnitude between J2 and J3 . • inertial regions that are basically located out of the six-arm star formed by the level sets of the Volterra kernel. For reasons similar to those exposed above, the contribution of these regions is negligible.

Let us consider again first the quasi-static contribution in the origin, where the integrand is approximated as Bx,0 (ω1 , ω2 ) = K2 (0, 0) Bf (ω1 , ω2 ). The integral of this function is



m3,b =

Bx,0 (ω1 , ω2 ) dω1 dω2 =


m3,f k3



After subtraction of this first contribution from the integrand, we are left with the estimation of

1Bx (ω1 , ω2 ) = Bx (ω1 , ω2 ) − Bx,0 (ω1 , ω2 ) [ ] 1 = K2 (ω1 , ω2 ) − 3 Bf (ω1 , ω2 ) k

which needs to be integrated along (ω1 , ω2 ). The next most important contributions are attributable to the quasi-static biresonant peaks (J2 ) and the triangular basins (J4 ). Extension of the integration procedure requires to select one of them, subtract


Considering (30) and α = ζ ω0 , 1Bx writes m3,f

1Bx =


K2 (ω0 (1 + ξ η1 ) , ω0 ξ η2 ) −

α2  ×B

1 + ξ η1 ξ η2









This function needs to be integrated but has no explicit solution as B essentially results from numerical computations. This integrand is therefore approximated by its Laurent series expansion for small ξ . The Laurent series expansion for ξ around 0 of the factor into the square brackets is K2 (ω0 (1 + ξ η1 ) , ω0 ξ η2 ) −


1 k3

1 4k3

(η1 − i) (η1 + η2 + i) ξ 2

+ o(ξ −1 ),


from which only the real part 4k3 ξ 2

η12 + η1 η2 + 1 + o(ξ −1 ) 2 η12 + η1 η2 + 1 + η22


is interesting since the contributions of the imaginary part cancel out throughout the whole domain of integration. The limit behavior of the last factor of 1Bx is governed by

 lim B

1 + ξ η1 ξ η2



ξ →0


 =B

1 ξ η2






Substitution of (54) and (55) into (52) yields an approximate integrand for small damping (ξ ≪ 1), and distinct peaks (ζ ≪ 1)

1 Bx =

η12 + η1 η2 + 1 Bf (ω0 , ξ η2 ω0 ) .  2  ξ η + η1 η2 + 1 2 + η2 1 2


4k3 2


This function is sketched in Fig. 5(b). It may be seen that the first term of the Laurent series (53) incorporates also parts of the triangular basins in the approximate integrand. For this reason, we may admit that a portion of J4 corresponding to the perimeter of the basins is actually included in J2 . The rest of J4 requires less attention due to the presumable rapid decrease of Bx along lines ω1 = cst and ω2 = cst. From this point onwards, we may thus already claim that the developments of the approximate integral will stop after the treatment of this contribution. Consideration of 1 Bx in (51) yields the approximate contribution related to one biresonant peak and parts of the connected triangular basins m3,r1 =


1Bx [ω0 (1 + ξ η1 ) , ω0 ξ η2 ] ω02 ξ 2 dη1 dη2 .


1 2 0 ord



R3 =

first approximation, the integral J2 on these small regions is therefore Bf (ω0 , 0)


in order to better emphasize the local behavior of the integral, see Fig. 5(a). The double integral that needs to be computed cuts down to

∫∫ ord(1);




Bf (ω0 , ξ η2 ω0 )


−∞ +∞

η12 + η1 η2 + 1 dη1 dη2 2 −∞ η12 + η1 η2 + 1 + η22 ∫ ξ ω03 +∞ Bf (ω0 , ω2 ) =π 3 dω 2 . 2 2 k −∞ (2ξ ω0 ) + ω2 ∫


V. Denoël / Engineering Structures 33 (2011) 2271–2283

Fig. 5. Sketch of the integrand after subtraction of the quasi-static component 1Bx as a function of the stretched co-ordinates (η1 , η2 ) introduced to focus on a particular quasi-static biresonance peak. The lower plots represent the asymptotic behavior of the integrand for small ξ ; it captures the shape of the peak as well as that of the connected triangular inertial bi-quasi-static region (triangular basin).

Because of the hexa-symmetry, the contributions of the six other quasi-static biresonant peaks are equal and the total biresonant contribution is m3,r = 6π

ξ ω03 k3

Bf (ω0 , ω2 )



(2ξ ω0 ) + ω 2

2 2

dω2 .


The consideration of the background contribution m3,b from (47) and this biresonant contribution yields an accurate background/biresonant (B/bR) estimation of the third central moment m3 = m3,b + m3,r .


For typical applications, it is not necessary to refine the computation with the consideration of other contributions like those coming from the six lower peaks. The approximate third order dynamic amplification factor 1 + m3,r /m3,b therefore writes

3 = 1 + 6π A

ξ ω03 m3,f

+∞ −∞

Bf (ω0 , ω2 ) dω2 . (2ξ ω0 )2 + ω22


Eq. (57) is the proposed approximation of the biresonant contribution to the third moment of the structural response. Application of (58) hinges on the same assumptions as those that are commonly adopted for the second order analysis and it

requires a simple integration instead of the heavier double integral (11). Furthermore, the integrand is expressed as the product of Bf (ω0 , ω2 ) and

φ (ω2 ) =


(2ξ ω0 )2 + ω22



a bell-shaped function that takes only significant values along a very short interval Ω = [−20ξ ω0 ; +20ξ ω0 ]. Even if there exists no explicit expression for this integral, the computation of the third order statistical response is greatly facilitated by use of the proposed B/bR decomposition. Retrospectively, it is interesting to underline that the introduction of the small parameter ζ was the key element in the developments. It was precisely aimed at reducing Bf to a function of only ω2 in the vicinity of the biresonant peak, see (56), and as a result allowed dropping by one the order of the integral. Application of the B/bR decomposition is however less immediate than the B/R because the bispectrum still needs to be integrated once. This is not necessarily a trivial task because the bispectrum is not given analytically. The practical implementation of the B/bR is discussed in the following section. We must acknowledge that the third central moment is a quantity with few physical meanings. This is compensated by the definition of the skewness coefficient which provides a tangible

V. Denoël / Engineering Structures 33 (2011) 2271–2283

insight on the asymmetry of the statistical distribution. It is formally defined as a function of the central moments of the response m2,x and m3,x , as

γ3,x =

m3,x 3/2 m2,x

= 

A3 m3,b

A2 m2,b

3/2 .



Because γ3,f = m3,b /m2,b corresponds to the skewness coefficient of the loading, we may therefore introduce Aγ3 , the dynamic amplification coefficient of the skewness, which represents the amount by which the skewness of the loading is increased/decreased

Aγ3 =

γ3,x A3 = 3/2 . γ3,f A2


ξ ω03

γ3 = A

3 A


1 + 6π m


= 

2 A


Bf (ω0 ,ω2 ) −∞ (2ξ ω0 )2 +ω2 dω2 2

 +∞

3/2 πω0 Sf (ω0 ) 2ξ m2,f



4.3. Practical implementation of the proposed formula In this section, we briefly study various methods to approximate the remaining integral in the B/bR decomposition. Local, global and mixed techniques are investigated. 4.3.1. Asymptotics of the approximate integral If ξ ω0 ≪ ζ ω0 , the function Bf (ω0 , ω2 ) may be considered to be constant along the domain Ω on which φ takes significant values. Under this assumption (73) reduces to ≪ζ ) (ξ A = 1 + 3π 2 ξ ω02 3

Bf (ω0 , 0) m3,f



This corresponds precisely to the case studied in Section 3.4, with a dynamic amplification factor independent from the damping coefficient. On the other hand, if ξ ω0 ≫ ζ ω0 , we may assume that the decrease of Bf (ω0 , ω2 ) is so rapid that it essentially takes place along Ω . A second asymptotic response is thus obtained by considering φ = φ(0) in (59), which reduces to (ξ ≫ζ )

3 A

= 1 + 3π

ω0 2ξ

4.3.2. Numerical integration A short brainwork is presented next about the (global) numerical integration of (57). This section is deliberately kept short as this is not meant to be the main objective of this paper. A standard integration procedure would require to truncate the integration domain to a finite interval Ω = [Ω1 , ΩN ], e.g. Ω = [−20ξ ω0 ; +20ξ ω0 ], and to sample the integrand for discrete values ω2k , k = 1, . . . , N with restrictions ω2k+1 > ω2k and also generally ω21 = Ω1 and ω2N = ΩN . Application of the rectangle rule consists in using uniformly distributed integration points. Its application to the estimation of (57) gives a numerical approximation of the dynamic amplification factor

(ω2 )


Invoking the central limit theorem, it is expected that the skewness of the response be smaller than that of the loading, which rather translates into a reduction factor. Estimation of this coefficient is of paramount importance in the applicability of the proposed method. Indeed, once statistical properties of the loading (the dissymmetry γ3,f here) are obtained by separate developments (see e.g. [36] for aerodynamic loading on bridge decks), the estimation of Aγ3 gives a straightforward estimation of the statistical properties of the response, with γ3,x = Aγ3 γ3,f , which in turn provides the required information for a better assessment of structural extreme values. The exact expression of Aγ3 is obtainable in academic examples only; nevertheless a variant of (62) obtained with the B/R and B/bR approximations may be used instead

+∞ −∞

Bf (ω0 , ω2 ) m3,f

dω2 .


Unfortunately, in practical application, ξ and ζ vary in operating ranges of similar orders of magnitude which makes the use of these approximate solutions rather vain.


  N Bf ω0 , ω2k ξ ω03 − 1ω2 = 1 + 6π m3,f k=1 (2ξ ω0 )2 + ω22 k


where 1ω2 = ω2k+1 − ω2k (k = 1, . . . , N − 1). Notice that the reduction of the domain of integration to Ω may be accompanied by an excessive absolute error on the integral. Debatable situations include the case where Bf (ω0 , ω2 ) is negligible along Ω but not elsewhere, conditions that are often met, e.g. ξ ≪ 1 and Bf (ω0 , 0) = 0. In this case, the integral consists of two main contributions (one along Ω and the other one out of Ω ) and the integral is necessarily small which translates into a small dynamic amplification factor. This warning notice may therefore be discussed by considering that only the dynamic amplification factor (73) is ultimately relevant. Precisely, errors on small biresonant components are meaningless to the dynamic amplification factor because biresonant and background components have to be combined. We may be concerned by the ability of the numerical procedure to cover accurately the asymptotic case ξ ≪ 1, which is usually encountered in typical applications. In this case, the integral in (57) reduces to the integration of φ (ω2 ), as seen in the former asymptotic developments. Actually, φ (ω2 ) is known as Runge’s function [37] and, in spite of a commodious analytical expression, it is rather troublesome to numerically integrate with precision. When integration is performed along a finite domain, a possible solution is to use a non uniform distribution of integration points [38]. In our case, because the domain of integration spans from −∞ to +∞, an adequate non uniform distribution may be obtained implicitly by introducing the change of variable

ω2 = 2ξ ω0 tan θ .


Integration of φ (ω2 ) becomes


φ dω 2 = −∞

+π /2

1 2ξ ω0

−π /2

dθ ,


i.e. integration of a simple constant function, which is exactly obtained with many simple integration procedures. In particular, we propose to use a uniform sampling with N integration points  along − π2 ; π2 and the rectangle rule (θ)

3 = 1 + 3π A

N ω02 −


Bf (ω0 , 2ξ ω0 tan θk ) 1θ



1 with θk = − π2 + Nk− π , k = 1, . . . , N. −1 The number N of integration points could be chosen as small as N = 10 or N = 20, in order to keep the function counts as low as possible but at the expense of an inaccurate dynamic amplification estimate.


V. Denoël / Engineering Structures 33 (2011) 2271–2283

m3,r = 2π

m3,f ξ ω03


4 3 uk


Su (ω0 ) Su (ω2 ) + Su (ω0 + ω2 ) Su (ω0 ) + Su (ω2 ) Su (ω0 + ω2 )

(2ξ ω0 )2 + ω22


dω2 .


Box II.

4.3.3. Approximate solution for decreasing bispectrum The advantage of a local method as (64) is to provide an estimation of the integral with a very limited number of function estimates, but with a worse accuracy than global methods as (66) and (69) which however require a larger computational effort. In this section, we propose an intermediate method that keeps a very low number of integration points but not strictly located at ω2 = 0, as in the asymptotic study. The idea is to replace the actual bispectrum Bf (ω0 , ω2 ) with a simple analytical expression B∗

 Bf (ω0 , ω2 ) =


a∗ = 

πξ ω0 Bf (ω0 ,0)




Introduction of (70) into (57) yields an explicit formulation for the biresonant contribution

ω02 B∗



a∗ + 2ξ ω0


and the dynamic amplification factor

(Bf )

3 A

= 1 + π2

ω02 B∗


2m3,f a∗ + 2ξ ω0



This method is efficient because it just requires 2 functions estimates but its application is definitely restricted to bispectrum with a local decrease similar to (70). 4.3.4. Dynamic amplification of the skewness Several approximation techniques were introduced in this section. In light of (63), other estimates of the dynamic amplification of the skewness can be obtained. These are written in a canonic form

(ψ) A γ3 =


3 A



2 A

where (ψ) symbolizes one of the aforementioned methods: Bf . (ξ ≪ ζ ) , (ξ ≫ ζ ), (ω2 ), (θ) ,  5. The particular case of quadratic form of a Gaussian process As a particular case, we study next how the approximate biresonant contribution simplifies when f is expressed as a squared Gaussian process, see Section 3.1. The introduction of (16) into (57) yields Eq. (75) given in Box II. On account that Su (ω0 ) ≪ Su (ω2 ) due to the assumed rapid decrease of Su , and because the major contribution to this integral comes from ω2 ∈ [−10ξ ω0 ; +10ξ ω0 ], this equation further

Su (ω0 )


+∞ −∞

Su (ω) /σu2

(2ξ ω0 )2 + ω2



Application of the second order B/R decomposition (9) yields

2 (ω0 , ξ ) = 1 + A

π ω0 Su (ω0 ) 2ξ σu2


so that the dynamic reduction factor on the skewness is

γ3 = A


Bf (ω0 ,πξ ω0 )

m3,r = π 2

3 (ω0 , ξ ) = 1 + 4π ξ ω03 A



where a∗ and B∗ are equivalence parameters to be determined. Because Bf (ω0 , ω2 ) is multiplied by Runge’s function φ (which takes significant values on Ω ), parameters a∗ and B∗ could be chosen so that the fitted function has the same ordinate than the original function at ω2 = 0 and ω2 = π ξ ω0 , i.e. B∗ = Bf (ω0 , 0) ;

simplifies and ultimately provides the following estimation of the dynamic amplification factor

1 + 4π ξ ω03

Su (ω)/σu2 −∞ (2ξ ω0 )2 +ω2 dω . 3/2 π ω0 Su (ω0 ) 2ξ σu2

Su (ω0 )

 +∞




6. Validation The proposed method as well as its practical implementation are illustrated. Two different loadings are considered; they are both expressed as squared Gaussian processes, see Section 3.1. The power spectral densities of the generating turbulences are

σu2 Su1 (ω) = π a2 + ω 2 a

√ and Su2 (ω) =

2 aω2 σu2

π a4 + ω 4



These two stochastic processes have a similar high-frequency decrease (in ω−2 ), but show a major difference in the vicinity of the origin. Indeed, u1 shows a certain low-frequency content whereas the very low-frequency content of u2 is null. Actually, u1 is nothing but the Ornstein–Uhlenbeck process studied in Section 3.2. Notice that the choice of a psd model decreasing as ω−2 is actually not crucial for the proposed method. The limitations are just the same as those of the classical B/R decomposition; so the asymptotic behavior of Sf (ω) for large ω is simply arbitrary, as long as it decreases fast enough to create disjunct peaks in the bispectrum. Exact expressions of the amplification factors A2 and A3 are given by (23)–(24) for u1 ; similar expressions may be obtained for u2 but are not provided for the sake of conciseness. The dynamic amplification of the skewness Aγ3 is given by (62). Both A3 and Aγ3 are represented by solid lines in Fig. 6. Fig. 6(a) refers to the first loading (u1 ); Fig. 6(b) refers to the second loading (u2 ). Plots are presented with a log-scale for the damping coefficient so as to observe objectively various orders of magnitude of the damping. Also, plots are presented for various values of a/ω0 . Typical values encountered in wind engineering would range around a/ω0 ≃ 0.01. Results obtained for larger ratio are also presented here in order to better assess the limits of validity of the proposed method. The dynamic amplification of the third moment of the response expresses the ratio of the total response to the quasi-static one. It is therefore logical to observe a decrease towards a unit value for large ξ (no dynamic amplification). The behavior for small ξ is however different for the two loadings: the first loading creates a dynamic amplification A3 ≃ 3 whilst the second loading provides no amplification at all for small damping. This is due to the absence of low-frequency content in u2 , as may be readily deduced from ≪ζ ) (ξ the first asymptotic scenario A , see (64). Notice also that a 3 dynamic amplification A3 ≃ 3 is not worrying at all, especially when to the amplification of the second moment A2 =  compared  O ξ −1 for ξ ≪ 1. Indeed in the low frequency range, A2 gains

V. Denoël / Engineering Structures 33 (2011) 2271–2283




Fig. 6. Dynamic amplification of the third central moment of the response A3 (left) and of the skewness coefficient Aγ3 (right). They correspond to two different loadings (a, b) expressed as the square of Gaussian processes with psds given in (79). Notice that A3 is the ratio of the total response to the quasi-static one (B + bR) /B, while Aγ3 is the ratio of the skewness of the response to that of the loading. Among various possible approximations of the exact results (solid line), the θ -method (black dots) provides very accurate results.

the upper hand and for this reason Aγ3 shows a somehow similar monotonic profile for both loadings. Limit behaviors Aγ3 = 0 and Aγ3 = 1 for small and large damping ratio respectively are easily explained by a totally resonant response (and thus symmetrically distributed) on the one hand, and totally quasi-static on the other hand, presenting therefore the same skewness coefficient as the loading. Between these two limits, Aγ3 is almost monotonic, as expected following common sense. Notice the existence of a short range where Aγ3 is slightly larger than unity for the first loading;

this indicates that the skewness coefficient of the response is in this case slightly larger than that of the loading. 3 given by The proposed approximate dynamic amplification A (59), as well as Aγ3 deduced therefrom, are represented by dashed lines. The correspondence with the exact  result is almost perfect and is as good as ξ ≪ 1 and/or ζ = O ωa ≪ 1. In this particular 0 illustration, the bispectrum of the loading is known analytically 3 is computed analytically. and the integration in the definition of A In realistic applications however, numerical implementations of (59), i.e. approximations of the approximate solution, as


V. Denoël / Engineering Structures 33 (2011) 2271–2283

suggested in Section 4.3, need to be used. For clarity of the plots, only the results obtained with the numerical integration procedures are presented and N = 20 is chosen without discussion. Results are represented by black and white dots. The basic ω2 -method (white dots) provides acceptable approximations for ξ ≪ 1 and ζ ≪ 1 but starts drifting away from the exact solutions for typical operating values of ξ and ζ . On the contrary, the θ -method (black dots) shows a very good agreement with the exact solution, up to very large damping coefficients. This suggests that this method is most likely to be optimized with a number of integration points smaller than N = 20. 7. Conclusions The third central moment of the response of a structure to a non-Gaussian random loading is formally obtained by the double integration of the corresponding bispectrum in the (ω1 , ω2 )frequency space. Because of the acute shape of the bispectrum around the biresonant peaks and the ridges connecting them together, which is as marked as the damping ratio is small, the integration is heavy to perform numerically. All the more, the estimation of the integrand is typically time consuming in the framework of large finite element models. A similar problem occurs for the computation of the variance of the response, where the power spectral density (psd) needs to be integrated along frequencies. In this second order context, a famous solution consists in decomposing the integral into two terms, the background and resonant contributions. The cheap estimation of these components is based on only the variance of the loading as well as the estimation of the psd of the force at the natural frequency. This paper discussed the possibility of developing a similar decomposition for the third statistical order, in order to stand in for the heavy numerical integration of the bispectrum. As a first attempt, we examined to possibility to roughly extend the second order philosophy by replacing the bispectrum of the force by a constant value, so that the integrand of the double integral cuts down to the second Volterra kernel, and that the integral may be computed once for all. This basic transposition is however not successful: the integral of the second Volterra kernel is too slightly damping-dependent, compared to the actual third order moment. This suggests that the particular shape of the bispectrum of the force needs to be considered carefully, in the vicinity of the biresonant peaks. The lack of strong dampingdependence in the integral of the Volterra kernel is explained as  follows. The transfer function H (ω) behaves as O ξ −1 in a region ℓ = O (ξ ω0 ) around the resonance peaks. In the second order analysis, H 2 (ω)is integrated along ℓ which yields a resonant term of order O ξ −1 . In the third order analysis however, the peaks of Volterra kernel -which span a region of order ℓ2 each- are integrated. Because they correspond to resonance for only two of   the three factors, the integral is of order O ξ 0 , and so not strongly damping dependent. We then investigated the possibility of simplifying the double integral by considering this problem as a perturbed integral with two small parameters, namely the damping ratio ξ and a parameter ζ related to the distinctness of the frequency contents of the turbulence and of the dynamic response. The method, inspired from a well-known method for integration of analytical functions, consists in identifying the zones in the domain of integration that contribute the most to the integral. In an iterative way, each zone is considered; the integrand is cleverly approximated in order to yield a closed-form expression for that contribution; the approximate integrand is subtracted from the current integrand and following zones are appraised. In this process, we first identified a quasi-static contribution around the origin of the

frequency space, with the same meaning as that of the classical background component. The second contribution to the integral comes from biresonant peaks. It is expressed as an integral, which confirms indeed that the particular shape of the bispectrum of the force has to be considered. In practical applications, the estimation of this integral requires a numerical treatment, at least because the bispectrum of the force requires itself a numerical estimation. Among the local, global or mixed procedures discussed in Section 4.3, an integration method -referred to as the θ -method- appears to be very efficient for any given shape of bispectrum. In particular, thanks to an adequate change of variable, the θ -method also features the interesting convenience to come down on an exact solution for a small damping ratio. Two simple examples have shown that 20 integration points, i.e. 20 estimates of the bispectrum of the force, are necessary to accurately estimate the skewness coefficient of the response. This number has to be compared to a thousand or more integration points necessary for the numerical computation of the double integral. In this view, the many hours of computation that are nowadays required for a third order analysis (large FE model, standard 1core processor) turn into several minutes of computation, and thus make this kind of analysis accessible to everyday practice. Acknowledgments We would like to acknowledge the support of the Belgian Fund for Scientific Research who provided a perfect sterilized environment that allowed to formulate this problem properly. Many thanks go to Prof. A. Pierce (Dept of Mathematics, U. of British Columbia) and Prof. E. Detournay (Dept of Civil Engineering, U. of Minnesota) for valuable discussions at the important milestones of these developments. References [1] Grigoriu M. Stochastic calculus: applications in science and engineering. Springer Verlag, Birkhäuser; 2002. [2] Holmes JD. Wind loading on structures. 2nd ed. London: SponPress; 2007. [3] Simiu E, Scanlan R. Wind effects on structures. 3rd ed. John Wiley and Sons; 1996. [4] Davenport AG. The application of statistical concepts to the wind loading of structures. Proc Inst Civil Eng 1961;19:449–72. [5] Davenport AG. The response of slender, line-like structures to a gusty wind. Proc Inst Civil Eng 1962;23:389–408. [6] Ito K. Stochastic differential equations in a differentiable manifold. Nagoya Math J 1950;1:35–47. [7] Chen L, Letchford CW. A deterministic-stochastic hybrid model of downbursts and its impact on a cantilevered structure. Eng Struct 2004;26(5):619–29. doi:10.1016/j.engstruct.2003.12.009. [8] Chen L, Letchford CW. Numerical simulation of extreme winds from thunderstorm downbursts. J Wind Eng Ind Aerodyn 2007;95(9–11):977–90. doi:10.1016/j.jweia.2007.01.021. [9] Ibrahim RA, Soundararajan A, Heo H. Stochastic response of nonlinear dynamic systems based on a non-Gaussian closure. Trans ASME J Appl Mech 1985; 52(4):965–70. [10] Di Paola M, Falsone G. Stochastic dynamics of nonlinear systems driven by non-normal delta-correlated processes. Trans ASME J Appl Mech 1993;60(1): 141–8. [11] Mirri D, Iuculano G, Traverso PA, Pasini G, Filicori F. Non-linear dynamic system modelling based on modified Volterra series approaches. Measurement 2003;33(1):9–21. [12] Namachchivaya NS, Kok D, Ariaratnam ST. Stability of noisy nonlinear autoparametric systems. Nonlinear Dynam 2007;47(1–3):143–65. [13] Benfratello S, Falsone G, Muscolino G. Influence of the quadratic term in the alongwind stochastic response of sdof structures. Eng Struct 1996;18(9): 685–95. [14] Benfratello S, Falsone G. Non-Gaussian approach for stochastic-analysis of offshore structures. J Eng Mech, ASCE 1995;121(11):1173–80. [15] Grigoriu M. Response of linear systems to quadratic Gaussian excitations. J Eng Mech, ASCE 1986;112(6):523–35. [16] Kareem A, Tognarelli MA, Gurley KR. Modeling and analysis of quadratic term in the wind effects on structures. J Wind Eng Ind Aerodyn 1998;74(6): 1101–10. [17] Lutes LD, Hu SLJ. Nonnormal stochastic response of linear-systems. J Eng Mech, ASCE 1986;112(2):127–41. [18] Muscolino G. Linear systems excited by polynomial forms of non-Gaussian filtered processes. Probabilist Eng Mech 1995;10(1):35–44.

V. Denoël / Engineering Structures 33 (2011) 2271–2283 [19] Soize C. Gust loading factors with nonlinear pressure terms. J Struct Div, ASCE 1978;104:991–1007. [20] Ashraf Ali M, Gould PL. On the resonant component of the response of single degree-of-freedom systems under random loading. Eng Struct 1985;7(4): 280–2. [21] Clough RW, Penzien J. Dynamics of structures. 2nd ed. New-York: McGrawHill; 1993. [22] Gurley KR, Tognarelli MA, Kareem A. Analysis and simulation tools for wind engineering. Probabilist Eng Mech 1997;12(1):9–31. [23] Denoël V. Application of stochastic analysis methods to the study of the effects of wind on civil engineering structures. Ph.D. thesis. 2005. [24] Swami A, Giannakis GB, Zhou G. Bibliography on higher-order statistics. Signal Process 1997;60(1):65–126. [25] Benfratello S, Di Paola M, Spanos PD. Stochastic response of mdof wind-excited structures by means of Volterra series approach. J Wind Eng Ind Aerodyn 1998; 74(6):1135–45. [26] Carassale L, Kareem A. Dynamic analysis of complex systems by Volterra approach. In: Computational stochastic mechanics. 2003. p. 107–17. [27] Carassale L, Kareem A. Modeling nonlinear systems by Volterra series. J Eng Mech, ASCE 2010;136:801–18. [28] Denoël V, Degée H. Influence of the non-linearity of the aerodynamic coefficients on the skewness of the buffeting drag force. Wind Struct 2006; 9(6):457–71.


[29] Gusella V, Materazzi AL. Non-Gaussian along-wind response analysis in time and frequency domains. Eng Struct 2000;22(1):49–57. [30] Uhlenbeck GE, Ornstein LS. On the theory of the Brownian motion. Phys Rev 1930;36(5):823. Copyright (C) 2010 The American Physical Society Please report any problems to [email protected] PR. [31] Grigoriu M, Ariaratnam ST. Response of linear systems to polynomials of Gaussian processes. Trans ASME J Appl Mech 1988;55(4):905–10. [32] Krantz SG. The residue theorem (Â Section 4.4.2). In: Handbook of complex variables. Boston (MA): Birkauser; 1999. [33] Bender CM, Orszag SA. Advanced mathematical methods for scientists and engineers. new ed. 1999. [34] Vanmarcke E. Properties of spectral moments with applications to random vibration. J Eng Mech Div, ASCE 1972;98(2):425–46. [35] Hinch EJ. Perturbation methods, vol. 1. Cambridge: Cambridge University Press; 1991. [36] Denoël V. Limit analysis of the statistics of quasi-steady non-linear aerodynamic forces for small turbulence intensities. Probabilist Eng Mech 2009; 24(4):552–64. doi:10.1016/j.probengmech.2009.04.001. [37] Runge C. Über empirische funktionen und die interpolation zwischen äquidistanten ordinaten. Z Math Phys 1901;46:224–43. [38] Berrut J-P, Trefethen LN. Barycentric lagrange interpolation. SIAM Rev 2004; 46(3):501–17.