High Efficiency Computation of the Variances of Structural Evolutionary Random Responses

For structures subjected to stationary or evolutionary white/colored random noise, their various response variances satisfy algebraic or differential Lyapunov equations. The solution of these Lyapunov equations used to be very difficult. A precise integration method is proposed in the present paper, which solves such Lyapunov equations accurately and very efficiently.


Introduction
Many engineering structures, e.g. bridges, tallbuildings, TV towers, offshore platforms, are very sensitive to the ambient random loads such as windgusts, earthquakes, irregular waves, etc. The analysis methodology for these problems have received much attention, and many contributions have been published (e.g. [1,2,8,9,14]). For most problems, the variance analysis is of major interest. However, if the structures are subjected to evolutionary random excitations, such analysis is quite difficult. Recently, Lin et al. [3][4][5][6][7] developed the pseudo excitation method (PEM) which has enabled the computation of such problems rather conveniently. Unfortunately, for structures subjected to wide band excitations, a great number of frequency points still have to be taken in order to obtain accurate variances, and the computation would still be rather time-consuming. Further improvement is still required.
In past few years, Zhong, X.N. [10] derived the closed solutions of the algebraic and differential Lyapunov equations based on the precise numerical computation of the exponential matrix exp([A]t) proposed by Zhong, W.X. [11][12][13]. In the present paper, these methods have been combined and extended to the solution of the above mentioned variances. Two examples are given which show that the proposed method is not only accurate, but also extremely efficient.

Differential Lyapunov equation
The equations of motion of an n DOF (degrees of freedom) structure subjected to an evolutionary random excitation {f (t)} can be expressed as Eq. (1) can be expressed as where (3) will vanish, and so its solution would be [3,4] [Φ(t, Usually, the initial conditions can be regarded as uncorrelated with the excitations, so that Eq. (5) leads to the covariance matrix That leads to the differential Lyapunov equation

Precise integration for Lyapunov equation
Consider the homogeneous differential Lyapunov equations The values [P 0 ] at initial time t 0 are known. It is readily verified that its solution is Next, let and is the particular solutions of the non-homogeneous equations (11). Therefore, the solutions of Eq. (11) with non-zero initial values should be In the following, only the case with [P 0 ] = [0] (when t 0 = 0) is considered. The non-zero initial-value solution can be separately solved and then be linearly added to the zero-initial-value solution in order to produce the complete solution. Obviously, at time t, therefore at time t + τ Similarly, the double-step solution (i.e. that at time 2t) can be obtained, being The formulas for some simple modulation functions are given in the following. It is always assumed that when t < 0, g(t) = 0. The response variances are going to be computed at time τ, 2τ, 3τ, . . .

For suddenly applied white noise
Similarly, at time 2t k , the double step solution is

For time dependent modulation function
Provided This matrix has the form at time t k + τ and the corresponding double-step form would be The integration formats for other modulation functions g(t), so long as they are not too complicated, can be similarly derived. Note that the above formulas are all exact. The key problem is how to ensure high precision on computers.

Precise integration of Lyapunov equation and its implementation on computers
In order to achieve precise integration on computers, the step size τ should be divided into 2 N equal intervals, N = 20 is usually adequate [10]. Take the values of [Φ(η)] and [P (η)] at time η = τ/s N as the initial values for the step-by-step integration. Because η is extremely small, therefore taking only the first 3 or 4 terms in the related Taylor's expansions would achieve extremely high precision. This is the essence of the precise integration approach. The first step in implementation is to execute the double step formulas for N times starting from the initial time t = 0, to produce their values at time τ . Secondly, the precise integration is executed with the same step size τ to generate the response variances at time 2τ , 3τ , . . . The implementation process for g(t) = √ t is given here in detail. For other modulation functions, the processes are quite similar, and will be referred to the users.
With N = 20, the errors due to neglecting the higher order terms in the Taylor's series are all of the order of O(η 5 ) = O(τ 5 /10 30 ), which is within the doubleprecision of available computers.
3. Repeat the following process for N times, to work out the variance matrix [P ] at time τ : The above process is for zero-initial value problems. The linear superposition should be executed for the contribution of non-zero initial conditions.

Numerical examples
Example 1. The equation of motion of a single degree of freedom system is Its RHS represents a suddenly applied filtered white noise excitation, i.e. g(t) = 1.0 (when t 0) and the PSD of the stationary random process x(t) is For simplicity, dimensionless parameters are used here. And it is assumed that ω = ω g = 1.0, ζ = ζ g = 0.5 and S 0 = 1.0. For this simple example, analytical solutions are available [9], i.e. the variances of y anḋ y are P 11 = 1.5π = 4.71238898 . . . and P 22 = π = 3.14159265 . . .  This example has also been computed in terms of the proposed scheme. Note that when the response variances at time τ have been obtained (in the computation, τ was arbitrarily taken as 0.0015707963); the doublestep formula, Eq. (24), was repeatedly used to compute the response variances at times 2τ, 4τ, 8τ, . . . until convergence is achieved. Through only 14 such passes, i.e. at t = 16384τ , the responses reach the stable values P 11 = 4.71238898 . . . and P 22 = 3.14159265 . . . Clearly, their first 9 effective digits are both identical to those of the analytical solutions.
Assume that the system is subjected to a suddenly applied filtered white noise excitation, i.e. when t 0, where the PSD of x(t) has the form of Eq. (37) with S 0 = 142.75, ω g = 19.07 and ζ g = 0.544. The integration step size was τ = 0.02 and 200 steps were computed. The displacement and velocity variance curves of the three masses are shown in Figs 2 and 3. Computation shows that when t = 4.0, the variances of such displacements are {0.5605, 1.807, 2.803}. This example was previously computed [7] where the pseudo excitation method was used; and when t = 4.0, the first four digits of these variances are still {0.5605 1.807 2.803}, exactly the same as those computed here. Some differences occur after their fourth digits because a limited frequency band ω ∈ [−200, 200] was used previously [7], while the variances based on the Lyapunov equation correspond to the integration frequency band ω ∈ (−∞, ∞).
Four different methods have been used [7] to compute the present example. The precision and CPU time used for all the five methods (including the present one) are listed in Table 1. It can be seen that previously the best way is the HPD-S method which used 226 s(CPU) to obtail the accurate result for σ 2 y3 . However, the method proposed in the present paper used only 16 s to produce the same result. All computations were executed on a 586 personal computer with main frequency 133 MHz.
For this example, if the modulation function is replaced by g(t) = √ t (when t 0), and repeat the above computations, the corresponding results are shown in Figs 4 and 5. Comparing Fig. 2 and Fig. 4 shows that when t 1, the responses due to g(t) = √ t are obviously smaller than those due to g(t) = 1.0. When t 2, however, inverse phenomenon takes place. That is because the responses due to g(t) = √ t increase much faster after t > 1 and so they demonstrate much stronger non-stationary property. The CPU used was 19 s, about 20% more than that for g(t) = 1.0.
Although only proportional damping matrix [C] was used in the above, the process would be all the same if any non-proportional damping matrix is used.

Conclusions
The solution of algebraic or differential Lyapunov equations is a basic problem in the random vibration field. This problem has long been considered as very difficult to resolve. Based on the previously developed precise integration method, i.e. its HPD-L and HPD-S forms, the authors now further propose a new precise integration scheme, which can solve such Lyapunov equations accurately and very efficiently. The excitations can be stationary or non-stationary (evolutionary), white or colored random noise, therefore the proposed method applies to a wide range of problems.