Energy dissipation in sloshing waves in a rolling rectangular tank Part I: Mathematical Theory Part II: Solution method and analysis of numerical technique Part III: Results and applications

Printed in Great Britain.

T4NI$CHE WilVEnsimi

Meketweg 2, 2628 co DeTft

r

Otb.7 _O$1.1818$J

Ocean Engng. Vol. 10. No 5, pp. 347-35S, 1983. 0029-8018/83 S3.00±.00

ENERGY DISSIPATION IN SLOSHING WAVES IN A

ROLLING RECTANGULAR TANK

- I. MATHEMATICAL

THEORY*

ZEKI DEMIRBILEKt

Research and Development Department, Conoco Inc., Ponca City, OK 74603, U.S.A. AbstractThe hydrodynamic aspects of the motion of a viscous fluid having a free surface in a

rolling tank have been investigated. In a sequence of three papers. an analytical technique together with a numerical solution method will be presented to describe the sloshing

phenomenon accurately and efficiently. This first paper introduces a linear theory of viscous liquid sloshing and formulates a boundary value problem subject to appropriate conditions. Viscosity is included in the problem formulation and its effects are properly accounted for. The second paper will describe a solution of the problem in function space by the truncation cf infinite series. Boundary conditions are satisfied through the usc of Fourier series expansions. However, the no-slip condition at the side walls can also be treated in a least-squares sense.

Among the results that will be reported in the third paper are the effects of viscosity on liquid sloshing phenomenon and the dependence of viscous dissipation on the Reynolds and Froude numbers. Furthermore, the influence of the tank aspect ratio on viscous dissipation has been explored. These results demonstrate some unknown features of the functional relationships that

exist between the dissipated energy and the Reynolds and Froude numbers. Similarly, the

dependence of the dissipated energy on the aspect ratio has been analytically studied. The results obtained agree with the physical laws for the range of parameters investigated.

1. INTRODUCTION

THE PROBLEM of determining the dynamic response of liquids in rigid containers is an engineering problem of great practical importance. In a partially full container, it is well known that the unrestrained free surface has a tendency to exhibit large excursions. This is true for even small excitations of the tank. The motion subsequent to this excitation is termed "sloshing".

Interest in the general problem of liquid sloshing has grown considerably in the past three decades, largely because of its occurrence in various engineering disciplines. This

fascinating subject has attracted the attention of many engineers, mathematicians and

other scientific workers. The importance of sloshing in connection with the developments

of large space vehicles, is well recognized. In spite of extensive work carried out by

numerous researchers for the spherical and cylindrical shape containers that are used in space applications, the related phenomenon of liquid sloshing in prismatic rectangular tanks has received scant attention to date. However, the problem is of growing concern to engineers, since rectangular tanks are in use in maritime applications, such as in LNG

tankers and on offshore structures. They are also considered to be used in orbiting satellite stations in the future. Thus, a mathematical model is needed to describe the

sloshing problem in a rectangular tank.

* This work is_{part of the dissertation by the author submitted to the Graduate college of Texas A&M} University, College Station, TX 77843, U.S.A.

Formerly Graduate Research Assistant, Department of Civil Engineering, Ocean En°incering Program, Texas A&M Untversity, College Station, TX 77843, U.S.A.

347

Pergamon Press Ltd.

An analytic study of the liquid motion in rigid containers has been subject to numerous investigations in the past. Early analytic analyses were crude and rather simplified until the 1950s, when man's interest in space exploration required a very accurate knowledge

of the effect of sloshing motion on the trajectory of the space vehicles. Consequently,

solution methods increasingly became more sophisticated. With the advent of high-speed computers, sloshing motion prediction techniques and computer programs have reached a fairly high level of sophistication. The majority of the present and past work is based on the simplified potential flow theory. Two classes of solutions evolved from this theory.

These are known as the linear and nonlinear solutions. No attempt will be made to repeat these, since there are many excellent surveys available (such as Lamb, 1945;

Stoker, 1957; Kinsman, 1965; Cooper, 1960; and Abramson, 1966). After the pioneering works of Stokes (1847), Levi-Cavita (1925), Struik (1926), and Penny and Price (1952), significant contributions to the nonlinear wave theories were made by Schwartz (1974),

Cokelet (1977), Holyer (1979) and Venezian (1979). Unfortunately, all of the

above-mentioned mathematical theories are again based on the potential flow

formulation and do not take into account the effect of viscosity.

A great deal of the recent nonlinear studies of liquid sloshing are founded on

Moiseev's (1958) theory. The practitioners of this method are Hutton (1963), Chu (1968) and Faltinsen (1974). Although the trends and basic behavior of the flow are predicted correctly by this method, it is not feasible to incorporate the viscosity into the problem formulation. There are as yet no explanations from all of the existing solutions as to the

effects of viscosity on sloshing. This is the question we shall be concerned with and

investigate analytically.

In the past, Graham and Rodriguez (1952) have developed a linear mathematical

model for rectangular tanks subject to roll oscillations. This study is based on the

potential flow theory and later was extended by Vasta et al. (1961) and Chu et al. (1968)

for ship-roll stabilization tanks. We note that antirolling tanks are in common use in

maritime applications as a ship-roll stabilization device which utilizes a passively acting, partially full tank system. In their analytical and experimental studies, Vasta et al. and

Chu et al. concluded that the essential features of fluid damping in the tank must be explored and understood before any significant improvement on the present design

methods can be foreseen. Thus, this paper is aimed to reveal some of the characteristics of fluid damping in a rectangular tank executing roll oscillations.

As far as the nonlinear theories of forced oscillations of liquid in rectangular

containers are concerned, there are only a few theories that are widely accepted. For

large depths, typical of nearly full LNG tanks, the nonlinear perturbation technique of

Faltinsen (1974) is applicable. This mathematical model is an extension of Moiseev's

theoretical work. In his analysis, Faltinsen considers the depth of the fluid to be either 0(1) or infinite and only the small-amplitude roll motions of the container are considered. The study is concerned with the frequency range near the lowest resonance frequency

and is not suitable for the shallow-depth case where hydraulic jumps are formed.

Furthermore, the viscosity is not included in the problem formulation, similar to the

other existing theories.

On the other hand, the shallow-depth case has been investigated by Chu and Ying

(1963), Chester (1968), Verhagen and Wijngaarden (1965) and Wijngaarden (1978). The solution is determined by using a perturbation method. The linear theory of acoustics,

which used. back know In term I the th the ti-artificj The nonliri and th viscou made. ViSCOSi Am Boussj to rept Howe formul Hou fluid. I-nonliru

tobeh

waves An by Boi. waves dampit finite d the vici combin The solid-ru exten& results Keuli finite-ar of wave bounda experirr Miles axisymn rotation viscous More 348 ZEKI DEMIRBILEK

Energy dissipation in sloshing waves - I 349 which is also known as the shallow-water wave theory in the area of hydrodynamics, is used. It predicts a single hydraulic jump of co1stant strength which moves periodically back and forth in the container. The viscous effects are neglected, even though it is well known that energy is dissipated across a jump.

In a recent paper, Faltinsen (1978) introduced, in the Euler equation, a fictitious small term to account for the viscous dissipation in order to resolve the discrepancy between the theory and experiments. This difference is attributed to the omission of viscosity in

the theory. However, in his boundary integral technique, Faltinsen found that the

artificial damping leads to difficulties.

The problem of liquid sloshing, in general, is a nonlinear phenomenon. Two types of nonlinear effects occur. One of these is at the free surface due to the large fluid motion and the other at the fluid-tank interface. In addition to these obvious nonlinearities, the viscous effects are believed to play a very important role. While much progress has been made regarding the nonlinearities, there has been relatively little work on the effect of viscosity in sloshing.

Among the works of early authors who treated the question of viscous dissipation are Boussinesq (1868), Basset (1888), Hough (1897) and Lamb (1945). We shall not attempt to reproduce these studies here, since they are mainly developed for progressive waves.

However, Hough's theory deserves special attention, since it uses the stream function

formulation similar to the present work, as we shall see later.

