On the Shaker Simulation of Wind-Induced Non-Gaussian Random Vibration

Gaussian signal is produced by ordinary random vibration controllers to test the products in the laboratory, while the field data is usually non-Gaussian. Two methodologies are presented in this paper for shaker simulation of wind-induced non-Gaussian vibration. The first methodology synthesizes the non-Gaussian signal offline and replicates it on the shaker in the TimeWaveform Replication (TWR) mode. A new synthesis method is used to model the non-Gaussian signal as a Gaussian signal multiplied by an amplitude modulation function (AMF). A case study is presented to show that the synthesized non-Gaussian signal has the same power spectral density (PSD), probability density function (PDF), and loading cycle distribution (LCD) as the field data.The second methodology derives a damage equivalent Gaussian signal from the non-Gaussian signal based on the fatigue damage spectrum (FDS) and the extreme response spectrum (ERS) and reproduces it on the shaker in the closed-loop frequency domain control mode. The PSD level and the duration time of the derived Gaussian signal can be manipulated for accelerated testing purpose. A case study is presented to show that the derived PSD matches the damage potential of the non-Gaussian environment for both fatigue and peak response.


Introduction
Random vibration testing is usually used to bring a test item to failure to identify weaknesses in the product or to verify if the product can survive a particular random vibration environment.Historically, random vibration controllers accomplish this goal by producing a power spectral density (PSD) that would expose the test item to the type of vibratory environment that the test item would experience in a real-world setting.However, a random process can be characterized by a PSD only when it is stationary and Gaussian distributed.In practice, non-Gaussian vibration is usually encountered, especially in the road transportation [1,2].The difference of the probability density function (PDF) between the measured acceleration and Gaussian distribution may lead to totally different accumulated fatigue damage [3].MIL-STD810 points out that we should always check if the field tested data is non-Gaussian and the testing hardware and software is appropriate [4].The Time Waveform Replication (TWR) is frequently referred to as a methodology for non-Gaussian testing [4].The basic idea of TWR is to reproduce a sequence of instantaneous values of the vibration process.Such a test may be non-Gaussian; however, this is only a replication of one particular measured record, not a simulation of a specified road type.Smallwood presented a zero memory nonlinear function (ZMNL) method to synthesize a given non-Gaussian data [5].Analytical solution was presented later to extend the original offline method to a method that can be used in the closed-loop frequency domain control mode.Iteration method was presented to address the PSD distortion problem introduced by the nonlinear transformation [6]. A. Steinwolf presented a phase selection method to simulate the road vehicle measured data [7].In particular, selected phase is transformed from random to deterministic in order to obtain a prescribed kurtosis.PSD and PDF are controlled independently.An analytical relation between kurtosis, amplitude, and phase at specific frequencies was presented later to make this method applicable in a closed-loop control implemented by the shake controller [8][9][10].These methods mentioned above are applicable when the non-Gaussianity is caused by significant shocks in vehicle vibration data.For wind-induced vibration data, however, the non-Gaussianity may be caused by other factors, for example, by the running root mean square.Rouillard and Sek presented a novel technique to synthesize non-Gaussian vibrations by generating a sequence of random Gaussian processes of varying RMS levels and durations [11].
Two methodologies are presented in this paper for shaker simulation of wind-induced vibration.The essential idea of the first methodology is to synthesize the non-Gaussian signals offline and replicate them on the shaker in the Time Waveform Replication mode.The non-Gaussian signal is modeled as a Gaussian signal multiplied by an amplitude modulation function (AMF).A two-parameter Weibull distribution is used to create the AMF.A case study is presented to show that the synthesized non-Gaussian signal has the same power spectral density (PSD), probability density function (PDF), and loading cycle distribution (LCD) as the field data.The rainflow cycle counting (RFCC) method [12] and Dirlik method [13] are used for the LCD calculation.
The fatigue damage spectrum (FDS) and the extreme response spectrum (ERS) have been used as a quantitative way to take the fatigue damage and overstress damage into consideration for field data [14][15][16][17][18][19][20].The FDS and ERS are used in the second methodology to derive a damage equivalent Gaussian signal (characterized by a PSD and time duration) from the non-Gaussian signal, so that they can be easily produced on the shaker in the closed-loop frequency domain control mode.Improvement of the derived PSD by taking its spectrum shape into consideration is also presented.The PSD level and the duration time of the derived Gaussian signal can be manipulated for accelerated testing purpose.A case study is presented to show that the derived PSD matches the damage potential of non-Gaussian environment for both fatigue and peak response.

