A method for the numerical simulation of unsteady gravity-capillary waves of finite amplitude

U.S.S.R. Comput.Maths.Math.Phys.,Vo1.29,No.3,pp.140-146,1989 Printed in Great Britain

0041-5553/89 $lO.OO+O.OO 01990 Pergamon Press plc

A METHOD FOR THE NUMERICAL SIMULATION OF UNSTEADY GRAVITY-CAPILLARY WAVES OF FINITE AMPLITUDE* V.R. KOGAN and V.V. KUZNETSOV

An algorithm is given for calculating unsteady gravity-capillarywaves on the surface of a fluid of infinite depth. From the computational stand-point, the method is based on fast inversion of the operator corresponding to the Dirichlet problem, and on the analytic apparatus of the theory of boundary-valueproblems. The dynamics of surface gravity-capillarywaves were studied in /l/. The difficulties of numerical simulation when calculating the surface tension forces were mentioned there. In the present paper we take a more modern form of the algorithm of /l/ in which we use the duality between the solution of the Dirichlet problem and the conformal mapping of the flow domain into the unit circle. This device is not new. Starting with the classical work of Nekrasov and Levi-Civita, it has been used by a number of research workers, see e.g., /2, 3/, when finding the profile and velocity of progressive waves. But steady waves were exclusively considered. The problem thus reduced to finding a single conformal mapping, and no question arose concerning the motion of points over the unit circle modefling the surface wave dynamics. In the case of unsteady surface waves, the conformal mapping G(z,t) depends on time. Its construction on each time layer is quite a complex computational problem, and its solution has only been made possible by the rapid solution of the Dirichlet problem by the method of boundary integral equations, and by using duality relations. By using variational theory of conformal mappings, we can obtain on the unit circle a family of non-autonomousdynamic systems which depends on the wave form and contains all the information about its dynamics. Thus the mapping G(z,t) is on the one hand a part of the computational algorithm, and on the other hand, enables the analytic apparatus of boundary value problems to be widely used in the study of unsteady surface waves, and in particular, in the study of their stability. 1. Formulation of the problem. We consider the plane-parallel irrotational flow, periodic in the horizontal coordinate, of a fluid of infinite depth. The fluid is assumed to be incompressible,inviscid, and of constant density, and to move under the action of gravity and surface tension forces. We are interested in the evolution of an initial disturbance, specified on the fluid surface. Since we assume that the flow is irrotational and periodic, it suffices to find the harmonic function 'p (potential1in the domain D, bounded by the wave period and the free surface. For each fixed t, the function cp(t, x, y) satisfies Laplace's equation inside D, and the non-linear Bernoulli equation on the free wave surface /l/. We use the conformal mapping c=e-':, z=x-tig, to map D, into the domain DC, bounded by a closed curve C in the j plane. We write b as <=fYe, and introduce local coordinates (s,0) into the neighbourhood of C, where s is the arc length of C and n is normal to C. The system of equations for q. r.8, defining the motion of a point of the free surface, is /l/

where K is the local curvature of C, and p., is the atmospheric pressure (throughout, p,=O). System (1.1) is written in dimensionless form, and \Ve is Weber's number. The system obviously has to be closed by the equation for the normal derivative d&n on C. We obtain this by passage to the limit, say in the Sokhotskii-Plemelformulae, and it has the form

ibZh.vychis~.Mat.mat.Fis.,29,6,844-852,1989 140

141 the integration is carried out with is the distance between points P,Q=C(t); Here, R(P, Q) respect to the coordinates of point P, parametrized either by the arc length s or by the "angle of observation" o. 2. A method for the numericat solution of system (1.1) and the Schmrtz operator. be the polar coordinates of the j-th particle on the contour c(t) at any Let (0,. 6) instant t, and let r&J be the value of the potential. The subscript j varies from 1 to N and indicates that a variable refers to the j-th particle, j being considered as the Lagrange coordinate for any t. System (1.1) then has the form (2.la)

(2.lb)

(2.lc)

We solve the Cauchy problem for system (2.1) by the Runge-Kutta method in conjunction with Adams' method. To evaluate the derivatives (d/as),, we introduce a new parameter j., which varies continuously and monotonically in the interval IGj,GN and is such that j.=j at integervalued points of the interval. The arc length S is then given by

[pg)‘+(ry~]”

+

(2.2)

.B

Notice that, on the wave surface, (&/aj.)#O,since the fluid is incompressible,so that the functions (do/as), (arias), (arpias) are well defined by the general relation

The functions evaluating

d’O/dj,’

v(i.),are approximated by cubic splines. in (2.2), and from the five-point difference relations

On

