RESPONSE STATISTICS OF OCEAN STRUCTURES TO NON-LINEAR HYDRODYNAMIC LOADING PART II: NON-GAUSSIAN OCEAN WAVES

Journal of Sound and Vibration (1996) 191(1), 107–128

RESPONSE STATISTICS OF OCEAN STRUCTURES TO NON-LINEAR HYDRODYNAMIC LOADING. PART II: NON-GAUSSIAN OCEAN WAVES N. M and R. A. I Department of Mechanical Engineering, Wayne State University, Detroit, MI 48202, U.S.A. (Received 7 November 1994, and in final form 21 February 1995) In this second part, the response statistics of ocean elastic systems to non-linear hydrodynamic loading represented by a non-Gaussian random process are considered. The non-Gaussian process is generated from a non-linear filter excited by a white noise process. The filter non-linearity is represented by a gradient of a potential function where an exact closed form solution of the stationary probability density exists. The filter non-linearity is selected such that its response probability density function exhibits the skewness feature of ocean waves. In order to estimate the filter power spectrum, a joint probability density function of the filter response co-ordinates at two different times is derived, using an orthogonal expansion technique together with the Fokker–Planck equation. This process yields a first order differential equation governing the time evolution of the filter response statistics. This equation forms an infinite hierarchy of coupled equations which are solved up to rank 10. The power spectra of the filter are obtained for different filter and excitation parameters. The response statistics of a linear elastic structure subjected to non-Gaussian non-linear hydrodynamic forces are determined using the stochastic averaging method and, alternatively, the orthogonal expansion method. 7 1996 Academic Press Limited

1. INTRODUCTION

The random response of non-linear single-degree-of-freedom systems has been extensively examined in the literature. These systems are usually modelled by a non-linear differential equation of the form q¨+2zvn q˙+vn2q+oC(q, q˙ , q¨ )=j(t),

(i)

where j(t) is an external excitation, C may include non-linear stiffness, damping, and inertia terms, and o is a small parameter. The Duffing oscillator is a special class of system (i) when C is cubic in q. The early study of a Duffing oscillator under random excitation goes back to the work of Lyon et al. [1]. They considered j(t) as a narrow band Gaussian process derived from wide-band excitation of a resonant filter centered at frequency v1 : j +2zf v1 j +v12j=W(t),

(ii)

where W(t) is white noise. The purposes of their study was to discover if the well-known jump phenomenon can occur under narrow-band random excitation. The response was studied by using the quasi-linearization method and by experimental testing. The results revealed multi-valued response characteristics which have the same general appearance as those for sinusoidal forcing, except that the peaks are much sharper. The same feature 107 0022–460X/96/110107+22 $18.00/0

7 1996 Academic Press Limited

108

.   . . 

has been observed by Lennox and Kuak [2]. They used a quasi-static method with small but finite bandwidth. However, their method does not allow one to investigate the influence of bandwidth. Fang and Dowell [3] analyzed the response of a Duffing oscillator to a narrow band random excitation by numerical simulation. Their results showed that multi-valued mean square responses can occur for mono-level excitations with very narrow bandwidths. As the bandwidth increases, the multi-valued responses give way to single-valued responses. Similar results were obtained using many different approaches, such as, the averaging method [4–7], stochastic linearization [8–11], the probabilistic linearization technique [12], numerical simulation [4–6, 9, 12–14], and experimental testing [1]. Zhu et al. [14] have shown that for a combination of system and excitation parameters all the statistics of the stationary response are unique and independent of the initial conditions. In a certain domain of the parameter space there are two possible motions in the stationary response, and jumps may occur. A jump is a transition from one probable motion to another or vice versa. Outside of this domain the stationary response is either nearly Gaussian or like a diffused limit cycle. As the parameters cross the boundary of the domain the qualitative behavior of the stationary response changes and it is a special form of bifurcation. An alternative approach has been proposed by Davies and Nandlall [15]. They modelled the excitation as the response of a lightly-damped second order filter to white noise. Thus the four-dimensional Fokker–Plank–Kolmogorov (FPK) equation can be written for the response. An approximate time-dependent solution based on a Gaussian closure scheme was obtained from the four-dimensional FPK equation. The result was smoothed time histories for the mean square displacement and velocity. When plotted in a phase plane, the solution showed two stable attractors and a saddle point for the case in which the excitation had very narrow bandwidth, in a certain frequency range. This behavior is related to the well-known jump phenomenon. As the bandwidth of the excitation increased, the phase plane was reduced to a single stable attractor (sink), indicating the elimination of the jump phenomenon. The FPK equation for the response of a Duffing oscillator to filtered excitation has been solved numerically by Kapitaniak [16] using the path-integral method. He presented the stationary probability density function of the stochastic process as bimodal. Davies and Liu [5] investigated the same system via an averaging method and numerical simulation. They showed that the probability density function obtained for narrow bandwidth excitation has multiple local maxima (multi-furcation) with random jumps. Also, for wider bandwidth excitation the probability density function shows only a single local maximum. A similar result was obtained by Iyengar [12] using numerical simulation. Koliopulos and Bishop [17] studied a method for the prediction of the response statistics of non-linear oscillators subjected to narrow band random noise based on a quasi-harmonic assumption. They showed that the method can easily incorporate the possible jumps of the system via the input–response amplitude curve. This curve defines critical input values which define the region of possible jumps. In this case the probability density of the response amplitude exhibits two peaks. Iyengar [18, 19] has presented a stochastic stability analysis of a non-linear hardening-type oscillator subjected to narrow band random excitation. He showed that only one statistical moment is stable and, hence, the multi-valued predictions of the equivalent linearization analysis do not correspond to the real behavior of the system. Koliopulos and Langley [20] improved the stability analysis for the same system. They showed that when jumps between competing response states occur, the choice of an appropriate value for the kurtosis of the response is crucial for a reliable estimation of the local statistical moments and for a more accurate stochastic stability analysis than the

 

109