Hough studied the two-dimensional, unsteady, laminar motion of an incompressible fluid. He introduced the stream function into the NavierStokes equations, neglected the nonlinear convective terms and assumed both the stream function and the fluid velocities to be harmonic in space and time. However, his derivations are good only for progressive waves over a flat bottom without vertical solid boundaries.

An attempt to account for the effects of solid boundaries on wave damping was made

by Boussinesq as early as 1868. His theory extends to both progressive and standing

waves and was used in more recent times by Keulegan (1959). Other earlier work on the damping of waves includes that of Biesel (1949), which deals with waves in a channel of finite depth but infinite width, and that of UrselI (1952), which concerns the dissipation in the vicinity of vertical walls when the depth is infinite. Also, Hunt (1952) has studied the combined effects of finite width and finite depth.

The damping of liquid free-surface oscillations in a circular, cylindrical tank by a flat, solid-ring baffle has been predicted theoretically by Miles (1958). The theory has been extended and modified by Bauer (1962) and is in good agreement with the experimental results but is applicable only to cylindrical tanks.

Keulegan (1959) studied the decay modulus in a rectangular closed basin for the

finite-amplitude waves. Briefly, in this analysis it is assumed that the entire loss of energy

of waves is localized in the boundary layers adjacent to the solid walls. Therefore, a

boundary layer type of solution is implemented and some satisfactory agreement with the experimental results can be observed under certain conditions.

Miles (1963) explored the free-surface oscillations of a liquid about a vertical axis of an

axisymmetric container for small Froude numbers. This work presents the effects of rotation on the free-surface oscillations of a rotating liquid but does not consider the

viscous effects.

350

-ZEKI DEMIRBILEK

surface of an infinitely deep viscous liquid is investigated by Prosperetti (1976). An

integro-differential equation of motion for the free surface is derived and solved.

However, for some time limits, his solution is significantly different from all of the

existing studies, except when the effect of viscosity is very small.

In addition to the analytical solutions cited, there have been other types of approaches

to represent the motion of liquids in containers. The finite-element method and

finite-difference schemes deserve recognition. A brief discussion of these methods is

given by Demirbilek (1982).

It appears from the foregoing that there are no well-established analytical methods to

analyze or predict the characteristics of the sloshing and, in particular, to reveal the important features of the viscous dissipation. Thus, it was decided to undertake a comprehensive theoretical investigation of the liquid sloshing in rectangular containers, with the emphasis focused on the dissipation.

2. MATHEMATICAL FORMULATION

2.1. Problem description

Consider a tank of rectangular shape partially filled with an incompressible and viscous liquid. The configuration variables are illustrated in Fig. 1, and all variables are defined later in this section. Assume that a constant gravitational acceleration is acting along the z-axis. The surface of the liquid then assumes a planar surface normal to this axis, which

we shall call a free surface or quiescent free surface. Superimposed on the constant

gravitational acceleration motion is the oscillatory motion of the tank system. It is this latter motion which causes perturbations or disturbances of the free surface.

In this analysis, it will be assumed that the tank is rigid with vertical walls, and the walls are high enough so that the fluid does not hit the tank topto produce a spraying phenomenon or generate air bubbles inside the domain of the flow field. We shall further assume that the motion of the fluid is two-dimensional in space and harmonic in time.

For the formulation of the problem, in order to include the nonlinearity that occurs at

the fluid-tank interface while avoiding the complex nonlinear boundary conditions, _a

coordinate system attached to the tank is used, as shown in Fig. 1. We note that the tank

coordinate system can be easily related to a Newtonian reference frame with the help of Fig. 1 (c). The tank coordinate system is also known as the "moving coordinate system."

In general, this noniriertial reference frame can have a translational speed and, at the same time, can be rotating about a fixed axis passing through its origin. It should be

pointed out that if the external boundary to the fluid is in motion, then it is often

advantageous to use a moving coordinate system relative to which the boundary is seen at rest. Had a Eulerian frame of reference been used in which the coordinate system is

fixed to the earth, then the tank walls would be moving through

space, yielding complicated boundary conditions. For this reason, it may be convenient to choose a

frame of reference so that the moving boundary is brought to rest.

The fundamental principle governing the equation of motion in _{an inertial frame of} reference is Newton's second law. In vectorial form, this law is identicalto

Dq

_{___VP_gVh*+v*V2q}

a Di p (2.1) FIG. in wli press paran fluid We systec "app (2.1). motio partic

Energy dissipation in sloshing waves - I 351

4

x (a) (ci

I

x (b)

Fio. 1 Tank geometry and coordinate systems. (a) Equilibrium position; (b) instantaneous position; (c) relation of tank coordinate systemLOfixed frame of reference.

in which a represents the absolute acceleration, q the fluid particle velocity vector, p the

pressure in the fluid and v'' the coefficient of kinematic viscosity of the fluid. The

parameter h* is the vertical distance and will be related to the instantaneous position ofa fluid particle as shown in Fig. 1 (c).

We emphasize that when Newton's second law is applied to a moving coordinate

system, its form is identical to that in an absolute frame, provided we include all of the "apparent accelerations" in addition to the absolute acceleration expressed in Equation

(2.1). The apparent accelerations are the results of the translational and

rotational motion of the moving coordinate system. These fictitious accelerations acting on a fluid particle include the following:

Apparent accelerations = translational or "drift" acceleration + centrifugal acceleration + Euler acceleration + Coriolis

352

in which we define these accelerations as

Translational or "drift" acceleration U dt Centrifugal acceleration = - 11 x (11 x r) dfI Euler acceleration = X

r

dt Coriolis acceleration = - 211 x q

where q is the fluid particle velocity vector in the noninertial coordinate system, 11 the

At the

angular velocity, U the translational velocity vector of the moving reference frame and r conditii the position vector relative to the origin of the moving coordinate system. Therefore, we surface have

dU d11

Xr

aapparent =

-dt 211 x q

-

dt

- 11 x (11 x r)

Adding (2.7) to the right-hand side of (2.1) yields the equation of motion for a fluid

particle in a moving coordinate system as

Dq 1

a==__Vp_gVh*+v*V2q_

dU

dt

2flxq

1Ix(flxr). (2.8)

for a tw

This last equation is the law of motion for an incompressible fluid with respect to a as noninertial frame of reference, and thus, Equation (2.8) is the equation of motion for the

present problem.

Conservation of mass is invariant under a coordinate transformation, and therefore, no modification of the continuity equation is necessary for the moving coordinate system.

E The velocities, however, are now relative to the moving frame of reference, and we can quat

write the equation ofcontinuity asunchanged. free-surf

differenc

V . q

= 0 (2.9) second

equi Ii bn

where the particle velocity vector q has the components u and w in the x- and

P With

z-directions, respectively, pressure

A boundary value problem is not complete without the proper boundary conditions. The specification of the boundary conditions plays a very important and often crucial role in theoretical work. In particular, the success of analytical work relies heavily on how

and what kind of boundary conditions ae applied.

In the present problem, the condition of no slip at a solid boundary prevails, since the This corr

fluid is considered to be viscous. This requires under co

i-..-. ZEKI DEMIRSILEK (2.3) (2.7) where i' the oth( bounda. other fe It can

w = 0 at x = ±a

u = 0 at z = h

Energy dissipation in sloshing waves - I

where h is the equilibrium water depth and 2a is the length of the tank as shown in Fig. 1. The normal velocity to a solid, impermeable boundary must be zero. Therefore, we have

u=Oatx=±a

(2.12)

w = 0 at z = h .

(2.13)

At the free surface, we have both the dynamic and kinematic type of boundary

conditions. The kinematic free-surface boundary condition states that a particle on the surface remains on the surface at any later time and can be expressed as

where r corresponds to the free-surface elevation above the equilibrium water depth. On the other hand, the effects of viscosity are manifested through the dynamic free-surface boundary conditions. A very detailed description of these conditions together with the other features of the problem formulation is given by Demirbilek (1982).

It can be shown that the normal stress component reduces to

(P - P0) + 2p v w = 0 at z =

(2.15)

for a two-dimensional flow. Similarly, the tangential stress component can be expressed as

v*(u + w) = 0 at z =

. (2.16)

Equations (2.15) and (2.16) replace the classical Bernoulli's equation as the dynamic free-surface boundary conditions. The first of these equations asserts that the pressure difference between the liquid and air balances the normal viscous shear stress, while the

second equation expresses the fact that the tangential sheer stress must vanish for

equilibrium to exist. We will replace the quantity (P - P0) in Equation (2.15) simply by P, with the understanding that the pressure in the fluid is referenced to the atmospheric pressure above it without any loss of generality. Thus, Equation (2.15) can be rewritten as

This completes the description of all boundary conditions neded to solve the problem

under consideration.

I a1

+u

=watz=

(2.14)

354 Ziu DEMIR1LEK

2.2. Streamfunction formulation

Two-dimensional flows can be conveniently characterized by introducing a mathema-tical artifice known as the stream function, J(x,z,t). This function relates the concept of streamlines to the principle of mass conservation. The equation of continuity in (2.9) is

identically satisfied if we define where

i the

q = (u,w) =

JxV1,

= (LjJ2,fl) . (2.18) conject

For a two-dimensional incompressible flow, the stream function formulation can be used

to express the combined continuity and momentum equation. Taking the curl of

Equation (2.8) and using (2.9) yields

where the subscript means partial differentiation with respect to the variable indicated. Equations (2. 19)(2.26) form a complete set of formulas representing the motion of a liquid in a rectangular tank executing roll oscillations.

Since the objective of this work is to obtain a linear or first-order solution, we shall

neglect any quadratic, i.e. nonlinear, terms which appear in the above equations.

Furthermore, the free-surface boundary conditions can be transferred to the equilibrium

position, z = 0, by expanding them into a Taylor series for the linearized solution

without any loss of generality. Also, the pressure term in Equation (2.26) can be

eliminated with the help of the momentum equations (see Appendix A). Moreover,

considering the response of the fluid only to harmonic oscillations with frequency o, then 0 and

ti

can be represented as

- 0 at x

=

±a

(2.20) = 0 at z = (2.21) = 0 at x = ±a (2.22) = 0 at z =

h

(2.23)

-

1I5zuix = at z = (2.24) v*(jI

-

_{) = 0 at z = (2.25)}

_P+2PV*I,,xz=0atZTI

(2.26) V21,,

-

4V2

+ i1jV2lj

2 + v*V44i

(2.19)

where = and 0 is the tank roll angle.

We note in passing that Equation (2.19) is a nonlinear partial differential equation of order four and results from the momentum equations, i.e. the NavierStokes equations.

It will be the objective of the second paper of this series to obtain a linear solution of

Equation (2.19).

Upon substituting Equation (2.18) into the boundary conditions of the problem, one could express all the conditions in terms of the stream function. This process yields

in whic depend the gov of a vis where The La respect

where 13 is the amplitude of roll excitation, i.o the angular frequency of the oscillation and

i the space-dependent part of

which, in general, is complex. It can be easily

conjectured that the free-surface elevation ii and pressure P have the form

-in which ii and P are complex and represent the parts of and P that have spatial

dependence only. Substitution of these relations yields the following linearized version of the governing equation of motion and its associated boundary conditions for the sloshing

of a viscous fluid in a rectangular tank:

where =

_{, K =}

2w13 ; K1 = +

od)