Theory
2.1.Non-Gaussian Random Signal.If a Gaussian process has a zero mean value and is ergodic, then the PDF of the instantaneous acceleration values () is given by the Gaussian distribution with zero mean: where   is standard deviation.For large durations , the variance  2  is given by where   () is the single-sided PSD.It shows that random Gaussian processes with zero mean can be completely described by the PSD function.One useful method for establishing how well a random process can be described by the Gaussian distribution is by

Response motion
Input motion Figure 1: SDOF systems for calculation of FDS.
computing the higher order moments of the process defined as (continuous and discrete form) When the mean value is zero, where  is the skewness and  is the kurtosis.For a truly Gaussian process, the skewness is 0 and the kurtosis is 3. Any deviation from these indicates that the process is non-Gaussian.From the perspective of simulating non-Gaussian vibration, kurtosis is a more important parameter than skewness, because it represents the probability of peak values in time history.

Response Spectra.
All structures experience fatigue as they are repeatedly exposed to adequate levels of stress.The total damage a product experiences in a particular time period can be calculated from field data and plotted for a specific range of frequencies.The FDS is essentially a plot that shows the responses of a number of single degree of freedom (SDOF) systems to a base input acceleration time history.Many SDOF systems that are tuned to a range of natural frequencies are assessed using the same input.The final plot, the FDS, shows the fatigue damage encountered for a particular SDOF system anywhere within the analyzed time, as shown in Figures 1 and 2.
It has been shown that the stress is roughly proportional to pseudovelocity [21].Thus, it is more reasonable to calculate the pseudovelocity FDS.For a SDOF system with a natural frequency   and a damping ratio , the output pseudovelocity  pv to an input acceleration   can be calculated using a ramp invariant digital filter method [16]: where   is the sampling frequency and  filter indicates filtering the input signal with a ramp invariant digital filter.
With the output  pv , the cumulative damage can be calculated in both time domain and frequency domain.In time domain, RFCC method is usually used for stress cycle counting.With cyclic numbers of different stress level, - curve and Miner's rule are combined to calculate fatigue damage [22] as follows: where   is the number of cycles exposure at   ( = 1, 2, . . ., ),  is the number of stress levels considered,   is fatigue life at stress   ,  is the fatigue exponent, and   is the total damage index calculated in time domain.
Assume that the stress is proportional to pseudovelocity as follows: The total damage can be calculated as In frequency domain, Rayleigh distribution of response stress maxima is assumed and used to calculate the stress cycles where  is the stress value of peaks and   is the RMS of the stress time history.The total damage can be calculated as where  +  is the number of positive maxima per unit time in the stress history, which can be replaced by the natural frequency   for a lightly damped SDOF system response and  is the total time of exposure to the stress environment.
Solving (10) with ( 9) leads to where Γ is the Gamma function and  pv is the RMS of pseudovelocity.
The RMS of pseudovelocity can be calculated using where () is the transmissibility of a SDOF system (pseudovelocity/acceleration).
If the input PSD is relatively flat in the half power band width of each SDOF system, then Mile's equation can be used to calculate  pv in a closed form approximately as follows: With ( 11) and ( 13), the PSD at each natural frequency can be derived from the FDS at the corresponding natural frequency: Note that if the flat spectrum assumption does not hold around some natural frequencies, the FDS calculated with (11) and ( 12) using the derived PSD in (14) will deviate from the FDS calculated with ( 11) and ( 13).In such case, an iteration procedure is needed to take the shape of spectrum into account; see Section 2.4.
The ERS shows the damage caused by vibration from an overstress viewpoint.Like the FDS, the ERS is essentially a plot that shows the responses of a number of SDOF systems to a base input acceleration time history.The only difference is that the ERS is generated by calculating the maximum response of a SDOF system to the input.The final plot, the ERS, shows the largest response encountered for a particular SDOF system anywhere within the analyzed time.For a SDOF system with a natural frequency   and a damping ratio , the ERS can also be calculated in both time domain and frequency domain [23] as follows: In this paper, the FDS and ERS are used to determine the damage equivalent PSD and test duration.
Since the essential idea of the new synthesis method is to model the non-Gaussian signal by a Gaussian signal multiplied by an AMF, the field data must be analyzed to see if it contains many shocks.Shocks must be removed if they are the main reason to cause the non-Gaussianity and processed with separate shock tests.To detect if there are shocks in the field data, the response of a number of SDOF systems to the field data is firstly calculated in time domain.Then, the response is divided into  consecutive segments and the maximum value is calculated for each segment.Finally, the ratio between this maximum value and the ERS computed from the PSD, which is calculated from the whole time history using Welch method, is used to show where the shocks are [24]: where  , indicates the maximum value calculated for the th segment and the th SDOF system and ERS  indicates the SRS calculated form PSD for the th SDOF system.