equivalent linearization method. Furthermore, a modification is proposed which takes into account the influence of certain non-linear terms in the variational equation for the case of strong nonlinearities. Under white noise random excitation, the well-known jump phenomenon in the response behavior will disappear. In addition, under Gaussian random excitation the response of the non-linear oscillator is certainly non-Gaussian. In this case one may take advantage of a shaping filter designed to simulate ocean waves. The output of such a filter can be used to describe the hydrodynamic forces acting on immersed ocean structures. The present study is a continuation of previous work [21] dealing with the random response of ocean systems to Gaussian waves. Most of the work is devoted to designing a non-linear shaping filter whose probabilistic characteristics emulate ocean waves. The random response of a linear elastic structure to a non-linear non-Gaussian hydrodynamic force is analyzed using the stochastic averaging technique. 2. STATEMENT OF THE PROBLEM

Consider a linear single-degree-of-freedom system subjected to non-linear hydrodynamic forces. The motion of such a system is usually described by the equation of motion

6

7

1 x¨+2z0 v0 x˙+v02x=(1/m0 ) rAcm U −rAcI x¨+ r =U−x˙ = (U−x)Ds cd , 2

(1)

where a dot denotes differentiation with respect to time t, x is the displacement of the structure in line with the flow, mo is the mass per unit length of the structure, and j0 and v0 are the damping factor and the natural frequency of the structure (in vacuum), respectively. The right-hand side of equation (1) represents Morison’s hydrodynamic force components (inertia and drag), where r is the fluid density, A is the area of the cross-section of the structure, U is the flow velocity (which is a random function of time), cm is the inertia coefficient, cI is the added mass coefficient, cd is the drag coefficients, and Ds is the width of the structure. The non-linear drag force acts in the direction of the flow U−x˙ , relative to the structure. The physical origin of the non-linear hydrodynamic drag term is due to the sum of the viscous drag and pressure drag produced by the relative velocity between the structure and the flow. Introducing the non-dimensional parameters t=Vt,

q=x/Ds ,

U/VDs ,

p=dq/dt,

V 2=m0 v02/(m0+rAcI ), equation (1) can be written in the form dq d2q +a +q=g =u0−p=(u0−p)+bp0 (t), dt 2 dt

p0 (t)0du0 /dt,

(2)

where a=2j0zm0 /(m0+rAcI ),

b=rAcm /(m0+rAcI ),

g=DcD r/2v0 zm0 (m0+rAcI ) The wave velocity function u0 (t) is assumed to be a non-Gaussian stationary zero-mean random processes which may be described as the output of the following second order non-linear shaping filter: du dP(u0 ) dB(t) d2u0 +z 0+ =zs , dt 2 dt du0 dt

(3)

.   . . 

110

where B(t) is a standard Brownian motion process, z and s are some positive parameters (damping and intensity) of the filter, and P(u0 ) is some non-linear smooth function of u0 . The probability density function of the response of equation (3) can be determined from the stationary solution of the Fokker–Planck equation of equation (3) and is given by the well-known Gibbs distribution f(u0 , p0 )=

C z2pD

6 $

1 p2 P(u0 )+ 0 D 2

exp −

%7

,

D=

s , 2z

(4)

where the normalization constant is defined as

>g

a

exp[−D−1P(u0 )] du0

C=1

(5)

−a

provided that the integral in equation (5) exists. It follows from equation (4) that we have normal distribution with respect to p0 (with zero expectation and covariance D) and non-normal distribution with respect to u0 . The non-linear function P(u0 ) will be selected in a way that the spectrum and probability distribution of u0 are similar to those measured for ocean wave processes. Note that system (3) has a clear mechanical meaning. It describes the behavior of a mechanical system with potential P(u0 ) under white noise stochastic excitation. The function P(u0 ) may be selected as a fourth-order polynomial, P(u0 )=a1 (u0+a2 )2[(u0−a3 )2+a4 ],

(6a)

where ai are some positive parameters. Most of the measured probability density functions of ocean waves exhibit a single peak. In order to guarantee the existence of one peak, the function P(u0 ) should have one well, and this is guaranteed if the following condition is satisfied: (a2+a3 )2Q8a4 , In this case the probability density function f(u0 )=C exp{−D−1a1 (u0+a2 )2[(u0−a3 )2+a4 ]}

(6b)

has only one peak. The function f(u0 ) is plotted in Figure 1 for different values of a2+a3 and a4 with D−1a1=1. The parameter a2 is chosen so that

g

a

u0 f(u0 ) du0=0,

−a

i.e., the stationary process u0 (t) has zero mean. Thus by choosing the appropriate values of a2 , a3 and a4 , we obtain probability density functions of u0 similar to those measured experimentally. Experimental measurements [22–24] are processed for waves of significant height. It was shown that the measured probability density functions of severe waves are essentially non-Gaussian with non-zero mean. Based on these measured wave probability density functions, one can select appropriate values of a2 , a3 and a4 such that relation (6b) adequately describes a realistic distribution. It is seen that some probability density plots in Figures 1(a) and 1(b) are very similar to those measured in references [22–24]. One should recall that the analysis given in reference [21] considered only linear Gaussian waves with zero-mean. The linear filters yield Gaussian waves but can be adjusted to produce power spectra identical to those measured experimentally, such as the Pierson–Moskowitz spectrum.

 

111

Figure 1. Dependence of the filter response probability density function on (a) the parameters a2+a3 , for a4=0·5 and a1 /D=1; (b) the parameters a4 for a2+a3=1 and a1 /D=1.

Two identical probability density plots do not guarantee the equivalence of their power spectral density functions. Accordingly one should take care to select the parameters ai which guarantee the equivalence of both probability and spectra density functions. The

.   . . 

112

estimation of the power spectral density of u0 (t) will be carried out using an approximate approach in the next section. 3. FILTER OUTPUT STATISTICS

The filter’s equation (3) may be written in the state vector form: dp0=−(s0+s1 u0+s2 u02+s3 u03+jp0 ) dt+zs dB,

du0=p0 dt,

(7)

where s0 , s1 , s2 , s3 are constants defined by a1 , a2 , a3 , a4 : s0=2a1 a2 (a32+a4−a2 a3 ),