we obtain the expression for the second derivatives 82rld.sz, d’O/ds’: +O((As)').

It has a higher order of approximation than the approximation by cubic splines, which, when evaluating the second derivatives, gives an error of O((As)*), This fact is important for finding the curvature. Let us briefly describe the scheme for solving Eq.(1.2), noting the main points. Let Q be a fixed point on the contour C(t) with coordinates (r,, O,), and let P be the running contour point. We measure the arc length from the point Q, i.e., s,=o. Obviously, InR(Q,P) has a singularity at Q=P, but the integral on the left of (1.2) can be regularized. We have

(2.3)

where imin(sN-6. 0.05.~~).Comparison of SL,sl.9 with sn. is necessary, since points on the contour can move irregularly during the computations. Denote the integrals on the right-hand side of (2.3) by I,,I?.IS,I;.IS respectively. The first three integrals have no singularities and are evaluated from the four-point Lagrange formulae with a local approximation error of O((As)'), where As is the maximum arc length between two adjacent points. The integral

I,

is evaluated

by summation

of integrals

of the type

142 A,.,

JT J[(,I,

+(s-s,)

‘I

(2 1,+...

(2.4)

Ins ds,

+ -qr(sn),]

and I, by summation of integrals of the type

(2.5)

The normal derivatives in (2.4), (2.5) were approximated by fourth-order Lagrange polynomials with centre at the point j. Each term in (2.41, (2.5) was then evaluated exactly using the expression aI.0

s

s"ln8ds =

9,

Using (2.3)-(2.51,we obtain the system of algebraic equations A.

(2 1=F.(cp).

(2.6)

” an ,

where A,, is a matrix whose elements depend only on the type of contour C(t), and P<(V) is determined by the value of the potential on C(t). System (2.6) was solved in /l/ by Gauss's method, i.e., it was reduced on the "forward run" to

(2) , =W,(cp)

Ra,

with lower triangular matrices R and Q, which depend on the coordinates of find on the "backward run":

C(t).

We

easily

av ( 272, 1 =R,A?,*F,(cp).

(2.5)

It follows dram Green's formula that the Dirichlet problem in D: is solved by finding @/an. into the unit Let us show that this fact can be used to obtain a conformal mapping of D; circle and included in the computational procedure. such that In the domain D: with boundary C(t) we take the analytic function g=g,Sig, g,I,,!l=-lnl:-:C,.I.

(2.8)

where %EC(I) and 5,, is the image of the point at infinity in the D, plane. Substituting (2.8) into the right-hand side of (2.7), we find by the backward run ag,lanlor,. From the Cauchy-Riemann formula dg,/dn=dg,ids

we have g2(s)=

+I.*.

In short, we know the boundary value of the analytic function g on boundary-valueof the function G(S.5,J=e'"(t-c,>)e'

(2.9)

C(t), and hence the

(p is a real constant). On continuing G(g,5.) analytically inside DC. we find G(C.50)=e'iY5 -f,t)@.

(2.10)

Obviously, InlGI,,,,=lnl~-S,,J+R,I,,,,=I). The algorithm is thus constructive method of obtaining the conformal napping of DL into the unit circle; in the terminology of /4/, G is the Schwartz operator. Since the matrices R and Q depend only on C(t). the construction of the conformal napping G additionally requires only one Gauss backward run. The entire procedure for finding @/an and G at each time step takes roughly 15 sec. of ES-1055 M computer time with A'-60. 3. Modernization of the algorithm based on the Schuartz operator. We put &=O in (2.10), and let G map D, into the unit circle Dy in the 5 plane. We introduce into Df the polar coordinates (p,?j).Obviously, p=re~~,~=e+gZ+~, where g,

143

is given by (2.9). Since G is conformal,

where c(t) is the image of c(t) under the mapping G. are the real and imaginary parts of the analytic function Since acp/N and &p& in the circle IPIG% they are connected on c(t) by the Hilbert transformation dWMl zn (3.1)

Let [a,sin(ka)+b,cos(ka)];

then /5/ [a,cos(k6)-b,sin(k0)1, i.e.,

a,cos(kQ-b,sin(ki3)1. *=a

(3.2)

Since g,(,,,,=-lnr,we easily see from (2.10) that IG'lo,,=d6/as+ag,/as, or, using the Cauchy-Riemann formula, JG'lc,,,-a6/as+dg,/dn. We used (3.2) to evaluate the normal derivative &p/&zlcc,,. Since 'p is harmonic in the belongs to the Holder class H(u) with LL=~. Then, see domain 9, the function &@a /5/, acp/ao belongs to the class H(l-a), where E is an infinitely small positive constant. In this case, therefore, the Hilbert transformation (3.1) is smoothing. As a result, the modernized algorithm is more stable and accounts than the algorithm of Sect.2. The computing time is increased by 5%. Table 1 shows comparative data on the mean curvature of the profile R, the total energy E, and the integral A of the normal derivative of the potential over the contour c(t), evaluated by the two algorithms. Here, R,E, and A are invariants of the irrotational flow: II Z = 2

E=Q+v+n

jK(e)dO=O, ”

=

COllSL,

(3.3)

where the wave kinetic and potential energies are respectively

while the energy due to surface tension forces is

rI=w,

A -

[ _p&-2n];

(3.4)

j,;dr=O.

It is clear from Table 1 that the invariants (3.31, (3.4) are much better preserved in the modernized algorithm than in that of Sect.2. The dynamics of gravity-capillarywaves with initial curvature 2a/h=0.06 were modelled numerically. The initial profile and potential were specified as

144 2an

cp.=-*exp(y.)sinr. h

y~=-TcoS5~

It was found from the numerical modelling that the ripple arises, not at the wave vertex, as was previously assumed in the physical literature /6/, but at the frontal slope close to the angle, and its appearance is connected with loss of symmetry of the profile. At the instant when the ripple appears, the highest space harmonics increase very rapidly in the profile and t=6.6 respectively. expansion. In Fig.l,a,b, we show the wave profiles at t=O.O Table 1

7:0.0

u.ti 1.2

1.8 2.4 3.0 3.6 4.2 A.8 5.4 6.0 ti.6

f Sect.

n

Algcjr.ithms -i-

r’ 4

-0.18.W” -0.48. IO-& -&81. IO-* -0.63. to-0.42.lk (1.42.10-S 0.25.10-’ 0 13 lo-,’ u.40.1lr 0.18.10~’ 0.21.10-” 0.56. lo-

0.69. itI-”

O.iO94682 0.45.IU-s&LO94663 U.b9. IO-” U.1094650 n.91. lo-~ 0.1094646 0.34.i0- Il.1094649 0.11 0.109465i 0.31. IO-” 0.1094678

IO-0

0.1 LlW5 0.20.10-~ 0.21.10-’

0.1094668

F

A

-

0.3i.10-*0.1094iOI 0.88.10-hO.109'1696

modernized

_-

_E

L

-0.23.10-8 -O.38.1O-B -0.15.10-6 -O.18.1O-8 -0.42. 1O-8 -0.35.10-7 0.12.10-e 0.23. IO-6 0.37. lo-”

O.i3.lO-B

0.1094649 0.12.10-5 0.1094624 0.24 10-a

I

E

I

0.18.10-80.1094698

0.11.10-” 0.20. IO-1 0.23. lo-8 0.33.10-E O.24.lO-o 0.30. IO-8 0.41’10-* 0.12.10-’ 0.49. IO-’ 0.82.10-’ 0.13~10-~

a

0.1094695 0.1094692 0.1094690 0.1094689 0.1094687 0.1094686

0.1094683 0.1094686 0.1094684 0.1094680 0.10946i5

b

Fig.2 Fig.1 4. Equations of the gravity-capillary wves on the unit circle. Let C(t!) and C(t,) be the shapes of the wave free surface on the 5 plane, corresponding to the instants t, and t? respectively, C(L) being obtained from C(t,) by a -__ shift along the system trajectories df/dt=dW/d~. where W(5.t)=cp(~,t)+i~(S,t)is the timedependent complex potential. Put T,=: C(&)-CCL). The transformation TIC and the Schwartz operator G generate a transformation of the circle in accordance with C(t) T,“-GzT,%-‘.

Since W(c, t) is analytic for every t, then TtC, Lagrange derivative in coordinate (r,i3) is

and hence T,‘,

are not analytic. The

It can easily be shown that (4.2)

Since G is conformal, and using (4.1), (4.2), we easily obtain the Lagrangian derivative in (I). ii) coordinates: (4.3)

Applying the operator Ll/Dt

to the variables

p,?, we obtain

DO/r!Jt=r’IG’{ %pir/fi.

(4.4)

145 (4.5)

In time At we have from (4.5):

in a time At under the This relation gives the radial deformation of the unit circle C(t) action of the vector field (4.41, (4.5). We map conformally the deformed contour into the unit circle. Here, as we know from the lead to rotation through the angle variational theory of conformal mappings /4/. small BP A% of points on the unit circle: 21 AB=$Ap(k)e@0

8-Q 2 df.

(4.6)

Finally. from (1.11. (4.41 and (4.61. we obtain the system of eauations on the unit circle _ C(t) for the dynamics of gravity-capillarywaves: *n (4.7)

(4.8)

Using (4.31, we can easily calculate the commutator of the operators D/Dt Applying it to the function v, we obtain

and

a/&i.

(4.9)

Differentiating (4.8) with respect to B, and using (4.9), we have [($)' +(~)']a~(~lG'l').

(4.10)

Eqs.(4.7) and (4.101, along with Hilbert's relation /5/

?.a

for known IG'( and r(g),respectively define the dynamics of periodic gravity-capillary waves on the surface of a fluid of infinite depth. The results of solving system (4.71, (4.8) with 2alh=0.06 are given in Fig.2, where we show the circle p=l, equipped by the vector field (4.7). For clarity, we have added to the field (4.7) the small vector field Dp/Dt=e. The essential differences between Fig.2, a (time t=O.O, no ripple has yet appeared), and Fig.2,b (time t=6.6, a capillary ripple has developed at the frontalslope),are obvious. When the ripple occurs, an unstable point splits into two unstable and one stable points. The situation is then repeated: the new unstable points again bifurcate, the nature of the bifurcation being dependent on the ratio between the curvature derivative and the steepness. Note that, during the dynamic behaviour of a gravity wave, there exist right up to collapse, just two singular points of the vector field D@Dt: unstable on the frontal slope, and stable on the rear slope. REFERENCES 1. KOGAN V.R. and CHURILOV YU.A., Numerical simulation of unsteady gravity-capillarywaves on the surface of a fluid of infinite depth, Izv. Akad. Nauk SSSR, Mekhan. Zhidkosti i Gasa, 22, 2, 124-129, 1987. 2. SCHWARTZ L.W. and VANDEN BROECK J.M., Numerical solution of the exact equation for capillary-gravity waves, J. Fluid Mech., 95, 1, 119-139, 1976. 3. CHEN B. and SAFFMAN P.G., Steady gravity-capillarywaves on deep water, Studies Appl. Math., 60, 3, 183-210, 1979, 4. LAVRENT'YEV M.A. and SHABAT B.V., Methods of the Theory of Functions of a Complex Variable, Moscow, 5.

Nauka,

1965.

MUSKHELISHVILIN.I., Singular Integral Equations, Nauka, Moscow, 1968.

USSR29:3-J

146 6. LONGUET-HIGGINSM.S., The generation of capillary waves by step gravity waves, J. Fluid Mech. , 16, 1, 138-159, 1963. 7. LONGUET-HIGGINSM.S. and COKELET E.D., The deformation of steep surface waves on water, 1, A numerical method of computation, Proc. Roy. Sot. A, 350, l-26, 1976.

Translated by D.E.B.

U.S.S.R. Comput.Maths.Math.Phys.,vo1.29,3,pp.146-X4,1989 Printed in Great Britain

0041-5553/89 $lO.OO+O.OO 01990 Pergamon Press plc

A SOLVABILITY CRITERION FOR THE BOUNDARY-VALUE PROBLEM FOR THE TRANSPORT EQUATION AND SOME COROLLARIES* V.N. LEVKIN

A criterion is formulated and proved for the existence of a non-negativesolution of thetransport equation on cones in K-spaces of real functions. The criterion is then used to derive constructive sufficient conditions for the existence of a non-negativesolution to the problem, a criterion for the existence of a bounded non-negative solution, and other results. Monotone sequences of upper solutions are obtained. Recent studies of the kinetic radiation transport equations /l-3/ have established the existence of solutions under certain conditions. In /l/ constructive sufficient conditions were established for the steady-state radiation transport problem to be solvable in a domain GcC?, with boundary conditions of external irradiation and with a prescribed law of reflection at the outer boundary of G; these conditions are appropriate if G is bounded, a plane-parallel layer or an infinite cylinder. In /3/ the existence and uniqueness of a solution of the unsteady (steady-state)radiation transport equations in U?" with isotropic scattering and the energy equations in a two-temperature approximationwere proved by a method of monotone sequences of upper and lower solutions, proposed in /2/, where the method was used to investigate the existence and uniqueness of the solution of the kinetic equations of a nuclear reactor with temperature feedback in a plane layer. The domain G in which the radiation transport is considered is assumed to be a bounded non-concave three-dimensionaldomain with a piecewise-smoothboundary. In this paper we establish a general criterion for the existence of non-negative solutions of boundary-valueproblems for the steady-state radiation transport equation in a bounded or unbounded domain Gc[RJ subject to external irradiation, where the solutions are to assumed to belong to cones in K-spaces of real functions. The criterion is then used derive constructive sufficient conditions for the existence of non-negative solutions in the K-spaces under consideration. It is shown that if the boundary-valueproblem for the transport equation is solvable, then at least one of the solutions may be obtained by successive approximations. Certain other propositions are stated. The results are used to describe certain methods to derive monotone sequences of upper and lower solutions of the problem satisfying the boundary Condition (2) of Sect.1. 1. Statement of the problem. We consider the solvability of the boundary-valueproblem for the transport equation /4/ A(D=B@+Q, where A=QV+T(~,f~),

"Zh.Vychisl.Mat.mat.Fiz.,29,6,853-866,1989

(Qtr.=+

(1.1)