Synthesis of Non-Gaussian Signal.
The procedure for synthesizing wind-induced non-Gaussian vibration signal consists of three phases.
(1) Check the characteristics of measured vibration signal.
Calculate the kurtosis of the measured signal to see if it is non-Gaussian signal.Then, calculate the PSD of the measured signal using Welch's method.Use (19) to detect shocks.
(2) Synthesize the non-Gaussian random vibration.Create a stationary Gaussian signal using the PSD calculated in phase 1. Calculate the probability distribution function  RMS of the running RMS of the field data.Use a two-parameter Weibull distribution to model the probability distribution of the running RMS.The cumulative distribution function of Weibull distribution can be expressed as where  is the random variable,  is the scale parameter, and  is the shape parameter.
To find the correct  and , the probability distribution function of the assumed Weibull distribution  Weibull and the running RMS  RMS are used to form an objective function using least squares method: where  are the values of the running RMS, depending on how many bins are used to calculate the probability distribution function of the running RMS,  indicates the th iteration, and std indicates standard deviation.Note that both  RMS and  Weibull are vectors.
The MATLAB function fminsearch is used to find the right  and , so that the standard deviation of error between these two probability distribution functions is minimized.
The initial AMF is essentially a vector with  elements following the estimated Weibull distribution: where  is the total number of elements in the initial AMF,  total is the total time of the field data, and   (correlation time) is the time interval (in second) that is used to construct the AMF.
The AMF is then resampled to the same length of the Gaussian signal created from the PSD.The non-Gaussian signal is synthesized by multiplying the AMF with the Gaussian signal.Using the kurtosis of the field data as an objective function, update the AMF from the Weibull distribution until the kurtosis of the field data is reached by the synthesized non-Gaussian signal.The synthesized non-Gaussian signal is scaled to have the same RMS as the field data.
(3) Compare PSD, PDF, and LCD of the synthesized non-Gaussian signal and the field data.

Derivation of Damage Equivalent Gaussian Signal.
To avoid overtesting or undertesting the specimen, both FDS and ERS of the damage equivalent Gaussian signal should be updated so that they match the FDS and ERS of the field data.The procedure of deriving the damage equivalent Gaussian signal from the non-Gaussian signal can be implemented as follows.
(1) Calculate the FDS from the field measured non-Gaussian signal in time domain, and resample the field data if needed.
As pointed out by Lalanne [25],  does not affect the PSD of derived equivalent Gaussian signal, and a smaller value of  leads to more conservative result.As recommended in MIL-STD 810G [4],  = 10 and  = 4 are used to calculate the pseudovelocity FDS from the field data.
(2) Calculate the ERS from the time history in time domain.
(3) Predetermine a test time duration and use FDS from step (1) and ( 14) to derive an initial PSD.
(5) Use an iteration procedure to update the derived PSD in step (3) if the error between two FDS in step (1) and step (4) is significant where (  ) indicates the PSD,  indicates the iteration number,  env (  ) indicates the FDS envelope of field data, and  ,num (  ) indicates the FDS calculated from the PSD in iteration number  using numerical integration as shown in (11) and ( 12).(6) Calculate the ERS with (updated) PSD in frequency domain using (16) and compare it with that calculated in step (2).
(7) Manipulate the level of PSD and the test time duration to match both FDS and ERS between derived Gaussian signal and measured non-Gaussian signal.for synthesizing the non-Gaussian vibration signal.The random acceleration signal induced by wind is taken from a test item installed on a mast, as shown in Figure 3, along with the running RMS.The test item is a dummy box containing electronics.The dummy box was equipped with B&K WB0179 triaxial accelerometer set and a wind speed measuring device.The box was a standard Ericsson Micro Radio Base equipment, manufactured in cast lightweight alloy with integrated heat sinks.The dimensions are about 530 * 400 * 185 mm and the weight was about 21 kg.The back side of the dummy box was mounted to the mast at a height of 50 meters.The running RMS is calculated using 10 s frames with 50% overlapping.This is a non-Gaussian signal with kurtosis equaling 9.4262.Equation ( 19) is used to detect shocks (see Figure 4).The ratio is shown in the third dimension by the color.As we can see from Figure 4, there are only few shocks superimposed on the underlying random vibration signal, which means that the non-Gaussianity is caused by the running RMS.The PSD of the field data is calculated using Welch method with 4096-point Hanning window and 50% overlapping.A Gaussian signal is generated using the same PSD.The PDF of both the field data and the synthesized Gaussian signal are shown in Figure 9.As we can see, the PDF of the field data is obviously different from that of the synthesized Gaussian signal.The PDF of running RMS is calculated using different number of bins and is shown in Figure 5. From Figure 5, we can see that despite of the numbers of bins used the PDF of the running RMS is similar to the PDF of a given Rayleigh distribution.To be more general, a two-parameter Weibull distribution is used to model the probability distribution function of the running RMS.