s1=2a1 (a22+a32+a4−4a2 a3 ),

s2=6a1 (a2−a3 ),

s3=4a1 .

In order to estimate the power spectral density function (which is the Fourier transform of the autocorrelation function) of the filter output, it is useful to determine the second-order joint probability density function f2 (u01 , p01 , u02 , p02 ; t1 , t2 ) of the process [u0 (t), p0 (t)] at t1 and t2 . f2 satisfies the partial differential equation 1 1 s 1 2f2 1f2 2 2 + ( p02 f2 )− [(s0+s1 u02+s2 u02 +s3 u02 +jp02 ]− 2 =0, 1t2 1u02 1p02 2 1p02

(8)

subject to the initial and normalized conditions, respectively: f2 (u01 , p01 , u02 , p02 ; t1 , t1 )=f1 (u01 , p01 ; t1 )d(u02−q01 )d( p02−p01 ),

gg

a

f2 (u01 , p01 , u02 , p02 ; t1 , t2 ) du02 dp02=f1 (u01 , p01 ; t1 ),

(9) (10)

−a

where f1 (u01 , p01 ; t1 ) is the first-order joint probability density function, which also satisfies the normalization condition and equation (8) when index 2 is replaced by 1. Unfortunately, an exact solution of equation (8) is not available. Instead, we will follow an approximate method by Pugachev and Sinitsyn [25] based on orthogonal expansion of the unknown probability density function. Further, we will use a modified form of orthogonal expansion method described by Moshchuk [26]. The method is based on the following representation of the unknown function f2 : f2 (u01 , p01 , u02 , p02 ; t1 , t2 )}

2 2 2 2 exp(−(u01 +u02 )/2g1−( p01 +p02 )/2g2 ) 2 (2p) g1 g2

a

× s cj 1 j 2 j 3 j 4 (t1 , t2 )Hj 1 (u01 , g1 )Hj 2 ( p01 , g1 )Hj 3 (u02 , g2 )Hj 4 ( p02 , g2 ),

(11)

j=0

where Hj are the Hermite polynomials defined by the relationship Hj ( y, g)=(−1)j

X 01

0 1

gj y2 d j y2 exp exp − ; j! 2g dy j 2g

gj are some positive parameters and cj 1 j 2 j 3 j 4 (t1 , t2 )=E[Hj 1 (u01 (t1 ), g1 )Hj 2 ( p01 (t1 ), g1 )Hj 3 (u02 (t2 ), g2 )Hj 4 ( p02 (t2 ), g2 )], The coefficients cj 1j 2j 3j 4 (t1 , t2 ) represent second-order statistical parameters of rank j1+j2+j3+j4 and are equivalent to a linear combination of second-order moments.

 

113

Substituting equation (11) into equation (8) gives the following first-order system of linear differential equations (t2et1 ):

X

1cj 1 , j 2 , j 3 , j 4 g = j3 2 (zj4+1cj 1 , j 2 , j 3−1, j 4+1+zj4cj 1 , j 2 , j 3−1, j 4−1 )−s0 1t2 g1 −s1

X

−s2 g1

X

j4 c g2 j 1 , j 2 , j 3 , j 4−1

g1 (zj3+1cj 1 , j 2 , j 3+1, j 4−1+zj3cj 1 , j 2 , j 3−1, j 4−1 ) g2

j4

X

j4 [z( j3+1)( j3+2)cj 1 , j 2 , j 3+2, j 4−1+(2j3+1)cj 1 , j 2 , j 3 , j 4−1 g2

+zj3 ( j3−1)cj 1 , j 2 , j 3−2, j 4−1 ]−s3 g1

X

j4

g1 [z( j3+1)( j3+2)j3+3)cj 1 , j 2 , j 3+3, j 4−1 g2

+3( j3+1)z( j3+1)cj 1 , j 2 , j 3+1, j 4−1+3j3 zj3cj 1 , j 2 , j 3−1, j 4−1 +zj3 ( j3−1)( j3−2)cj 1 , j 2 , j 3−3, j 4−1 ]+

0

1

s −j zj4 ( j4−1)cj 1 , j 2 , j 3 , j 4−2 2g2

−jj4 cj 1 , j 2 , j 3 , j 4 ;

(12)

with the initial conditions cj 1 , j 2 , j 3 , j 4 (t1 , t1 )=

gg

a

f1 (u01 , p01 ; t1 )Hj 1 (u01 )Hj 2 ( p01 )Hj 3 (u01 )Hj 4 ( p01 ) du01 dp01 ,

(13)

−a

which follow from equation (9). Equation (12) constitutes a set of an infinite hierarchy of coupled equations. Since we are interested in estimating the autocorrelation function of u0 (t), we can decompose equation (12) into independent subsystems by fixing j1=1 and j2=0. It should be noted that 1 1 c1010= E[u0 (t1 )u0 (t2 )]= Ru 0 (t1 , t2 ), g1 g1

(14)

where Ru 0 (t1 , t2 ) is the correlation function of the process u0 (t). Thus, for c1,0, j 3 , j 4 we have an infinite set of equations generated from equation (12). For the steady state case, f1 does not depend on t1 (and is given by relation (4)) and f2 depends on t=t2−t1 . It is convenient now to set g1=Est [u02], and g2=Est [p02]=s/(2j), where Est means stationary expected value. In this case the initial conditions given by equation (13) take the form

6

c1,0, j 3 ,0 (t1 , t1 )=C zj3+1

+zj3

g

a

Hj 3+1 (x, g1 ) exp[−D−1P(x)] dx

−a

g

a

7