(0

o = i 13

eio)t

= Ii

e' V2

+ ivV4

i4i = K

Ii = 0 at z = h

v

= 0 at z = h

LTJz = 0 at x = ±a v

= 0 at x = ±a

v(i)=0atz=0

LlJz + iv(3i +

'1') +

2

Energy dissipation in sloshing waves - 1 355

(2.27) (2.28) a2 + a4 a4

+2

22+4

ax az az Iixx = K1 at z = 0 (2.31) (2.32) (2.33) (2.34) (2.35) (2.36) (2.37) (2.39) The Laplacian and biharmonic operators which appear in Equation (2.31) are defined, respectively, as =

e"

(2.29)

P = i P e)t

(2.30)

at z = 0

(2.38) (2.40) (2.41)

356 ZEKI DEMIRBILEK

This completes the problem formulation. To present a solution technique for the

above boundary value problem will be the next objective of this work. For more detailed information on the preceding derivations, the reader is referred to Demirbilek (1982).

3. SUMMARY

The sloshing motion of a viscous fluid inside a rectangular container has been

formulated. Viscosity is included in the present formulation whereas the other existing approaches neglect this essential property. It is impossible to study the dissipated energy analytically without the presence of viscosity. Furthermore, the above novel formulation permits one to directly recover the necessary equations for an inviscid problem by simply setting v = 0 in these final equations. For instance, Su and Demirbilek's (1981) inviscid formulation could easily be obtained by this process.

A method of solution for this linearized boundary value problem will be discussed in a subsequent paper.

REFERENCES

ABRAMSON, H.N. (Ed.) 1966. The dynamic behavior of liquids in containers. NASA SP 106.

BASSET A.B. 1888. A Treatise on Hydrodynamics. Deighton, London. Also, Dover, New York, 1961.

BAyER, H. 1962. The damping factor provided by flat annular ring baffles for free surface oscillations.

NASA-MSF, MTP-AERO-62.

BIESEL, F. 1949. Calcul de lamortissement d'une houle dans un liquide visqueux de profondeur finie. Houille blanche 4, 630.

BoussINEsQ, J. 1868. Memoir on the influence of friction in the regular motion of fluids. J. Math. 2, 377. CHESTER, W. 1968. Resonant oscillations of water waves. Proc. R. Soc. 306A, 5.

CHU, W.H. 1968. Subharmonic oscillations in an arbitrary tank resulting from axial excitation. I. app!. Mech. 35, 143.

Ciju, B.T. and YING, S.J. 1963. Thermally driven nonlinear oscillations in a pipe with traveling shock waves. Physics Fluids 6, 1625.

CHU, W.H., DALZELL, J.F. and M0DISETrE, J.E. 1968. Theoretical and experimental study of ship-roll stabilization tanks. J. Ship Res. 12, 165.

COKELET, ED. 1977. Steep gravity waves in water of arbitrary uniform depth. Phi!. Trans. R. Soc. 286A, 183. COOPER, R.M. 1960. Dynamics of liquids in moving containers. I. Am. Rockei Soc. 30, 725.

DEMIRUILEK, Z. 1982. A linear theory of viscous liquid sloshing. Ph.D. thesis, Texas A&M University, College Station, Texas.

FALTINSEN, O.M. 1974. A nonlinear theory of sloshing in rectangular tanks. J. Ship Res. 18, 224.

FALTINSEN, O.M. 1978. A numerical nonlinear method of sloshing in tanks with two-dimensional flow. I. Ship Res. 22, 193.

GRAHAM, E.W. and RODRIGUEZ, AM. 1952. The characteristics of fuel motion which affect airplane dynamics. 1. appl. Mech. 19, 381.

HOLYER, J.Y. 1979. Large amplitude progressive interfacial waves. J. Fluid Mech. 93, 433.

HOUGH, M.A. 1897. On the influence of viscosity on waves and currents. Proc. Lond. Math. Soc. 28, 264. Hur, J.N. 1952. Amortissement par viscosite de Ia houle sur un fond incline dans un canal de largeur finie.

Houille blanche 7, 836.

HuTr0N, R.E. 1963. An investigation of resonant, nonlinear free surface oscillations of a fluid. Ph.D. thesis, Univ. of California at Los Angeles.

KEULEGAN, G.R.L. 1959. Energy dissipation in standing waves in rectangular basins. 1. Fluid Mech. 6, 33. KINSMAN, B. 1965. Wind Waves. Prentice-Hall, Englewood Cliffs.

LAMB, H. 1945. Hydrodynamics, 6th Edo. Dover, New York.

LEVI-CAVITA, T. 1925. Determination ngoureuse des ondes permanentes d'ampleur finie. Moth. Annln 93.264. MILES, J.W. 1958. Ring damping of free surface oscillations in circular tanks. J. app!. Mech. 25, 274. MILES, J.W. 1963. Free-surface oscillations in a slowly rotating liquid. J. Fluid Mech. 18, 187.

MoIsEEv, N.N. 1958. On the theory of nonlinear vibrations of a liquid of finite volume. App!. Math. Mech. N.Y. 22, 612.

