Stochastic modelling of dynamic properties of nonlinear water waves

Stochastic modelling of dynamic properties of nonlinear water waves N.C.G.

Markatos*

Department of Chemical Engineering, Imperial College of Science and Technology, London SW7, UK (Received 29 October 1976; revised 10 March 1978)

This paper describes the stochastic modelling of water waves and presents predictions of their structure. The work is based on the analysis of more than 10~ data points, obtained in the form of wave records during experiments in a channeP. Both turbulent wind-driven waves and artificially generated waves superimposed on the former, were studied. The wave records were considered as random time series of events; and amplitudes, surface elevations and other characteristics were considered as stochastic processes, characterisable by their moments which were calculated from the wave records by means of a computer program. Despite the twodimensional character of the flow, the suggested distributions are not inconsistent with oceanic data, provided that the latter are treated as two-dimensional.

Introduction The results of the analysis of the experimental data are summarized below. More details and a complete analysis of a wave record may be found elsewhere 1. (1) The wave series is stationary, where stationarity is defined in terms of the number of events (waves) up to time t, {Nt}. (2) The wave process is not a Poisson process, i.e. it is not completely random. There are trends in the series of waves, which are functions of the time parameter t, and not functions of the serial number j of the wave. (3) There is a strong correlation between time intervals {X j}. Thus, the process is not a renewal process either. (4) Despite the lack of stationarity in the interval process {Xj}, the estimated spectra and correlograms indicate the existence of at least one cyclic trend. The same harmonic term was found at the same frequency, in the spectra of the counting process {N,}. Investigations of the underlying mechanism included * Present address: CHAM Ltd., Bakery House, 40 High Street, London SW19 5AU, UK

0307-904X/78/040227-12 $02.00 © 1978 IPC Business Press

the consideration of various models, e.g. semi-Markov models, Wold's Markov chains, pooled outputs, and simple series onto which a completely random series was super imposed. The best representation was achieved through a model generated by the superimposition of a regular series, B cos(co0t + 0o), and a completely random, stationary one, t/(t), which can be Gaussian.

Statistical distributions for surface elevation of a random non linear wavy interface In the linear theory of wind-generated water waves, the distribution of the surface elevation ~(t), is Gaussian. Analysis of the present channel flows proved that for air velocities greater than 0.75 m/s or for artificial waves superimposed on wind-driven waves, a considerable skewness exists so that the above distribution is no longer Gaussian. The Fourier transform of the characteristic function, in terms of the correlation functions, cannot be calculated directly, and we used the so-called quasimoment functions2; then, the n-dimensional probability densities, W., are given by:

Appl. Math. Modelling, 1978, Vol 2, December

227