Hj 3−1 (x, g1 ) exp[−D−1P(x) dx ;

−a

c1,0, j 3 , j 4 (t1 , t1 )=0

for j4$0.

(15)

114

.   . . 

Generally it is not possible to find an exact solution of the infinite hierarchy of equations (12) for c1,0, j 3 , j 4 . We can therefore confine our attention to the expansion coefficients c1,0, j 3 , j 4 up to and including some rank N of the subset j3 and j4 i.e. 0Ej3+j4EN (coefficients of a higher rank are ignored). Solving this truncated system with initial conditions (15), we obtain an approximation for the autocorrelation function Ru 0 (t) by numerical integration. The spectral density function Su 0 (v) is the Fourier transform of Ru 0 (t) which is estimated numerically at discrete values of v: Su 0 (v)=

1 p

g

a

Ru 0 (t) cos vt dt.

(16)

0

In the present study, computations showed that N=10 provides satisfactory convergence, and it was found that additional equations (N=12 and 14) do not significantly improve the results. For N=10, 65 equations of the statistical parameters c1,0, j 3 , j 4 are generated by developing a computer program that changes the indices of the j3 and j4 . The numerical integration is performed by using the fourth-order Runge–Kutta method with integration step Dt=0·001. The estimated value c1010 is then used in relation (14) to evaluate the autocorrelation function Ru 0 (t). The Fourier transform (16) is performed by fixing v, and numerical integration is carried out within a time region starting from zero up to a value which yields convergence. The value of v is changed in steps of 0·01 starting from zero up to 3. Figures 2 and 3 show power spectra Su 0 (v) of the non-linear filter (3) for different values of a1 and j, respectively. It is seen that when a1 increases the peak of the spectrum is reduced and moves to the right. At the same time the width

Figure 2. Dependence of the filter power spectra on a1 for a2+a3=1, a4=2, j=0·1 and s=0·2.

 

115

Figure 3. Dependence of the filter power spectra on j for a1=0·02, a2+a3=1, a4=2 and s=0·2.

of the spectrum remains almost the same. When the damping coefficient j increases, the effective bandwidth of the filter increases as in the linear case. Figure 4 shows Su 0 (v) for different values of the white noise intensity s. In contrast to the linear filter, it is seen that the peak of the non-linear spectrum moves to the right when s increases. The reason for this shift is that the natural frequency of the filter is amplitude-dependent, a fact which is well known in deterministic non-linear vibration theory. Careful inspection of the filter spectra, and comparison with the linear filter spectra, reveal some asymmetry which is not common for linear oscillators. The mean square response of the filter is shown in Figure 5 as a function of the white noise intensity (curve 1). For three selected points in Figure 5, the center frequency of the non-linear filter is estimated and the corresponding linear response is plotted by three straight lines. Both linear and non-linear responses intersect at points which correspond almost exactly to the center frequency of the non-linear filter. It is interesting to note that each point on the non-linear response curve can be predicted using the linear relationship E[u02]=s/(2vc2+j 2 )j, where vc is the center frequency of the filter for a particular excitation level, and corresponds to the peak of the power spectrum density. Figure 5 also shows another mean square response (curve 2) for an auxiliary filter whose linear and non-linear parts do not vary as the excitation level changes. For relatively low excitation levels, both the linear and auxiliary non-linear filters (for filter natural frequency vf=zs0=0·8) have almost the same mean square responses. Note that vc and vf are not equal. As the excitation level increases it is seen that the non-linearity affects the filter response.

.   . . 

116

Figure 4. Dependence of the filter power spectra on s for a1=0·02, a2+a3=1, a4=2 and j=0·2.

4. SYSTEM RESPONSE STATISTICS

The system response to the excitation generated from the non-linear filter outlined in the previous section is studied by using the stochastic averaging method of Stratonovich [27] and Khasminskii [28]. The parameters a, g and b 2 of equations (2) are assumed small and of the same order e. Introducing the polar co-ordinates I and 8 defined by the transformation q=I cos 8,

dq/dt=−I sin 8,

8=c+t,

(17)

equations (2) can be rewritten in the standard form of stochastic averaging: dI/dt=−aI sin2 8−g sin 8=u0+I sin 8=(u0+I sin 8)−b sin 8p0 (t), dc/dt=−a sin 8 cos 8−gI−1 cos 8=u0+I sin 8=(u0+I sin 8)−I−1b cos 8p0 (t).

(18)

Note that the state co-ordinates I and c are slowly varying with time t, while u0 and p0 exhibit fast variation. According to Khasminskii’s limit theorem [28], it follows that the solution of equations (2) converges weakly to a diffusion process described by the following Ito stochastic differential equations:

7 X

6

1 c dI= − aI+ −gU(I) dt+ 2 4I

c 1 dB (t), dc= 2 1 I

X

c dB (t), 2 2

(19)

where U(I)=

1 2p

g g 2p

a

d8

0

−a

=u0+I sin 8=(u0+I sin 8)f(u0 ) du0,

c=b 2v,

v=2pSp 0 (1),

(20)

 

117

Figure 5. Mean square response of linear and non-linear filters as function of excitation level: 1, non-linear filter with coefficient which vary with s such that the mean value E[u0 ]=0; 2, nonlinear filter with constant coefficients, where E[u0 ]$0; linear filter responses with different frequency vf=0·8, 0·97, and 1·08.

and B1 (t) and B2 (t) are two independent Brownian motion processes, each of unit intensity. As g}e, the averaging of the non-linear drag term (term with g) in equations (18) is taken with respect to 8 and with respect to the stationary distribution of the ergodic process u0 (t), keeping I and c frozen. In equation (20), Sp 0 (v) is the power spectral density of the process p0 (t) and f(u0 ) is the stationary probability density function of u0 . As Sp 0 (v)=v 2Su 0 (v), it follows that Sp 0 (1)=Su 0 (1). It is not difficult to show that [21] U(I)=I 2[=u0 =+h(u0 /I)],

W=

g

a

(W)f(u0 ) du0 ,

(21)

−a

where h(y)=(2/3p)[3z1−y2+3y sin−1 y−(1−y 2 )3/2]−ye0, h(y)=h(−y),

for 1eye0;

and h(y)=0 for =y =q1.

Introducing the Hamiltonian function H=I 2/2, the first equation (19) takes a form which includes the Ito correction term after the transformation,

$

%

c dH= −aH+ −8gH 2n(H) dt+zcH dB1 , 2

(22)

where a=a+2g=u0 =,

n(H)=

g

1

h(y)f(yz2H) dyq0.

0