PENNY, W.G. and PRICE, A.T. 1952. Finite periodicstationary gravity waves in a perfect liquid. Phil. Trans. R. Soc. 244A, 254.

PROSPERETT1, A. 1976. Viscous effects on small amplitude surface waves. Physics Fluids 19, 195.

SCHWARTZ, L.W. 1974. Computer extension and analytic continuation of Stokes' expansion for gravity waves. J. Fluid Mech. 62, 553. STOKER, STOKES, STRUIK, finie. Su, T.C. nonrc URSELL, VASTA, 1. Trans. VENEZIA, Hono VERHAGE 22,73 WUNG,t, Ooste The pr equatio It is n Thus, af to x an Note th moment of mom Solving function The v referrin where 0 Furthe located motion

I

APPENDIX A

The pressure term which appears in Equation (2.26) can be obtained from the momentum equations. The following is a brief description of this objective.

It is noted that Equation (2.6) is applied at a fixed value of z but must hold at every value of x. Thus, after eliminating the time dependence and linearization, one can differentiate it with respect to x and combine it with Equation (2.24). This yields

p ci)

(A.1) Note that the first term in Equation (A.1) is explicitly defined by the x-component of the momentum equations. After we linearize Equation (2.8), which is a statement for the conservation of momentum in vector form, it reduces to

1

ii,

= - - P gh.' + v V2u -

ii

-

z (A.2)

p 1

w,=--Pgh+v'V2w iJ3.lx

(A.3)

p

Solving for - i/p P from Equation (A.2) and expressing the velocity in terms of the stream function yields

- -- P = -

+ gh + v

V2I. + J at z = 0 (A.4)

The vertical distance h* can be related to the instantaneous position of a fluid particle by referring to Fig. 1(c). It can be shown that

= z cos 0 - x sin 0 + j U3 cit (A.5)

where 0 is the tank roll angle. The angle 0 is related to the angular velocity of oscillations fl by

dO 0

fl=jfl=j

_dt

Furthermore, the excitation functions lJ and 113 can be obtained as follows: if the center of roll is located at a distance d below the equilibrium position of the tank and the angular velocity of roll motion is fi, then the linear velocity of the tank trajectory with respect to an inertial system is

(A.6)

Energy dissipation in sloshing waves - I 357 STOKER, J.J. 1957. Water Waves. Interscience, New York.

STOKES, G.G. 1847. On the theory of oscillatory waves. Trans. Camb. phil. Soc. 8, 441.

Sntui,c, D.J. 1926. Determination rigoureuse des ondes irrotationelles pCriodiques dans un canal a profondeur firiie. Math. Annln 95, 595.

Su, T.C. and DEMIRSILEK, Z. 19S1. A nonlinear analysis of liquid sloshing in a rectangular tank subject to nonresonant roll excitations. Rept. No. COE-223, Texas A&M Univ., College Station,Texas.

URSELL, F. 1952. Edge waves on a sloping beach. Proc. R. Soc. London 214A, 79.

VASTA, J., GLDDINCS, A.J., TAPLIN, A. and ST!LWELL, J.J. 1961. Roll stabilization by means of passive tanks. Trans. Soc. flat. Archil. mar. Engrs, N.Y. 69, 411.

VENEZ1AN, G. 1979. Hydrodynamic wave-current and pressure. Technical Report No. 46, Univ. of Hawaii, Honolulu, Hawaii.

Va1utGEN, J.H.G. and WUNGAARDEN, V.L. 1965. Nonlinear oscillations of fluid in a container. J. Fluid Mech. 22, 737.

WUNGAARDEN, V.L. 1978. Nonlinear Acoustics. In Symp. on. Appi. Math. Edited by Hermans A.J.and Oosterveld M.W.C. Deift, The Netherlands.

g3

+wdatz=Q

Equaton (A.11) is the final form for the normal stress condition at the free surface. Similar to the remaining boundary conditions, it is now in the proper form to be used.

Ocean Printed j THIS PA. prob!err these e maybe 1.1. M In vie which a of satisf given b3 simple where 4 expresse providec In termS p

Equations (A.1) and (A.1O) can be combined to yield

+ iv(V2. +

_{+ -s-}

_{= K1 at z = 0}

Ki=B(---+d).

(A.1O) (A.11) (A.12) Al fo so ob pr ph Co in' de V(t) = . (A.7)

Therefore, the accelerations of the tank in the x- and z-dircctions, respectively,are

&=.

c'

(A.8)

°

d

--fl2d.

(A.9)

After inserting the necessary parameters into the right-hand side of Equation (A.4) and factoring out the time-dependent part, one obtains

= -

- iv V2

i. Ocean Engng, Vol. 10, No 5, pp. 359-374, 1983. 0029-8018/83 $3.00+.00

Printed in Great Britain. Pergamon Press Ltd.

ENERGY DISSIPATION IN SLOSHING WAVES IN A

ROLLING RECTANGULAR TANK

- II. SOLUTION

METHOD AND ANALYSIS OF NUMERICAL TECHNIQUE

ZEKI DEMIRBILEK

Research and Development Department, Conoco Inc., Ponca City, OK 74603, U.S.A. AbstractIn Part I of these series of papers, the complete problem formulation in a linearized form was presented. In order to provide the engineer with an in-depth knowledge about the exact solution of the problem, it is natural and essential to start with a linear solution. l'his will be the objective of Part H. together with an exposition to the analysis of numerical technique utilized. A truncated infinite Fourier series-type solution is adopted for the linearized boundary value

problem. It is shown that such a solution is mathematically consistent and represents the phenomenon properly by satisfying all of the field equations and the imposed boundary conditions. The dependence of the Fourier coefficients on the truncation limit has been

investigated. The best lower and upper "cutoff limits" for the truncation of an infinite series are determined. An error analysis of the solution technique is performed.

1. INTRODUCTION

Ti-us PAPER iS devoted to the description of a mathematical solution of the boundary value problem summarized in Equations (2.31)(2.3S) by Demirbilek (1983). We shall refer to

these equations without attempting to reproduce here. Therefore, the present material

may be considered as a continuation of this earlier work.

Li.

Method of solution

In view of the existence of symmetry in the i-direction, we shall seek solutions of i which are even in x. These solutions must be sufficiently general so that they are capable

of satisfying all of the boundary conditions. Although the governing Equation (2.31)

given by Demirbilek (1983) is nonhomogeneous, it can be made homogeneous with a

simple transformation. Let

= + 4)

where 4) corresponds to the homogenous solution of the partial differential equation

expressed in (2.31). From Equation (2.31), one can show that 4) must satisfy

V24)+ivV44)=O

(1.2)

provided that the particular solution, tJ, is chosen as

= Kx2

. (1.3)

In terms of 4), the final problem formulation can be outlined as follows:

V24)+ivV44)=O

(1.4)

360

- (4

+ Kx) at z = 0

where

K2 -

f3(2 d

- g)

(1.12)

This final formulation of the problem will be used henceforth. Next, we shall introduce a direct and versatile solution technique to the governing Equation (1.4).

It

is noted that as the tank is set in motion, the fluid motion produced will be

two-dimensional and wavelike. Each single wave is described by a certain wave number and frequency associated with it. For a progressive wave, the wave number is defined in the direction of wave propagation. In a plane flow, the waves produced inside the tank,

in general, have different wave numbers in two different coordinate directions.

Therefore, since represents the waves produced inside the tank, then it can have a wave either in the form of e coskx or e cosXz. The first type of wave is a possibility provided that k and cc are related to each other in a manner uniquely determined by the governing equation of motion. It can be shown that this relationship is imbedded in an indical equation as

r+ivr2O

(1.13)

where

r=a2k2

(1.14)

We note that the indical Equation (1.13) can identically be satisfied with two

possibilities. The first of these is

r=Oora2k2

(1.15)

and the second possibility requires

2 £

r= or-a =

IC +

-V V ZEKI DEMRBTLEK

= Kx at x = ±a

(1.5)

= 0 at x = ±a

(1.6)

&=Kxatz=h

(1.7)

4=Oatz=h

(1.8) (1.9) + i v w2(3

+ 4) g'xx = K2 at z = 0

(1.10) (1.16) where Equati Simi where and X throug We Equati to x de by sim can be Fourie such e (1.4) f Since accordi superp in whic numbe

