- Email: [email protected]
Response statistics of ocean structures to non-linear hydrodynamic loading
Journal of Sound and Vibration (1995) 184(4), 681–701
RESPONSE STATISTICS OF OCEAN STRUCTURES TO NON-LINEAR HYDRODYNAMIC LOADING, PART I: GAUSSIAN OCEAN WAVES N. M R. A. I Department of Mechanical Engineering, Wayne State University, Detroit, Michigan 48202, U.S.A.
R. K Department of Mathematics, Wayne State University, Detroit, Michigan 48202, U.S.A. (Received 7 March 1994, and in final form 21 June 1994) This paper presents an analytical approach based on the stochastic averaging of the energy envelope to treat the dynamic behavior of single-degree-of-freedom elastic ocean structures. Such systems are usually subjected to a narrow-band random process which may be modelled as the output of a shaping filter. The response co-ordinates of the system and filter experience slow and fast variations with respect to time, respectively. For this class of problems, the method of stochastic averaging is used to establish stochastic Ito differential equations for the process of slow variation. This approach has not previously been used for non-linear mechanical problems. Three different shaping filters (including those possessing a Pierson–Moskowitz spectrum) are employed to model Gaussian random sea waves. The response statistics of the structure are estimated in terms of the excitation spectrum for different levels of non-linear hydrodynamic drag force. It is found that the hydrodynamic drag reduces the system response energy and consequently suppresses the motion of the structure. In addition, the response probability density deviates from normality as the non-linear hydrodynamic drag parameter increases. The case of non-Gaussian random sea waves will be considered in Part II. 7 1995 Academic Press Limited
1. INTRODUCTION
The non-linear dynamics of ocean structures exposed to severe wave conditions pose complex problems for engineers and mathematicians. The dynamic hydroelastic problem of ocean systems involves the interaction of inertia forces, elastic forces and hydrodynamic forces. Basically, the hydrodynamic forces induced on structural elements are non-linear and a combination of inertial and drag contributions [1]. Each force is time-dependent and depends on the geometrical properties of the structure, fluid properties describing the flow field and some ‘‘unknown constants’’ which are determined experimentally. Furthermore, under large deflection these structures experience geometrical non-linearities. The hydrodynamic forces due to sea waves are random, and under extreme environmental conditions are non-Gaussian. The basic problems of these types of systems include the response statistics of non-linear dynamical systems to non-Gaussian excitations, level crossings, local extremes of the response and associated reliability problems. Experimental measurements of the response of laboratory models under wave action have indicated a significant deviation from a Gaussian probability distribution, 681 0022–460X/95/290681 + 21 $12.00/0
7 1995 Academic Press Limited
682
. .
particularly at a high level of response [2]. This behavior is taken to be the result of the non-linearity inherent in the drag component of Morison’s equation, and leads to an increased probability of exceeding a given load or response when compared with a Gaussian distribution at high levels of response. The kurtosis of the basic process serves as a useful measure of the non-linearity. The results of Tickell [2] are based on a long-crested or uni-directional surface elevation whose horizontal and vertical velocities of flow are non-zero mean Gaussian variables. The dynamic behavior of steel-framed offshore structures, including stochastic analysis in the frequency domain and simulation in the time domain, was examined within the framework of the linear theory of small oscillations [3]. Approximate analytical methods such as the equivalent linearization technique [4, 5], perturbation analysis [6], stochastic averaging [6] and Gaussian and non-Gaussian closures [7] have been applied for structural elements described by Morison’s equation [1]. Digital simulation techniques have been used to determine the response statistics of offshore structural elements [8–10]. The main results of these studies indicate that the non-linear hydrodynamic damping has the effect of increasing the rate of decay of the response amplitude. At high wind velocities the response is quasi-static exhibiting a wide-band random process characteristic [8]. A comprehensive study [11] reviewed some recent studies of the non-linear collapse and fatigue analysis of jacket structures, dynamic analysis of jack-up structures, springing and slow-drift response of a tension leg platform, and hysteretic foundation behavior of a concrete gravity platform. The mathematical modelling of ocean systems depends on whether the system is treated as a rigid body or as a continuous elastic structure. The models discussed in the literature involve only a finite number of independent co-ordinates which are described by a set of ordinary differential equations. These include ship roll and coupled pitch roll motions, floating buoys, moored semi-submersibles, rigid gravity offshore platforms with co-ordinates describing sliding and rocking motions, and jacket template platforms with co-ordinates to describe the motion of discrete masses lumped at node points. However, for line components such as rather long beams, pipelines and cables, alternative continuous system models may provide more precise descriptions of component motion. Since a partial differential equation is used to characterize the motion of a continuous line component, a solution may be more involved mathematically than for a corresponding lumped system. For beams, the bending stiffness and longitudinal tension are incorporated in the model. Cables, on the other hand, resist tension but their bending and torsional stiffnesses are usually neglected. These structural elements are usually encountered in pipelines, steel ropes and chains used to stay buoys, floating platforms, compliant towers and ships. The inclusion of large deformation in the mathematical modelling will result in non-linear partial differential equations which include non-homogeneous random forcing terms and/or parametric excitation terms. The dynamic behavior of ocean systems has been treated by using one of four approaches mentioned earlier. These approaches have been used based on the assumption that the waves are represented by a white noise or a wide-band random process. In some cases the non-linear hydrodynamic drag force was replaced by an equivalent linear force. It is obvious that the reported results are not accurate and not realistic since the spectral density of ocean waves occupies a narrow frequency band and the non-linear hydrodynamic drag term requires special treatment. The present paper provides an alternative approach based on the stochastic average concept with a special treatment to Morison’s non-linear drag term. The approach has not been previously used for non-linear mechanical problems. It is convenient for the investigation of ocean systems subjected to non-linear hydryodynamic loading. In addition, the fluid velocity is represented by a Gaussian random function of time which can be interpreted as the response of a certain
683
linear shaping filter to a broadband random excitation. The case of non-Gaussian ocean waves will be considered in Part II. 2. STATEMENT OF THE PROBLEM
The effect of ocean waves on elastic structures introduces hydrodynamic forces of two types, namely inertia and drag forces. The inertia force is the sum of two components. The first is due to a buoyancy force acting on a structure in the fluid due to a pressure gradient generated from the flow acceleration. The buoyancy force is equal to the mass of the fluid displaced by the structure multiplied by the acceleration of the flow U , where U is the flow velocity and a dot denotes differentiation with respect to time t. The second component of the inertia force is due to the added mass which is proportional to the relative acceleration between the structure and the fluid accounts for the flow entrained by the structure and is given by the expression cI rA(U − X ), where cI is the added mass coefficient, A is the area of the cross-section of the structure, r is the fluid density and X is acceleration of the structure. The drag force is the sum of the viscous and pressure drags produced by the relative velocity between the structure and the flow. This type of hydrodynamic drag is proportional to the square of the relative velocity. In general, most hydrodynamic modelling has been based on Morison’s equation [1], which contains inertia and non-linear drag forces only. The inclusion of the hydrodynamic forces in the mathematical modeling will lead to a non-linear partial differential equation. Using the Galerkin approach one may generate a set of coupled ordinary non-linear differential equations. In order to demonstrate the proposed stochastic averaging approach, the analysis will be restricted to single mode excitation. A simple model of a single-degree-of-freedom elastic ocean structure is shown in Figure 1. This system has been extensively considered by several authors who used different approximate techniques. The differences between the present analysis and the previous work by others will be discussed. The mathematical model of this system is [12] X + 2z0 v0 X + v02 X =
1 {rAcm U + rAcI X + 12 r=U − X =(U − X )DcD }, m0
(1)
where X is the displacement of the structure in line with the flow, m0 is the mass per unit length of the structure, and z0 and v0 are the damping factor and the natural frequency of the structure (in vacuum), respectively. The right-hand side of equation (1) represents the Morison’s hydrodynamic force components (inertia and drag), where cm is an inertia coefficient of the flow, U is the flow velocity (which is a random function of time), cD is
Figure 1. A schematic diagram of a single-degree-of-freedom oscillator in a fluid flow.
. .
684
the drag coefficient and D is the width of the structure. The non-linear drag force acts in direction of the relative flow velocity U − X . Introducing the non-dimensional parameters t = Vt,
q = X/D,
u0 = U/VD,
p = dq/dt,
where V = m0 v /(m0 + rAcI ), equation (1) can be written in the state vector form 2
2 0
dq = p, dt
dp du = −q − ap + b 0 + g=u0 − p=(u0 − p), dt dt
(2)
where a = 2z0 zm0 /(m0 + rAcI ), b = rAcm /(m0 + rAcI ) and g = DcD r/2v0 zm0 (m0 + rAcI ). The non-linear term =u0 − p=(u0 − p) has received special treatment in the literature based on approximate representation. For example, Dalzell [13] replaced this term by the polynomial expansion =u0 − p=(u0 − p) 1 s Ck k
(u0 − p)k , D0k − 2
k = 1, 3, 5, . . . ,
(3)
where the Ck are constant coefficients to be determined by using a least-squares fit and D0 q 0 is the maximum relative velocity defined by the range −D0 Q u0 − p Q D0 . Another treatment adopted by others (see, for example, [6]) is given by the representation =u0 − p=(u0 − p) 1 =u0 =u0 − 2=u0 =p + p 2 sgn (u0 ),
(=u0 = w =p=).
(4)
The two approximations (3) and (4) will be used in the present analysis for comparison with results obtained for the expression of the non-linear hydrodynamic drag g=u0 − p=(u0 − p). In the analysis of reference [6], the hydrodynamic drag coefficient g was taken to be of order ze, while in the present analysis it is taken to be of order e. The two types of ordering will yield different results when one deals with Ito stochastic calculus. The difference between the two is discussed in section 6. Furthermore, the expressions =u0 =u0 and =u0 = were replaced in reference [6] by Hermite polynomials which are similar to a certain extent to Dalzell’s expansion (3). 3. FIRST ORDER FILTER
The wave velocity function u0 (t) is assumed to be a narrow-band random process. It is possible to model this process as the response of a shaping filter to a wide-band random process. In this section, this process is represented as the output of the first order shaping filter du0 + du0 = V(t), dt
(5a)
where d is a positive filter parameter which governs the filter frequency bandwidth and V(t) is a scalar wide-band random process whose correlation time is much smaller than any characteristic period of the system oscillations. It is well known that the output of this filter can be approximated by the solution of the Ito stochastic differential equation (SDE) du0 = −du0 dt + zn dB(t),
(5b)
where V(t) has been replaced by the formal derivative of a standard Brownian motion B(t), i.e., V(t) = zn dB(t)/dt, and n is the intensity of the random process V(t).
685
Equation (5b) has a well-known stationary solution (Ornstein–Uhlenbeck process) with the following power spectral density Su 0 (v) and autocorrelation function Ru 0 (t): Su 0 (v) =
n 1 , 2p d 2 + v 2
Ru 0 (t) =
n −d=t= e . 2d
(6)
The function Su 0 (v) is plotted in Figure 2 (solid curve) for d = 0·6, n = 1. The expression of Su 0 (v) reveals that, for small values of d, the spectrum of the stationary solution of (5b) is mainly concentrated about a narrow frequency band close to the origin. Blevins [6] reported that most offshore structures have a fundamental natural frequency around 1 Hz, while the energy content of ocean waves is concentrated at frequencies of less than 0·2 Hz. The filter output stationary probability density function f(u0 ) and the stationary mean of the absolute value of =u0 = are, respectively, f(u0 ) =
0 1
1 u2 exp − 0 , 2D1 z2pD1
E[=u0 =] = z2D1 /p,
(7)
where D1 = n/(2d). Equations (2) and (5b) constitute a set of Ito SDE with the corresponding initial conditions. The statistical analysis of these equations will be carried out using the stochastic averaging of the response energy envelope. 4. STOCHASTIC AVERAGING OF THE RESPONSE ENERGY
The stochastic averaging method of the energy envelope was originally introduced by Stratonovich [14]. We will use the generalization of this method [15] (see also Appendix A). The parameters a, g, d 2 and b 2 are assumed small and of the same order. Introducing the polar co-ordinates I and 8 defined by the transformation q = I cos 8,
p = −I sin 8,
8 = c + t,
(8)
Figure 2. Power spectra of different shaping filters: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; — – —, fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, qk = 1·5, n = 2/17.
. .
686
equations (2) can be rewritten in the standard form dI = {−aI sin2 8 + bdu0 sin 8 − g sin 8=u0 + I sin 8=(u0 + I sin 8) + (2I)−1b 2n cos2 8} dt − bzn sin 8 dB(t), dc = {−a sin 8 cos 8 + I−1bdu0 cos 8 − gI−1 cos 8=u0 + I sin 8=(u0 + I sin 8) − (I)−2b 2n sin 8 cos 8} dt − I−1bzn cos 8 dB(t),
(9)
where the Ito correction term is included in these equations. It should be noticed that the state co-ordinates I and c are slowly varying with time t, while u0 and 8 exhibit fast variation. In this case, the averaging of the drift vector and diffusion matrix of equations (9) will be taken with respect to 8 and with respect to the stationary distribution of the ergodic process u0 (t), keeping I and c frozen. This process gives the averaged equations (see Appendix A). dI = {−12 aI + c/4I − gI 2x¯ (I)} dt + zc/2 dB1 (t), dc = (1/I)zc/2 dB2 (t),
(10)
a x(x)f(u0 ) du0 , x = u0 /I, c = b 2n and B1 (t) and B2 (t) are two independent where x¯ (I) = f−a Brownian motion processes, each of unit intensity. The function x is some deterministic function of the new wave velocity parameter x evaluated for two different domains (see Appendix B).
x(x) =
1 2p
g
2p
sin 8=x + sin 8=(x + sin 8) d8,
(11a)
0
or
F g 1 (4 cos u + 4u sin u − 4 cos3 u), 3 2p G f
x(x) = G ,
for =x= e 1 (11b) for =x= E 1,
where (x = sin u).
The function x(x) is plotted in Figure 3(a). It satisfies the two inequalities x(x) e =x=
and x e 4/3p.
(12)
The function x(x) may also be written in the form (see also Figure 3(b)) x(x) = =x= + h(x),
Figure 3. The dependence of (a) the function x(x) on x = u0 /I, and (b) the function h(x) on x = u0 /I.
687
where h(x) = (2/3p)[3z1 − x + 3x sin x − (1 − x ) ] − x, for 1 e x e 0, and h(x) = h(−x), h(x) = 0 for =x= q 1. The component h is only significant over the range 1 e x e 0, as shown in Figure 3(b). Statements (11) are accurate and provide better representation of the non-linear Morison’s drag term in the averaging process. They facilitate the analysis of the averaged system. To this end let us introduce the Hamiltonian H = I 2/2. In terms of H, the first equation (10) takes the form which includes the Ito correction term after the transformation. 2
−1
2 3/2
$
%
c dH = −(a + 2gE=u0 =)H + − 2gHz2Hh¯ (H) dt + zcH dB1 , 2 or
$
%
c dH = −aH + − U(H) dt + zcH dB1, 2
(13a)
a = a + 2gE[=u0 =] = a + 2gzn/dp,
(13b)
g
(13c)
where
U(H) = 2gHz2Hh¯ (H) = 2gHz2H
a
h(u0 /z2H)f(u0 ) du0 = 8gH 2y(H),
−a
y(H) =
g
1
h(x)f(xz2H) dx.
(13d)
0
Notice that the integration limits are changed from 2a to 21 (or from 0 to +1 due to symmetry) because h(x) is non-zero only over the domain −1 Q x Q +1, and both h(x) and f(xz2H) are even functions with respect to x. The phase angle c(t) only constitutes a diffusion component, as the drift coefficient vanishes. It follows from relations (11) and (12) that the non-linear term g=u0 − p=(u0 − p) in the expression for the drag force ‘‘in average’’ has a dissipative nature (−x Q 0). The stationary probability density of the response process H of system (13) can be derived by using the Fokker–Planck equation in an explicit form:
$ $
exp[−2(a/c)H] exp −(16g/c) f(H) =
g6 a
% %7
hn(h) dh
0
exp[−2(a/c)H] exp −(16g/c)
0
g g
H
H
hn(h) dh
0
.
(14)
dH
The following inequality may be established based on equation (13) after using equation (12):
$
%
0
1
c c dt + zcH dB1 , dH = −aH + − U(H) dt + zcH dB1 E −aH + 2 2
(15)
where the right side does not include the contribution of h(H). The corresponding level of H may be governed by the differential equation (16) dH' = (−aH' + c/2) dt + zcH' dB1 . The corresponding energy level forms an ‘‘upper bound’’ for the original averaged equation (13). It is clear from equations (16) and (13a) that the non-linear hydrodynamic drag force enhances the damping coefficient. Equation (16) will be referred to as the ‘‘simplified’’ model.
. .
688
The function h(x) (or the drift term U(H)) indicates the difference between Morison’s non-linear term and its approximate representation (4) as adopted in reference [6]. It is worth noting that if representation (4) is used it will lead to the approximate averaged equation (16). Later we will discuss the conditions under which this difference becomes considerable. In fact, the simplified model is the result of some equivalent linearization, when Morison’s non-linear drag term g=u0 − p=(u0 − p) is replaced by −2gE[=u0 =]p, as shown by Spanos and Chen [4] for the deterministic case. The effect of h(x) on the system is negligible if =x= e 1 and it becomes significant if =x= Q 1. For the simplified model the response energy level is governed by the Gibbs distribution, i.e., f (H') =
0
1
2a 2a exp − H' . c c
(17)
The corresponding density for the amplitude I = zq 2 + p 2 is known to follow Rayleigh’s distribution f(I) = (2aI/c) exp(−aI 2/c). If the original co-ordinates are used, the response will be governed by a Gaussian distribution. Equation (16) can now be analyzed for the statistical behavior of the response energy function H'. Instead of solving for the probability density function, it is convenient to establish the corresponding equations for the n-dimensional characteristic functions Fn (l1 , . . . , ln ; t1 , . . . , tn ), where the li are arbitrary dummy real parameters. These characteristic functions are represented by a set of linear partial differential equations of the first order and may be integrated in a closed form as outlined in Appendix C. The stationary characteristic functions of system (16) for n = 1, 2 are obtained in the closed form: 1 F1 (l) = , F2 (l1 , l2 ; t2 − t1 ) = 1/{1 − id(l1 + l2 ) − l1 l2 d 2[1 − e−a(t2 − t2 )]}, (18) 1 − idl where d = c/(2a). The autocorrelation function RH' (t) and the power spectral density SH' (v) associated with the stationary process H'(t) are obtained from the solution F2 , i.e., RH' (t) =
b
1 2F 2 = d 2{e−a=t= + 1}, 1(il1 ) 1(il2 ) l = 0
SH' (v) =
i
d2 a . p a2 + v2
(19)
Expressions (19) are of the same type as those of the excitation u0 (t). The stationary mean value of H' can be obtained from the first relation (18) and is found to be c/(2a). It is interesting to notice that, for relatively large values of n, the stationary mean value is proportional to zn, i.e., E[H'] 2 (b 2/4g)zpdn . One may check that, for the one-dimensional equation (16), the border H = 0 is not accessible for the process H(t), and any perturbation will bring the system energy back to its original value H(t) (see Appendix D). Equally important is to note that the existence of the stationary solution of the SDEs (2) and (5b) (for the general case) can be established by applying the generator L( · ) (see, for example, reference [16]): L( · ) = f(X, t)
$
%
7
1( · ) 1 1 1T + 2 Tr ( · ) [g(X, t)ngT(X, t)] , 1X 1X 1X
(20)
on any response function, such as the Lyapunov function, based on the general Ito stochastic differential equation dX = f(X, t) dt + g(X, t) dB(t), (21) where X is a state vector, e.g., for the present system X = {q, p, u0 }T, n is the covariance matrix of the vector dB. For simplicity we consider the case when the co-ordinate q is
689
bounded. It is convenient to shift the co-ordinate p by z = p − u0 . Accordingly, the governing equations of the system and filter take the form dq = {z + u0 } dt, dz = {−q − az − g=z=z − r¯ u0 } dt + (b − 1)zn dB(t),
(22)
du0 = −du0 dt + zn dB(t), where r¯ = a + bd − d. Introducing the response function u = (A1 z 2 + A2 u02 )/2 as a Lyapunov function (where A1 and A2 are positive constants), and applying equation (20) gives
0
1 0
1
r¯ r¯ 2A1 2 n n u0 − A 2 d − u0 + A1 (b − 1)2 + A2 . 2a 2 2 4a 2
Lu = −A1 zq − gA1 z 2=z= − aA1 z +
(23) One of the main conditions for the existence of a stationary solution (see reference [16], p. 90) is that the supremum of the right side of equation (23) decreases as both =u0 = and =z= increase, i.e., r¯ 2 sup Lu(z, u0 ):−a as R:a, if A2 d q A1 . (24) 2 2 4a z +u qR 0 The result of equation (24) indicates that the stationary solution exists and forms a stable and unique attractor. Note that the approximation given by relation (3) will lead to the averaged equation dH = {−(a + 2gc˜1 )H + b 2n/2 − g[c˜3 (2H)2 + c˜5 (2H)3 + · · ·]} dt + bznH dB1 .
(25)
Each of the coefficients c˜i represents some series, so one should keep sufficiently many terms in order to achieve high accuracy. Notice that the difference between equation (25) and the exact equation (13) is in the expressions of the drift coefficient. 5. SECOND ORDER SHAPING FILTER
5.1. - vf $ 1 The first order shaping filter characterizes the wave velocity whose peak spectrum occurs at zero frequency and its bandwidth is governed by the system parameters. However, if the peak spectrum is shifted from zero frequency, as in the case of any real ocean system, one should use a higher order filter. In this section, we consider the case when the wave velocity u0 (t) is represented by the output of the shaping filter du0 = p0 dt,
dp0 = {−vf2 u0 − 2vf zp0 } dt + zn dB(t),
(26)
where vf and z are the natural frequency and damping factor of the filter, respectively. Filter (26) has a well-known stationary solution, with the following power spectral density and autocorrelation function, respectively: Su 0 (v) =
Ru 0 (t) = where v02 = vf2 (1 − z 2).
0
n 1 , 2p (vf2 − v 2)2 + (2vvf z)2
1
n vz cos v0 t + f sin v0 =t= exp(−zvf =t=), 4vf3 z v0
(27)
. .
690
The function Su 0 (v) is plotted by the dashed curve in Figure 2 for z = 0·4, vf = 0·5 and n = 0·1. The output stationary probability density function for the filter (26) is f(u0 , p0 ) =
6 0
1 p2 u2 exp −12 0 + 0 D1 D2 2pzD1 D2
17
,
where D2 = n/4vf z and D1 = D2 /vf2 = n/4vf3 z. In order to derive the response statistics of systems (2) and (26), the parameters a, g, z 2 and b 2 with be assumed to be of the same order of smallness. In this case, these equations must be transformed into a form suitable for averaging, in the non-resonance case (i.e., vf $ 1). Introducing the transformation 1 0 K u0L K 0 1 G p0G=G 0 b/(1 − vf2 ) G qG G 0 k p l k−bvf2 /(1 − vf2 )
0 0 1 0
0LKu0 L 0GGp0 G 0GGq1 G 1lkp1 l
(28)
into equations (2) gives dq1 = {p1 + 2vf zb p0 } dt − b zn dB, dp1 = {−q1 − ap1 + ab vf2 u0 + g=d u0 − p1 =(d u0 − p1 ) dt,
(29)
where b = b/(1 − vf2 ), d = [1 + (b − 1)vf2 ]/(1 − vf2 ), vf $ 1, and the filter state equations (26) remain the same. Transformation (28) eliminates the linear coupling from equation (2), where the unperturbed conservative parts of equations (2) and (26) take a Hamilton form. Introducing the slow variables I and c, defined by the co-ordinate transformation q1 = I cos 8,
p1 = −I sin 8,
8 = c + t,
(30)
equations (29) can be rewritten in the standard form dI={−aI sin2 8 − avf2 b u0 sin 8 + 2zvf b p0 cos 8 − g sin 8=d u0 + I sin 8=(d u0 + I sin 8) + (2I)−1b 2n sin2 8} dt − b cos 8zn dB(t), dc = {−2I−1zvf b p0 sin 8 − a sin 8 cos 8 − I−1ab vf2 u0 cos 8 − gI−1 cos 8=d u0 + I sin 8=(d u0 + I sin 8) + (I)−2b 2nsin 8 cos 8} dt + I−1b sin 8zn dB(t). (31) Averaging gives a similar set of averaged SDE (10), where c = b 2n, x(I) =
gg
a
x(d u0 I−1 )f(u0 , p0 ) du0 dp0 .
(32)
−a
The equation for the energy, H = I 2/2, takes the form of equation (13), with the following expression for the parameter a: a = a + 2gE[=d u0 =] = a + 2gd zn/2pzvf3 ,
(33a)
U(H) = 8gd −1H 2n(H),
(33b)
and function U(H):
691
where n(H) = (1/z2pD1 )
g
1
h(x) exp(−Hx 2/d 2D1 ) dx.
(33c)
0
The phase angle c(t) only constitutes a diffusion component as the drift coefficient vanishes. It follows from relations (11) and (12) that the non-linear term g=u0 − p=(u0 − p) in the drag force expression has a dissipative nature (−x Q 0) and reduces the stationary expected value of H. The stationary probability density of the response process of system (13) can be derived by using the Fokker–Planck equation in an explicit form:
$ $
g g
exp(−2aH/c) exp −(16g/d c) f(H) =
g6 a
% %7
hn(h) dh
0
H
exp(−2aH/c) exp −(16g/d c)
0
H
hn(h) dh
0
.
(34)
dH
This expression is identical to the probability density function obtained for the first order filter case. However, the two results are affected by the filter constants d, d, vf and z, which are included in the parameters a and c. Notice also that the function n(h) is different in both cases as given by relations (13c) and (33c). Following the analysis of the previous section, it is not difficult to show that the non-linear drag force on average enhances the damping and we also may establish the inequality (15). Again, equality (16) establishes the upper bound of the response energy level of the structure and one may determine all finite-dimensional non-stationary characteristic functions of the energy H, as outlined in Appendix C. In conclusion, it is possible to write the expression for the mean value of the energy in terms of the initial variables q, p: E[12 (q 2 + p 2)] = E[H] + 12 b 2(1 + vf2 )D2 . 5.2. : vf 1 1 Consider the case of the shaping filter (26) when the frequency vf is in the neighborhood of 1. Therefore, the detuning (vf − 1) is a small value. Introducing the linear transformation (which eliminates the linear coupling from the system and filter equations when g = z = 0, n = 0)
K qL K 1 G p G=G −a/2 G u0G G 0 k p0l k 0
0 V 0 0
a1 −a2 vf2 1 0
a2LK q1 L a1GG p1 G 0GG u0 G 1lk p0 l
(35)
into equations (2) gives dq1 = {−aq1 /2 + Vp1 + 2za2 p0 } dt − (a2 /vf )zn dB, dp1 = {−Vq1 − ap1 /2 + 2za3 vf p0 + g/V=aq1 /2 − Vp1 + Q=(aq1 /2 − Vp1 + Q)} dt − a3 zn dB,
(36)
where V = z1 − a 2/4, a1 = D−1abvf2 , a2 = D−1b(1 − vf2 ), D = (1 − vf2 )2 + a 2vf2 , a3 = ab(1 + vf2 )/2DV and Q = (1 + a2 vf2 )u0 − a1 p0 . It is assumed that a, g, vf − 1, z 2 and b 2/3 are of the same order of smallness. Introducing the slow variables I, and c, defined by q = I cos 8,
p = −I sin 8,
8 = c + Vt,
. .
692
and averaging, the following averaged SDE are obtained (see Appendix A): x(I)} dt + zc/2 dB1 (t) dI = {−12 aI + c/4I − gI 2 dc = −(ga/2V)I 2 x(I) dt + (1/I)zc/2 dB2 (t)
(37)
a x(QI−1 )f(u0 , p0 ) du0 dp0 , c = n(a22 /vf2 + a32 ), and dB1 (t) and dB2 (t) are where x(I) = ff−a two independent Brownian motion processes, each of unit intensity. In terms of H, equations (37) take the form
6
c dH = −aH + − 2gHz2Hh¯ (H)} dt + zcH dB1 (t), 2 dc = {−(ga/V)(zH/2E[=Q=] + Hh¯ (H))} dt + 12 zc/H dB2 (t),
(38)
where a = a + 2gE=Q=, E=Q= = z2DQ /p, h¯ (H) = Eh(Q/z2H) and DQ is the covariance of the variable Q which can be derived from the filter equations. It is seen from equation (38) that the phase angle c(t) constitutes drift and diffusion components in contrast to the non-resonance case. The stationary probability density of the response process H can now be derived by using the Fokker–Planck equation in an explicit form. It has the form (14), where the expression for the n(H), given by equation (13c) as well as the probability density (7), are of the same form, and the co-ordinate u0 should be replaced by Q (which is normally distributed with zero mean and covariance DQ ). 6. FOURTH ORDER SHAPING FILTER FOR THE WAVE EXCITATION
For better fitting of the known Pierson–Moskowitz spectrum [17], following Spanos [18], a fourth-order shaping filter may be considered: d3u0 d2u0 du0 d3V d4u0 + a0 u0 = 3 , 4 + a3 3 + a2 2 + a1 dt dt dt dt dt
(39)
where a3 = c + c¯ , a2 = k + k + cc¯ , a1 = ck + c¯k, a0 = kk , c, c¯ , k and k are some constants satisfying the stability condition a1 (a2 a3 − a1 ) − a0 a32 q 0.
(40)
If condition (40) holds then Ito’s SDE (39) has a stable stationary solution with the power spectral density Su 0 (v) =
nv 4 . [(v − k) + (cv) )][(v 2 − k)2 + (c¯v)2] 2
2
2
(41)
The parameters a1 and a3 in equation (39) are damping coefficients. When a1 = a3 = 0 and V = 0, the system (39) is conservative with the two natural frequencies v1 = za2 /2 − za22 /4 − a0 ,
v2 = za2/2 + za22 /4 − a0 .
Figure 2 includes two typical spectra plotted according to relation (41) for two different constant coefficient sets. The dotted curve is plotted for the case of c = 0·17, c¯ = 0·476, k = 0·147 and k = 0·375, and for a white excitation level n = 1/130. The spectrum shown by short and long dashes is plotted for the case of c = 0·338, c¯ = 0·946, k = 0·586 and k = 1·5, and for the white excitation level n = 2/17. It is seen that this spectrum is closer to the resonance frequency than the previous one.
693
Equation (39) can be written as a set of four linear SDE: dx1 = x2 dt + zn dB,
dx2 = x3 dt − a3 zn dB,
dx3 = x4 dt + (a32 − a2 )zn dB, (42)
dx4 = {−a0 x1 − a1 x2 − a2 x3 − a3 x4 } dt + (−a33 + 2a2 a3 − a1 )zn dB, where x1 (t) = u0 (t). For further analysis, it is convenient to apply the following linear transformation of the variables xj which eliminates the linear coupling from equation (42) when a1 = a3 = 0, V = 0:
Ky1L Kv22 0 Gy2G = G 0 v22 Gy3G Gv12 0 ky4l k 0 v12
1 0 1 0
0 LKx1L 1 GGx2G , 0 GGx3G 1 lkx4l
provided that v1 $ v2 .
(43)
Now the system and filter equations take the form dq = p dt, dp = {−q − ap + bD−1(dy1 /dt − dy3 /dt) + g=D−1(y1 − y3 ) − p=(D−1(y1 − y3 ) − p)} dt, dy1 = y2 dt + (a32 − v12 )zn dB, dy2 = {−v12 y1 − h1 y2 + h2 y4 } dt − g2 zn dB, dy3 = y4 dt + (a32 − v22 )zn dB, dy4 = {−v22 y3 − h1 y2 + h2 y4 } dt − g1 zn dB,
(44)
where D = v22 − v12 , h1 = D−1(a1 − a3 v12 ), h2 = D−1(a1 − a3 v22 ), g1 = a1 + a33 + v12 a3 − 2a2 a3 , and g2 = a1 + a33 + v22 a3 − 2a2 a3 . The linear transformation of the variables q and p (which is valid in the non-resonance case v1 $ 1, v2 $ 1), q1 = q − b1 y2 + b2 y4 , b1 = b/D(1 − v12 ),
p1 = p + b1 v12 y1 − b2 v22 y3 , b2 = b/D(1 − v22 ),
(45)
eliminates the coupling from equation (44) when the damping vanishes (i.e., a1 = a3 = a = 0) and V = 0. In this case the governing equations of q1 and p1 are dq1 = {p1 − d1 y2 + d2 y4 } dt + e1 zn dB, dp1 = {−q1 − ap1 + ab1 v12 y1 − ab2 v22 y3 + g=D−1(y1 − y3 ) − p1 + b1 v12 y1 − b2 v22 y3 = × (D−1(y1 − y3 ) − p1 + b1 v12 y1 − b2 v22 y3 )} dt + e2 zn dB,
(46)
where e1 = −b2 g1 , e2 = b + b1 v12 (a32 − v12 ) − b2 v22 (a32 − v22 ), d1 = bh1 /(1 − v12 )(1 − v22 ) and d2 = bh2 /(1 − v12 )(1 − v22 ). Introducing the slow variables I and c defined by equation (30) and averaging, we obtain equations (10) in which c = n(e12 + e22 ), d1 =
x¯ (I) = Ex(QI−1 ),
1 + (b − 1)v12 , D(1 − v12 )
d2 =
Q = d1 y1 − d2 y3 , 1 + (b − 1)v22 . D(1 − v22 )
(47)
694
. .
In terms of H the first equation takes the form (13) by setting a = a + 2gE[=Q=], and U(H) is given by equation (13b), where q0 must be replaced by Q. The mean square of Q can be estimated from the last four equations (44) by using, for example, the moments method. The stationary probability density function f(Q) is governed by equation (7) where D1 must be replaced by the covariance of Q. The analysis will lead to results similar to those of equations (13) and (14).
7. DISCUSSION OF RESULTS AND CONCLUSIONS
The analytical results derived in the previous two sections provide complete information regarding the statistics of a linear single-degree-of-freedom system subjected to a non-linear random hydrodynamic force generated by three different filters. With reference to the first order shaping filter, the spectrum is governed mainly by the parameter d and becomes almost flat for large values of d. For the case of the second order filter, the filter damping coefficient governs the wave bandwidth. As the filter damping ratio decreases, the band width of the wave velocity becomes more narrow. The fourth order filter emulates to a certain extent the Pierson–Moskowitz spectrum which represents the actual spectrum of ocean waves. The spectra of these three filters are shown in Figure 1. These filters are subjected to a wide-band random excitation with different power spectrum levels. It is seen that as the order of the shaping filter increases it requires a different level of the excitation, depending on its physical parameters. The probability density functions of the system response to the output of these three filters are obtained in a closed form and are plotted in Figure 4 for the system parameters a = 0·01, b = 0·25 and g = 0·01. The characteristics of the density curves are primarily
Figure 4. Probability density functions of the system response energy to the output of three filters for system parameters a = 0·01, b = 0·25 and g = 0·01: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
695
governed by the spectral density level of the white noise input to each filter and how close the peak of the spectrum is to the system’s natural frequency. In general, the system response is non-Gaussian due to the non-linear hydrodynamic drag. If the approximation of Eatock Taylor and Rajagopalan as given by relation (4) is used, the present analysis will give Gaussian probability density for every shaping filter. However, according to the stochastic averaging adopted by Eatock Taylor and Rajogopalan [6], the response density is not Gaussian because of two main factors. The first is that they replaced the hydrodynamic drag force by a Hermite polynomial and the second is that they considered g to be of the order ze, which results in both drift and diffusion terms in the Ito equation of the governing equations of motion. In the present analysis, the coefficient g is taken to be of the order e in order to take into account the complete effect of non-linear hydrodynamic drag without approximation. The advantage of this is that the diffusion coefficient due to the non-linear drag term will vanish and this avoids several analytical complications. It is found that the probability density curves, obtained by using the approximation (4) and by using the exact drag expression g=u0 − p=(u0 − p), are close for small values of g and the difference becomes larger for larger values of g, as shown in Figure 5. This difference is a measure of the deviation of the response from normality. It is seen that the first order filter (solid curve), whose center frequency is zero and shifted away from the system natural frequency, requires a very high level of filter excitation spectrum, i.e., n = 1. In view of this large level of n, the probability density of the system response spreads over a relatively higher response energy level. For the second order filter (dashed curve), the probability density function spreads over higher response levels as the filter natural frequency approaches the system natural frequency. Fourth order filters (dotted and dash-dot curves) yield different probability density curves depending on the filter frequency
Figure 5. Probability density functions of the system response energy to the output of a second order filter for z = 0·4, vf = 0·5, n = 0·2, a = 0·01, b = 0·25 and two different values of the non-linear drag coefficient g: ----, using exact drag expression g(u0 − p)=u0 − p=; – – – –, using the approximate drag expression given by equation (4).
696
. .
Figure 6. The dependence of the expected value of the response energy on the white noise excitation level for system parameters a = 0·01, b = 0·25 and g = 0·02: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
and the value of n. It is seen that for small n and lower center frequency, the response probability density has a peak which is concentrated around the zero response level. On the other hand, as both the filter center frequency and excitation level n increase, the probability of higher response level increases.
Figure 7. The dependence of the expected value of the response energy on the non-linear hydrodynamic drag coefficient for system parameters a = 0·01, b = 0·25 and g = 0·01: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
697
Figure 8. The normalized mean value of the system response energy (with respect to the expected value of the linear system response energy) versus the white noise excitation level for system parameters a = 0·01, b = 0·25 and g = 0·01: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; — – —, fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
The dependence of the expected value of the response energy level E[H] on the excitation level n is shown in Figure 6 for a non-linear hydrodynamic drag coefficient g = 0·02. For the three filters, the response mean value increases monotonically with n in a non-linear
Figure 9. The dependence of the mean value of the system response energy on the filter center frequency for the case of a second-order shaping filter, in the neighborhood of resonance for system parameters a = 0·01, b = 0·25 and g = 0·01.
698
. .
fashion and is of the form zn. The dependence of E[H] on the non-linear hydrodynamic drag coefficient g is shown in Figure 7 for the three filters. It is seen that the non-linear drag has a remarkable effect on suppressing the system response. The normalized response energy level with respect to the corresponding linear value is plotted as a function of the excitation level n in Figure 8. The deviation of the normalized response from the value 1 measures the influence of the hydrodynamic drag. Again, comparison with results by Eatock Taylor and Rajagopalan reveals that their results are more conservative. Figure 9 shows the expected value of the response energy level for the case of the second order shaping filter in the neighborhood of resonance, vf = 1, for the system parameters a = 0·01, b = 0·25, g = 0·01. Unlike other cases, the response exhibits peak values only at resonance. The analysis of this paper is based on three different shaping filters. The Pierson– Moskowitz (PM) spectrum [17] is more realistic and Part II of this investigation will consider the PM spectrum and a shaping filter which takes into account the non-normality of sea waves. Furthermore, the system inherent non-linearities have to be included in the future analysis. These factors are important in order to examine important issues such as level crossings, local extremes of the response and other related reliability problems. ACKNOWLEDGMENTS
This research is supported by a grant from Office of Naval Research under grant No. N000149310936, Dr Julia Abrahams is the Scientific Officer. REFERENCES 1. J. R. M, M. P. O’B, J. W. J and S. A. S 1950 Petroleum Transactions 189, 149–154. The forces exerted by surface waves on piles. 2. R. G. T 1977 The Structural Engineer 55(5), 209–222. Continuous random wave loading on structural members. 3. R. E T 1975 Off-shore Structures. 125–132. London: ICE. Structural dynamics of offshore platforms. 4. P.-T. D. S and T. W. C 1980 International Journal of Non-linear Mechanics 15, 115–126. Response of a dynamical system to flow-induced load. 5. A. W. L 1986 Transaction of the American Society of Civil Engineers, Journal of Structural Engineering 112, 2416–2429. Nonlinear structural response in random waves. 6. R. E T and A. R 1982 Journal of Sound and Vibration 83, 401–431. Dynamics of offshore structures: parts 1 and 2. 7. S. K. T and J. M. N 1992 Transaction of the American Society of Civil Engineers, Journal of Engineering Mechanics 118, 942–960. Parametric and external excitation of marine risers. 8. M. S, C. Y and R. V 1977 Transaction of the American Society of Civil Engineers, Journal of Structural Mechanics 5(2), 135–146. Dynamic analysis of fixed offshore structures subjected to wind generated waves. 9. P.-T. D. S and J. E. H 1981 Transaction of the American Society of Mechanical Engineers, Journal of Energy Resources Technology 103, 43–249. Linear prediction theory for digital simulation of sea waves. 10. N. K. L and W. A. H 1984 Transaction of the American Society of Mechanical Engineers, Journal of Energy Resources Technology 106, 466–472. Time series simulations of wide-band spectra for fatigue tests of offshore structures. 11. P. B, I. L, S. R. W, T. M and D. K 1992 in Nonlinear Stochastic Mechanics, 83–101. Berlin: Springer. Application of nonlinear stochastic mechanics in offshore engineering. 12. R. D. B 1977 Flow-induced Vibration. New York: Van Nostrand Reinhold. 13. J. F. D 1978 Journal of Ship Research 22(3), 178–185. A note on the form of ship roll damping.
699
14. R. L. S 1963 Topics in the Theory of Random Noise. Vol. 1. New York: Gordon & Breach. 15. R. Z. K 1968 Cybernetics (Kybernetika Cislo 3, Rocnik) 3, 260–279. On averaging principle for Ito stochastic differential equations (in Russian). 16. R. Z. K 1980 Stochastic Stability of Differential Equations. Alphen: Sijthoff and Noordhoff. 17. W. J. P and L. M 1964 Journal of Geophysical Research 69, 5181–5190. A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. 18. P.-T. D. S 1983 Transation of the American Society of Mechanical Engineers, Journal of Energy Resources Technology 105, 300–309. ARMA algorithms for ocean wave modeling. 19. V. S. P and I. N. S 1987 Stochastic Differential Systems. New York: John Wiley. 20. W. E. W 1980 Partial Differential Equations. Oxford: Clarendon Press. 21. Y. V. P and Y. A. R 1969 Probability Theory: Basic Concepts, Limit Theorems, Random Processes. New York: Springer-Verlag. APPENDIX A: STATEMENT OF AVERAGING PRINCIPLE [15]
Let {X (t), Y (e)(t)}T be a Markov vector in El with X e(t) giving as a slow random process, while Y e(t) experiences fast variation. The two processes may be described by the following system of Ito SDE: (e)
dX (e)(t) = eA1 (X, Y) dt + zes(X, Y) dB(t),
dY (e)(t) = A2 (X, Y) dt + 8(X, Y) dB(t),
with dim X = l1 , dim Y = l2 and dim B = l3 , and subject to the initial conditions X (e)(0) = x0 ,
Y (e)(0) = y0 .
Here, A1 (x, y) and A2 (x, y) are given functions of the dimensions l1 and l2 , respectively, and s(x, y) and f(x, y) are given matrix functions of the dimension l1 × l3 and l2 × l3 , respectively. The elements of the vector B are independent standard Brownian processes. Assume that the vector functions A1 , A2 , s and f satisfy the Lipschitz condition with respect to x and y. Since the process X is slowly varying with respect to time, one may establish the Ito SDE for the process Y (x,y)(t) by freezing the process X: dY (x,y)(t) = A2 (x, Y (x,y) ) dt + 8(x, Y (x,y) ) dB(t),
Y (x,y)(0) = y.
If the following limits exist A 1 (x) = lim E t:a
$g
D (x) = lim E t:a
1 t
$g
t+t
1 t
t+t
%
A1 (x, Y (x,y)(s)) ds ,
t
%
s(x, Y (x,y)(s))s T(x, Y (x,y)(s)) ds ,
t
then under some assumptions [15] the process X (e)(t) weakly converges, on the time interval 0 E t E T/e as e:0, to the Markov process X (0)(t) governed by the Ito SDE dX (0)(t) = eA 1 (X (0) ) dt + zes¯ (X (0) ) dB'(t),
s¯ s¯ T = D .
APPENDIX B: EVALUATING FUNCTION x (11a)
Consider the function x(x) =
1 2p
g
2p
sin 8=x + sin 8=(x + sin 8) d8.
(B1)
0
We split the x-axis into the four domains: (1) x E −1,
(2) −1 E x E 0,
(3) 0 E x E 1
and
(4) x e 1,
. .
700
and evaluate integral (B.1) for each domain as follows: x(x) = −
(1)
1 2p
g
2p
sin 8(x + sin 8)2 d8 = −x.
0
(2) Introducing auxiliary variable u $ [0, p/2], as sin u = −x, we have x(x) =
=
1 2p
$g
p−u
sin 8(x + sin 8)2 d8 −
u
g
2p + u
%
sin 8(x + sin 8)2 d8
p−u
1 {4 cos u + 4 sin u − 43 cos3 u}. 2p
(3) Introducing auxiliary variable u $ [0, p/2], as sin u = x, we have x(x) =
= (4)
1 2p
$g
p+u
sin 8(x + sin 8)2 df −
g
2p − u
%
sin 8(x + sin 8)2 d8
p+u
−u
1 {4 cos u + 4u sin u − 43 cos3 u}. 2p x(x) =
1 2p
g
2p
sin 8(x + sin 8)2 d8 = x.
0
APPENDIX C: EVOLUTION OF CHARACTERISTIC FUNCTION OF ONE-DIMENSIONAL SYSTEM
Consider the one-dimensional Ito SDE (C.1) dX(t) = f(X, t) dt + g(X, t)zn dB, X(t0 ) = x0 , where B(t) is a vector of the independent standard Brownian motions; f(x, t) and g(x, t) are given functions of the x $ R 1. For definiteness we consider the case when both f(x, t) and g(x, t) are given as polynomials in x. The evolution of the first order characteristic function F1 (l; t) = E[eilX(t)] of the random process X(t) at time t, where l is a dummy arbitrary real parameter, is described by the partial differential equation
6 0 1
0 1 0 17
1 1 1 1F1 (l; t) = ilf , t − 12 lg , t ng , t l F1 (l; t). 1t i1l i1l i1l
(C.2)
The second and nth orders characteristic functions of X(t) are, respectively, F2 (l1 , l2 ; t1 , t2 ) = E[eil1 X(t1 ) + il2 X(t2 )], Fn (l1 , . . . , ln ; t1 , . . . , tn ) = E[eil1 X(t1 ) + ··· + iln X(tn )]. The partial differential equation which describes the evolution of Fn is [19]
6 0
1
0
10
1 7
1 1 1 1Fn (l1 , . . . , ln ; t1 , . . . , tn ) = ilnT f , t − 12 lnT g , t ng ,t l F . 1tn i1ln n i1ln n i1ln n n n T
(C.3) It is not difficult to show that the differential equation for the characteristic function of the system described by equation (16) is 1F c 1Fn = ln (−a + 12 icnln ) n + i ln Fn , 2 1tn 1ln
(C.4)
701
with appropriate initial conditions F1 (l1 ; t0 ) = F0 (l1 ), Fn (l1 , . . . , ln ; t1 , . . . , tn − 1 , tn − 1 ) = Fn − 1 (l1 , . . . , ln − 2 , ln + 1 + ln ; t1 , . . . , tn − 1 ). The solution of equation (C.4) is obtained by using the method of characteristics (see, for example, reference [20]) and is given in the form F1 (l1 ; t1 ) = F0 (l1 C−1(l1 , t − t0 ) e−a(t − t0 ) )C−1(l1 , t − t0 ), Fn = Fn − 1 (l1 , . . . , ln − 2 , ln − 1 + ln C−1(ln , tn − tn − 1 ) e−a(tn − tn − 1 ); t1 , . . . , tn − 1 ) × C−1(ln , tn − tn − 1 ), where C(l, t) = 1 − (ilc/2a) [1 − exp(−at)], and F0 is estimated from the initial distribution of the response. For example, if the initial distribution is given by the Gibbs function p0 (x) = d−1 e−x/d , the corresponding characteristic function is F0 (l) = 1/(1 − idl). APPENDIX D: NON-ACCESSIBILITY OF THE BORDER FOR THE ENERGY-AVERAGED EQUATION (16)
Consider the one-dimensional Ito stochastic differential equation (C.1) on the half-line [0, a) with continuous drift and diffusion coefficients and unit intensity white noise. The left border x = 0 is accessible if and only if the function R(x) [21], R(x) = R1 (x)
g
x
0
$g
dy , g (y)R1 (y)
x
R1 (x) = exp −
2
0
%
dy , g 2(y)R1 (y)
is integrable in the neighborhood of zero. For equation (16) we have f = −ax + c/2,
g = zcxn.
In this case, R(x) =
0 1
1 a exp x , x c
and now it is clear that the function
0 1g >6 0 17
1 a exp x x c
x
dy
0
c exp
a y c
is not integrable in the neighborhood of zero, so for the process x(t), the left border x = 0 is not accessible.
RESPONSE STATISTICS OF OCEAN STRUCTURES TO NON-LINEAR HYDRODYNAMIC LOADING, PART I: GAUSSIAN OCEAN WAVES N. M R. A. I Department of Mechanical Engineering, Wayne State University, Detroit, Michigan 48202, U.S.A.
R. K Department of Mathematics, Wayne State University, Detroit, Michigan 48202, U.S.A. (Received 7 March 1994, and in final form 21 June 1994) This paper presents an analytical approach based on the stochastic averaging of the energy envelope to treat the dynamic behavior of single-degree-of-freedom elastic ocean structures. Such systems are usually subjected to a narrow-band random process which may be modelled as the output of a shaping filter. The response co-ordinates of the system and filter experience slow and fast variations with respect to time, respectively. For this class of problems, the method of stochastic averaging is used to establish stochastic Ito differential equations for the process of slow variation. This approach has not previously been used for non-linear mechanical problems. Three different shaping filters (including those possessing a Pierson–Moskowitz spectrum) are employed to model Gaussian random sea waves. The response statistics of the structure are estimated in terms of the excitation spectrum for different levels of non-linear hydrodynamic drag force. It is found that the hydrodynamic drag reduces the system response energy and consequently suppresses the motion of the structure. In addition, the response probability density deviates from normality as the non-linear hydrodynamic drag parameter increases. The case of non-Gaussian random sea waves will be considered in Part II. 7 1995 Academic Press Limited
1. INTRODUCTION
The non-linear dynamics of ocean structures exposed to severe wave conditions pose complex problems for engineers and mathematicians. The dynamic hydroelastic problem of ocean systems involves the interaction of inertia forces, elastic forces and hydrodynamic forces. Basically, the hydrodynamic forces induced on structural elements are non-linear and a combination of inertial and drag contributions [1]. Each force is time-dependent and depends on the geometrical properties of the structure, fluid properties describing the flow field and some ‘‘unknown constants’’ which are determined experimentally. Furthermore, under large deflection these structures experience geometrical non-linearities. The hydrodynamic forces due to sea waves are random, and under extreme environmental conditions are non-Gaussian. The basic problems of these types of systems include the response statistics of non-linear dynamical systems to non-Gaussian excitations, level crossings, local extremes of the response and associated reliability problems. Experimental measurements of the response of laboratory models under wave action have indicated a significant deviation from a Gaussian probability distribution, 681 0022–460X/95/290681 + 21 $12.00/0
7 1995 Academic Press Limited
682
. .
particularly at a high level of response [2]. This behavior is taken to be the result of the non-linearity inherent in the drag component of Morison’s equation, and leads to an increased probability of exceeding a given load or response when compared with a Gaussian distribution at high levels of response. The kurtosis of the basic process serves as a useful measure of the non-linearity. The results of Tickell [2] are based on a long-crested or uni-directional surface elevation whose horizontal and vertical velocities of flow are non-zero mean Gaussian variables. The dynamic behavior of steel-framed offshore structures, including stochastic analysis in the frequency domain and simulation in the time domain, was examined within the framework of the linear theory of small oscillations [3]. Approximate analytical methods such as the equivalent linearization technique [4, 5], perturbation analysis [6], stochastic averaging [6] and Gaussian and non-Gaussian closures [7] have been applied for structural elements described by Morison’s equation [1]. Digital simulation techniques have been used to determine the response statistics of offshore structural elements [8–10]. The main results of these studies indicate that the non-linear hydrodynamic damping has the effect of increasing the rate of decay of the response amplitude. At high wind velocities the response is quasi-static exhibiting a wide-band random process characteristic [8]. A comprehensive study [11] reviewed some recent studies of the non-linear collapse and fatigue analysis of jacket structures, dynamic analysis of jack-up structures, springing and slow-drift response of a tension leg platform, and hysteretic foundation behavior of a concrete gravity platform. The mathematical modelling of ocean systems depends on whether the system is treated as a rigid body or as a continuous elastic structure. The models discussed in the literature involve only a finite number of independent co-ordinates which are described by a set of ordinary differential equations. These include ship roll and coupled pitch roll motions, floating buoys, moored semi-submersibles, rigid gravity offshore platforms with co-ordinates describing sliding and rocking motions, and jacket template platforms with co-ordinates to describe the motion of discrete masses lumped at node points. However, for line components such as rather long beams, pipelines and cables, alternative continuous system models may provide more precise descriptions of component motion. Since a partial differential equation is used to characterize the motion of a continuous line component, a solution may be more involved mathematically than for a corresponding lumped system. For beams, the bending stiffness and longitudinal tension are incorporated in the model. Cables, on the other hand, resist tension but their bending and torsional stiffnesses are usually neglected. These structural elements are usually encountered in pipelines, steel ropes and chains used to stay buoys, floating platforms, compliant towers and ships. The inclusion of large deformation in the mathematical modelling will result in non-linear partial differential equations which include non-homogeneous random forcing terms and/or parametric excitation terms. The dynamic behavior of ocean systems has been treated by using one of four approaches mentioned earlier. These approaches have been used based on the assumption that the waves are represented by a white noise or a wide-band random process. In some cases the non-linear hydrodynamic drag force was replaced by an equivalent linear force. It is obvious that the reported results are not accurate and not realistic since the spectral density of ocean waves occupies a narrow frequency band and the non-linear hydrodynamic drag term requires special treatment. The present paper provides an alternative approach based on the stochastic average concept with a special treatment to Morison’s non-linear drag term. The approach has not been previously used for non-linear mechanical problems. It is convenient for the investigation of ocean systems subjected to non-linear hydryodynamic loading. In addition, the fluid velocity is represented by a Gaussian random function of time which can be interpreted as the response of a certain
683
linear shaping filter to a broadband random excitation. The case of non-Gaussian ocean waves will be considered in Part II. 2. STATEMENT OF THE PROBLEM
The effect of ocean waves on elastic structures introduces hydrodynamic forces of two types, namely inertia and drag forces. The inertia force is the sum of two components. The first is due to a buoyancy force acting on a structure in the fluid due to a pressure gradient generated from the flow acceleration. The buoyancy force is equal to the mass of the fluid displaced by the structure multiplied by the acceleration of the flow U , where U is the flow velocity and a dot denotes differentiation with respect to time t. The second component of the inertia force is due to the added mass which is proportional to the relative acceleration between the structure and the fluid accounts for the flow entrained by the structure and is given by the expression cI rA(U − X ), where cI is the added mass coefficient, A is the area of the cross-section of the structure, r is the fluid density and X is acceleration of the structure. The drag force is the sum of the viscous and pressure drags produced by the relative velocity between the structure and the flow. This type of hydrodynamic drag is proportional to the square of the relative velocity. In general, most hydrodynamic modelling has been based on Morison’s equation [1], which contains inertia and non-linear drag forces only. The inclusion of the hydrodynamic forces in the mathematical modeling will lead to a non-linear partial differential equation. Using the Galerkin approach one may generate a set of coupled ordinary non-linear differential equations. In order to demonstrate the proposed stochastic averaging approach, the analysis will be restricted to single mode excitation. A simple model of a single-degree-of-freedom elastic ocean structure is shown in Figure 1. This system has been extensively considered by several authors who used different approximate techniques. The differences between the present analysis and the previous work by others will be discussed. The mathematical model of this system is [12] X + 2z0 v0 X + v02 X =
1 {rAcm U + rAcI X + 12 r=U − X =(U − X )DcD }, m0
(1)
where X is the displacement of the structure in line with the flow, m0 is the mass per unit length of the structure, and z0 and v0 are the damping factor and the natural frequency of the structure (in vacuum), respectively. The right-hand side of equation (1) represents the Morison’s hydrodynamic force components (inertia and drag), where cm is an inertia coefficient of the flow, U is the flow velocity (which is a random function of time), cD is
Figure 1. A schematic diagram of a single-degree-of-freedom oscillator in a fluid flow.
. .
684
the drag coefficient and D is the width of the structure. The non-linear drag force acts in direction of the relative flow velocity U − X . Introducing the non-dimensional parameters t = Vt,
q = X/D,
u0 = U/VD,
p = dq/dt,
where V = m0 v /(m0 + rAcI ), equation (1) can be written in the state vector form 2
2 0
dq = p, dt
dp du = −q − ap + b 0 + g=u0 − p=(u0 − p), dt dt
(2)
where a = 2z0 zm0 /(m0 + rAcI ), b = rAcm /(m0 + rAcI ) and g = DcD r/2v0 zm0 (m0 + rAcI ). The non-linear term =u0 − p=(u0 − p) has received special treatment in the literature based on approximate representation. For example, Dalzell [13] replaced this term by the polynomial expansion =u0 − p=(u0 − p) 1 s Ck k
(u0 − p)k , D0k − 2
k = 1, 3, 5, . . . ,
(3)
where the Ck are constant coefficients to be determined by using a least-squares fit and D0 q 0 is the maximum relative velocity defined by the range −D0 Q u0 − p Q D0 . Another treatment adopted by others (see, for example, [6]) is given by the representation =u0 − p=(u0 − p) 1 =u0 =u0 − 2=u0 =p + p 2 sgn (u0 ),
(=u0 = w =p=).
(4)
The two approximations (3) and (4) will be used in the present analysis for comparison with results obtained for the expression of the non-linear hydrodynamic drag g=u0 − p=(u0 − p). In the analysis of reference [6], the hydrodynamic drag coefficient g was taken to be of order ze, while in the present analysis it is taken to be of order e. The two types of ordering will yield different results when one deals with Ito stochastic calculus. The difference between the two is discussed in section 6. Furthermore, the expressions =u0 =u0 and =u0 = were replaced in reference [6] by Hermite polynomials which are similar to a certain extent to Dalzell’s expansion (3). 3. FIRST ORDER FILTER
The wave velocity function u0 (t) is assumed to be a narrow-band random process. It is possible to model this process as the response of a shaping filter to a wide-band random process. In this section, this process is represented as the output of the first order shaping filter du0 + du0 = V(t), dt
(5a)
where d is a positive filter parameter which governs the filter frequency bandwidth and V(t) is a scalar wide-band random process whose correlation time is much smaller than any characteristic period of the system oscillations. It is well known that the output of this filter can be approximated by the solution of the Ito stochastic differential equation (SDE) du0 = −du0 dt + zn dB(t),
(5b)
where V(t) has been replaced by the formal derivative of a standard Brownian motion B(t), i.e., V(t) = zn dB(t)/dt, and n is the intensity of the random process V(t).
685
Equation (5b) has a well-known stationary solution (Ornstein–Uhlenbeck process) with the following power spectral density Su 0 (v) and autocorrelation function Ru 0 (t): Su 0 (v) =
n 1 , 2p d 2 + v 2
Ru 0 (t) =
n −d=t= e . 2d
(6)
The function Su 0 (v) is plotted in Figure 2 (solid curve) for d = 0·6, n = 1. The expression of Su 0 (v) reveals that, for small values of d, the spectrum of the stationary solution of (5b) is mainly concentrated about a narrow frequency band close to the origin. Blevins [6] reported that most offshore structures have a fundamental natural frequency around 1 Hz, while the energy content of ocean waves is concentrated at frequencies of less than 0·2 Hz. The filter output stationary probability density function f(u0 ) and the stationary mean of the absolute value of =u0 = are, respectively, f(u0 ) =
0 1
1 u2 exp − 0 , 2D1 z2pD1
E[=u0 =] = z2D1 /p,
(7)
where D1 = n/(2d). Equations (2) and (5b) constitute a set of Ito SDE with the corresponding initial conditions. The statistical analysis of these equations will be carried out using the stochastic averaging of the response energy envelope. 4. STOCHASTIC AVERAGING OF THE RESPONSE ENERGY
The stochastic averaging method of the energy envelope was originally introduced by Stratonovich [14]. We will use the generalization of this method [15] (see also Appendix A). The parameters a, g, d 2 and b 2 are assumed small and of the same order. Introducing the polar co-ordinates I and 8 defined by the transformation q = I cos 8,
p = −I sin 8,
8 = c + t,
(8)
Figure 2. Power spectra of different shaping filters: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; — – —, fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, qk = 1·5, n = 2/17.
. .
686
equations (2) can be rewritten in the standard form dI = {−aI sin2 8 + bdu0 sin 8 − g sin 8=u0 + I sin 8=(u0 + I sin 8) + (2I)−1b 2n cos2 8} dt − bzn sin 8 dB(t), dc = {−a sin 8 cos 8 + I−1bdu0 cos 8 − gI−1 cos 8=u0 + I sin 8=(u0 + I sin 8) − (I)−2b 2n sin 8 cos 8} dt − I−1bzn cos 8 dB(t),
(9)
where the Ito correction term is included in these equations. It should be noticed that the state co-ordinates I and c are slowly varying with time t, while u0 and 8 exhibit fast variation. In this case, the averaging of the drift vector and diffusion matrix of equations (9) will be taken with respect to 8 and with respect to the stationary distribution of the ergodic process u0 (t), keeping I and c frozen. This process gives the averaged equations (see Appendix A). dI = {−12 aI + c/4I − gI 2x¯ (I)} dt + zc/2 dB1 (t), dc = (1/I)zc/2 dB2 (t),
(10)
a x(x)f(u0 ) du0 , x = u0 /I, c = b 2n and B1 (t) and B2 (t) are two independent where x¯ (I) = f−a Brownian motion processes, each of unit intensity. The function x is some deterministic function of the new wave velocity parameter x evaluated for two different domains (see Appendix B).
x(x) =
1 2p
g
2p
sin 8=x + sin 8=(x + sin 8) d8,
(11a)
0
or
F g 1 (4 cos u + 4u sin u − 4 cos3 u), 3 2p G f
x(x) = G ,
for =x= e 1 (11b) for =x= E 1,
where (x = sin u).
The function x(x) is plotted in Figure 3(a). It satisfies the two inequalities x(x) e =x=
and x e 4/3p.
(12)
The function x(x) may also be written in the form (see also Figure 3(b)) x(x) = =x= + h(x),
Figure 3. The dependence of (a) the function x(x) on x = u0 /I, and (b) the function h(x) on x = u0 /I.
687
where h(x) = (2/3p)[3z1 − x + 3x sin x − (1 − x ) ] − x, for 1 e x e 0, and h(x) = h(−x), h(x) = 0 for =x= q 1. The component h is only significant over the range 1 e x e 0, as shown in Figure 3(b). Statements (11) are accurate and provide better representation of the non-linear Morison’s drag term in the averaging process. They facilitate the analysis of the averaged system. To this end let us introduce the Hamiltonian H = I 2/2. In terms of H, the first equation (10) takes the form which includes the Ito correction term after the transformation. 2
−1
2 3/2
$
%
c dH = −(a + 2gE=u0 =)H + − 2gHz2Hh¯ (H) dt + zcH dB1 , 2 or
$
%
c dH = −aH + − U(H) dt + zcH dB1, 2
(13a)
a = a + 2gE[=u0 =] = a + 2gzn/dp,
(13b)
g
(13c)
where
U(H) = 2gHz2Hh¯ (H) = 2gHz2H
a
h(u0 /z2H)f(u0 ) du0 = 8gH 2y(H),
−a
y(H) =
g
1
h(x)f(xz2H) dx.
(13d)
0
Notice that the integration limits are changed from 2a to 21 (or from 0 to +1 due to symmetry) because h(x) is non-zero only over the domain −1 Q x Q +1, and both h(x) and f(xz2H) are even functions with respect to x. The phase angle c(t) only constitutes a diffusion component, as the drift coefficient vanishes. It follows from relations (11) and (12) that the non-linear term g=u0 − p=(u0 − p) in the expression for the drag force ‘‘in average’’ has a dissipative nature (−x Q 0). The stationary probability density of the response process H of system (13) can be derived by using the Fokker–Planck equation in an explicit form:
$ $
exp[−2(a/c)H] exp −(16g/c) f(H) =
g6 a
% %7
hn(h) dh
0
exp[−2(a/c)H] exp −(16g/c)
0
g g
H
H
hn(h) dh
0
.
(14)
dH
The following inequality may be established based on equation (13) after using equation (12):
$
%
0
1
c c dt + zcH dB1 , dH = −aH + − U(H) dt + zcH dB1 E −aH + 2 2
(15)
where the right side does not include the contribution of h(H). The corresponding level of H may be governed by the differential equation (16) dH' = (−aH' + c/2) dt + zcH' dB1 . The corresponding energy level forms an ‘‘upper bound’’ for the original averaged equation (13). It is clear from equations (16) and (13a) that the non-linear hydrodynamic drag force enhances the damping coefficient. Equation (16) will be referred to as the ‘‘simplified’’ model.
. .
688
The function h(x) (or the drift term U(H)) indicates the difference between Morison’s non-linear term and its approximate representation (4) as adopted in reference [6]. It is worth noting that if representation (4) is used it will lead to the approximate averaged equation (16). Later we will discuss the conditions under which this difference becomes considerable. In fact, the simplified model is the result of some equivalent linearization, when Morison’s non-linear drag term g=u0 − p=(u0 − p) is replaced by −2gE[=u0 =]p, as shown by Spanos and Chen [4] for the deterministic case. The effect of h(x) on the system is negligible if =x= e 1 and it becomes significant if =x= Q 1. For the simplified model the response energy level is governed by the Gibbs distribution, i.e., f (H') =
0
1
2a 2a exp − H' . c c
(17)
The corresponding density for the amplitude I = zq 2 + p 2 is known to follow Rayleigh’s distribution f(I) = (2aI/c) exp(−aI 2/c). If the original co-ordinates are used, the response will be governed by a Gaussian distribution. Equation (16) can now be analyzed for the statistical behavior of the response energy function H'. Instead of solving for the probability density function, it is convenient to establish the corresponding equations for the n-dimensional characteristic functions Fn (l1 , . . . , ln ; t1 , . . . , tn ), where the li are arbitrary dummy real parameters. These characteristic functions are represented by a set of linear partial differential equations of the first order and may be integrated in a closed form as outlined in Appendix C. The stationary characteristic functions of system (16) for n = 1, 2 are obtained in the closed form: 1 F1 (l) = , F2 (l1 , l2 ; t2 − t1 ) = 1/{1 − id(l1 + l2 ) − l1 l2 d 2[1 − e−a(t2 − t2 )]}, (18) 1 − idl where d = c/(2a). The autocorrelation function RH' (t) and the power spectral density SH' (v) associated with the stationary process H'(t) are obtained from the solution F2 , i.e., RH' (t) =
b
1 2F 2 = d 2{e−a=t= + 1}, 1(il1 ) 1(il2 ) l = 0
SH' (v) =
i
d2 a . p a2 + v2
(19)
Expressions (19) are of the same type as those of the excitation u0 (t). The stationary mean value of H' can be obtained from the first relation (18) and is found to be c/(2a). It is interesting to notice that, for relatively large values of n, the stationary mean value is proportional to zn, i.e., E[H'] 2 (b 2/4g)zpdn . One may check that, for the one-dimensional equation (16), the border H = 0 is not accessible for the process H(t), and any perturbation will bring the system energy back to its original value H(t) (see Appendix D). Equally important is to note that the existence of the stationary solution of the SDEs (2) and (5b) (for the general case) can be established by applying the generator L( · ) (see, for example, reference [16]): L( · ) = f(X, t)
$
%
7
1( · ) 1 1 1T + 2 Tr ( · ) [g(X, t)ngT(X, t)] , 1X 1X 1X
(20)
on any response function, such as the Lyapunov function, based on the general Ito stochastic differential equation dX = f(X, t) dt + g(X, t) dB(t), (21) where X is a state vector, e.g., for the present system X = {q, p, u0 }T, n is the covariance matrix of the vector dB. For simplicity we consider the case when the co-ordinate q is
689
bounded. It is convenient to shift the co-ordinate p by z = p − u0 . Accordingly, the governing equations of the system and filter take the form dq = {z + u0 } dt, dz = {−q − az − g=z=z − r¯ u0 } dt + (b − 1)zn dB(t),
(22)
du0 = −du0 dt + zn dB(t), where r¯ = a + bd − d. Introducing the response function u = (A1 z 2 + A2 u02 )/2 as a Lyapunov function (where A1 and A2 are positive constants), and applying equation (20) gives
0
1 0
1
r¯ r¯ 2A1 2 n n u0 − A 2 d − u0 + A1 (b − 1)2 + A2 . 2a 2 2 4a 2
Lu = −A1 zq − gA1 z 2=z= − aA1 z +
(23) One of the main conditions for the existence of a stationary solution (see reference [16], p. 90) is that the supremum of the right side of equation (23) decreases as both =u0 = and =z= increase, i.e., r¯ 2 sup Lu(z, u0 ):−a as R:a, if A2 d q A1 . (24) 2 2 4a z +u qR 0 The result of equation (24) indicates that the stationary solution exists and forms a stable and unique attractor. Note that the approximation given by relation (3) will lead to the averaged equation dH = {−(a + 2gc˜1 )H + b 2n/2 − g[c˜3 (2H)2 + c˜5 (2H)3 + · · ·]} dt + bznH dB1 .
(25)
Each of the coefficients c˜i represents some series, so one should keep sufficiently many terms in order to achieve high accuracy. Notice that the difference between equation (25) and the exact equation (13) is in the expressions of the drift coefficient. 5. SECOND ORDER SHAPING FILTER
5.1. - vf $ 1 The first order shaping filter characterizes the wave velocity whose peak spectrum occurs at zero frequency and its bandwidth is governed by the system parameters. However, if the peak spectrum is shifted from zero frequency, as in the case of any real ocean system, one should use a higher order filter. In this section, we consider the case when the wave velocity u0 (t) is represented by the output of the shaping filter du0 = p0 dt,
dp0 = {−vf2 u0 − 2vf zp0 } dt + zn dB(t),
(26)
where vf and z are the natural frequency and damping factor of the filter, respectively. Filter (26) has a well-known stationary solution, with the following power spectral density and autocorrelation function, respectively: Su 0 (v) =
Ru 0 (t) = where v02 = vf2 (1 − z 2).
0
n 1 , 2p (vf2 − v 2)2 + (2vvf z)2
1
n vz cos v0 t + f sin v0 =t= exp(−zvf =t=), 4vf3 z v0
(27)
. .
690
The function Su 0 (v) is plotted by the dashed curve in Figure 2 for z = 0·4, vf = 0·5 and n = 0·1. The output stationary probability density function for the filter (26) is f(u0 , p0 ) =
6 0
1 p2 u2 exp −12 0 + 0 D1 D2 2pzD1 D2
17
,
where D2 = n/4vf z and D1 = D2 /vf2 = n/4vf3 z. In order to derive the response statistics of systems (2) and (26), the parameters a, g, z 2 and b 2 with be assumed to be of the same order of smallness. In this case, these equations must be transformed into a form suitable for averaging, in the non-resonance case (i.e., vf $ 1). Introducing the transformation 1 0 K u0L K 0 1 G p0G=G 0 b/(1 − vf2 ) G qG G 0 k p l k−bvf2 /(1 − vf2 )
0 0 1 0
0LKu0 L 0GGp0 G 0GGq1 G 1lkp1 l
(28)
into equations (2) gives dq1 = {p1 + 2vf zb p0 } dt − b zn dB, dp1 = {−q1 − ap1 + ab vf2 u0 + g=d u0 − p1 =(d u0 − p1 ) dt,
(29)
where b = b/(1 − vf2 ), d = [1 + (b − 1)vf2 ]/(1 − vf2 ), vf $ 1, and the filter state equations (26) remain the same. Transformation (28) eliminates the linear coupling from equation (2), where the unperturbed conservative parts of equations (2) and (26) take a Hamilton form. Introducing the slow variables I and c, defined by the co-ordinate transformation q1 = I cos 8,
p1 = −I sin 8,
8 = c + t,
(30)
equations (29) can be rewritten in the standard form dI={−aI sin2 8 − avf2 b u0 sin 8 + 2zvf b p0 cos 8 − g sin 8=d u0 + I sin 8=(d u0 + I sin 8) + (2I)−1b 2n sin2 8} dt − b cos 8zn dB(t), dc = {−2I−1zvf b p0 sin 8 − a sin 8 cos 8 − I−1ab vf2 u0 cos 8 − gI−1 cos 8=d u0 + I sin 8=(d u0 + I sin 8) + (I)−2b 2nsin 8 cos 8} dt + I−1b sin 8zn dB(t). (31) Averaging gives a similar set of averaged SDE (10), where c = b 2n, x(I) =
gg
a
x(d u0 I−1 )f(u0 , p0 ) du0 dp0 .
(32)
−a
The equation for the energy, H = I 2/2, takes the form of equation (13), with the following expression for the parameter a: a = a + 2gE[=d u0 =] = a + 2gd zn/2pzvf3 ,
(33a)
U(H) = 8gd −1H 2n(H),
(33b)
and function U(H):
691
where n(H) = (1/z2pD1 )
g
1
h(x) exp(−Hx 2/d 2D1 ) dx.
(33c)
0
The phase angle c(t) only constitutes a diffusion component as the drift coefficient vanishes. It follows from relations (11) and (12) that the non-linear term g=u0 − p=(u0 − p) in the drag force expression has a dissipative nature (−x Q 0) and reduces the stationary expected value of H. The stationary probability density of the response process of system (13) can be derived by using the Fokker–Planck equation in an explicit form:
$ $
g g
exp(−2aH/c) exp −(16g/d c) f(H) =
g6 a
% %7
hn(h) dh
0
H
exp(−2aH/c) exp −(16g/d c)
0
H
hn(h) dh
0
.
(34)
dH
This expression is identical to the probability density function obtained for the first order filter case. However, the two results are affected by the filter constants d, d, vf and z, which are included in the parameters a and c. Notice also that the function n(h) is different in both cases as given by relations (13c) and (33c). Following the analysis of the previous section, it is not difficult to show that the non-linear drag force on average enhances the damping and we also may establish the inequality (15). Again, equality (16) establishes the upper bound of the response energy level of the structure and one may determine all finite-dimensional non-stationary characteristic functions of the energy H, as outlined in Appendix C. In conclusion, it is possible to write the expression for the mean value of the energy in terms of the initial variables q, p: E[12 (q 2 + p 2)] = E[H] + 12 b 2(1 + vf2 )D2 . 5.2. : vf 1 1 Consider the case of the shaping filter (26) when the frequency vf is in the neighborhood of 1. Therefore, the detuning (vf − 1) is a small value. Introducing the linear transformation (which eliminates the linear coupling from the system and filter equations when g = z = 0, n = 0)
K qL K 1 G p G=G −a/2 G u0G G 0 k p0l k 0
0 V 0 0
a1 −a2 vf2 1 0
a2LK q1 L a1GG p1 G 0GG u0 G 1lk p0 l
(35)
into equations (2) gives dq1 = {−aq1 /2 + Vp1 + 2za2 p0 } dt − (a2 /vf )zn dB, dp1 = {−Vq1 − ap1 /2 + 2za3 vf p0 + g/V=aq1 /2 − Vp1 + Q=(aq1 /2 − Vp1 + Q)} dt − a3 zn dB,
(36)
where V = z1 − a 2/4, a1 = D−1abvf2 , a2 = D−1b(1 − vf2 ), D = (1 − vf2 )2 + a 2vf2 , a3 = ab(1 + vf2 )/2DV and Q = (1 + a2 vf2 )u0 − a1 p0 . It is assumed that a, g, vf − 1, z 2 and b 2/3 are of the same order of smallness. Introducing the slow variables I, and c, defined by q = I cos 8,
p = −I sin 8,
8 = c + Vt,
. .
692
and averaging, the following averaged SDE are obtained (see Appendix A): x(I)} dt + zc/2 dB1 (t) dI = {−12 aI + c/4I − gI 2 dc = −(ga/2V)I 2 x(I) dt + (1/I)zc/2 dB2 (t)
(37)
a x(QI−1 )f(u0 , p0 ) du0 dp0 , c = n(a22 /vf2 + a32 ), and dB1 (t) and dB2 (t) are where x(I) = ff−a two independent Brownian motion processes, each of unit intensity. In terms of H, equations (37) take the form
6
c dH = −aH + − 2gHz2Hh¯ (H)} dt + zcH dB1 (t), 2 dc = {−(ga/V)(zH/2E[=Q=] + Hh¯ (H))} dt + 12 zc/H dB2 (t),
(38)
where a = a + 2gE=Q=, E=Q= = z2DQ /p, h¯ (H) = Eh(Q/z2H) and DQ is the covariance of the variable Q which can be derived from the filter equations. It is seen from equation (38) that the phase angle c(t) constitutes drift and diffusion components in contrast to the non-resonance case. The stationary probability density of the response process H can now be derived by using the Fokker–Planck equation in an explicit form. It has the form (14), where the expression for the n(H), given by equation (13c) as well as the probability density (7), are of the same form, and the co-ordinate u0 should be replaced by Q (which is normally distributed with zero mean and covariance DQ ). 6. FOURTH ORDER SHAPING FILTER FOR THE WAVE EXCITATION
For better fitting of the known Pierson–Moskowitz spectrum [17], following Spanos [18], a fourth-order shaping filter may be considered: d3u0 d2u0 du0 d3V d4u0 + a0 u0 = 3 , 4 + a3 3 + a2 2 + a1 dt dt dt dt dt
(39)
where a3 = c + c¯ , a2 = k + k + cc¯ , a1 = ck + c¯k, a0 = kk , c, c¯ , k and k are some constants satisfying the stability condition a1 (a2 a3 − a1 ) − a0 a32 q 0.
(40)
If condition (40) holds then Ito’s SDE (39) has a stable stationary solution with the power spectral density Su 0 (v) =
nv 4 . [(v − k) + (cv) )][(v 2 − k)2 + (c¯v)2] 2
2
2
(41)
The parameters a1 and a3 in equation (39) are damping coefficients. When a1 = a3 = 0 and V = 0, the system (39) is conservative with the two natural frequencies v1 = za2 /2 − za22 /4 − a0 ,
v2 = za2/2 + za22 /4 − a0 .
Figure 2 includes two typical spectra plotted according to relation (41) for two different constant coefficient sets. The dotted curve is plotted for the case of c = 0·17, c¯ = 0·476, k = 0·147 and k = 0·375, and for a white excitation level n = 1/130. The spectrum shown by short and long dashes is plotted for the case of c = 0·338, c¯ = 0·946, k = 0·586 and k = 1·5, and for the white excitation level n = 2/17. It is seen that this spectrum is closer to the resonance frequency than the previous one.
693
Equation (39) can be written as a set of four linear SDE: dx1 = x2 dt + zn dB,
dx2 = x3 dt − a3 zn dB,
dx3 = x4 dt + (a32 − a2 )zn dB, (42)
dx4 = {−a0 x1 − a1 x2 − a2 x3 − a3 x4 } dt + (−a33 + 2a2 a3 − a1 )zn dB, where x1 (t) = u0 (t). For further analysis, it is convenient to apply the following linear transformation of the variables xj which eliminates the linear coupling from equation (42) when a1 = a3 = 0, V = 0:
Ky1L Kv22 0 Gy2G = G 0 v22 Gy3G Gv12 0 ky4l k 0 v12
1 0 1 0
0 LKx1L 1 GGx2G , 0 GGx3G 1 lkx4l
provided that v1 $ v2 .
(43)
Now the system and filter equations take the form dq = p dt, dp = {−q − ap + bD−1(dy1 /dt − dy3 /dt) + g=D−1(y1 − y3 ) − p=(D−1(y1 − y3 ) − p)} dt, dy1 = y2 dt + (a32 − v12 )zn dB, dy2 = {−v12 y1 − h1 y2 + h2 y4 } dt − g2 zn dB, dy3 = y4 dt + (a32 − v22 )zn dB, dy4 = {−v22 y3 − h1 y2 + h2 y4 } dt − g1 zn dB,
(44)
where D = v22 − v12 , h1 = D−1(a1 − a3 v12 ), h2 = D−1(a1 − a3 v22 ), g1 = a1 + a33 + v12 a3 − 2a2 a3 , and g2 = a1 + a33 + v22 a3 − 2a2 a3 . The linear transformation of the variables q and p (which is valid in the non-resonance case v1 $ 1, v2 $ 1), q1 = q − b1 y2 + b2 y4 , b1 = b/D(1 − v12 ),
p1 = p + b1 v12 y1 − b2 v22 y3 , b2 = b/D(1 − v22 ),
(45)
eliminates the coupling from equation (44) when the damping vanishes (i.e., a1 = a3 = a = 0) and V = 0. In this case the governing equations of q1 and p1 are dq1 = {p1 − d1 y2 + d2 y4 } dt + e1 zn dB, dp1 = {−q1 − ap1 + ab1 v12 y1 − ab2 v22 y3 + g=D−1(y1 − y3 ) − p1 + b1 v12 y1 − b2 v22 y3 = × (D−1(y1 − y3 ) − p1 + b1 v12 y1 − b2 v22 y3 )} dt + e2 zn dB,
(46)
where e1 = −b2 g1 , e2 = b + b1 v12 (a32 − v12 ) − b2 v22 (a32 − v22 ), d1 = bh1 /(1 − v12 )(1 − v22 ) and d2 = bh2 /(1 − v12 )(1 − v22 ). Introducing the slow variables I and c defined by equation (30) and averaging, we obtain equations (10) in which c = n(e12 + e22 ), d1 =
x¯ (I) = Ex(QI−1 ),
1 + (b − 1)v12 , D(1 − v12 )
d2 =
Q = d1 y1 − d2 y3 , 1 + (b − 1)v22 . D(1 − v22 )
(47)
694
. .
In terms of H the first equation takes the form (13) by setting a = a + 2gE[=Q=], and U(H) is given by equation (13b), where q0 must be replaced by Q. The mean square of Q can be estimated from the last four equations (44) by using, for example, the moments method. The stationary probability density function f(Q) is governed by equation (7) where D1 must be replaced by the covariance of Q. The analysis will lead to results similar to those of equations (13) and (14).
7. DISCUSSION OF RESULTS AND CONCLUSIONS
The analytical results derived in the previous two sections provide complete information regarding the statistics of a linear single-degree-of-freedom system subjected to a non-linear random hydrodynamic force generated by three different filters. With reference to the first order shaping filter, the spectrum is governed mainly by the parameter d and becomes almost flat for large values of d. For the case of the second order filter, the filter damping coefficient governs the wave bandwidth. As the filter damping ratio decreases, the band width of the wave velocity becomes more narrow. The fourth order filter emulates to a certain extent the Pierson–Moskowitz spectrum which represents the actual spectrum of ocean waves. The spectra of these three filters are shown in Figure 1. These filters are subjected to a wide-band random excitation with different power spectrum levels. It is seen that as the order of the shaping filter increases it requires a different level of the excitation, depending on its physical parameters. The probability density functions of the system response to the output of these three filters are obtained in a closed form and are plotted in Figure 4 for the system parameters a = 0·01, b = 0·25 and g = 0·01. The characteristics of the density curves are primarily
Figure 4. Probability density functions of the system response energy to the output of three filters for system parameters a = 0·01, b = 0·25 and g = 0·01: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
695
governed by the spectral density level of the white noise input to each filter and how close the peak of the spectrum is to the system’s natural frequency. In general, the system response is non-Gaussian due to the non-linear hydrodynamic drag. If the approximation of Eatock Taylor and Rajagopalan as given by relation (4) is used, the present analysis will give Gaussian probability density for every shaping filter. However, according to the stochastic averaging adopted by Eatock Taylor and Rajogopalan [6], the response density is not Gaussian because of two main factors. The first is that they replaced the hydrodynamic drag force by a Hermite polynomial and the second is that they considered g to be of the order ze, which results in both drift and diffusion terms in the Ito equation of the governing equations of motion. In the present analysis, the coefficient g is taken to be of the order e in order to take into account the complete effect of non-linear hydrodynamic drag without approximation. The advantage of this is that the diffusion coefficient due to the non-linear drag term will vanish and this avoids several analytical complications. It is found that the probability density curves, obtained by using the approximation (4) and by using the exact drag expression g=u0 − p=(u0 − p), are close for small values of g and the difference becomes larger for larger values of g, as shown in Figure 5. This difference is a measure of the deviation of the response from normality. It is seen that the first order filter (solid curve), whose center frequency is zero and shifted away from the system natural frequency, requires a very high level of filter excitation spectrum, i.e., n = 1. In view of this large level of n, the probability density of the system response spreads over a relatively higher response energy level. For the second order filter (dashed curve), the probability density function spreads over higher response levels as the filter natural frequency approaches the system natural frequency. Fourth order filters (dotted and dash-dot curves) yield different probability density curves depending on the filter frequency
Figure 5. Probability density functions of the system response energy to the output of a second order filter for z = 0·4, vf = 0·5, n = 0·2, a = 0·01, b = 0·25 and two different values of the non-linear drag coefficient g: ----, using exact drag expression g(u0 − p)=u0 − p=; – – – –, using the approximate drag expression given by equation (4).
696
. .
Figure 6. The dependence of the expected value of the response energy on the white noise excitation level for system parameters a = 0·01, b = 0·25 and g = 0·02: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
and the value of n. It is seen that for small n and lower center frequency, the response probability density has a peak which is concentrated around the zero response level. On the other hand, as both the filter center frequency and excitation level n increase, the probability of higher response level increases.
Figure 7. The dependence of the expected value of the response energy on the non-linear hydrodynamic drag coefficient for system parameters a = 0·01, b = 0·25 and g = 0·01: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
697
Figure 8. The normalized mean value of the system response energy (with respect to the expected value of the linear system response energy) versus the white noise excitation level for system parameters a = 0·01, b = 0·25 and g = 0·01: ----, first order shaping filter, d = 0·6, n = 1; – – – –, second order shaping filter, z = 0·4, vf = 0·5, n = 0·1; · · · ·, fourth order shaping filter, c = 0·17, c¯ = 0·476, k = 0·147, k = 0·375, n = 1/130; — – —, fourth order shaping filter, c = 0·338, c¯ = 0·946, k = 0·586, k = 1·5, n = 2/17.
The dependence of the expected value of the response energy level E[H] on the excitation level n is shown in Figure 6 for a non-linear hydrodynamic drag coefficient g = 0·02. For the three filters, the response mean value increases monotonically with n in a non-linear
Figure 9. The dependence of the mean value of the system response energy on the filter center frequency for the case of a second-order shaping filter, in the neighborhood of resonance for system parameters a = 0·01, b = 0·25 and g = 0·01.
698
. .
fashion and is of the form zn. The dependence of E[H] on the non-linear hydrodynamic drag coefficient g is shown in Figure 7 for the three filters. It is seen that the non-linear drag has a remarkable effect on suppressing the system response. The normalized response energy level with respect to the corresponding linear value is plotted as a function of the excitation level n in Figure 8. The deviation of the normalized response from the value 1 measures the influence of the hydrodynamic drag. Again, comparison with results by Eatock Taylor and Rajagopalan reveals that their results are more conservative. Figure 9 shows the expected value of the response energy level for the case of the second order shaping filter in the neighborhood of resonance, vf = 1, for the system parameters a = 0·01, b = 0·25, g = 0·01. Unlike other cases, the response exhibits peak values only at resonance. The analysis of this paper is based on three different shaping filters. The Pierson– Moskowitz (PM) spectrum [17] is more realistic and Part II of this investigation will consider the PM spectrum and a shaping filter which takes into account the non-normality of sea waves. Furthermore, the system inherent non-linearities have to be included in the future analysis. These factors are important in order to examine important issues such as level crossings, local extremes of the response and other related reliability problems. ACKNOWLEDGMENTS
This research is supported by a grant from Office of Naval Research under grant No. N000149310936, Dr Julia Abrahams is the Scientific Officer. REFERENCES 1. J. R. M, M. P. O’B, J. W. J and S. A. S 1950 Petroleum Transactions 189, 149–154. The forces exerted by surface waves on piles. 2. R. G. T 1977 The Structural Engineer 55(5), 209–222. Continuous random wave loading on structural members. 3. R. E T 1975 Off-shore Structures. 125–132. London: ICE. Structural dynamics of offshore platforms. 4. P.-T. D. S and T. W. C 1980 International Journal of Non-linear Mechanics 15, 115–126. Response of a dynamical system to flow-induced load. 5. A. W. L 1986 Transaction of the American Society of Civil Engineers, Journal of Structural Engineering 112, 2416–2429. Nonlinear structural response in random waves. 6. R. E T and A. R 1982 Journal of Sound and Vibration 83, 401–431. Dynamics of offshore structures: parts 1 and 2. 7. S. K. T and J. M. N 1992 Transaction of the American Society of Civil Engineers, Journal of Engineering Mechanics 118, 942–960. Parametric and external excitation of marine risers. 8. M. S, C. Y and R. V 1977 Transaction of the American Society of Civil Engineers, Journal of Structural Mechanics 5(2), 135–146. Dynamic analysis of fixed offshore structures subjected to wind generated waves. 9. P.-T. D. S and J. E. H 1981 Transaction of the American Society of Mechanical Engineers, Journal of Energy Resources Technology 103, 43–249. Linear prediction theory for digital simulation of sea waves. 10. N. K. L and W. A. H 1984 Transaction of the American Society of Mechanical Engineers, Journal of Energy Resources Technology 106, 466–472. Time series simulations of wide-band spectra for fatigue tests of offshore structures. 11. P. B, I. L, S. R. W, T. M and D. K 1992 in Nonlinear Stochastic Mechanics, 83–101. Berlin: Springer. Application of nonlinear stochastic mechanics in offshore engineering. 12. R. D. B 1977 Flow-induced Vibration. New York: Van Nostrand Reinhold. 13. J. F. D 1978 Journal of Ship Research 22(3), 178–185. A note on the form of ship roll damping.
699
14. R. L. S 1963 Topics in the Theory of Random Noise. Vol. 1. New York: Gordon & Breach. 15. R. Z. K 1968 Cybernetics (Kybernetika Cislo 3, Rocnik) 3, 260–279. On averaging principle for Ito stochastic differential equations (in Russian). 16. R. Z. K 1980 Stochastic Stability of Differential Equations. Alphen: Sijthoff and Noordhoff. 17. W. J. P and L. M 1964 Journal of Geophysical Research 69, 5181–5190. A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. 18. P.-T. D. S 1983 Transation of the American Society of Mechanical Engineers, Journal of Energy Resources Technology 105, 300–309. ARMA algorithms for ocean wave modeling. 19. V. S. P and I. N. S 1987 Stochastic Differential Systems. New York: John Wiley. 20. W. E. W 1980 Partial Differential Equations. Oxford: Clarendon Press. 21. Y. V. P and Y. A. R 1969 Probability Theory: Basic Concepts, Limit Theorems, Random Processes. New York: Springer-Verlag. APPENDIX A: STATEMENT OF AVERAGING PRINCIPLE [15]
Let {X (t), Y (e)(t)}T be a Markov vector in El with X e(t) giving as a slow random process, while Y e(t) experiences fast variation. The two processes may be described by the following system of Ito SDE: (e)
dX (e)(t) = eA1 (X, Y) dt + zes(X, Y) dB(t),
dY (e)(t) = A2 (X, Y) dt + 8(X, Y) dB(t),
with dim X = l1 , dim Y = l2 and dim B = l3 , and subject to the initial conditions X (e)(0) = x0 ,
Y (e)(0) = y0 .
Here, A1 (x, y) and A2 (x, y) are given functions of the dimensions l1 and l2 , respectively, and s(x, y) and f(x, y) are given matrix functions of the dimension l1 × l3 and l2 × l3 , respectively. The elements of the vector B are independent standard Brownian processes. Assume that the vector functions A1 , A2 , s and f satisfy the Lipschitz condition with respect to x and y. Since the process X is slowly varying with respect to time, one may establish the Ito SDE for the process Y (x,y)(t) by freezing the process X: dY (x,y)(t) = A2 (x, Y (x,y) ) dt + 8(x, Y (x,y) ) dB(t),
Y (x,y)(0) = y.
If the following limits exist A 1 (x) = lim E t:a
$g
D (x) = lim E t:a
1 t
$g
t+t
1 t
t+t
%
A1 (x, Y (x,y)(s)) ds ,
t
%
s(x, Y (x,y)(s))s T(x, Y (x,y)(s)) ds ,
t
then under some assumptions [15] the process X (e)(t) weakly converges, on the time interval 0 E t E T/e as e:0, to the Markov process X (0)(t) governed by the Ito SDE dX (0)(t) = eA 1 (X (0) ) dt + zes¯ (X (0) ) dB'(t),
s¯ s¯ T = D .
APPENDIX B: EVALUATING FUNCTION x (11a)
Consider the function x(x) =
1 2p
g
2p
sin 8=x + sin 8=(x + sin 8) d8.
(B1)
0
We split the x-axis into the four domains: (1) x E −1,
(2) −1 E x E 0,
(3) 0 E x E 1
and
(4) x e 1,
. .
700
and evaluate integral (B.1) for each domain as follows: x(x) = −
(1)
1 2p
g
2p
sin 8(x + sin 8)2 d8 = −x.
0
(2) Introducing auxiliary variable u $ [0, p/2], as sin u = −x, we have x(x) =
=
1 2p
$g
p−u
sin 8(x + sin 8)2 d8 −
u
g
2p + u
%
sin 8(x + sin 8)2 d8
p−u
1 {4 cos u + 4 sin u − 43 cos3 u}. 2p
(3) Introducing auxiliary variable u $ [0, p/2], as sin u = x, we have x(x) =
= (4)
1 2p
$g
p+u
sin 8(x + sin 8)2 df −
g
2p − u
%
sin 8(x + sin 8)2 d8
p+u
−u
1 {4 cos u + 4u sin u − 43 cos3 u}. 2p x(x) =
1 2p
g
2p
sin 8(x + sin 8)2 d8 = x.
0
APPENDIX C: EVOLUTION OF CHARACTERISTIC FUNCTION OF ONE-DIMENSIONAL SYSTEM
Consider the one-dimensional Ito SDE (C.1) dX(t) = f(X, t) dt + g(X, t)zn dB, X(t0 ) = x0 , where B(t) is a vector of the independent standard Brownian motions; f(x, t) and g(x, t) are given functions of the x $ R 1. For definiteness we consider the case when both f(x, t) and g(x, t) are given as polynomials in x. The evolution of the first order characteristic function F1 (l; t) = E[eilX(t)] of the random process X(t) at time t, where l is a dummy arbitrary real parameter, is described by the partial differential equation
6 0 1
0 1 0 17
1 1 1 1F1 (l; t) = ilf , t − 12 lg , t ng , t l F1 (l; t). 1t i1l i1l i1l
(C.2)
The second and nth orders characteristic functions of X(t) are, respectively, F2 (l1 , l2 ; t1 , t2 ) = E[eil1 X(t1 ) + il2 X(t2 )], Fn (l1 , . . . , ln ; t1 , . . . , tn ) = E[eil1 X(t1 ) + ··· + iln X(tn )]. The partial differential equation which describes the evolution of Fn is [19]
6 0
1
0
10
1 7
1 1 1 1Fn (l1 , . . . , ln ; t1 , . . . , tn ) = ilnT f , t − 12 lnT g , t ng ,t l F . 1tn i1ln n i1ln n i1ln n n n T
(C.3) It is not difficult to show that the differential equation for the characteristic function of the system described by equation (16) is 1F c 1Fn = ln (−a + 12 icnln ) n + i ln Fn , 2 1tn 1ln
(C.4)
701
with appropriate initial conditions F1 (l1 ; t0 ) = F0 (l1 ), Fn (l1 , . . . , ln ; t1 , . . . , tn − 1 , tn − 1 ) = Fn − 1 (l1 , . . . , ln − 2 , ln + 1 + ln ; t1 , . . . , tn − 1 ). The solution of equation (C.4) is obtained by using the method of characteristics (see, for example, reference [20]) and is given in the form F1 (l1 ; t1 ) = F0 (l1 C−1(l1 , t − t0 ) e−a(t − t0 ) )C−1(l1 , t − t0 ), Fn = Fn − 1 (l1 , . . . , ln − 2 , ln − 1 + ln C−1(ln , tn − tn − 1 ) e−a(tn − tn − 1 ); t1 , . . . , tn − 1 ) × C−1(ln , tn − tn − 1 ), where C(l, t) = 1 − (ilc/2a) [1 − exp(−at)], and F0 is estimated from the initial distribution of the response. For example, if the initial distribution is given by the Gibbs function p0 (x) = d−1 e−x/d , the corresponding characteristic function is F0 (l) = 1/(1 − idl). APPENDIX D: NON-ACCESSIBILITY OF THE BORDER FOR THE ENERGY-AVERAGED EQUATION (16)
Consider the one-dimensional Ito stochastic differential equation (C.1) on the half-line [0, a) with continuous drift and diffusion coefficients and unit intensity white noise. The left border x = 0 is accessible if and only if the function R(x) [21], R(x) = R1 (x)
g
x
0
$g
dy , g (y)R1 (y)
x
R1 (x) = exp −
2
0
%
dy , g 2(y)R1 (y)
is integrable in the neighborhood of zero. For equation (16) we have f = −ax + c/2,
g = zcxn.
In this case, R(x) =
0 1
1 a exp x , x c
and now it is clear that the function
0 1g >6 0 17
1 a exp x x c
x
dy
0
c exp
a y c
is not integrable in the neighborhood of zero, so for the process x(t), the left border x = 0 is not accessible.