Note that the integration limits are changed from 2a to 21 (or from 0 to +1 due to symmetry) because h(y) is non-zero only over the domain −1QyQ+1, and both h(y) and f(yz2H) are even functions with respect to y.

.   . . 

118

The phase angle c(t) only constitutes a diffusion component as the drift coefficient vanishes. It follows from equations (19) and (21) that the non-linear term g =u0−p =(u0−p) in the expression of the drag force on average has a dissipative nature (Uq0). It is difficult to show that the solution of the stationary FPK equation of system (22) is given by the expression

$ $

exp[−2(a/c)H] exp −(16g/c) f(H)=

g6 a

exp[−2(a/c)H] exp −(16g/c)

0

g g

H

% %7

hn(h) dh

0 H

0

hn(h) dh

.

(23)

dH

The following inequality for the response stationary mean value may be established based on equation (22): (E[H])stE

c pb 2Su 0 (1) , = 2a a+2g=u0 =

(24)

where the right-hand side does not include the contribution of h(H). It is seen that the response statistics of the system depend on the value of the filter spectral density Su 0 at v=1 (that is, the natural frequency of the system) and on the stationary distribution of u0 . The dependence of the system response mean value E[H] on the white noise intensity is plotted in Figure 6, for system parameters a/g=1, b 2/g=4 and for filter parameters a2+a3=1, a4=2, and a1=0·002. The value of a2 is not fixed in this figure. However, it is selected such that the mean value of the filter response is always zero. Table 1 shows the values of a2 corresponding to white noise intensity s. The corresponding values of the filter spectra and =u0 = are also included in the table. The dependence of E[H] on the non-linear hydrodynamic drag coefficient g is shown in Figure 7. It is seen that the non-linear drag force suppresses the system response. The probability density of the response energy H is estimated numerically from equation (23) and is plotted in Figure 8 for three different values of the non-linear drag parameter

Figure 6. Dependence of system mean energy response on the white noise excitation level for filter parameters a1=0·02, a2+a3=1, a4=2 and j=0·1, and system parameters a/g=1 and b 2/g=4.

 

119

T 1 Statics of the non-linear filter s

a2

Sq 0 (1)

=q0 =

0·08 0·1 0·2 0·3 0·4 0·6 0·8 1·0 1·2 1·5 2·0 2·5 4·0

0·368 0·380 0·412 0·427 0·436 0·447 0·454 0·459 0·462 0·466 0·471 0·473 0·479

0·037 0·056 0·251 0·553 1·129 2·827 3·782 4·733 5·066 5·0 4·607 4·322 3·740

0·917 0·976 1·181 1·319 1·423 1·584 1·708 1·810 1·897 2·010 2·164 2·291 2·575

g/a=0, 0·1 and 0·5. It is clear that the non-linear drag results in non-normality of the system response. The corresponding normal density shown by dotted curves is obtained by using the well-known Gibbs function f(H)=

0

1

1 H exp − , E[H] E[H]

(25)

where the mean value E[H] is taken from the corresponding non-linear estimate. 5. ORTHOGONAL EXPANSION OF THE SYSTEM-FILTER EQUATIONS

System response statistics obtained by applying stochastic averaging are accurate when a, g and b 2 are of order e when e is sufficiently small. If these assumptions are not valid, the results obained by the averaging method are no longer accurate. Alternatively, the response of system (2) can be analyzed by the normal approximation method, moments

Figure 7. Dependence of system mean energy response on the non-linear hydrodynamic drag coefficient g for filter parameters a1=0·02, a2+a3=1, a4=2, s=0·6 and j=0·1, and system parameters a=1, and b=2.

.   . . 

120

Figure 8. Probability density function of the system response energy level for different system parameters and for filter parameters a1=0·02, a2=0·447 a2+a3=1, a4=2, s=0·6 and j=0·1: ——, non-linear drag; – – –, normal density.

method, cumulants method, and any other method of orthogonal expansion. In this section, one version of the orthogonal expansion technique [24] will be adopted. The method is based on orthogonal expansions of the unknown probability density function in a suitable functional space and is similar to that used in section 3. Introducing a new variable p*=dq/dt−u0 , the system and filter equations can be rewritten as a set of the Ito stochastic differential equations: dq=( p*+u0 ) dt, du0=p0 dt,

dp*=[−q−a( p*+u0 )−gp*=p*=+(b−1)p0 ] dt, dp0=−(s0+s1 u0+s2 u02+s3 s03+jp0 ) dt+zs dB.

(26)

The joint probability density function f1 (q, p*, u0 , p0 , t) of the process [q(t), p*(t), u0 (t), p0 (t)] satisfies the FPK equation and can be written in terms of orthogonal expansion similar to the one given by relation (11): a

f1 (q, p*, u0 , p0 , t)}m s cj 1 j 2 j 3 j 4 (t)Hj 1 (q, g1 )Hj 2 ( p*, g2 )Hj 3 (u0 , g3 )Hj 4 ( p0 , g4 ), j=0

m=

exp(−q 2/2g1−p*2/2g2−u02/2g3−p*2/2g4 ) , (2p)2zg1 g2 g3 g4

(27)

where Hj are Hermite polynomials defined by equation (11a), gj are some positive parameters, and cj 1 j 2 j 3 j 4 (t)=E[Hj 1 (q(t), g1 )Hj 2 ( p*(t), g2 )Hj 3 (u0 (t), g3 )Hj 4 ( p0 (t), g4 )].

(28)

 

121

Series (27) converges if and only if

g

a

f 12m−1 dq dp* du0 dp0Qa.

−a

Coefficients cj 1 j 2 j 3 j 4 (t) satisfy a set of infinite hierarchy of linear differential equations which can be obtained from equation (28) by using the Ito formula of complex functions:

$ $0 0

dcj 1 j 2 j 3 j 4 (t) d =E (H H H H ) dt dt j 1 j 2 j 3 j 4 =E

×

%

1 1

0 0

1 1

1Hj 1 dq 1Hj 2 dp* H H H +Hj 1 Hj 3 Hj 4+Hj 1 Hj 2 1q dt j 2 j 3 j 4 1p* dt

0 1%