a and

Ther of view comput avoide normal be esse the fina

I

and X is the wave number in the z-direction and p. is a complex parameter related to X through Equation (1.19).

We note in passing that the reason for discarding the sinhAx and sinhix terms in

Equation (1.18) is due to the symmetry property of the boundary conditions with respect to x dependence. Furthermore, it is obvious that one can have another possible solution by simply replacing cosXz term with sinXz in Equation (1.18). This possibility, however, can be shown that is not needed because of the so-called half-range expansions of the Fourier series. The reader is referred to Churchill (1941) for some intrinsic properties of such expansions. In conclusion, we shall construct a more general solution of Equation (1.4) from a combination of the two possibilities given in Equations (1.17) and (1.18). Since the system is linear, there is no loss of generality in adding these two solutions, according to the principle of superposition. Therefore, a more general solution of4) is a

superposition of the two families of solution given above. This can be written as

4) = (A,, sinhk,,z + B,, coshk,,z + C,, sinha,,z +

D,, cosha,,z) cosk,,x + (E,,, coshXmx +

m

Fm coshp.,,,x) cosXmz (1.20)

in which A,,, B,,, C,,, D,,, E, and Fm are some complex constant coefficients. The wave numbers in the x- and z-directions are k,, and Xm, respectively. The complex parameters a,, and p.,,, are defined in Equations (1.16) and (1.19).

There is a need to normalize each term of Equation (1.20) from a computational point of view. As the arguments of the hyperbolic functions reach a moderate value, serious computational difficulties are to be expected. These numerical problems can be partly avoided by dividing each term of the two infinite sums in Equation (1.20) with a proper normalization factor. However, despite the use of normalization factors, at times it may be essential to introduce the asymptotic expansions for hyperbolic functions. Therefore, the final form of4)in view of the above considerations becomes

Energy dissipation in sloshing waves - II 361

where k is the wave number in the x-direction and a is a complex quantity defined in

Equation (1.16). Then, with a wave of the first type, we could have

4) = (A1 sinhkz + B1 conshkz + C1 sinhaz + D1 coshaz) coskx . (1.17)

Similarly, it is straightforward to show that the second-type wave as a possibility yields

4) = (E1 coshXx + F1 coshp) cosXz (1.18)

where

p.2 = 2+ (1.19)

I

I

and

The answer to the question of whether a given problem is well posed or not can be

obtained by implementing its boundary conditions. We start with the vanishing of normal velocity conditions at the side walls given by Equation (1.6), which requires

>Z(z) cosk,,a -

X,,,(a) >m 5IflXrnZ = 0 at any z. This can be accomplished simply by enforcing

(1.24)

where

cosk,1a = 0 (1.25)

Xm(a) = 0 . (1.26)

We note that Equation (1.25) defines the wave numbers in the x-direction while

Equation (1.26) yields a relationship between the constant coefficients Em and F,,,

respectively, as

Similarly, no flow can penetrate the bottom of the tank. Equation (1.7) is an

expression of this statement, and it demands that

Z,,(h) k,, sink,,x + >X,,(x) CO5hXmh = Kx

(1.29)

at any x. It is preferable to select Equation (1.29) as the relation for defining the wave numbers in the z-direction. This could be achieved if one chooses

Insej Equ velocit' for all .z Theref k

(2n - 1)

n = 1, 2, 3, ...

. (1.27) form a 2a one to

Fm = Em

(1.28) fUflctio functio This re multipi 362 ZEKI D1IRBILEK = >Z,,(z) cosk,x + >X,,,(x) cosX,,,z where we defined (1.21) which 1

Z,,(z) = (A,, sinhk,1z + B,, coshk,,z) +

coshk,,h 1 (1.22) The re D,, (C,, sinha,,z + cosha,,z) cosha,,h and coshX,,,x

X,,,(x) = Em + F,,, coshli,,,x (1.23) Howe_intervE

I

I

Energy dissipation in s!oshing waves - II 363

COSXmh = 0 (1.30)

which renders

in = 1, 2, 3, ...

(1.31)

The relation (1.29), therefore, simplifies to

However, the Fourier series expansion of the functionf(x) = x with period 4a within the

interval 2a

x 2a is

(2nz-1)ir

m

2h

8a(-1)''

(2i - 1)irx

(2n - 1)21T2 2a

Inserting (1.33) into (1.32) after utilizing (1.27) and rearranging, we get

A,, tanhk,,h - B,, + C,, tanhe,2h -

= -s,,

where

2K(-1)

-

ak,,3 >Z,, (h) k,, sink,7x = Kx . (1.32)

f coskcosk,,th =(

a if,i = m

0 otherwise

a

-(1.33) (1.34) (1.35)

Equation (1.8) implies that on the bottom solid wall which is at rest, the tangential

velocity must be zero (the "no-slip condition"). This gives

>Z,(h) cosk,,x + >Xm(X) Xm sinX,,1h = 0 (1.36)

for all x. A quick glance at Equation (1.23) indicates that X,,,(x) is an even function of x.. Therefore, it can be expanded into a Fourier series in terms of cosk,,x. Since these bases

form a complete set of orthogonal functions, then the orthogonality principle permits one to tackle the relation given in Equation (1.36). We shall remember that the set of

functions cosk,,x is said to be orthogonal with respect to the interval (a,a) if these

functions satisfy the relations

(1.37)

This remarkable feature of the orthogonality principle will be used henceforth. Thus

i

Y5(n) + D,, Y6(n) = Y7(n) - >.10(n,rn) Em V _m

A,, Y8(n) + B,, Y9(n) + C,, Y10(n) + D,, Y11(n) = i

Y12(n) - - 'y(n,m) Em

(1.40)

V

>Yi3(m,n) A,, + Y14(m,n) B,, + Y15(m,n) C,, +

D,, Y16(m,n) Em - Co(m) (1.41)

where for the description of the above functions, Appendix A can be consulted. An

elegant mathematical proof is given by Demirbilek (1982) to verify that a least-squares approach for the no-slip condition at the side walls yields identical results to the Fourier series treatment.

The expressions derived from the boundary conditions, namely, Equations (1.34),

(1.38) and (1.39)(1.41) constitute a linear system in the complex domain whose solution will determine the only remaining unknowns of the present study. These are the constant coefficients A,,, B,,, C,,, D,, and Em. The following section briefly explores the main

features of the numerical analysis of this system.

2. THE ANALYSIS OF NUMERICAL TECHNIQUE

2.1. Introduction

The mathematical model developed in this investigation has been programmed for execution on a digital computer. The analytical expressions obtained in the previous

sections are solved by a numerical technique which will briefly be discussed next.

Numerical methods are commonly based on approximations. Clearly, all of the

numerical schemes contain some inherent errors, since a digital computer can perform the arithmetic to only a certain degree of accuracy with a limited number of operations. Therefore, even if an error-free mathematical model could be developed, it could not, in general, be solved exactly on a digital computer.

We are concerned here with the accurate solution of a system of simultaneous linear algebraic equations. Early observations and computations prpved that the matrix is very

close to being singular and thus render the system on the verge of becoming

(1.39) ill-corid pivotin Compu double-the syst without literatu. 1967). j obtaine It is e in ordei index n MM tet system is faced values c the cho such as would Ii A clai from th( B, C,,, they mu the class through Fourii nonperi determi general, Fourier of overs Fourier times wi curious oscillatic

'It shc

converg to a fini accuracy Section everywh magnituc coefficie magnitw Howevei the 364 ZEKL DEMIRBILEK

Y1(n) - B,, Y2(n) + C,, Y3(n) - D Y4(n)

--t-->i3(n,m)Em

(1.38)

yin

in which all of the functions appearing in the above expression are defined in Appendix A. It should further be pointed out that a similar procedure can be followed to satisfy the remaining boundary conditions. It suffices to give the results, and the reader can find all

of the essential details in Demirbilek (1982) if necessary. Be this as it may, it can be

I

Energy dissipation in sloshing waves - H 365

ill-conditioned. Hence, it was decided to use the Gaussian elimination with partial

pivoting technique for solving the ill-conditioned system of equations. The analysis and

computations are based on the complex domain matrix operations using a

double-precision arithmetic. In general, the adopted technique gives either a solution of the system or alternatively indicates that the system is too ill-conditioned to be solved without working with a higher-precision arithmetic. This method is well known in the

literature, and no attempt is necessary to explain its features (see Forsythe and Moler,

1967). An iterative scheme suggested by Reid (1981) was also used to verify the results

obtained by the first method. Both methods, in general, were in good agreement.

It is essential to truncate the infinite series appearing in our linear system of equations

in order to obtain the desired values of the unknown coefficients. Let the series with

index n be truncated at NN terms, while the other series is truncated with the index m at

MM terms. Depending on the magnitudes of the truncation terms NN and MM, the

system could easily be made as small or as large as we wish. In view of this, one naturally is faced with the following questions: What is the criterion for deciding on the optimum values of NN and MM? What effects do the truncation terms have on the solution? Are the choices of NN and MM related to some other physical parameters of the problem, such as the Reynolds or Froude numbers? These are only some of the questions that we would like to address and investigate in the following paragraphs.

A clarification needs to be made in connection with the Fourier coefficients. It is noted from the form of the solution given in Equation (1.21) that the constant coefficients A,

B,2, C,,, D, Em and F,,, alone should not be considered as the Fourier coefficients. Instead, they must be grouped together with their respective normalization factors so as to fit into the classical definition of the Fourier coefficients. It is important to keep this fact in mind

throughout the remainder of this work.

Fourier series, in essence, is an approximation of the exact value of a periodic or

nonperiodic function. In the Fourier series representation of a function, it is of interest to

determine how well the first few terms of expansion represent the given function. In

general, the accuracy of the representation improves with an increase in the number of Fourier coefficients taken. However, if the function has a discontinuity, then the amount

of overshoot at a jump or discontinuity does not approach zero, no matter how many Fourier coefficients are used to better its representation. In fact, the approximation at

times winds itself around the true value in the form of high-frequency oscillations. This

curious behavior is known as the "noise" or the "Gibbs phenomenon" or the "Fejer

oscillations" (Lanczos, 1964, 1966).

It should be remembered that the Fourier series

is an infinite series, but its

convergence properties permit reasonable engineering approximations with truncation

to a finite number of terms. How large this number of terms must be for a desired accuracy depends on the nature of the function being approximated. Lanczos (1961,

Section 2.6) points out that for the trigonometric Fourier series, if the function is

everywhere Continuous and its m-th derivative is discontinuous at some point, the magnitudes of the Fourier coefficients decrease at a rate

_z±), where

