- Email: [email protected]
Designing random vibration tests
Pergamon
Inr J Non-Linear
A4echanic.s. Vol. 31, No. 5, pp. 563-572, 1996 Copyright 0 1996 Elsevier Science Ltd Pnnted m Great Britam. All rights reserved 0020-7462/96 $15.00 + 0.00
SOOZO-7462(96)00021-2
DESIGNING
RANDOM VIBRATION TESTS Ludomir
Politechnika
Rzeszowska,
M. Laudafiski
ul. Iranka-Osmeckiego
15 m 4, 35 106 Rzeszow,
Poland
Abstract-The
paper provides information about two concepts of digital simulation of samples of the stationary Gaussian stochastic processes possessing multi-modal spectra. They have been applied in order to imitate dynamic loads arising on an airplane undergoing gusty flying conditions. Therefore the particular spectra typical for an airplane gust study were involved, reflecting also elastic properties of the flying vehicle. The presented details are devoted to the problem of solving the system of algebraic non-linear equations describing the desired linear filter. At this stage the given results can be also applied in studying earthquakes, modelling gusty winds for civil engineering and other purposes. Copyright ICI 1996 Elsevier Science Ltd.
Keywords: random
vibration,
1. INTRODUCTION-A
digital
simulation,
airplane
SURVEY
ON
loads
THE
STATE
OF
THE
ART
In general there are two competitive approaches to the problem of designing random vibration tests. The first one is based upon the method of simulation related to the earliest investigation in applying stationary random processes to the engineering problems which was initiated exactly half a century ago by Rice [l], and which is known as the harmonic method. The second approach-known as thejliltration method-probably originated independently in the West (USA) and in the East (USSR) and conveniently is associated with the names of Wiener and Kchinchin accordingly. Both approaches were subsequently and extensively developed throughout the decades, therefore we are able to mention here only a few papers [l-9] belonging to this wide scope of investigations which approximately 20 years ago reached a mature stage of experimental application. Seemingly the method of harmonic (based upon specialized processors of FFT) became a leading method in order to feed the specially designed actuators shaking not only particular components of a design, but also the entire structures. One such application is aeronautical engineering, and from this particular field we display some details in Fig. 1, see [ll]. Leaving out the details given in this figure we again turn towards the general description of the methods used to simulate the samples of the stochastic processes necessary to manage such an experiment. Papers by Tunik [2,3] provide the reader with many details related to the harmonic method employed in experimental projects imitating random vibrations as they were developed in USSR. It is interesting to note that ref. [2] states that if the number of harmonic components used in the digital simulation of the desired samples of the stationary stochastic processes becomes of order N = 10-20, then the probability distributions defined upon them are from the practical point of view identical with the Gaussian distributions. On the other hand Yang proved in ref. [4] that the number of harmonics which secure the 2-a interval of the Gaussian distribution is of the order of N = 500 components. This leaves open the question about the desired accuracy of the samples simulated in this way. But there is a point which can be undertaken in favour of the methods of filtration: the samples generated in this way may omit the problem of the periodicity which characterizes the samples generated by the method of harmonics. Also in simulating two-dimensional random field samples we have here two possibilities: either to go through the procedure of the appropriate cross-correlations, or to use the samples of the two-dimensional white noise random field (see Fig. 2, repeated from ref. [lo]). Its suitable cross-sections may provide all the one-dimensional samples of white noise, preserving required time lags and cross-correlations. It is especially useful and straightforward if we have in mind the so-called eight point model of an airplane [8]. 563
564
L. M. Laudakki
Fig. 1. Random
vibration
Fig. 2. Two-dimensional
laboratory
test layout
sample of white noise.
We would like to mention here also the latest works by Spanos [S, 61 where the recursive filter technique was extensively developed. It seems that Spanos successfully overcame the difficulties related to the recursive filters and their usage in the simulation technique. The approach of using the recurstve filter for simulation purposes seemingly can be dated back as far as to the paper by Sragovitch 171. The author of this paper tried to use both the Pade and Tchebyshev polynomials as means to determine both components of the recursive filter before it can be applied to the eventual simulation procedure. Due to the method developed by Sragovitch, this intermediate stage requires expression of a spectral density function of the target process by the rational fraction composed by two polynomial expansions. Unfortunately we failed in such an approach, not being able to secure numerically stable simulation procedures. Simply speaking, due to this method each value of the simulated sample became the result of the subtraction of two numbers: the problem lay in the fact that these two numbers were by two orders bigger than the simulated one. Therefore to obtain the sample values within a reasonable accuracy, the accuracy in evaluating both terms of
Designing
random
vibration
tests
565
.
4;
30
50
A F=O-50
Fig. 3. Spectral
N=lO
density
IP=l
M-18
L=200
of the gust response
U=200
DELF=O.OSO
of the PZL gravity.
zsc
M-18 agricultural
W-W
A/P at its centre
-J
of
the recursive filter was beyond the scope of the practical possibilities. This negative experience with recursive filtration pushed us towards the method of deriving the nonrecursive filters presented in the paper. A typical spectrum for airplanes/gliders looks like the one presented in Fig. 3, showing the response calculated at the centre of gravity for the PZL M-18 (a Polish agricultural airplane), which flies at a speed of 100 km h- ’ horizontally in vertical gusts characterized by the scale of turbulence L = 50 m (which means rather severe conditions). We see for the particular airframe typical peaks associated with the sensitive frequencies reflecting its resonances at the particular point chosen for it. We call such spectra multi-modal spectra, and the processes associated with them multi-modal processes, accordingly. For the Gaussian stochastic processes possessing such spectral densities we developed special technique of simulation of the appropriate sample functions. It seemingly fills a gap in the wide class of publications related to this subject (for other references see also ref. [S]). Below we present the essence of two numerical methods aiming at such a purpose. They both have a common origin and we begin by showing this origin by stating the numerical problem which has to be solved. Then, we show two algorithms solving this problem. In the end we give a summary of the latest numerical results obtained by us in this field.
2. THE
FILTER
EQUATION
We begin with the definition of the correlation function which is related to the properties of the linear systems while examining their input-output relations. We define the correlation function for the output stochastic process {y(t)), which left linear devices called filters being fed by the input stochastic process {x(t)}:
K(z) = 6
2 k(k) x(t k=O
- k) ; n=O
k(n)x(t
+ z-n).
1
(1)
There is now an interesting innovation: we shall look at the output process with the intention of recovering the impulse characteristic of the filter k(t) by assuming that the input
L. M. Laudahski
566
was a white noise process. Which after all simplifications K(z) = 1
k(n)k(n
+ z)
with
leads to the final form:
r = 0, 1,2,.
. . , N.
n=O
To see the problem we shall face in this way exactly it is recommended (2) into its explicit form: K(0) = k(O)k(O) + k(l)k(l)
to expand
equation
+ k(2)k(2) + . . . + k(N)k(N)
K(1) = k(O)k(l) + k(l)k(2) + k(2)k(3) + . . K(2) = k(O)k(2) + k(l)k(3) + k(2)k(4) + .
K(N - 2) = k(O)k(N - 2) + k(l)k(N
+ k(N - l)k(N) . + k(N - 2)k(N)
- 1) + k(2)k(N)
K(N - 1) = k(O)k(N - 1) + k(l)k(N) K(N) = k(O)k(N).
(3)
Equations (3) represent a system of N algebraic equations of order two each. The equations are linearly independent, i.e. their Jacobian is non-zero. Therefore in general the system has 2N solutions, say c(~, a2,. . , ulN, which can be both real and complex. Each solution of equations (3) preserving the form k’(l), k’(2), . , k’(N) composes the filter, although it is not clear at the beginning whether each solution of equations (3) possesses physical meaning, i.e. can be considered as the physically realized filter. There are some basic properties of the solutions which have to be mentioned. If hi (1) hi (2), . . . , k’(N) is the solution of the system (3), then k’(N), k’(N - l), . . . , hi (1) also becomes the solution. Also, solutions of the form -k’(l), -hi (2), . . . , -k’(N) and -k’(N), -k’(N - l), . . ,-k’(1) satisfy equations (3). If we call the solution k’(l), k’ (2), . . , k’(N) the basic solution, we can state that in general there will always be N/2 basic solutions, either real or complex. Each system (3) has its normalized correspondent obtained by dividing (3) side-by-side with value K(1). It is easy to prove that normalized solutions can be obtained dividing ordinary solutions
by m.
3. TWO
POINTS
FILTER
Let us consider a special case of equation (3) i.e. situation when N = 2. For this case one can easily find the analytical solution. Some details of the procedure are shown below. The system of equations (3) now has the form: K(1) = k(l)k(l)
+ k(2)k(2)
K(2) = k(l)k(2).
(4)
Let us denote: D1 = JK(1) With this notation
+ 2K(2)
we shall get analytical
and
D2 = JK(1)
solutions
- 2K(2).
(5)
of (4) in the form:
k(l) = f (01 + D,) k(2) = f (DI - D2).
(6)
It is clear now that once: K(1) < 2K(2) value D2 becomes imaginary, therefore the corresponding solution complex. Such a solution evidently has no physical meaning.
(7) of equation
(4) becomes
Designing
random
It is also easy to see that with notation solutions:
tests
given by equation
[-Ml), -wJ)l,
C&3 wn
CW),WI>
vibration
561
(6) the system
(4) has four
C-m), -WI.
(8)
Let us acknowledge by the way that the condition (7) assessing existence of the real solutions of equation (4) can be physically interpreted as a case with sufficiently rapidly decreasing correlation function K(r) with respect to increasing value of its argument.
4. NUMERICAL
TEST
Testing and debugging stages require reliable and suitable numerical examples. In this case a special value has examples offering solutions expressed by integers. We enclose below two such examples. For the case N = 2 we use: K(1) = 34 having
the basic solutions
and
K(2) = 15
(9)
as below: k(l) = 5,
h(2) = 3.
(IO)
For the case N = 4 we use: K(1) = 30, which possesses
K(2) = 20,
K(3) = 11,
a single exact basic solution
of the form:
k(2) = 2,
k(3) = 3,
h(1) = 1,
and another basic solution, although this time approximate, rounded to two decimal digits after comma: h(l)r
1.65,
k(2) z 1.58,
5. ALGORITHM
k(3) E 4.35,
K(4) = 4
h(4) = 4
(11)
(12)
which we give here in a form
k(4)%2.42.
(13)
MESORS
In the beginning the algebraic problem stated by equation (3) is reformulated transformed into a geometrical one. To do that we define the norm: I/K(z) - I?(z) 11= =f K(z) - ‘f k(n)k(n + z) ’ r=O n=O
and
(14)
with the different notations for the source process { Y(t)} and the simulated one called the target process {p(t)}. Now some formal requirement concerning the accuracy can be specified by demanding: 111((T)- I?(z) /I d 6. (15) Therefore N + 1 values k(n) can be understood as the arguments of an N + l-dimensional function defined by equation (14). With the condition (15) we seek the minima of this multi-dimensional geometrical object. It was done numerically by resorting to the Svejgaard algorithm, the main idea of which is based on a gradient searching approach and was developed following the Algal implementation given in refs [12, 151. It is also shown in Fig. 4 where one can find some indications as to how the iterations proceed for a chosen step in searching for the minimum for a two-dimensional problem used as an illustrative example.
6. ALGORITHM
ANNA
The algorithm solving the system (3) that we are going to present now does it in a direct way. The essence in directly solving system (3) lies in finding such a solution LX” E B?’ for
L. M. Laudafiski
568
Lake filter
‘.,,
‘..
‘,.,, ,,,./’
h(l)
,:’
Fig. 4. An idea of a single iteration search by Mesors
which mapping F E 2” x 2’ becomes zero. Non-linear character application of an iterative method based upon an algorithm:
of mapping
xk+ ’ = xk - [J- ‘lx”)] F(xk)
F suggests the
(16)
where J is the Jacobian of the mapping matrix F, xk the kth approximation Practical implementation of the algorithm (16) which we call Anna following theorem [13].
of the vector x. is based on the
Theorem Let the function F(x) become differentiable according to Frechet within the neighbourhood K(cw,p) of such a point x that K(U) = 0. Moreover the first derivative of F(x) is a continuous function at the same point c( and non-singular. With these circumstances the point OLbecomes an attractive point of the Newton iterative method: X k+l
[J-‘(xk)]F(xk).
= xk -
Strong assumptions imposed upon the function F(x+ontinuity, differentiability and non-singularity of its derivative at the point a-guarantee the local convergence of the method. The system (3) satisfies all the requirements mentioned. Methods based on the above Theorem having the Jacobian given numerically are called quasi-Newtonian methods. One particular example of such a method offers the hybrid method of Powell [14]. Here the Jacobian is calculated by the method of finite differences. In order to improve the convergence this method supposes moreover: /IF(xk+’ )1/z > Xk+l
To find pk, which is necessary to determine system of equations given below:
=
IIF(xk)llz
Xk + pk.
(18)
the next iteration,
J(xk)pk = - F(xk). The idea of the algorithm explained as:
(17)
Anna implementing
iterations
xk+l must be solved by the
(19) with respect to k can be briefly
Designing random vibration tests
569
1st step: If F(xk) = 0 =+ STOP 2nd step: Calculate
pk by solving
equation
(19)
3rd step: If I/F(xk + p”) I/2 < 11 F(xk) I/2 the step is accepted X
k+l
=
xk
+
pk
then k=zk + 1 and return
resulting
in:
to step 1 otherwise:
4th step: Calculations
7. ON NUMERICAL
are interrupted.
RESULTS,
CONCLUSIONS
The algorithm Mesors was implemented into a program written first in Fortran IV and widely used in numerous calculations by resorting to the electronic computers of the generations IBM 360, PDP 1 l/70 and CDC 6600. They imposed hard restrictions upon the volume of the problem which did not reach at that time above N = 32. The arrival and fast development of personal computers gave rise to new implementations by using Fortran 77 and by using special compilers like the NDP Fortran. Both implementations, the earlier from the mid seventies and the actual one, were successful in the sense that they lead to right solutions-which became completely evident lastly by testing this algorithm with the examples described above in paragraphs 3 and 4 developed only now. To complete successful computations for N = 64 by using an IBM PC 486-50 requires about 3-5 minutes. Increase in time goes approximately in such a way that doubling the volume N leads to about 10 times longer computations. Some other information about the computer time can be read in Fig. 5, summarizing our numerous investigations in this field. It is worth mentioning that the nature of the calculations within the airplane gust response studies needs even values N = 512 or may sometimes be as big as N = 1024. Completing the solution to the system (3) by using the Svejgaard algorithm became a difficult problem of reaching the desired accuracy for the results derived in this way. In particular it concerns those distanced (latest) components of the approximate numerical solution which are usually hundreds of times smaller than the greatest initial values. We guess that the accuracy in their estimation may become dramatically low. The question is: whether they have to be derived exactly? Or in other words-which level of errors in their estimations can be accepted, and which not? And for these questions we have not had a satisfying answer until now. It seems reasonable to demonstrate finally an example of the results derived in this way. We use this opportunity to show the results which are related to the spectrum displayed by Fig. 3. They describe the airplane gust response study. First in Fig. 6 we present two filters which were obtained for the rigid airplane approximation and for the elastic airplane approximation. Then we show in Fig. 7 two sample records derived for those two filters. They follow the time history of the vertical acceleration computed for an airplane centre of gravity. The algorithm Anna was developed at the end of 1993 and our numerical experience with its application to the particular aeronautical purposes is just at its beginning (for a concise description see ref. [16]). There is no convincing way to compare both algorithms quite literally. Seemingly the direct solution of Anna goes about 10 times faster than geometrical approach which follows Mesors. In both cases, the choice of the initial point is crucial for the successful solution. It is so far only possible to say that Mesors almost always produces a solution of equations (3) disregarding the particular choice of the initial point. There are some doubts-as we said before-about the accuracy of these solutions. Nevertheless
L. M. Laudahki
570
1OOmin
/NDP
/
486
lmin
20
O.lmin
Fig. 5. Computer
100
-..-._- 80
60
40
time to solve five iterations
by Mesnrs.
(a> 0.4 0.3
N=2
0.2 z =
0.1 I\
(b)
-0.1 -
-o.20
I 0.1
I 0.2
I 0.4
I 0.3
t Fig. 6. Non-recursive
filters.
I 0.5
I 0.6
120
N
Designing
random
vibration
tests
571
(a) 0.3 -
= ;
0.7 0.6 0.5 0.4 0.3 0.2 0.1 -0.‘: -0.2 -0.3 -0.4 -OS 0.5
1.0
1.5 t
Fig. 7. Records
of simulated
samples.
throughout lengthy numerical practice there were only a few cases of divergent behaviour observed. On the other hand, Anna solves equations (3) only when the choice of the initial point is luckily done, so there were only few cases that the solution was obtained. The solutions produced by Anna go significantly faster than solution by using Meson and their accuracy is high. Both algorithms still remain as a potential field for gaining further experience and results. It would have been useful to compare the spectrum yielded by the computer filter to the target spectrum, as shown for instance in Fig. 3 (as it was suggested by both the reviewers of the paper). Such a comparison remains one of the future problems to be solved. To assess the accuracy of the final product i.e. to state whether the simulated sample of the stochastic process under attention conforms to some requirements related to its geometry, numerical characteristics such as zeros, maxima and inflection points have been used. They were counted directly from the simulated sample (such as that shown in Fig. 7) and independently calculated from the theoretical formulae as given by Rice [l]. The agreement between these two approaches was quite satisfying, also corresponding to the cases shown in Fig. 3, possessing the multi-modal spectral densities. Such a point of view seems to be natural, bearing in mind the application of the simulated samples to run random vibrations tests in order to predict the fatigue life of the airplanes. The idea was discussed at length in ref. [S]. Acknowledgement~The author would like to thank A. A. Tunik for his kind attention to the problems undertaken in the paper and remarks on the harmonic method. The author also would like to thank P. D. Spanos for the information about the ARMA method given during the discussion in the course of the colloquium. The author wishes to thank both reviewers of this paper for their interesting remarks which led to the removal of some imperfections and helped the author in the critical evaluation of the paper. The results presented in this paper, derived between 1992 and 1993, were financed by the Polish Scientific Research Council KBN grant no. 770039203. REFERENCES 1. S. 0. Rice, Mathematical analysis of random noise. Bell System Tech. J. 23, 282-332 (1944); 24,466156 (1945). 2. A. A. Tunik, An Automatic Control of the Vibration Tests. The Energy, Moscow (1978) (in Russian). 3. A. A. Tunik, On algorithms of the digital simulation of the random multi-dimensional processes with given spectral matrices. Cybernet. Comput. Tech. J. Ukrain. Acad. Sci. 28, 98-107 (1975) (in Russian). 4. J.-N. Yang, On the normality and accuracy of simulated random processes. JSV 26,417428 (1973). 5. M. P. Mignolet and P. D. Spanos, Recursive simulation of stationary multivariate random processes-Part I. Trans. ASME J. Appl. Mech. 54, 6744680 (1987). 6. P. D. Spanos and M. P. Mignolet, Recursive simulation of stationary multivariate random processes-Part II. Trans. ASME J. Appl. Mech. 54, 681-687 (1987). 7. V. G. Sragovitch, Modelling a class of the random processes. J. Comput. Math. Math. Phys. 3, 5866593 (1963) (in Russian). 8. L. M. Laudafiski, An Introduction to Stochastic Dynamics of Gliders. Rzeszow Technical University Press, Rzeszow (1978) (in Polish). 9. M. Shinozuka and C.-M. Jan, Digital simulation of random processes and its simulation. JSV 25, 1 (1972).
572
L. M. Laudanski
10. P. A. Robinson, The modelling of turbulence and downbursts for flight simulators. UTIAS Report No. 339, Downsview, Toronto (1991). 11. Testing of Materials and Structuresfor the Aerospace Industry. Special Publication no. 2217/2e, Schenck, Darmstadt (1982). 12. W. M. Turski, Usaye of Digital Computers in Scientijc Engineering Centres. PWN, Warsaw (1973) (in Polish). 13. 2. Fortuna, B. Macugow and J. Wysowski, Numerical Methods. WNT, Warsaw (1982) (in Polish). 14. M. J. D. Powell, A hybrid method for nonlinear equations. Methods for Nonlinear Algebraic Equations, P. Rabinowitz (ed.). Gordon and Breach (1970). 15. L. Laudaiiski and K. Kunysz, On digital simulation of Gaussian loads. In Proc. 8th Int. Conf. Non-Linear Oscillations, Prague, pp. 4299434 (1978). 16. A. Kucaba-Pietal and L. Laudatiski, Modelling stationary Gaussian loads. In Proc. XXXIV Symp. Modelling in Mechanics, Wisla, Zeszyty Naukowe Politechniki Slaskiej, Seria: Mechanika, Vol. 121, pp. 173-182 (1995).
Inr J Non-Linear
A4echanic.s. Vol. 31, No. 5, pp. 563-572, 1996 Copyright 0 1996 Elsevier Science Ltd Pnnted m Great Britam. All rights reserved 0020-7462/96 $15.00 + 0.00
SOOZO-7462(96)00021-2
DESIGNING
RANDOM VIBRATION TESTS Ludomir
Politechnika
Rzeszowska,
M. Laudafiski
ul. Iranka-Osmeckiego
15 m 4, 35 106 Rzeszow,
Poland
Abstract-The
paper provides information about two concepts of digital simulation of samples of the stationary Gaussian stochastic processes possessing multi-modal spectra. They have been applied in order to imitate dynamic loads arising on an airplane undergoing gusty flying conditions. Therefore the particular spectra typical for an airplane gust study were involved, reflecting also elastic properties of the flying vehicle. The presented details are devoted to the problem of solving the system of algebraic non-linear equations describing the desired linear filter. At this stage the given results can be also applied in studying earthquakes, modelling gusty winds for civil engineering and other purposes. Copyright ICI 1996 Elsevier Science Ltd.
Keywords: random
vibration,
1. INTRODUCTION-A
digital
simulation,
airplane
SURVEY
ON
loads
THE
STATE
OF
THE
ART
In general there are two competitive approaches to the problem of designing random vibration tests. The first one is based upon the method of simulation related to the earliest investigation in applying stationary random processes to the engineering problems which was initiated exactly half a century ago by Rice [l], and which is known as the harmonic method. The second approach-known as thejliltration method-probably originated independently in the West (USA) and in the East (USSR) and conveniently is associated with the names of Wiener and Kchinchin accordingly. Both approaches were subsequently and extensively developed throughout the decades, therefore we are able to mention here only a few papers [l-9] belonging to this wide scope of investigations which approximately 20 years ago reached a mature stage of experimental application. Seemingly the method of harmonic (based upon specialized processors of FFT) became a leading method in order to feed the specially designed actuators shaking not only particular components of a design, but also the entire structures. One such application is aeronautical engineering, and from this particular field we display some details in Fig. 1, see [ll]. Leaving out the details given in this figure we again turn towards the general description of the methods used to simulate the samples of the stochastic processes necessary to manage such an experiment. Papers by Tunik [2,3] provide the reader with many details related to the harmonic method employed in experimental projects imitating random vibrations as they were developed in USSR. It is interesting to note that ref. [2] states that if the number of harmonic components used in the digital simulation of the desired samples of the stationary stochastic processes becomes of order N = 10-20, then the probability distributions defined upon them are from the practical point of view identical with the Gaussian distributions. On the other hand Yang proved in ref. [4] that the number of harmonics which secure the 2-a interval of the Gaussian distribution is of the order of N = 500 components. This leaves open the question about the desired accuracy of the samples simulated in this way. But there is a point which can be undertaken in favour of the methods of filtration: the samples generated in this way may omit the problem of the periodicity which characterizes the samples generated by the method of harmonics. Also in simulating two-dimensional random field samples we have here two possibilities: either to go through the procedure of the appropriate cross-correlations, or to use the samples of the two-dimensional white noise random field (see Fig. 2, repeated from ref. [lo]). Its suitable cross-sections may provide all the one-dimensional samples of white noise, preserving required time lags and cross-correlations. It is especially useful and straightforward if we have in mind the so-called eight point model of an airplane [8]. 563
564
L. M. Laudakki
Fig. 1. Random
vibration
Fig. 2. Two-dimensional
laboratory
test layout
sample of white noise.
We would like to mention here also the latest works by Spanos [S, 61 where the recursive filter technique was extensively developed. It seems that Spanos successfully overcame the difficulties related to the recursive filters and their usage in the simulation technique. The approach of using the recurstve filter for simulation purposes seemingly can be dated back as far as to the paper by Sragovitch 171. The author of this paper tried to use both the Pade and Tchebyshev polynomials as means to determine both components of the recursive filter before it can be applied to the eventual simulation procedure. Due to the method developed by Sragovitch, this intermediate stage requires expression of a spectral density function of the target process by the rational fraction composed by two polynomial expansions. Unfortunately we failed in such an approach, not being able to secure numerically stable simulation procedures. Simply speaking, due to this method each value of the simulated sample became the result of the subtraction of two numbers: the problem lay in the fact that these two numbers were by two orders bigger than the simulated one. Therefore to obtain the sample values within a reasonable accuracy, the accuracy in evaluating both terms of
Designing
random
vibration
tests
565
.
4;
30
50
A F=O-50
Fig. 3. Spectral
N=lO
density
IP=l
M-18
L=200
of the gust response
U=200
DELF=O.OSO
of the PZL gravity.
zsc
M-18 agricultural
W-W
A/P at its centre
-J
of
the recursive filter was beyond the scope of the practical possibilities. This negative experience with recursive filtration pushed us towards the method of deriving the nonrecursive filters presented in the paper. A typical spectrum for airplanes/gliders looks like the one presented in Fig. 3, showing the response calculated at the centre of gravity for the PZL M-18 (a Polish agricultural airplane), which flies at a speed of 100 km h- ’ horizontally in vertical gusts characterized by the scale of turbulence L = 50 m (which means rather severe conditions). We see for the particular airframe typical peaks associated with the sensitive frequencies reflecting its resonances at the particular point chosen for it. We call such spectra multi-modal spectra, and the processes associated with them multi-modal processes, accordingly. For the Gaussian stochastic processes possessing such spectral densities we developed special technique of simulation of the appropriate sample functions. It seemingly fills a gap in the wide class of publications related to this subject (for other references see also ref. [S]). Below we present the essence of two numerical methods aiming at such a purpose. They both have a common origin and we begin by showing this origin by stating the numerical problem which has to be solved. Then, we show two algorithms solving this problem. In the end we give a summary of the latest numerical results obtained by us in this field.
2. THE
FILTER
EQUATION
We begin with the definition of the correlation function which is related to the properties of the linear systems while examining their input-output relations. We define the correlation function for the output stochastic process {y(t)), which left linear devices called filters being fed by the input stochastic process {x(t)}:
K(z) = 6
2 k(k) x(t k=O
- k) ; n=O
k(n)x(t
+ z-n).
1
(1)
There is now an interesting innovation: we shall look at the output process with the intention of recovering the impulse characteristic of the filter k(t) by assuming that the input
L. M. Laudahski
566
was a white noise process. Which after all simplifications K(z) = 1
k(n)k(n
+ z)
with
leads to the final form:
r = 0, 1,2,.
. . , N.
n=O
To see the problem we shall face in this way exactly it is recommended (2) into its explicit form: K(0) = k(O)k(O) + k(l)k(l)
to expand
equation
+ k(2)k(2) + . . . + k(N)k(N)
K(1) = k(O)k(l) + k(l)k(2) + k(2)k(3) + . . K(2) = k(O)k(2) + k(l)k(3) + k(2)k(4) + .
K(N - 2) = k(O)k(N - 2) + k(l)k(N
+ k(N - l)k(N) . + k(N - 2)k(N)
- 1) + k(2)k(N)
K(N - 1) = k(O)k(N - 1) + k(l)k(N) K(N) = k(O)k(N).
(3)
Equations (3) represent a system of N algebraic equations of order two each. The equations are linearly independent, i.e. their Jacobian is non-zero. Therefore in general the system has 2N solutions, say c(~, a2,. . , ulN, which can be both real and complex. Each solution of equations (3) preserving the form k’(l), k’(2), . , k’(N) composes the filter, although it is not clear at the beginning whether each solution of equations (3) possesses physical meaning, i.e. can be considered as the physically realized filter. There are some basic properties of the solutions which have to be mentioned. If hi (1) hi (2), . . . , k’(N) is the solution of the system (3), then k’(N), k’(N - l), . . . , hi (1) also becomes the solution. Also, solutions of the form -k’(l), -hi (2), . . . , -k’(N) and -k’(N), -k’(N - l), . . ,-k’(1) satisfy equations (3). If we call the solution k’(l), k’ (2), . . , k’(N) the basic solution, we can state that in general there will always be N/2 basic solutions, either real or complex. Each system (3) has its normalized correspondent obtained by dividing (3) side-by-side with value K(1). It is easy to prove that normalized solutions can be obtained dividing ordinary solutions
by m.
3. TWO
POINTS
FILTER
Let us consider a special case of equation (3) i.e. situation when N = 2. For this case one can easily find the analytical solution. Some details of the procedure are shown below. The system of equations (3) now has the form: K(1) = k(l)k(l)
+ k(2)k(2)
K(2) = k(l)k(2).
(4)
Let us denote: D1 = JK(1) With this notation
+ 2K(2)
we shall get analytical
and
D2 = JK(1)
solutions
- 2K(2).
(5)
of (4) in the form:
k(l) = f (01 + D,) k(2) = f (DI - D2).
(6)
It is clear now that once: K(1) < 2K(2) value D2 becomes imaginary, therefore the corresponding solution complex. Such a solution evidently has no physical meaning.
(7) of equation
(4) becomes
Designing
random
It is also easy to see that with notation solutions:
tests
given by equation
[-Ml), -wJ)l,
C&3 wn
CW),WI>
vibration
561
(6) the system
(4) has four
C-m), -WI.
(8)
Let us acknowledge by the way that the condition (7) assessing existence of the real solutions of equation (4) can be physically interpreted as a case with sufficiently rapidly decreasing correlation function K(r) with respect to increasing value of its argument.
4. NUMERICAL
TEST
Testing and debugging stages require reliable and suitable numerical examples. In this case a special value has examples offering solutions expressed by integers. We enclose below two such examples. For the case N = 2 we use: K(1) = 34 having
the basic solutions
and
K(2) = 15
(9)
as below: k(l) = 5,
h(2) = 3.
(IO)
For the case N = 4 we use: K(1) = 30, which possesses
K(2) = 20,
K(3) = 11,
a single exact basic solution
of the form:
k(2) = 2,
k(3) = 3,
h(1) = 1,
and another basic solution, although this time approximate, rounded to two decimal digits after comma: h(l)r
1.65,
k(2) z 1.58,
5. ALGORITHM
k(3) E 4.35,
K(4) = 4
h(4) = 4
(11)
(12)
which we give here in a form
k(4)%2.42.
(13)
MESORS
In the beginning the algebraic problem stated by equation (3) is reformulated transformed into a geometrical one. To do that we define the norm: I/K(z) - I?(z) 11= =f K(z) - ‘f k(n)k(n + z) ’ r=O n=O
and
(14)
with the different notations for the source process { Y(t)} and the simulated one called the target process {p(t)}. Now some formal requirement concerning the accuracy can be specified by demanding: 111((T)- I?(z) /I d 6. (15) Therefore N + 1 values k(n) can be understood as the arguments of an N + l-dimensional function defined by equation (14). With the condition (15) we seek the minima of this multi-dimensional geometrical object. It was done numerically by resorting to the Svejgaard algorithm, the main idea of which is based on a gradient searching approach and was developed following the Algal implementation given in refs [12, 151. It is also shown in Fig. 4 where one can find some indications as to how the iterations proceed for a chosen step in searching for the minimum for a two-dimensional problem used as an illustrative example.
6. ALGORITHM
ANNA
The algorithm solving the system (3) that we are going to present now does it in a direct way. The essence in directly solving system (3) lies in finding such a solution LX” E B?’ for
L. M. Laudafiski
568
Lake filter
‘.,,
‘..
‘,.,, ,,,./’
h(l)
,:’
Fig. 4. An idea of a single iteration search by Mesors
which mapping F E 2” x 2’ becomes zero. Non-linear character application of an iterative method based upon an algorithm:
of mapping
xk+ ’ = xk - [J- ‘lx”)] F(xk)
F suggests the
(16)
where J is the Jacobian of the mapping matrix F, xk the kth approximation Practical implementation of the algorithm (16) which we call Anna following theorem [13].
of the vector x. is based on the
Theorem Let the function F(x) become differentiable according to Frechet within the neighbourhood K(cw,p) of such a point x that K(U) = 0. Moreover the first derivative of F(x) is a continuous function at the same point c( and non-singular. With these circumstances the point OLbecomes an attractive point of the Newton iterative method: X k+l
[J-‘(xk)]F(xk).
= xk -
Strong assumptions imposed upon the function F(x+ontinuity, differentiability and non-singularity of its derivative at the point a-guarantee the local convergence of the method. The system (3) satisfies all the requirements mentioned. Methods based on the above Theorem having the Jacobian given numerically are called quasi-Newtonian methods. One particular example of such a method offers the hybrid method of Powell [14]. Here the Jacobian is calculated by the method of finite differences. In order to improve the convergence this method supposes moreover: /IF(xk+’ )1/z > Xk+l
To find pk, which is necessary to determine system of equations given below:
=
IIF(xk)llz
Xk + pk.
(18)
the next iteration,
J(xk)pk = - F(xk). The idea of the algorithm explained as:
(17)
Anna implementing
iterations
xk+l must be solved by the
(19) with respect to k can be briefly
Designing random vibration tests
569
1st step: If F(xk) = 0 =+ STOP 2nd step: Calculate
pk by solving
equation
(19)
3rd step: If I/F(xk + p”) I/2 < 11 F(xk) I/2 the step is accepted X
k+l
=
xk
+
pk
then k=zk + 1 and return
resulting
in:
to step 1 otherwise:
4th step: Calculations
7. ON NUMERICAL
are interrupted.
RESULTS,
CONCLUSIONS
The algorithm Mesors was implemented into a program written first in Fortran IV and widely used in numerous calculations by resorting to the electronic computers of the generations IBM 360, PDP 1 l/70 and CDC 6600. They imposed hard restrictions upon the volume of the problem which did not reach at that time above N = 32. The arrival and fast development of personal computers gave rise to new implementations by using Fortran 77 and by using special compilers like the NDP Fortran. Both implementations, the earlier from the mid seventies and the actual one, were successful in the sense that they lead to right solutions-which became completely evident lastly by testing this algorithm with the examples described above in paragraphs 3 and 4 developed only now. To complete successful computations for N = 64 by using an IBM PC 486-50 requires about 3-5 minutes. Increase in time goes approximately in such a way that doubling the volume N leads to about 10 times longer computations. Some other information about the computer time can be read in Fig. 5, summarizing our numerous investigations in this field. It is worth mentioning that the nature of the calculations within the airplane gust response studies needs even values N = 512 or may sometimes be as big as N = 1024. Completing the solution to the system (3) by using the Svejgaard algorithm became a difficult problem of reaching the desired accuracy for the results derived in this way. In particular it concerns those distanced (latest) components of the approximate numerical solution which are usually hundreds of times smaller than the greatest initial values. We guess that the accuracy in their estimation may become dramatically low. The question is: whether they have to be derived exactly? Or in other words-which level of errors in their estimations can be accepted, and which not? And for these questions we have not had a satisfying answer until now. It seems reasonable to demonstrate finally an example of the results derived in this way. We use this opportunity to show the results which are related to the spectrum displayed by Fig. 3. They describe the airplane gust response study. First in Fig. 6 we present two filters which were obtained for the rigid airplane approximation and for the elastic airplane approximation. Then we show in Fig. 7 two sample records derived for those two filters. They follow the time history of the vertical acceleration computed for an airplane centre of gravity. The algorithm Anna was developed at the end of 1993 and our numerical experience with its application to the particular aeronautical purposes is just at its beginning (for a concise description see ref. [16]). There is no convincing way to compare both algorithms quite literally. Seemingly the direct solution of Anna goes about 10 times faster than geometrical approach which follows Mesors. In both cases, the choice of the initial point is crucial for the successful solution. It is so far only possible to say that Mesors almost always produces a solution of equations (3) disregarding the particular choice of the initial point. There are some doubts-as we said before-about the accuracy of these solutions. Nevertheless
L. M. Laudahki
570
1OOmin
/NDP
/
486
lmin
20
O.lmin
Fig. 5. Computer
100
-..-._- 80
60
40
time to solve five iterations
by Mesnrs.
(a> 0.4 0.3
N=2
0.2 z =
0.1 I\
(b)
-0.1 -
-o.20
I 0.1
I 0.2
I 0.4
I 0.3
t Fig. 6. Non-recursive
filters.
I 0.5
I 0.6
120
N
Designing
random
vibration
tests
571
(a) 0.3 -
= ;
0.7 0.6 0.5 0.4 0.3 0.2 0.1 -0.‘: -0.2 -0.3 -0.4 -OS 0.5
1.0
1.5 t
Fig. 7. Records
of simulated
samples.
throughout lengthy numerical practice there were only a few cases of divergent behaviour observed. On the other hand, Anna solves equations (3) only when the choice of the initial point is luckily done, so there were only few cases that the solution was obtained. The solutions produced by Anna go significantly faster than solution by using Meson and their accuracy is high. Both algorithms still remain as a potential field for gaining further experience and results. It would have been useful to compare the spectrum yielded by the computer filter to the target spectrum, as shown for instance in Fig. 3 (as it was suggested by both the reviewers of the paper). Such a comparison remains one of the future problems to be solved. To assess the accuracy of the final product i.e. to state whether the simulated sample of the stochastic process under attention conforms to some requirements related to its geometry, numerical characteristics such as zeros, maxima and inflection points have been used. They were counted directly from the simulated sample (such as that shown in Fig. 7) and independently calculated from the theoretical formulae as given by Rice [l]. The agreement between these two approaches was quite satisfying, also corresponding to the cases shown in Fig. 3, possessing the multi-modal spectral densities. Such a point of view seems to be natural, bearing in mind the application of the simulated samples to run random vibrations tests in order to predict the fatigue life of the airplanes. The idea was discussed at length in ref. [S]. Acknowledgement~The author would like to thank A. A. Tunik for his kind attention to the problems undertaken in the paper and remarks on the harmonic method. The author also would like to thank P. D. Spanos for the information about the ARMA method given during the discussion in the course of the colloquium. The author wishes to thank both reviewers of this paper for their interesting remarks which led to the removal of some imperfections and helped the author in the critical evaluation of the paper. The results presented in this paper, derived between 1992 and 1993, were financed by the Polish Scientific Research Council KBN grant no. 770039203. REFERENCES 1. S. 0. Rice, Mathematical analysis of random noise. Bell System Tech. J. 23, 282-332 (1944); 24,466156 (1945). 2. A. A. Tunik, An Automatic Control of the Vibration Tests. The Energy, Moscow (1978) (in Russian). 3. A. A. Tunik, On algorithms of the digital simulation of the random multi-dimensional processes with given spectral matrices. Cybernet. Comput. Tech. J. Ukrain. Acad. Sci. 28, 98-107 (1975) (in Russian). 4. J.-N. Yang, On the normality and accuracy of simulated random processes. JSV 26,417428 (1973). 5. M. P. Mignolet and P. D. Spanos, Recursive simulation of stationary multivariate random processes-Part I. Trans. ASME J. Appl. Mech. 54, 6744680 (1987). 6. P. D. Spanos and M. P. Mignolet, Recursive simulation of stationary multivariate random processes-Part II. Trans. ASME J. Appl. Mech. 54, 681-687 (1987). 7. V. G. Sragovitch, Modelling a class of the random processes. J. Comput. Math. Math. Phys. 3, 5866593 (1963) (in Russian). 8. L. M. Laudafiski, An Introduction to Stochastic Dynamics of Gliders. Rzeszow Technical University Press, Rzeszow (1978) (in Polish). 9. M. Shinozuka and C.-M. Jan, Digital simulation of random processes and its simulation. JSV 25, 1 (1972).
572
L. M. Laudanski
10. P. A. Robinson, The modelling of turbulence and downbursts for flight simulators. UTIAS Report No. 339, Downsview, Toronto (1991). 11. Testing of Materials and Structuresfor the Aerospace Industry. Special Publication no. 2217/2e, Schenck, Darmstadt (1982). 12. W. M. Turski, Usaye of Digital Computers in Scientijc Engineering Centres. PWN, Warsaw (1973) (in Polish). 13. 2. Fortuna, B. Macugow and J. Wysowski, Numerical Methods. WNT, Warsaw (1982) (in Polish). 14. M. J. D. Powell, A hybrid method for nonlinear equations. Methods for Nonlinear Algebraic Equations, P. Rabinowitz (ed.). Gordon and Breach (1970). 15. L. Laudaiiski and K. Kunysz, On digital simulation of Gaussian loads. In Proc. 8th Int. Conf. Non-Linear Oscillations, Prague, pp. 4299434 (1978). 16. A. Kucaba-Pietal and L. Laudatiski, Modelling stationary Gaussian loads. In Proc. XXXIV Symp. Modelling in Mechanics, Wisla, Zeszyty Naukowe Politechniki Slaskiej, Seria: Mechanika, Vol. 121, pp. 173-182 (1995).