1Hj 3 du0 1Hj 4 dp0 s 1Hj24 Hj 4+Hj 1 Hj 2 Hj 3 + H j 1 Hj 2 H j 3 2 1u0 dt 1p0 dt 1p02

.

(29)

The time derivatives dq/dt, dp*/dt, du0 /dt, dp0 /dt must be substituted according to equation (26). Using the following properties of the Hermite polynomials:

X

1Hj ( y, gi ) j = H ( y, gi ), 1y gi j−1

yHj ( y, gi )=z( j+1)giHj+1 ( y, gi )+zjgiHj−1 ( y, gi ), (30)

it is not difficult to express the right-hand side of equation (29) in terms of the coefficients cj 1 j 2 j 3 j 4 (t). There is only one hydrodynamic term gp*=p*= which requires special attention. A detailed treatment of the component of dcj 1 j 2 j 3 j 4 (t)/dt is given:

$

E p*=p*=

%

1 (H H H H ) =zj2 /g2E[Hj 1 (Hj 2−1 p*=p* =)Hj 3 Hj 4 1p* j 1 j 2 j 3 j 4 =zj2 /g2

g

+a

f1 (q, p*, u0 , p0 ,t)Hj 1 (Hj 2−1 p*=p*=)

−a

×Hj 3 Hj 4 dq dp* du0 dp0 =zj2 /g2

g

a

+a

m

−a

s

ck 1 k 2 k 3 k 4 Hk 1 Hk 2 Hk 3 Hk 4 H j 1

k 1 ,k 2 ,k 3 ,k 4=0

×(Hj 2−1 p*=p*=)Hj 3 Hj 4 ] dq dp* du0 dp0 a

=zj2 /g2 s cj 1 k 2 j 3 j 4 k 2=0

g

+a

−a

exp(−p*2/2g2 ) z2pg2 a

×Hk 2 Hj 2−1 p*=p*= dp*=zj2 /g2 s cj 1 k2 j3 j4 k 2=0

×

6g

+a

0

exp(−p*2/2g2 ) z2pg2

Hk 2 [zj2 ( j2+1)Hj 2+1+(2j2−1)

×Hj 2−1+z( j2−1)( j2−2)Hj 2−3 ] dp*

.   . . 

122

−

g

0

exp(−p*2/2g2 )

−a

z2pg2

Hk 2 [zj2 ( j2+1)Hj 2+1

+(2j2−1)Hj 2−1+z( j2−1)( j2−2)Hj 2−3 ] dp*

=

7

X

j2 a s c k j j [zj2 ( j2+1)ak 2 , j 2+1 g2 k 2=0 j 1 2 3 4

+(2j2−1)ak 2 j 2−1+z( j2−1)( j2−2)ak 2 , j 2−3 ], where aij0

g

a

0

exp(−p*2/2g2 ) z2pg2

Hi ( p*)Hj ( p*) dp*−

g

0

exp(−p*2/2g2 )

−a

z2pg2

Hi ( p*)Hj ( p*) dp*

and {gi }={E[q 2 ], E[p*2 ], E[u02], E[p02]}T; T denotes the transpose. The following recursive formula can be obtained for aij : aij=

X

i H (0, g2 )Hj−1 (0, g2 ) a +2 i , j i−1, j−1 z2pj

which is useful for numerical analysis. Finally, the equation for cj 1 j 2 j 3 j 4 take the following form:

X X X

g2 (zj2+1cj 1−1, j 2+1, j 3 , j 4+zj2cj 1−1, j 2−1, j 3 j 4 ) g1

cj 1 j 2 j 3 j 4 /dt= j1

+ j1

g3 (zj3+1cj 1−1, j 2 , j 3+1, j 4+zj3cj 1−1, j 2 , j 3−1, j 4 ) g1

− j2

g1 (zj1+1cj 1+1, j 2−1, j 3 j 4+zj1cj 1−1, j 2−1, j 3 , j 4 ) g2

−a( j2 cj 1 , j 2 , j 3 , j 4+zj2 ( j2−1)cj 1 , j 2−2, j 3 , j 4 )

X

+(b−1) j2

X 0 1 X X

+ j3

+

g4 (zj4+1cj 1 , j 2−1, j 3 , j 4+1+zj4cj 1 , j 2−1, j 3 j 4−1 ) g2

g4 (zj4+1cj 1 , j 2 , j 3−1, j 4+1+zj4cj 1 , j 2 , j 3−1, j 4−1 ) g3

s −j zj4 ( j4−1)cj 1 , j 2 , j 3 , j 4−2−jj4 cj 1 , j 2 , j 3 , j 4 2g4

−a j2

−s0

g3 (zj3+1cj 1 , j 2−1, j 3+1, j 4+zj3cj 1 , j 2−1, j 3−1, j 4 ) g2

j4 c −s g4 j 1 , j 2 , j 3 , j 4−1 1

X

j4

g3 (zj3+1cj 1 , j 2 , j 3+1, j 4−1 g4

  +zj3 cj 1 , j 2 , j 3−1, j 4−1 )−s2 g3

123

X

j4 [z( j3+1)( j3+2)cj 1 , j 2 , j 3+2, j 4−1 g4

+(2j3+1)cj 1 , j 2 , j 3 , j 4−1+zj3 ( j3−1)cj 1 , j 2 , j 3−2, j 4−1 ] −s3 g3

X

j4

g3 [zj3+1)( j3+2)( j3+3)cj 1 , j 2 , j 3+3, j 4−1 g4

+3( j3+1)zj3+1cj 1 , j 2 , j 3+1, j 4−1+3j3 zj3cj 1 , j 2 , j 3−1, j 4−1 a

+zj3 ( j3−1)( j3−2)cj 1 , j 2 , j 3−3, j 4−1 ]−gzj2 g2 s cj 1 ,k, j 3 , j 4 [zj2 ( j2+1)ak, j 2+1 k=0

+(2j2−1)ak, j 2+z( j2−1)( j2−2)ak, j 2−3 ].

(31)