ii is the

coefficient number. In the absence of the noise, after some value of n, say NN, the magnitudes of the coefficients are small enough so that the series may be truncated.

However, in the presence of the noise, the coefficients, say a,,, for n > NN depart from the m+fl relationship and may no longer be negligible, nor do they have a tendency to

(2.1)

366 Zri DEMIRBILEK

diminish. These coefficients represent the nonconvergent part of the Fourier series caused by the nonanalytical behavior of noise, i.e. the Gibbs oscillations. The point at

which a departs from the is usually referred to as the "cutoff frequency"

discussed by Lanczos (1964, pp. 336 and 337). The following discusses this and another method for determining the cutoff frequency. A study of the FourierBessel coefficients using Lanczos' technique was undertaken by Riley (1969).

Two principles were used in determining the cutoff values. The first is based on the above mathematical theory of Lanczos, and the second is the analysis of the magnitude of the error involved in fulfilling each boundary condition of the problem.

Using Lanczos' idea, the dependence of the Fourier coefficients on the coefficient number can be shown clearly in a loglog plot. We note that the (n"') relationship

indicates a straight line with a negative slope on the loglog plot. Therefore, any

deviation from the straight line behaviour can be considered as an indication of the noise phenomenon. This can occur both at some small values of the coefficient number as well

as at

its large values. The first

is assumed to determine the minimum number of

truncation terms, whereas the second is assumed to predict the upper bound of the truncation value. The points at which the Fourier coefficients start to depart from a straight line behavior were decided to be he representatives of the cutoff values.

It

is useful and essential that we shall now introduce the Reynolds and Froude

numbers, which are the two most important parameters appearing in the study of

free-surface flows. For liquid sloshing in a rectangular tank, we shall define these

parameters as a2w Reynolds number = r = _v* aco2 Froude number = F = (2.2) g

where a is the half-tank length, w the frequency of excitation, v* the coefficient f kinematic viscosity of the fluid and g the gravitational acceleration. These dimensionless numbers play a very important role in the sloshing phenomenon. Therefore, it was not

surprising to observe the dependency of the cutoff values (NN and MM for the two series FiG. 1. truncated) on the Reynolds number for several examples studied. Intuitively, only a few

terms are needed to properly represent the flow field when Reynolds number is small:

On the contrary, when Reynolds number is large, the fluid displays more movement. There is a relatively more visible motion of the fluid particles near the bottom of the

container. Thus, one may need more terms in the solution for an accurate description of 2.2. the flow field. It should further be stressed that in view of Prandtl's (1904) boundary layer Base concept, the boundary layer developed near the solid walls is considerably thinner in the coeffici case of large Reynolds numbers. Therefore, it is essential to have a sufficient number of parame

Fourier terms in the solution in order to get a better resolution of the boundary layer Figu itself. Many numerical test runs have proved this expectation to be true. It should also be these p pointed out that no significant effect of Froude number on the Fourier coefficients was value o

I I 6.0 5.0 4.0 3.0 8 '-_..-' 2.0 0' 0 I.0 10

Energy dissipation in sloshing waves

- H

5 367 -2.0 (a) _(b) 2.1 -2 -2.2 a 0. -2.3 a 2.4 -U

\0

2$ -0 0

--

" 26 -0 _-2.7 .2 0 ₀ -2.8 C S -2.9 S _-3.0 -3 I 0 I I 2 2 _{0.2 0.4 0.6 0.8 I.0} '7 _{Log In)} (c) 5 L (dl 0 5-S 2 U

'

0 I0 -1.0 0.2 0.4 0.6 100 2 5 101 S ₀₂ Log In) A

Fic. 1. _{Fourier coefficient A,, vs coefficient}

number for hI2a = 0.50 and (a) R _{0.01. F 0.001; (b) R = 10.0,}

F= 0.01; (c) R = 10.0, F = 0.03; (d) R = 10.0, F= 0.601.

2.2. Results0/ithe behavior of Fourier coefficients

Based on the above observations, the

_{magnitude of the Fourier coefficients} vs

coefficient number is used to determine the truncation limits. The results of this

parametric study are depicted in Figs 1(a)-(d) and 2(a)-(d).