Case Study
Iteration method described in Section 2.3 (phase 2) is used to determine the  and , as shown in Figure 6.
As the parameters of the Weibull distribution are determined, the AMF can be created and the non-Gaussian signal can be synthesized using the iteration method presented in Section 2.3 (phase 2).The synthesized non-Gaussian signal is shown in Figure 7, together with the field data (the time duration of the synthesized non-Gaussian signal is adjustable).
The PSD and PDF of the field data and the synthesized non-Gaussian signal are very similar (see Figures 8 and 9), which proves the effectiveness of the new method.
For comparison, the cycle counting is performed on the field data, the synthesized Gaussian, and non-Gaussian signal, using the RFCC and Dirlik method.The results are shown in Figure 10.As we can see, the Dirlik method is equivalent to the RFCC for Gaussian signal.Big error can be introduced by the Dirlik method for non-Gaussian signal.The RFCC results for the field data and synthesized non-Gaussian signal are very similar.
Care must be taken when creating the AMF.The AMF can be generated by using different correlation time and being resampled.The AMFs using correlation time of 0.1 s, 10 s, and 100 s are shown in Figure 11.The corresponding synthesized    time signals are shown in Figure 12.From Figure 12, we can see that inappropriate choice of the correlation time (0.1 s and 100 s) leads to dissimilar time histories to the field data.
The PSD are calculated for the field data and the synthesized non-Gaussian signals using different correlation time (see Figure 13).From Figure 13, we can see that the PSD of the synthesized non-Gaussian signal with a correlation time of 0.1 s differs from others at low frequencies (between 1 Hz and 10 Hz).This is caused by new significant frequency components introduced by the AMF when it is created.This indicates that small correlation time may lead to inaccurate PSD.The PDF are calculated for the field data, synthesized Gaussian, and non-Gaussian signals using different correlation time (see Figure 14).From Figure 14, we can see that the PDF of the synthesized non-Gaussian signal using 100 s correlation time shows different characteristics from other non-Gaussian signals.The middle of this PDF is shaper and the tails are narrower.This indicates that large correlation time may lead to inaccurate PDF.

Derivation of Damage Equivalent Gaussian Signal.
In this case study, the same non-Gaussian signal from Section 3.1 is used to validate the procedure of generating damage equivalent Gaussian signal.Figure 15 shows the pseudovelocity FDS calculated from the field data.
Using the same time duration as the field data (1200 s), an equivalent Gaussian PSD is derived from the FDS using (14).To account for the shape of the PSD, the FDS is calculated again using (11) and ( 12) and compared with the FDS that is used to derive the PSD.Using the derived PSD as an initial value, ( 23) is used to update the derived PSD.The FDS

Figure 2 :
Figure 2: Calculation procedure of the FDS.

Figure 5 :
Figure 5: PDF of the running RMS using different number of bins and Rayleigh distribution.
Std of error between probability distribution functions

Figure 9 :
Figure 9: The PDF of the field data and synthesized Gaussian and non-Gaussian signal.

Figure 10 :
Figure 10: Number of cycles versus data ranges for RFCC and Dirlik.