All statistical moments of the system response may be expressed in terms of cj 1 j 2 j 3 j 4 . For example, expressions for E[q 2 ], E[p*2 ], E[dq/dt)2 ] and E[H] are E[q 2 ]=g1 (z2c2000+1), E

$0 1 % dq dt

E[p*2 ]=g2 (z2c0200+1),

2

=g2 (z2c0200+1)+g3 (z2c0020+1)+2zg2 g3 c0110 ,

1 E[H]= [g1 (z2c2000+1)+g2 (z2c0200+1)+g3 (z2c0020+1)+2zg2 g3 c0110 ]. 2

(32)

The probability density function of the response displacement q is a

f(q, t)= s cj 1 000 (t)Hj 1 (q).

(33)

j 1=1

It is not possible to obtain an exact solution of the infinity hierarchy of equations (31) for the coefficients cj 1 j 2 j 3 j 4 . In this case attention is focused only on coefficients up to some rank N of the subset j1 , j2 , j3 , j4 , i.e. 0Ej1+j2+j3+j4EN. Solving this truncated system with appropriate initial conditions, we obtain approximate values for the statistical parameters cj 1 j 2 j 3 j 4 and the probability density function f1 . The order of approximation N can easily be changed in the numerical algorithm. In the present case, computations showed that N=10 provides satisfactory convergence, and that additional equations (N=12 and 14) do not significantly improve the results. For N=10, 1000 equations of the statistical parameters cj 1 j 2 j 3 j 4 are generated by developing a computer program that changes the indices of j1 , j2 , j3 and j4 . Numerical integration is performed by using the fourth-order Runge–Kutta scheme with integration step Dt=0·001. The estimated value of cj 1 j 2 j 3 j 4 is then used in relations (32) and (33) to evaluate the statistical parameters E[q 2 ], E[p*2 ], E[dq/dt)2 ], E[H] and the probability density function of the response displacement q. The dependence of the system response mean value E[H] on the non-linear hydrodynamic drag coefficient g is shown in Figure 9 by a solid line for the following system and filter parameters: a=0·04, b=0·4, a1=0·02, a2+a3=1, a2=0·447, a4=2, s=0·6, j=0·1. The corresponding dependence obtained from the averaging method is plotted by a dotted line. It is seen that the two curves are very close for small values of g and then deviate as g increases. The mean value of the system energy response as determined from the averaging method is higher than the one obtained from orthogonal expansion. The dependence of different system statistical parameters on the non-linear hydrodynamic drag coefficient g is shown in Figure 10. The probability density of the

124

.   . . 

Figure 9. Dependence of the system mean energy response on the non-linear hydrodynamic drag coefficient g for filter parameters a1=0·02, a2=0·447, a2+a3=1, a4=2, s=0·6 and j=0·1, and system parameters a=0·04, b=0·4: ——, orthogonal expansion method; – – –, averaging method.

response co-ordinate q is plotted in Figure 11 for three values of the hydrodynamic coefficient g=0·0, 0·1 and 1·0. The response corresponding to zero hydrodynamic drag is included and reveals the influence of the wave non-normality on the linear oscillator response. The system and filter parameters used in this figure are a=1·0, b=2·0, a1=0·02, a2+a3=1, a2=0·447, a4=2, s=0·6 and j=0·1. The corresponding normal densities are shown by dotted curves. Covariances for the normal distribution are taken from the corresponding non-linear estimate. It is seen that a significant deviation from Gaussianity takes place when the non-linear drag coefficient increases. Furthermore, the hydrodynamic non-linear drag reduces the extreme values of the response, as indicated by

Figure 10. Dependence of different system and filter statistical parameters on the non-linear hydrodynamic drag coefficient g for filter parameters a1=0·02, a2=0·447, a2+a3=1, a4=2, s=0·6 and j=0·1, and system parameters a=1, b=2.

 

125

Figure 11. Probability density function of the response co-ordinate q for system and filter parameters a=1, b=2, a1=0·02, a2=0·447, a2+a3=1, a4=2, s=0·6 and j=0·1: ——, probability density function for two values of g=0·0, 0·1 and 1·0; – – – –, corresponding Gaussian probability density function.

the tails of the probability density plots. The wave non-linearity alone (in the absence of hydrodynamic drag) does not result in a significant deviation from the Gaussian distribution. However, it results in an associated increase of the extreme response amplitude which is slightly higher than the one given by the Gaussian plot. In order to verify the accuracy of the approximate scheme outlined in section 3, we consider the first-order probability density function f1 (u0 , p0 , t). In this case the unknown function f1 is written in the form f1 (u0 , p0 , t)}

exp(−u02/2g1−p02/2g2 ) 2pzg1 g2

a

s cj 3 j 4 (t)Hj 3 (u0 , g1 )Hj 4 ( p0 , g2 ),

(34)

j 3 , j 4=0

where cj 3 j 4 (t)=E[Hj 3 (u0 , g1 )Hj 4 ( p0 , g2 )]. Coefficients cj 3 , j 4 (t) represent first-order statistical parameters of rank j3+j4 (they are equivalent to a linear combination of first-order moments) and satisfy and infinite hierarchy given by equations (12), where indices j1 and j2 are dropped. For the stationary case, after truncation we obtain a linear algebraic system of equations which can easily be solved numerically. The exact probability density function is estimated using relation (6b). This exact solution is compared with the approximate truncated solution for N=8 and the two results are shown in Figure 12. It is clear that the approximate solution (dotted line) is very close to the exact one (solid line). 6. CONCLUSIONS

Non-Gaussian random ocean waves are represented by the response of a shaping non-linear filter to a Gaussian white noise random process. The filter non-linearity is selected such that its probabilistic behavior emulates the probabilistic measurements of ocean waves. In one case the filter response statistics are estimated by using the method of orthogonal expansion of multi-dimensional random processes, while the system response statistics are estimated by using the method of stochastic averaging.

126

.   . . 

Figure 12. Probability density functions of the filter response for parameters a1=0·02, a2=0·447, a2+a3=1, a4=2, s=0·1 and j=0·1; ——, exact; – – – –, approximate.