Stochastic modelling of nonlinear water waves. N. C. G. Markatos Wn(C1 .....

~n'] t l .....

{£1

=

1 +

For n = 1 the polynomials (equation 5) become the ordinary Hermite polynomials, Hi, leading to:

t,,)

S--~.

b~(t . . . . . .

a..,to=l

S=3

t~,)

W(n)=

(1)

x w,,G (~, ..... ~,,)

where:

W,G (Ca,-.., ~,,) = (2rc)-"/2[IA I-

l +s~3~.-~.s[

((n) and a are the mean value and the standard deviation of n(t), respectively) which is simply an Edgeworth 3 series. Keeping cumulants, ki, up to the fourth, the a i r water elevations probability density becomes:

1/2]

W(n) = . 2/~_--~exp = 1 ~*"(~* - m , ) ( ~ - mm)

x exp - ~ , ,

(2) x

is the probability density of a Gaussian process with the same correlation functions as the non-Gaussian process, b~ are the quasi-moment functions, ~ are random variables, m~ are their mean values, [A[ is the determinant of the correlation matrix and ~t~mis its inverse, so that:

~ a4,mK2(tm, t~) = 6¢~ = S1 m

=

X

L

~

0

for for

tk=2 4) 4:2

~= =,.=1

(3)

_

H4(C) = C4 - 6C2 + 3,

10~3 + 15~

n6(~ )

=

~ 6 __ 1 5 C 4

n = CiCi + CuCiCi + CukCiCjCk + "'"

el,m : A - 1

~X-~) e -y[~I

-F " "

-t- 4 5 ~ 2 - - 15 formula derived for ocean waves by Longuet-Higgins 4 in a different way. In the second approximation, where the skewness (SK) is appreciable, but its square and the kurtosis (kurt) can be neglected, the mean value of the elevations n is proved to be shifted by an amount ( n ) = kl. If:

Then with the coefficient matrix ~t¢. we can associate the generalized Hermite polynomials defined as:

H,p...,,[x] = eY[=I(- ~ - ~ ) . . - ( -

n6(~)

(8)

H3(C) = C3 - 3~, Hs(~) = ¢5

where: Xm : Cm -- mm

[ 1 + ~SKH3(~) 1 , + ~g~r,n~(~)

where:

OldpraXdPxm

X , = CO -- m ,

Ca

+ -~SK

where 6,~a is the Dirac function and K2 are the correlation functions. Suppose we are given the quadratic form, in n variables X,:

y[x] = ~

W6(n) (7)

(4)

where the C's are constants and the C~ are independent random variables, symmetrically distributed about 0 with variance a 2, then we have (assuming that the C's are symmetric in their suffices): ( n ) = Cija 2 q- 3Clijjal2 ffj2

If we now write,

Longuet-Higgins 4 showed that kl = C,a2; therefore the shift is due to the quadratic terms involved. In this case the probability density is multiplied by the factor:

F, = F~,[x] = --aY[X] = ~ ~ , . x . we obtain: //,[x-I = &

(5) H4,m[x] = FeFm -- o:,,~ HC,,,,~[x] = feFmFv - ot~,,,fv -- c%vFm- o~moF, etc. so that equation (1) becomes: Wn(~l . . . . . =

C., tl,---,

{ 1 + s~=3 ~.w 1 . . . . .f~,

t.)

1 bs(t~,'",t,~)

H [~(t)-Kl(t)]} .....

x W~(¢, . . . . . ¢,,)

The analysis of the experimental results proved that only the first few terms in the above expansion are required for an adequate representation; thus, the calculation of the matrix ate., becomes relatively easy.

One-dimensional probability density and its various moments

228

Appl. Math. Modelling, 1978, Vol 2, December

(6)

and skewness is introduced. In the third approximation, the mean is unshifted, but a non-zero kurtosis appears given by: k4/k 2. The four first central moments, calculated through the defining integrals are 0, a 2, a 3 SK 3a 4 + K.,t °', respectively. The first four moments about the origin are similarly calculated. The above results are in good agreement with those of Kinsman s over short fetches. The present expression differs from Kinsman's by the term (2//0.2)- s/2{exp( --½ n2/tr2)} 1/72SK2H6 which becomes negligible, when SK is small.

Two-dimensional probability density and its various moments Using equation (6) and the polynomials (equation 5) for n = 2, we conclude after some algebra that:

Stochastic modelling of nonlinear water waves: N. C. G. Markatos W2(nl,n2)

=

{

1

bom)H.m )

1 + S=3

• l+ra=S

)

x WG2(nl,nz)

(9)

where: b0,, ) = b, +,. (t, ..... t , , t2 . . . . . t2) l times m times

Therefore:

which is an expansion in orthogonal polynomials. In the second approximation, the mean of the distribution is shifted to: (<~1),(~2))

#,o=#ox,

(~)

=~

,ff

F(iu)F(iu 0 < exp{i(u~ + u,¢,)} > dudu,

L L

(klo,k01), k,o=ko,

= (#10,U01)

where ~(t) is the signal and nl(t) is the output. We can simply apply Rice's method by considering the transformation ndt) = g(¢) = ¢(t) for 141 < Z. Then since g(¢) vanishes for ~ less than some fixed number, we can choose equation (11) to be the Fourier transform, and we have to consider only real values of u.

=

=41 f :4

'

Ju-~u2eXp~-~o',(u

where/z,, are the bivariant moments, and the abscissa is multiplied by the factor:

+ 2R,(z)uu, +

L L

x Jo(B.x/u 2 + u~2 + 2uu~cosogoZ) dudu,

(12)

1

{1 + g[bt30)Ht30) + 3b(21)Ht21) + 3b(12)H(12) + b(03)H(o3~ } introducing various kinds of skewness, specified by the parameters b,j} above. For a process of mean ( n ) and variance a 2 the modification of the derived distribution is achieved simply by using the Gaussian bivariant distribution for an unnormalized variable and by calculating the Hermite polynomials H{.,). The above results can be applied to evaluate the joint distribution of two random variables other than the elevations as, for example, the two components of slope of a random air-water surface. It is proved that in the joint distribution of surface slopes, the third-order cumulants are, at least, of order a 6, proving that the distribution of slopes is not symmetrical about the mean. This symmetry would certainly follow if, for example, it were possible to reverse the direction of time (so that forward slopes become rear slopes and vice versa). In fact this cannot happen, even for free undamped waves, as it has been shown that there exists a flow transfer of energy from one part of the spectrum to another 6'7. The third-order quasi-moments bok and the Hermite polynomials are easily calculated through known relationships.

Correlation functions and power spectra of the wavy surface elevations The model, which explains the data reasonably, and has a sound theoretical justification, represents the air-water interface as an oscillator, under the influence of a regular force Bcos(co0t + 0o)* and a stationary random force n(t). Therefore: ~(t) = n(t) + Bcos(coot + 0o)

if F(iu)ei"~du L

* B is an amplitudeand mo a characteristicfrequency

L

x

__ 2

exp/__fu~/

L

x J o ( B J u 2 + U~2 + 2uu~ cos tOoZ) -~

~.[-~

fdUJk(BU)exp( ( - 1 ) .+k n!

a2

× - - e k R '

'

-12

~ coskcooZ

(13)

The correlation function k.(z) for ~(0 is obtained from equation (13) after subtracting the term corresponding to n = 0, k = 0 (the d.c. component). Calculating the integrals involved in equation (13) (taking care of two of those which have a pole at the origin) we obtain: 2

oO

oO

|

cri ~ ~' "__l__R.h2 _ g.(z) = rtu .~o k~=On! 1 k.~kCOSkoooz

(14)

where: eo = 1 and

(10)

with n(t) being a Gaussian process with correlation function a2Rl(z), and initial phase, 00, being completely random and statistically independent of ~(t). To find the correlation function of ~(t) we used the characteristic function method given by Rice s'9. The basic idea is to represent the transfer characteristics of a nonlinear device by a contour integral:

ql = g(O = ~

where Jr(x) is the Bessel function of the first kind and order K, and i the imaginary unit. Extending the addition theorem for Bessel functions to the modified Bessel functions (see Appendix), and expanding the exponentials we obtain:

ek=2

0 h~. =

for

for

k/>l K odd and n even k, n either both odd (15)

2Qk.(B)

or both even 0

for

k even and n odd

+ n - - 1)

(11)

Qk.(B) =

F( k 2 /a2~,_2/2 BI2 /

°xP(-5) F(k + l)

Appl. Math. Modelling, 1978, Vol 2, December

229

Stochastic modelling of nonlinear water waves."N. C. G. Markatos (16)

X M 1 / 2 ( n_ 2),k/2

F(Z) is the G a m m a function and M~,u{Z } is one form of the Whittaker function (reference 10), which is directly related to the degenerate hypergeometric function ~b(a; b;~). The amplitude of the ko~o/2n component of ¢ is calculated to be:

B2kF(k _ ½)

2tk- 5/2)~0"(12k + 3)F(2k

7961

"~

5.985

-~

to 4"00~

+ 1)

¢ k-

2033

1.

2,2k + 1;

(17)

Using

0057

Pl(f) =

t

S[_~____

0

2"136

5":272 7.c)09

10"545

13:181

to

and (18)

+ oo

Figure 1 S p e c t r u m of w a v y surface elevations. ~ = 0.3; coo = 4.60c/s; a~= = 3.89; z = 0.5

[p.(f) =1 f SEn;f_ x]P.-,(x)dx 4"397

- oo

where S(X; g) is the spectral density, we calculate the continuous portion of the spectrum to be: oo 2a~ ~ ~ ~k = 7_.., Z_, --7 h2K,, Sd¢;f] --~-,,=~ ~=on!

3382

n+k> 1

2"366

x [P,,(f - k¢oo) +

P.(f + kmo)]

(19)

(k, n both either even or odd). The terms (2, 0) and (1, 1) give the greatest contribution to the series, and the remaining terms contribute less and less as n and k increase. Assuming that n(t) has the very frequently occurring correlation function al2 exp( - blz]) cos ~OoZ. Then:

sd¢: q = V-/,.o2

n3

+

L

/~2+(f

1 2(-°o)2 )

2~0o)2

aul

-04 1 12re3 16ct2 + ( f -- 2(.00)2 +

x

~10/,,2 O(u 1 t~04

48rc 3

o-1~h22 _

4090)2

+

2re---W-

4c(2 + ( f _ 4¢Oo)2 + ...

(20)

where:

h,, - B --x/~e-(82/4"P[I°( L , -~

x//~ -

h22=

e-

~ - ~ ¢[ 3, he, = ~

230

(92/4a2) [

5

0

~

0-335 - 0680

0

I

0210

0421

0632

0"843

1054

where l,(x) is the modified Bessel function of the first kind, of order n with imaginary argument. Typical correlation and spectral functions for wavy surface elevations are given in Figures 1 and 2. Their analysis proves that:

)

~I0/~2

( 16C(2 ( f,

3

Correlation function of wavy surface, c( = 0.3; (% = 4.60 c/s; a~ = 3.89; z = 05; time scaled with respect to period

a41ah21 try°r1 i. 2 [ 1 } + ha(a2 + f2) + 8--~-"o4/16a2 + f 2 /

x

-

Figure2

+ a~ / 2

+

1

%

=]

a4ah2'[

0c

- 'l( -ffff~ ) ]

B2

B2

I~" 1; -2--~a12)

A p p l . M a t h . M o d e l l i n g , 1 9 7 8 , V o l 2, D e c e m b e r

(21)

(1) The correlation and spectral functions for the Gaussian component are independent of z (z = (B2)/(2a2)), all the other parameters being the same. (2) For different variances (a2), the correlation and spectral functions for the Gaussian component have identical shapes, and magnitudes directly proportional to the variances, all the other parameters being the same. (3) For different variances (a~), the correlation and spectral functions tor the wavy surface elevations have identical shapes, and magnitudes directly proportional to the square of the variances, all the other parameters being the same. (4) The correlation and spectral functions for the wavy surface elevations have the same shapes, but different magnitudes for different z, all the other parameters being the same. The shapes of the various functions vary with the damping factor ~. With ~ as a free parameter, we can easily fit any correlation and spectral functions observed during the experiments.

Stochastic modelling of nonlinear water waves: N. C. G. Markatos D e s c r i p t i o n o f wave profile

0733

The Fourier series representation of wavy interfaces is not the best representation. Not only do we need several Fourier components for an accurate representation, thus making the functional form of the surface inconvenient for further investigation, but also, this representation is not unique, in the sense that the different Fourier components have different importance under varying measuring conditions. Further, the series associated with a sequence of disturbances occurring at random 9:

~044C E

0293 o.,

a0

l(t) = ~- + ° =Z,

[

2rent

ta"c°s-- -÷

27rot b, sin ~ - - )

(22)

regards the coefficients a, and b, as independent random variables, distributed about zero according to a normal law. Einstein and Hope 11 have discussed the normal distribution of the above coefficients, when equation (22) is used to represent black-body radiation. They concluded that these coefficients are independent. However, it can be proved that they are not independent in the general case of a random stationary process, with the only exception being that of the 'chaos function'. The alternative to the Fourier series representation, which is theoretically more sound, and has the same accuracy with the great advantage of using fewer terms, is the Edgeworth expansion in terms of the orthogonal Hermite polynomials, of order r in y, H,,

~(y) = ~, hiH,(y)a(y)

(23)

i=0

where y is measured in a coordinate system travelling at the wave propagation velocity: 1

a(y) = ~ e x p (

x/2z

- ½x2)

(24)

and

i-2. i! Hi(y)= Z (-½)" .yi-2, ,,= o n !(i - 2n) !

(25)

o J -4", )18

. -2-673

.

.

-1"328

. 0-017

1'362

2-708

X

Figure 3 Typical wave profile and its representations. (1), representation by a six term Fourier series; (2), representation by a 3 term Fourier series; (3), exoerimental profile; (4), representation by a 3 term Edgeworth series

Statistical distributions for the w a v e amplitudes

The amplitude distribution of a non-Gaussian narrowband process is a modification of the Rayleigh one, given by: ( A2+B2 ) f(A) = A/tr 2 exp 20.2 Io(AB/0. 2) (27) where f(A) is the probability density function of the amplitude, A, a is the standard deviation and B may be considered as a fictitious 'additional' contribution to the amplitude. Fitting the experimental data by equation (27) proved very satisfactory and it is from this that equation (10) was devised to represent the air-water interface; because equation (10) leads to equation (27) for the amplitudes. It may be shown that the amplitude function can be replaced by a onedimensional Markov process. Figures 4 and 5 present, as functions of Z = BA/0. 2, the wave amplitude probability density and distribution functions. The numbers on the curves correspond (increasing magnitude) to the values z - 0.10, 0.50, 1.00, 2.00, 3.00 and 4.00. 0291

+~

h, = -~.

0

~(y)H,(y) dy = r! 0233

-oo

x I /#1_

r[2] z r[4] 1 .... 1 2~.T#,_ z + ~,2.v/~,-4 (26)

where #,~ is the moment of order r about the point a. The 0th term of equation (23) is a Gaussian and successive terms represent modifications to this shape. With the above expansion, we fitted all the experimental profiles by using only two and three terms. Figure 3 presents an experimental wave profile and its theoretical approximation by Fourier and Edgeworth series. The advantages of the latter are obvious.

0-175 0"116 0058

0

r 30

3- 61

s-1'92 6. ;23

8- ,54

Figure 4 Probability density of wave amplitudes for various Z(a 2 = 3.89)

Appl. Math. Modelling, 1978, Vol 2, December

231

Stochastic modelling of nonlinear water waves." N. C. G. Markatos

800 I'000

particularly for the cases of waves generated by the simultaneous action of wind and of a mechanical wavemaker. The mean value for a Rayleigh distribution, proposed by Longuet-Higgins 4 and used by several other authors to describe the amplitude distribution of the sea surface is: {A> = ,,f-~0.. Thus, the deviation of the present distributioffs mean from the above is the factor e- =~b(3; 1; z), and is increasing with z. According to Rice s, for the Rayleigh distribution (A)/0.o = 0.886, a0 being the r.m.s, amplitude (a0 = (A 2) = 2a 2 + B2). Figure 6 presents the present experimental data (least-square fitted to (A>/ao = 0.926), together with Longuet-Higgins predictions #, and our theoretical values averaged:

. . . . . .

080C

/

0-600

iiii 0850

2"550

4"350

5c)50

7"650

9"350

Figure 5 Distribution function of wave amplitudes for various Z(o2 = 3.89)

Moments of one-dimensional probability density function The mean wave height ( A ) is given by: oo

( A ) = f Af(A) dA tl 0

B2 iF(1 +½) ,/3

= x/~aexp - - ~ 2 ] ~ # [ ~ ;

B2)

1;f~G2

( A 2 ) -- 20"24 -- 1 ; 1 ; -

(29)

(Aa> = 2a/Ee34~

(30)

2' 1" -

and in general: n

q~ --~; 1;

- ~-~-~2

(31)

The nth central moment of the distribution is easily calculated from its defining integral; for example: < ( A - )~> = 2 a ~ ~

0.886e-Z~b(3; 1;z))/x/~ + z)

The theoretical values of (A)/ao for a range of z values (limited by the limitations of the present experimental facilities, as far as changing z is concerned) are given in Table 1. As z ~ O, ( A)/ao tends towards its Rayleigh value. The average of these values leads to (A)/ao = 0.9148, which is very close to the experimental value. For the present model, the standard deviation is given by:

(28)

Similarly, using some properties of the Whittaker and the degenerate hypergeometric functions 1° we derive the second and third moments about the origin:

(A"> = 2"/2o-"F 1 +

((A)/ao

S(A) = ((A -- (A>)2> 1/2 = 0.570" 1 + - ~/

20-2

and is approximated by S(A)/ao = 0.4037, the coefficient of variation is simply C = 0.441. Longuet-Higgins 4 derived for the Rayleigh model, x/((A -- (A))2)/(A 2) = 0.453. The experimental results give the value 0.383, and the arithmetic mean of the exact values predicted by the present model is 0.399. Values of the most interesting parameters of the amplitude distribution are given in Tables 1-3, and are plotted in Figures 7-10. Inspection of the values for the central moments reveals that their variation is large, and that no simplified formula can be devised; they must be calculated by the exact equation

- 1 ; 1; 8-80C

_

4(-~b 2

1

1;-

(32)

2'

7-06C

Using known properties of the degenerate hypergeometric function, we obtain an interesting result: (A2> = 20.2 + B 2

A

5320

V 3580

and ( A 4 ) = 80 -4 + 8G2B 2 + B 4

(33)

Expressions for the coefficient of variation, skewness and kurtosis are also easily calculated. The values of the degenerate hypergeometric functions, needed in the calculations, were given by a computer program which calculates the Whittaker functions in the complex domain. The above model of the distribution of wave heights describes the bulk of the measurement well,

232

3.~7.~/

Appl. Math. Modelling, 1978, Vol 2, December

1'840 o-loo 01'07

1-894

3' 8 0

5 66

7252

9039

frO Figure 6 Variation of mean wave amplitude with ao. (1), theoretical averaged; (2). Rayleigh model; (3), experimental values

Stochastic modelling

o f n o n f i n e a r w a t e r waves." N . C. G. M a r k a t o s

Table I Values of some interesting parameters of amplitude distribution function, for various values of argument Z

(A>/~/(A2>

Z 0.10 0.20 0.30 0.40 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 Average

((A - (A>)2>'/=/a

0.8865 0.8877 0.8895 0.8916 0.8938 0.9062 0.9176 0.9275 0.9356 0.9422 0.9477 0.9522 0.9148

((A - (A>)=>312/@ ~/((A - (A>)=)/~/

0.6863 0.7132 0.7366 0.7576 0.7766 0.8455 0.8884 0.9155 0.9339 0.9476 0.9571 0.9653

0.3232 0.3628 0.3997 0.4349 0.4684 0.6045 0.7014 0.7675 0.8146 0.8510 0.8769 0.8995

Table 2 Values of coefficients of variation, skewness and kurtosis of wave amplitude distribution for various values of Z

0.4627 0.4604 0.4568 0.4527 0.4484 0.4227 0.3973 0.3737 0.3530 0.3350 0.3190 0.3052 0.3989

0-613

O405

Z

C

SK

Kurt

0.10 0.20 0.30 0.40 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00

0.5219 0.5186 0.5135 0.5078 0.5016 0.4665 0.4329 0.4029 0.3772 0.3555 0.3366 0.3205

0.6135 0.5912 0.5688 0.5397 0.5068 0.3722 0.2653 0.1969 0.1456 0.0976 0.0655 0.0267

0.2486 0.2058 0.1387 0.0869 0.0513 - 0.0945 -0.1121 -0.1239 -0.1018 - 0.0121 0.0508 0.1997

"t

031E

~o O171 ,j 0023

-O123

0100 Table 3 Values for first four moments around origin of amplitude distribution, for various values of Z

Z

(A>/a

/o=

/o -~

/a =

0.10 0.20 0.30 0.40 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00

1.3148 1.3752 1.4343 1.4919 1.5481 1.8124 2.0519 2.2719 2.4754 2.6649 2.8432 3.0113

2.20 2.40 2.60 2.80 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00

4.3299 4.9144 5.5133 6.1251 6.7496 10.0665 13.6857 17.5918 21.7652 26.1891 30.8560 35.7495

9.64 11.36 13.16 15.04 17.00 28.00 41.00 56.00 73.00 92.00 113.00 136.00

0.965

1

1

0"809

2227

2.c)36

3"645

Z Figure 8 Variation of coefficients of variation. (1), skewness; (2), kurtosis; (3); with Z

_---------------3

O'443

J OO31 A ~-- -O-381 A V 4:-O793 V

f

-1205 -I.611

0-836

1"518

0"100

i

O" 09

i

1"518 2"227 Z

i

2.936

i

3"645

Figure 9 Variation of 2rid, 3rd and 4th central moments with Z (log scale)

0-708

% 6

to 0.580 O451 0"323 0"100.

i

0"809

1"5'18 2":~27

i

2'c)36 3"645

Z Figure 7 with Z

(1), variation of S/G with Z; (2), variation of S3/G 3

The coefficient of skewness varies inversely to the variation of z. This coefficient tends to zero for large values of z, as is predicted by the model. The amplitude distribution presents a positive skewness. The second and fourth central moments increase continuously with z, as fractions of 0.2 and tr4, respectively. The mode of the amplitude distribution must be found from the solution of an equation of the form:

A p p l . M a t h . Modelling, 1978, V o l 2, December

233

Stochastic m o d e l l i n g o f nonlinear water waves. N. C. G. Markatos 2 133

/

4

0.995

/ ~

1-730

A V

i

......_l_~

/

0.95"2

1-327

0.919

0-924 / O521

~

o.u8

/

O-881

I 0-843

0.1OO

I

I

0809

I

1'518

I

2"227 Z

I

2-936

3.645 0805

Figure 10 Variation of first 4 moments about origin, with Z. (log scale, numbers correspond to order)

(zl # 0)

Io(Zl)/lo(Zl) = -- 1/zl + zl~r2/B 2

(34)

Results of the numerical solution are given in Table 4 and are plotted in Figure 11.

I

05

1954

3-409

4863 Z

63'18

Figure 11 Variation of mode re(A) with Z

where IAI is the determinant of the covariance matrix, and stationarity of n(t) implies: ( Y l Y I , ) = (Y2Y2~) = R('c)

Multi-dimensional probability density

(YtY2r) = - (Y2Yl,) = S(z)

The problem is very difficult for non-Gaussian oscillations. We used a procedure 12 for a normal process, which we extended to the non-Gaussian amplitude oscillations, n(t) in equation (10) can be written as: (Yl cos mot - Y2 sin ~oot)

say ~

(37)

Z, -- ~ a..y,, the covariance matrix/~ of the Z's is n=l

(35)

1 f ( Y l , Y 2 , Y l ~ , Y z d - 4zr2lA]1/2

related to that of the y's (A) by p = A A A 1, where A 1 is the transpose of A, which is the matrix of this transformation (Z -- Ay). Transforming from the above random variables to the wavy surface variables (by substituting equation (35) into equation (10) and comparing with ¢(t) = A(t)cos(eJot + ~(t)) and observing that the Jacobian of the transtormation is simply AA~, we obtain the following four-dimensional distribution: AAz f ( A , ~b,At, ~r) - 4~z21A[1/2

exp - 2lAlX/2[a2(y2 + y~ + y2 + y2r

x exp

- 2R(z) x (YxYlr + Y2Y2~) (36)

J

Table 4

I

21A~/~[lo-(A

+ 2B 2 + A~) - 2BER(z)

+ {(2R(z)ArB - 2ArBa 2) cos(0r - 0o)

-- 2S(z)(ylY2r - Y2Yl~)]

+ 2S(z)A~B sin(0r - 0o)} + { ( 2 R ( z ) A B - 2ABa2)cos(0 - 0o)

Numerical solution of equation for mode of wave amplitude distribution

Z,

I°(Z' ) Io(Z, )

Z = B=/2a=

x~m(A)/x/(A =)

m(A)/a

m(A)/(A)

0.10 0.50 0.70 0.90 1.10 1.30 2.00 2.50 3.40 4.40 5.50 6.40 7.50 8.50

0.0499 0.2425 0.3301 0.4098 0.4806" 0.5426 0.6977 0.7650 0.8357 0.8773 0.9038 0.9181 0.9307 0.9391

0.0049 0.1114 0.1990 0.2958 0.3957 0.4954 0.8348 1.0729 1.5046 1.9916 2.5331 2.9783 3.5242 4.0214

1.0000 1.0043 1.0133 1.0277 1.0465 1.0678 1.1425 1.1853 1.2384 1.2745 1.2999 1.3146 1.3281 1.3375

1.0024 1.0590 1.1094 1.1699 1.2364 1.3059 1.5477 1.7066 1.9599 2.2045 2.4436 2.6222 2.8250 2.9971

--* 0.8053 0.8067 0.81 56 0.8287 0.8435 --* 0.9415 0.9551 0.9703 0.9871 0.9839 0.9936 0.9953

*These values of Z are not found in experimental runs, and no been calculated

234

)

It is proved that if n Gaussian variables y. are transformed into n variables Z. by the transformation

Yl,Y2 being linearly related to n(t), are also Gaussian with zero means. Multiplying equation (35) by the similar expression at time t + z, averaging and imposing the condition that n be stationary we obtain the joint probability density of the four random variables Yt, Y2, Ylr, Y2t ;

x

7772

Appl. Math. Modelling, 1978, Vol 2, December

m(A)/(A)

ratio has

t

Stochastic modelling of nonlinear water waves."N. C. G. Markatos -

2S(z)AB sin(~k - 0o)}

- { ( 2 R ( z ) A A , cos(qJ~ -- ~) + 2S(z)AA, sin(ql, - ¢)}] } (38) for A, A,/> 0 and 0 ~< W, W, ~< 2n and zero otherwise. For ~(t) to be a stationary process, the initial phase 0o must be completely random, f(Oo) = ½r~; and because of its independence from the other variables, we get:

f(A, ~k,A,, ~,, 0o) = f(Oo) .f(A, ~k,A,, 0.)

(39)

The joint probability density function of A and A, was obtained by integrating equation (39) over all values of ~O,~b~ and 0o. The integration is given in the Appendix. The result is:

f(A, h,) = AAJIAI ~12 x exp

{,

21£[1/210-2(A2 + 2B 2 + A~)

R,,(z) = Ro(r) + RI(r) cos COot + Rz(r) cos 2O9oZ + -...

(40)

where Io, Ik, are the modified Bessel functions of the first kind of order 0 and k respectively, and K, L and Ro are defined in the Appendix. In case of a fixed initial phase, f(A, AO has the same form as equation (40) except for a constant, and for the modulus of the complex number involved which degenerates from Ro to 1.

Correlation functions and power spectra of wave amplitude The correlation function is the average value of the function~

~ n,k=0 n+k>0

Knk(r)COSkz

(44)

which vanishes as r - > oo, since Knk contains q~ terms multiplied by the cosine terms and ~b - > 0 as Z - - > o(3. K n k is a function of the correlation functions of the processes Yl, Y2, Yl~,Y2T (defined in the previous section) which vanish as z - > oo since the processes are random and not periodic, and the whole combination of them giving gnk is a decreasing function of z, leading to zero continuous portion of the correlation function as z - > oo. With:

K,,k = L,,k-a"'l~l

(45)

which is a very frequently occurring relation in many physical phenomena, and which satisfies our boundary conditions, we obtain: ov

S[A - (A) : f] = F ( A ' A O = { O A~

(43)

The number of the spectral components lying in bands near the frequency k~oo(k = 0, 1,2...) depends on the damping exerted by the functions Ro, R1 ..... Using R,(z) = R~(z) + Rc(z) where R~(z) is that part that does not vanish as z -~ ~ and Re(z) is the correlation function of the 'continuous' portion of the spectrum, substituting into the definition of the spectrum, and integrating, we get the exact expression for the continuous part of the spectrum of the wave amplitudes. We have: Rc(z) =

--2R(z)B2]}IIo(K).lo(L).Io(Ro) + 2k=1 ~ Ik(K)Ik(L)Ik(Ro)COSkt]

which is, of course, the same as the mean value of the amplitude and that, second, there are no periodic components in the wave amplitude function since all the cos terms vanish as r - > ~ . For limited range of the time interval z, we may write:

allwhenb°thA'A~>0other A's

E

n,k = 0 n+k>0

f2

4ankLnk(f2

--

k2

2 + k 2 + a.k 2 2

2

-- ank) + 4ankf

2

(46) which is given by: +oo

+~x)

: dA f dA,F(A, AOf(A,A 0 -oo

(41)

-oo

where a.k are coefficients and f is the frequency. If a.k << k then the spectral density is concentrated in a relatively narrow band f - k = O(a.k) and we get:

S[A - ( A ) ; f ] -~ ~, .,k=0 (f

Carrying out the integrations we finish up with:

R(z) = ~

1

f

B2( 0-2

exp~-

1

R(z))

n+k>0

"C~o~ " ~ -

/'I

r--" _ B2/2#2

q~(3; 1;2~2 )

- -

(48)

TC n,k = O ank n+k>0

(42)

where the definitions of "W, q and ~ are given in the Appendix. This correlation function is not periodic as, superficially, could be suggested by its form. It is easy to see that first the amplitude of d.c. component of wave height A(t) is given by ~ax/rte

(47)

The time up to which there is appreciable correlation is:

}

x ~ ~ . 1 .,ittok,2.+k+a(q) .=Ok=otkn!tn + k)! × 'Yv'in+kcos kt

2a"kL"k

-77+ a.k2

A necessary condition for the above representation (equation 45) is the correlation functions of ya,y2 to be functions of I~l and not of any other power of [z[. We showed Yl,Y2 to be Gaussian. The simplest process which is sufficient-to give equation (45) is the 'exponentially correlated process' which has a 'stationary' correlation function e -bIll and spectral density 4b/(f 2 + b 2) and is simultaneously stationary, Gaussian and Markovian. A typical correlation function for the amplitude is given in Figure 12 and its spectrum in Figure 13.

Appl. Math. Modelling, 1978, Vol 2, December 235

Stochastic modelling of nonlinear water waves. N. C. G. Markatos I000

expansion for the amplitude time derivatives and applying the Bunimovich 12 equation, we derived a complicated equation which for (dA(t)/dt) = 0 simplifies to:

0-792

h2+B 2] ' h B ' I

N(h)=~2exp {

0.583 q: 0 3 7 5

x I D - 6 ( 0 ) { ~ a l Kurt + --~-axSK3452 + 5SK + 3 0 o -2}

0-167

_..//I 0041

0

Figure 12

1318

2.636

, 3-954 l:

+ 8x/-2alD_5(O)SK 1 (50)

"I

5"273

where Dp(z) are the parabolic cylinder functions.

6"591

Normalized correlation function of w a v e amplitude

function

Probability density of wavy surface slope at moment of inception of a peak It is easily proved that the distribution of this variable is a modification of the Rayleigh distribution, towards which it tends as a first approximation.

2 306 t.849

Distribution of wave peak durations and area occupied by waves over fixed level h

1-392

To calculate those distributions we expanded the elevation function ¢(t) to a Taylor series and regarded the second derivative at the time of inception of a peak as non-random. Under the above assumption and after some algebra the results are easily obtained.

O

m 0.935 0.477

Wave period 0020 0-300

L

i

I

2-936

5"573

8-209

i

10"845

i

13-482

f

Figure 13 Spectrum of wave amplitude function Statistical analysis of other important wave characteristics

Analysis of some further important wave characteristics is needed to complete the model for the description of a random, nonlinear wavy interface. The results follow:

Mean number of wavy surface elevations greater than a fixed level h Using the two-dimensional probability density for the surface elevations (equation 9) and their time derivatives together with the Bunimovich 12 equation we get: o1

f

1,2] j~

N(h) = ~-~exp~ - s n

The wave profile is skewed towards the air velocity direction. The intervals between troughs and their preceding crests (Type A, say) are smaller than the intervals between these crests and their preceding troughs (Type B interval). The experimental correlation coefficients prove that the fluctuations in the interval lengths can make an A interval to become greater than a B interval, L intervals distant, or a B interval to become smaller than an A interval, L apart, but this L is never less than 6; and it is on average, 11. There exist neither longer term cycles of change within the sample nor non-stationarity. The frequency goes above and below the mean value alternating, and so there is a tight regulation on a cycle-to-cycle basis. There is no evidence of a frequency control in which A's length is raised whilst B's length is lowered. Figure 14 presents the percentage that the A type interval is smaller than the B type interval; and also the ratio of the length of the ith A interval, ai, to the interval between two consecutive crests I~ versus wave amplitude. Evidently, there is an increase in

x[l+l[bi3o)h(h2-3)+abi21,(h2-1)~

,49,

1300

where o'1 is the standard deviation of the random process of the surface-time slope and b's are the quasimoments.

°/°

1500

o11-00

j.J . . . . . . .

~-

•~. "-'~/"

500

236

AppI. Math. Modelling, 1978, Vol 2, December

Figure 14

'

047 "~'~'-'~'~""~

I

o

There is no correlation between the wave amplitude function and its derivative; and also we assumed their statistical independence. Then, using a Gram-Charlier

049

o.48 o

~'9-00 700

Mean number of wave amplitudes, greater or equal to a fixed level h

A

./

[

I

~

B

346

i

2.00 4.60 660 860 10.00'12.60,45 Percentage by which A Wpe intervals are smaller than

B Wpe (curve A), and ei/li ratio (curve B) versus mean wave amplitude

Stochastic modelling of nonlinear water waves: N. C. G. Markatos asymmetry with increasing wave amplitude: up to ( A ) ,-- 3.2ram this increase is at a rather low rate; then it becomes greater and it becomes lower again after ( A ) = 5.8 ram. The problem of finding the probability density function of the time intervals between the points at which the wavy surface crosses the mean level is extremely difficult. An effort by Rice 9 gave only very conditional results. T o approximate It we used a perturbation method similar to one reported by Stratonovich 2 to derive the jitter of a relay. If the process S(t) creates the equally spaced points t,, then S(to) = 0 where to is the time at which the curve has a zero crossing. In fact a n o t h e r random process exists (see section on correlation functions and power spectra of wavy surface deviations) which leads to a time shift from to to ti such that S(ti) + tt(ti) = 0

(51)

The perturbation method consists of multiplying t/(t) by a small parameter, which finally is set equal to 1 and, searching for a solution of the form t i = t o + ~ t t + gzt 2 +

.

.

.

such that: (52)

S(ti) - S(to) = X ( t i ) = - ~rl(ti)

(the weakness of q(t)justifies these successive approximations). Expanding X and t / i n t o Taylor series, equating coefficients of identical powers of e, solving the resulting equations and neglecting terms of order 3 and higher, we get the equation giving the time shifts, depending on the derivatives of S(t) and t/(t) at to. Using the statistical independence of tt, t/' (q(t) is Gaussian), expanding the characteristic function in a series with respect to t p f / a S ' - S " t l 2 / 2 a S 2 ( a being the standard deviation of t/(t)), averaging and Fourier transforming, we finish up with: W(v) =

aS" 2 3 ~ e11 + 2 - ~ ( z - "c ) ~/2~

~2/2)

)~

(53)

For the higher velocities, without the action of the wave-maker, the Rayleigh model fits the data only within ~ 1 2 ~ , and the accuracy improves if we use the proposed distribution function, searching for values of a 2 and B 2 which minimize the deviations. F o r these higher velocities, the wavy surface m a y be considered to consist of a periodic process with constant amplitude, on which a completely r a n d o m Gaussian process is superimposed. The resultant process is no longer Gaussian. The same always holds true when the wave-maker is in action. If the experimental errors were of the order of 1 0 ~ , they could present the Rayleigh model as less adequate than the one proposed here. Extensive analysis of the errors proved that their effect is well below 5 %. It is reasonable, therefore, to believe that the present conclusions are valid. There is also theoretical justification. The Rayleigh distribution for the amplitudes is derived from a Gaussian model for the surface elevations. F o r the higher air velocities, the surface elevations are definitely not Gaussian, since their distribution possesses both skewness and kurtosis. Hence, from the theoretical point of view, a Rayleigh distribution is not justified for the higher air velocities. The best representation of the wave profile is achieved by an expansion, starting with the normal function and continuing with certain of its derivatives. The second approximation introduces the skewness, and the third the kurtosis. The theoretically correct procedure is to allow the third term to contain both the fourth and the sixth derivatives. The reason being that the coefficients of the terms in the fourth and sixth derivatives are of the same order of magnitude for large samples. The coefficients are expressed in terms of the central moments of second-, third- and fourthorder. We prefer not to choose functions with more than four parameters. Moments of high order, as it is easily seen from the definition and the method of computing them, depend to a great extent on the 'boundaries' of the distribution, which may not be indicative of the distribution.

where z = S ' ( ( t i - t o ) / a ( n ) ) . This distribution is not symmetrical. References

Conclusions The waves studied were both wind-driven waves and artificially generated waves superimposed on winddriven waves. The wavy structure is in most cases covered by small capillary waves. T h e latter have been integrated over the large waves for the purposes of this study. The wave amplitude, the surface elevations and other wave characteristics were considered as stochastic processes, characterizable by their moments. There is a strong two-dimensional character of the flow. However, the results are not inconsistent with the ocean data, thus encouraging l a b o r a t o r y studies for oceanographic purposes. For the lower air velocities, up to 0.8 m/s, and without the action of the artificial wave-maker, the Rayleigh model fits the d a t a for the wave amplitude well; and the surface elevations follow a Gaussian law.

1 2

Markatos, N. C. G. Trans. Instn. Chem. Engrs. 1976, 54, 184 Stratonovich, R. L. 'Topics in the Theory of Random Noise" (Translated by R. A. Silverman), Science Publishers, Rev. English Edn. 1967 3 Edgeworth, F. Y. Trans. Camb. Phil. Soc. 1906, 20, 36 and 20, 113 4 Longuet-Higgins,M. S. J. Fluid Mech. 1963, 16, 138 and 17, 459 5 Kinsman,B. Chesapeake Bay Inst. Tech. Rep. no 19, 1960 6 Phillips, O. M. J. Fluid Mech. 1960, 9, 193 and 11, 143 7 Hasselmann,K. Proc. Conf. Ocean Wave Spectra, Easton, USA, 1961 8 Rice,S. O. (1944) Bell System Tech. J. 1944, 23, 282 9 Rice,S. O. Bell System Tech. J. 1945, 24, 46; and 'Selected Papers on Noise and Stochastic Processes'. (Ed. N. Wax), Dover Publications, New York, 1954 10 Rysbik,I. M. and Gradstein, I. S. Table of Integrals, Series and Products, 4th edn. Translation edit. by A. Jeffrey, Academic Press, New York, 1965 11 Einstein,A. and Hope, D. Ann. d. Physik. 1910, 33, 1095 12 Kuznetsov,P. I. et al. In, Non-linear Transformations of Stochastic Processes. (Eds. J. Wise and D. C. Cooper), Pergamon, Oxford, 1965

Appl. Math. Modelling, 1978, Vol 2, December

237

Stochastic modelling of nonlinear water waves. N. C. G. Markatos

Appendix Integrating equation (38) with respect to 0o from 0 to 2n we obtain:

dimensional space.) Taking B~ = J~, v to be an integer and passing to the limit as v ~ 0 we get:

J o(mR) = J o(mS)J o(mp)

2~

+ 2 ~ ( - 1)kJk(mS)Jk(mp)COSk¢

f f(A, O, A,, 0,, 0o) d0o o

Setting m = i and using:

2n

Je(iz) = iele(z)

= f f ~ / f exp(P cos 0o + q sin 0o) d0o

(A8)

we have:

0

(A1)

= f#Jln{Io(x/~) × 2}

where, if, ..1/, P and q are easily defined by comparing (A1) with equation (38) and: c = p2 + qZ (A2) Let 2n

2n

2n

lo(R) = Io(S)Io(p) + 2 ~ Ik(s)Ik(p)COSkq5

o ~ = f d O f d~b~f dOof(A,O,A~,O~,Oo, 0

Then: 2~

2n

= f dO2nfflo(K)Io(L).2nlo(Ro) + f x 2~

= 8n3fflo(K)Io(L)Io(Ro)

2S(z)[R(z) -

0-2] - 0-2]2 _ S2(z)

(A4)

2~

AA~ fR(z)[{R(z) -

+

where

+ [--~-[(R(z) - 0-2)2 + S(z)2]

r S(z)[{R(z)-

{ AA~

-R(z){(R(z) -

2re

--

Equation (A10) leads directly to equation (40) Concerning equation (42) the definitions of q, ~¢r and ~ follow:

- 2S(z)R(r)[R(r) - 0-~ - - 0-2~ 2 ..]_ S 2 ( T )

0-2]

/'2

J

[R(r)

-

S-~

Re

COS t =

0-2)2 _ S2(r)} -] + 2SB(z)[R(z) _ 0-2] [R(z) - o 2] 2 + S2(z) cos

+

[--~T~--~ +

]

R0 = (p2 + p 2 ) 1 / 2 0 ~ t ~ 2n

+ 4S2(z)[R(z) - 0-212}1/2 cos ~) exp IAII/2

x

a2} 2 -- S2(z)] -- 2S(z)R(z)[R(z)

AA, [ e~ = ~A--ST~L

2A~AB2

+ __[A----i--{([R(z) - 0-2]2 _ $2(z))2

B

q = [A-~7~[(R(z)__- 0-2)2 ..1_S(z)] 1/2

sin

L([R(z){(R(r)-a2) 2

2n

1

=fdq/fdctff2nlo(x/K2+L2-2KLcoscO

~r2 =

I-~

(AID

-t- S 2 ( ' c ) } - - 2 S 2 ( ' c ) 0 - 2 ] 2

[{g(*) - 0-2}2 + sa(~)]2

+

0

x exp{P1 cosct + P2 sine} (A5) where K, L, P1 and P2 are easily defined. The well-known 'summation theorem' for Bessel functions is not valid here since L < 0. To avoid this condition we write the alternative relation for the Bessel functions:

B~(mR) R~ =

2~m-°F(v)

B~+k(mp) p~

[S(z){(R(r) - 0-2)2 S2(z)} /7 -- 2S(z)R(z)[ R(z) - 0-2]]z _

+

~ (v + k) Jv+k(m0-)

A2"+k+2exp

{

o2A2 ~Jk(qA)dA

21AI.2j

(A6)

where Cy,(cosq6) is the ultraspherical or Gegenbauer polynomial of degree k and order v (connected with the theory of spherical wave functions in (2v + 2) -

Appl. Math. Modelling, 1978, Vol 2, December

= 0

S~

• C~,(cos ¢)

[{R(z) -- 0-2}2 + $2(z)]2 (A12t

~k,z,,+k+3(q) k=O

238

0-2}2 + SZ(z)] _ 2S2(z)0-2]

0

A2B 2

x

(A10)

K=I

=

0

lk(K)1k(L)lk(Ro)coskt

+ 16n3ff ~

we obtain:

0

Ik(K)Ik(L)[rcR-klk(Ro). R k • 2 cos kt]

k=l

By making the transformation:

2n

dO2nff

0

(A3)

0

c~ = 0~ - ~ - arc tgrR( r)t

(A9)

k=l

0

o

(A7)

k=l

X

qkF(k + n + ~) 1( o2 I,,+k+3/2 r(k + l) 2~+ ~ /

(A13)