Figure 1(a) is a plot of AI/coshkh_{vs n in log-log scale showing the relation between} these parameters. It indicates the probable cut-off values of n, called NN. for _{a fixed} value of the Reynolds and Froude numbers and an aspect ratio, h/2a. _{Lanczos' theory}

368 0 U , 0 ZEKI DEMIRBILEK (b) (dl

Fin. 2. Fourier coefficients A,, and E,,, vs coefficient numbers for h/2a = 0.50 and (a) R = 2200.0, F = 0.001;

(b) R = 10.0, F = 0.03; (c) R = 10.0, F = 0.001; (d) R 2200.00, F 0.001.

the Fourier coefficient and coefficient number. Any deviation from this straight line is an indication of the noise or Gibbs phenomenon, and the point of its occurrence should be considered as a cutoff value for the number of terms required to truncate the series. It is

noted that in Fig. 1(a), the straight line ends at about n = 4 and oscillations start

occurring afterwards. This value of n is taken as the cutoff value, NN, for the given

nondimensional parameters. It was observed that in Some cases, due to large variations in the magnitude of the Fourier coefficients, a log-log plot would be impractical. Thus,

instead, one first needs to take the logarithm of the quantities and then plot these in a

regular cartesian scale. Figs 1(b), (c) and 2(b)-(d) are plotted accordingly. Obviously, all

the plots of the Fourier coefficients vs coefficient number do not exhibit a very clear

straight straight points. lines, in cannot checked best va1 It caii Lanczo to check a series The f that the compar more te. value of through of more the sens conditio displaye Lanczos solution need fo Reynoic determii Table 1. observec 6 -1.0 0 _(a) -3.0 .% ..-.

-

'S.. 'So 0% S. 15' - 5,. -7.0 -9.0 -11.0 IÔ5 I I I I -130 I I 00 2 G n 2 4 6 0 _0.2 0.4 0.6 log Cm) 0 7.0 (C)

.\

-30 S -9.0 'S -4,0 - 'S._"So -11.0 'C 'S 'S 'S. -6.0 _S. 3.0 "O

-0

-15.0 -lao -17.0

--l2.0 -19.0 0 -14.0 -2I.0 C -160 I I I I -23.0 I I 0 0.2 0.4 0.6 08 1.0 I2 0 0.2 0.4 0.6 0.8 1.0 12 0.8 I.0 12

I

Energy dissipation in sloshing waves - II 369

straight line trend. In fact, Figs 2(b)(d) seem to indicate that it is possible to draw two

straight lines in each figure. These two straight lines form an envelope of the computed points. It was interesting and quite surprising to observe that the intersection of the two lines, in general, served as a "good" estimate for the cutoff values. This point obviously

cannot be substantiated on the basis of Lanczos' theory. Such cases in particular were

checked by a second method to see whether the cutoff limits determined were in fact the best values.

It can be difficult to explain the behavior of the Fourier coefficients on the basis of

Lanczos' theory when two straight lines appear in each plot. In such dases, it was decided to check and verify the results obtained by Lanczos' method. The implication of truncating a series at a finite number of terms served as a cornerstone for the second method.

The fundamental principle in truncating a series at a finite number of terms requires

that the magnitudes of the indivual terms following the truncation limit be smaller

compared to the partial sum of the terms up to that limit. In other words, the addition of more terms beyond the truncation limit should not have a considerable influence on the

value of the partial sum of the truncated series. Thus, if the cutoff values determined

through Lanczos' theory are correct, then the results must be insensitive to the inclusion of more terms. With a gradual increase in the number of terms beyond the cutoff values, the sensitivity of the results can be investigated. This check was done on the boundary

conditions, velocity field and the dissipation function. If no significant variations are displayed, then it Was decided that the criterion of determining the cutoff values by Lanczos' theory is accurate and justified. Otherwise, the process is repeated until the

solution is not influenced by the trunôation limit. Computations showed that there was a

need for such a trial only in limited cases. These cases usually corresponded to the Reynolds number greater than iO. Based on these observations, the truncation limits determined for the ranges of Reynolds numbers at different aspect ratios are listed in Table 1. Furthermore, no such pronounced dependence of the truncation limits was

observed in connection with the Froude numbers tested.

TABLE 1. TRUNCATION CRITERIA FOR DIFFERENT REYNOLDS NUMBERS AT VARIOUS ASPECT RATIOS

Aspect ratio (h/2a) Reynolds number (R) NN MM 0.25 10_2 _{R 9 x 10' 4 5}

R102

4 10 0.5 10_2 R 9 x 10 4 5

10°R102

5 5

2x102R103

5 10 R>103 5 15 1.0 10 R

9 x 10'

4 5

Ri0°

4 10 2.0

R102

2 5

70 ZEKI DEMIRBILEK

The truncation values, NN and MM, for the two series in Table 1 show a noticeable instead i

dependence on the Reynolds numbers for a fixed value of the tank aspect ratio.

able to

However, this dependency seems to be irregular and does not follow a clear pattern, as boundal the aspect ratio varies. Therefore, Table 1 can be used as a guide for this purpose. This error cri

table is constructed based on the studies of the Fourier coefficients AI/coshk,7h and

I En,/OShXrn in Figs 1 and 2. Similar vork was also carried out for the other Fourier 2.4. Th

coefficients. However, this did not reveal any different or additional information. Some roundof difficulf It was decided that an error of 10-6 or less in the boundary conditions is satisfactory It is for the solution to be acceptable. Table 2 displays a typical case of the computed errors in point o the boundary conditions after a decision on the truncation terms is reached. Clearly, the discrepa

orders of magnitude of these errors are quite small and well within the range of

numeric

acceptance. Usually, the largest errors resulted as Reynolds number approached either a circumv

value of less than 10

or exceed 5 x iO3. For the Reynolds number within these

descript

indicated values, errors in the boundary conditions were not greater than 10_6.

First

However, outside these limits, the order of magnitude of the error terms was about iO4. our sol

accorda 2.3. Verifi cation of analytical solution

limits o sense, i they ca comput limitati Secon frequen either s simplest progra reductio troubles before process in terms gradient tionally greatly analytic can caus of magn

TABLE 2. COMPUTED ERRORS IN THE BOUNDARY CONDITIONS FOR R = 0.625, F = 0.001 AND h/2a

= 0.25

In fact, during the computations it was observed that in some cases, not only did the

errors increase outside the Reynolds numbers mentioned above, but also, the system of

linear equations became so ill-conditioned that a singular behavior resulted. Clearly, no The solution could be reached under these circumstances. It can therefore be concluded that investig

the reliability of the numerical computations for Reynolds numbers less than 10 or satisfied

greater than 5 x 10 is uncertain. This is not a limitation of the theoretical solution, but _{the sens}

dissipati

n BC 1 (RE) BC I (Im) BC 2 (Re) BC 2 (Tm)

1

-2.15 x io

-2.17 x i0'3

1.19 x i0'

4.22 x iO

2

-1.53 x iO'

-1.35 x iO 1.30 x iO-' 1.12 x

3 -3.86 x -7.11 x 10_Is

-4.17 x 10b

-2.18 x

4 -3.22 x 10_16 -7.11 x I0

2.28 x 106

4.78 x iO

n BC 3 (Re) BC 3 (Im) BC 4 (Re) BV 4 (Im)

1

-9.71 x 10's

1.62 x iO-' -8.60 x i0

8.95 x iO'

2 -3.52 x 10_16.

-5.73 X 10'

-9.68 x 1016

3.97 x iO'

3

-4.00 x 106

-7.04 X iO'

-9.19 x iO'

_{-4.53 X i0} 4 6.85 x )Ø_17 4.11 x 1.77 x 10_16

3.94 x 10's

m BC 5 (Re) BC 5 (Tm) 1 2.89 x

-3.89 x iO'4

2

6.94 x i0'3

_{8.66 x iO4} 3 -1.21 x -7.00 x 4 9.91 x -3.35 x 5 -1.48 x 10 -1.36 x

Energy dissipation in stoshing waves - II 371 instead it is an artificial imposition due to the numerical algorithm used. The author was able to test the program for Reynolds number up to 4 x iO with significant errors in the

boundary conditions. The results presented in this study are well within the bounds of

error criterion set for the acceptability of the solution. 2.4. Roundoff and cancellation errors

Some important precautions that were taken to reduce or eliminate the troublesome roundoff and cancellation type of errors deserve few comments. These are inevitable

difficulties challenging today's numerical analysts and engineers.

It is not surprising that expressions which are completely equivalent from an analytical

point of view may turn out to be quite different numerically. The reason for such a discrepancy is due to improper handling of the numbers in digital computers. The

numerical analysis literature is full of the methods and suggestions about how to

circumvent this serious difficulty. We shall not repeat these but, instead, give a brief

description of the remedies exercised in this study.

First of all, we tacitly normalized every hyperbolic function appearing in the form of

our solution. Hyperbolic functions of real or complex arguments are normalized in

accordance with their arguments so that these functions can be kept bounded within the

limits of the capabilities of the computer system used. Furthermore, this process, in a sense, is a way to control the magnitudes of the values of hyperbolic functions so that they can be computed correctly. Without normalization, it is believed that a general

computer programming of the present problem would be unfeasible with the numerical limitations of the available computer systems.

Secondly, the mathematical formulas presented in Part I and Part II of this v'ork were

frequently rearranged before they were coded. This requires that every expression, either separate or combined with the others, must be simplified analytically to the simplest final form possible. These final and very simple forms are not only easy to

program but also free of the cancellation errors. Although the simplification or analytical reduction process could be very time-consuming, but its superior benefits far outweigh its

troubles. Therefore, most of the equations programmed were treated in this manner

before they were programmed. In particular, the author has witnessed a necessity of this

process in computing the viscous dissipation. We should point out that the viscous

dissipation consists of two terms, as shown in Appendix B. The first part of it is expressed

in terms of the vorticity, whereas the second part is related to the sum of the velocity

gradients. Although other forms for the dissipation function were also derived, computa-tionally this form seemed to be superior to the others. It can be shown that this form may

greatly be simplified after a lengthy but tedious rearranging. However, in spite of the

analytical manipulations, there still will be some terms in the remaining final form that can cause large cancellation errors. Under such circumstances, the terms of similar order of magnitude were grouped together in order to minimize the roundoff errors.

3. SUMMARY

The solution of a fourth-order linear partial differential equation has been

investigated. The governing equation and all of the boundary conditions are fully

satisfied by the assumed solution. It is clear that the solution obtained is approximate in the sense that the infinite series have to be truncated at some point. However, a detailed f f ) t r t

372 ZEKI DEMIRBILEK

study on the behaviour of the Fourier coefficients is performed in order to determine the best truncation cutoff limits. The dependence of the Fourier coefficients is investigated as

a function of the Reynolds and Froude numbers. An error analysis of the numerical

solution has been carried out. These parametric studies revealed information to

determine the limitations and range of applicability of the present solution.

The variation of the dissipated energy in terms of the Reynolds and Froude numbers

and the aspect ratio will be the subject of a subsequent paper.

REFERENCES

CHURCHILL, R.V. 1941. Fourier Series and Boundary Value Problems. McGraw-Hill, New York.

DEMIRBILEK, Z. 1982. A linear theory of viscous liquid sloshing. Ph.D. thesis, Texas A&M University, College Station, Texas.

DEMIRBILEK, Z. 1983. Energy dissipation in sloshing waves in a rectangular tank. - I. Mathematical

theory Ocean Engng 10, 347-358.

FORSYTHE, G. and Moler, C.B. 1967. Computer Solution of Linear Algebraic Systems. Prentice-Hall, Englcwood Cliffs.

LANCZOS, C. 1961. Linear Differential Operators. Van Nostrand, New York.

LANCZOS,C. 1964. Applied Analysis. Prentice-Hall, Englewood Cliffs.

LANCZOS,C. 1966. Discourse on Fourier Series. Hafner, NewYork.

LANDAU,L.D. andLIFSHITZ, E.M. 1959. Fluid Mechanics. Pergamon, Oxford.

PRANDTL, L. 1904. Uber tlussikeitsbewegung bei sehr kleiner reiburig. Proc. Third mt. Math. Congress, p. 484. RAYLEIGH, L. 1873. Some general theorems relating to vibration. Proc. Lond. math. Soc. 1, 357.

REID, R.O. 1981. Private Communications.

RILEY, J.R. 1969. Finite FourierBessel expansions. Master thesis, University of Houston, Texas.

APPENDIX A

The following is a list of the functions which appeared in Equation (1.38)(1.41):

T(nz)

= (Xm tanhXma - Pin tanhp.ma) (A.!)

_Ka(_1)m (A.2) Co(tn)

-x_m T(m) Xm( 1) 3(n,m) 2_{+ k2) (Xm2 + k2)} (A.3) (Pin I y(n,m)

-(Xm2 + k2)

(2

+ k2) (A.4) O(n,m) -

-

(i

Xm \ 1

-

k ) (p.2+ k2) (Xm2 + k2) (A.5) Yo(n) = 1 + _(A.6)

Y1(n)

= a(l)'

_(A.7)

Y2(n) = Y1(n) tanhkh (A.8)

Y3(n) = Y1(n) (A.9)

n

Y4(n) = Y3(n) tanhah (A.1O) I

I The p as the co containe th th th th vo It is liquid sb that is dissipati accordin dynamic We no in terms simplifyi

ie as al to rs ) )