Another case considered the method of orthogonal expansion for both the filter and system equations. The probability density function of the response is described by the Fokker–Planck equation at two different time intervals whose difference measures the correlation time for stationary processes. A general first-order differential equation describing the response statistics is derived and is found to form an infinite hierarchy. 65 equations are generated from this equation which include statistics up to order 10, while all higher order statistics are ignored. The filter response spectra are then estimated for different filter and excitation parameters. It is found that the filter non-linear and excitation level have strong effects on the spectra peaks and on their location with respect to the frequency axis. The system response to such non-Gaussian ocean waves is then analyzed using the stochastic averaging method. It is shown that the response probability density becomes significantly non-normal as the non-linear drag force increases. The results obtained by the method of orthogonal expansion for the combination of the system and filter equations are less conservative than those given by stochastic averaging. The deviation of the results of the two methods is more pronounced when system non-linearity increases. The proposed non-linear filter described by equation (7) gives more reliable results than those modelled using the Gram–Charlier series or Edgworth series distributions. Ochi [24] has indicated that these distribution models can lead to negative probability density at the tails. The proposed filter preserves the statistical properties of non-Gaussian waves and approximates the distribution to those measured experimentally by properly selecting the coefficients of the filter potential function P(u0 ).

 

127

The analysis presented in this paper has been applied to single-degree-of-freedom oscillators. The authors are currently investigating more practical problems modelled by more degrees-of-freedom, such as immersed cables in the ocean. This type of problem involves non-linear modal interaction due to internal resonance conditions among the interesting modes. The accuracy of the results is also verified by using Monte Carlo simulation. ACKNOWLEDGMENT

This research is supported by a grant from ONR, No. N000149310936. Dr. Julia Abrahams is the Scientific Officer. REFERENCES 1. R. H. L, M. H and C. B. H 1961 The Journal of the Acoustical Society of America 33, 1404–1411. Response of hard-spring oscillator to narrow band excitation. 2. W. C. L and Y. C. K 1976 ASME Journal of Applied Mechanics 43, 340–344. Narrow band excitation of a nonlinear oscillator. 3. T. F and E. H. D 1987 International Journal of Non-linear Mechanics 22(3), 267–274. Numerical simulations of jump phenomena in stable Duffing systems. 4. H. G. D and S. R 1988 Journal of Sound and Vibration 126(2), 195–208. Random superharmonic and subharmonic response: multiple time scaling of a Duffing oscillator. 5. H. G. D and Q. L 1990 Journal of Sound and Vibration 139, 1–8. The response envelop probability density function of a Duffing oscillator with random narrow band excitation. 6. H. G. D and Q. L 1992 International Journal of Non-linear Mechanics 27(5), 805–816. On the narrow band random response distribution function of a non-linear oscillator. 7. J. B. R and P-T. D. S 1986 International Journal of Non-linear Mechanics 21(2), 111–134. Stochastic averaging: an approximate method of solving random vibration problems. 8. M. F. D 1971 Mechanics of Solids 6, 142–146. Oscillations of a system with nonlinear cubic characteristics under narrow band random excitation. 9. R. N. I 1989 Structural Safety 6, 177–185. Response of nonlinear systems to narrow band excitation. 10. J. B. R and P.-T. D. S 1990 Random Vibrations and Statistical Linearization. New York: John Wiley. 11. J. B. R 1991, International Journal of Non-linear Mechanics 26, 945–959. Multiple solutions generated by statistical linearization and their physical significance. 12. R. N. I 1992 Nonlinear Stochastic Mechanics, IUTAM Symposium, Turin. (N. Bellomo and F. Casciati, Editors), Springer-Verlag. Approximate analysis of non-linear systems under narrow band random inputs. 13. B. F. S and L. A. B 1993 Nonlinear Dynamics 4(4), 357–372. On the numerical solution of the Fokker–Planck equation for nonlinear stochastic systems. 14. W. Q. Z, M. Q. L and Q. T. W 1993 Journal of Sound and Vibration 165(2), 285–304. Stochastic jump and bifurcation of a Duffing oscillator under narrow-band excitation. 15. H. G. D and D. N 1986 Journal of Sound and Vibration 104(2), 277–283. Phase plane for narrow band random excitation of a Duffing oscillator. 16. T. K 1985 Journal of Sound and Vibration 102(3), 440–441. Stochastic response with bifurcation on nonlinear Duffing’s oscillators. 17. P. K. K and S. R. B 1993 Nonlinear Dynamics 4, 279–288. Quasi-harmonic analysis of the behavior of hardening Duffing oscillator subjected to filtered white noise. 18. R. N. I 1986 Journal of Sound and Vibration 126, 255–263. Stochastic response and stability of the Duffing oscillator under narrow band excitation. 19. R. N. I 1988 International Journal of Non-linear Mechanics 23, 385–391. Higher order linearization in non-linear random vibration. 20. P. K. K and R. S. L 1993 International Journal of Non-linear Mechanics 28(2), 145–155. Improved stability analysis of the response of a Duffing oscillator under filtered white noise.

128

.   . . 

21. N. K. M, R. A. I and R. Z. K 1995 Journal of Sound and Vibration 184, 681–701. Response statistics of ocean structures hydrodynamic loading. Gaussian ocean waves. Part I. 22. E. M. B 1980 Journal of Applied Ocean Research 2(2), 63–74. Nonlinear effects on the statistical model of shallow-water wind waves. 23. M. K. O and W. C. W 1984 19th Conference of Coated Engineering. Non-Gaussian characteristics of coastal waves. 24. M. K. O 1986 Probabilistic Engineering Mechanics 1(1), 28–39. Non-Gaussian random processes in ocean engineering. 25. V. S. P and I. N. S 1987 Stochastic Differential Systems. New York: John Wiley. 26. N. K. M 1994 Automation and Remote Control 1, 72–87. An approximate method for determining finite-dimensional distributions in stochastic differential systems. 27. R. L. S 1963 Topics in the Theory of Random Noise, Vol. 1 . New York: Gordon & Breach. 28. R. Z. K 1966 Theory of Probability and its Applications 11(3), 444–462. A limit theorem for the solutions of differential equations with random right-hand sides.