-Energy dissipation in sloshing waves - II 373

Y5(n) = Yj(n)sechk,1h (A.11)

Y6(n) = Y1(n) Yo(n) sechnh (A.12)

K

Y7(n) = (A.13)

n

cu2(1 - 2vk2

Ys(n) = Yi(n) _gk2 sechkh (A.14)

iv2

Y9(n) 2Y1(n) i-- sechah (A.15)

Yjo(n) = Yj(n) sechkh (A.16)

Y11(n) = Yi(n) ---- sechah (A.17)

K1

-

gK _(A.18)

Y12(n)

= gk4

(-1)'k(k sechkh -

(_1m_IXm_tankh)

(\19)

Y13(in, iz)

= T(m)(K2 + k2)

Y4(m n)

_-

(_1)_lknXm(_1)m_l (A.20)

T(m)(X,2 + k2)

Y15(,n,n)

(-1Y'k(an

secha,,h - (_1)m_lXm tanhah)

T(in) (X,2 + a2) (A.21)

(1)1knXm( 1)

Y16(rn,n)

-T(m) (Km2 + an2) (A.22)

APPENDIX B

The principles of conservation of mass and momentum are supplemented by a third one known as the conservation of energy. This fundamental law relatcs the rate of change of the kinetic energy contained in a volume of fluid, -V, to the sum of the following four terms:

the rate at which the external forces do work on the element of fluid in -V the rate at which the stresses do work on the surface S boundingV-. the rate at which the energy flows out ofV'across S, and

the rate at which the energy is dissipated, irreversibly, by viscosity in each element of volume of the fluid.

It is the major objective of this study to investigate the fourth term mentioned in (4) for the liquid sloshing problem. Customarily, the viscous dissipation is denoted by 4), and it can be shown that 4)is, in fact, equivalent to Rayleigh's (1873) 2F function. It suffices to say that the viscous

dissipation function, 4), in tensor notation can be expressed as

(2

(B.1)

\äx1 äx,J

according to Landau and Lifshitz (1959) where in Equation (B.1) p. represents the coefficient of dynamic viscosity of the fluid.

We note that for a plane two-dimensional flow, one can express part of the dissipation function 4)

in terms of vorticity. Adding the continuity equation to the right-hand side of Equation (B.1) and simplifying yields

T'''

-

. -':--' 3 ge I al II,

374 ZEKIDEMIRBILEK

= *

4p.(uw1 +

w)

(B.2)

where u _{and w are the components of the fluid particle velocity in the x- and z-directions,}

respectively, and is the vorticity, defined as

We shall refer to the expression of4) given in Equation (B.2) henceforth. Likewise, during the

computational algorithm, 4) is determined through this equation, and the total dissipation is

obtained by simply integrating4)over the volume of the occupied fluid.

In the problem studied, the viscous dissipation can be assumed to depend on the following quantities:

the size of the tank, a, the fill depth, h, the center of the roll, d, the gravitational acceleration, g, the frequency of excitation, o,, and

the viscosity and density of the fluid, and p.

Moreover, the viscous dissipation is proportional to the square of the amplitude of roll excitation. This dependence can be expressed in a mathematical form as

4)=4)(a,h,d,g,w,p*,p,3)

It can be shown that in nondimensional form, this functional relationship reduces to

4)

f2

w2a d h

pu3a532 - ' g ' a ' a

or

ITj =f(ir2,

i,

ir4, irs)

where v is the coefficient of kinematic viscosity of the fluid. We note that the parameterir1

represents the dimensionless viscous dissipation, ir2 the well-known Reynolds number and ir3 the Froude number. The last two iiterms are the nondimensional center of roll and the aspect ratio, respectively. We shall assume that the tank is rolling about the center of the bottom of the tank henceforth. This implies that d = h or r4 = ir5, which eliminates the effect of d on 4), rendering

only one aspect ratio to work with. Consequently, in computing 4), we need only to consider the

remaining three most important parameters, i.e. Reynolds number, Froude number and the tank aspect ratio. Furthermore, since 4)is proportional to the viscosity in Equation (B.2), then it is

appropriate to define

(B.7) so that the effect of the viscosity is directly incorporated into the viscous dissipation.

= uz - wx (B.3) Oceani Printed A cc er ar in Fi di e w ju fo le bc THE D contain point oi of natui being concent

the lo

rectangi where I gravitati fundam inviscid Viscosit value. Labor slosh ing the freq depth, a!