Department of Civil Engineering Stanford University Stanford, California
COMPUTER STUDIES OF FINITE-AMpLITUDE WATER WAVES
by
Robert K.-C. Chan Robert L. Street
and
Theodor Strelkoff
Technical Report No. 104
This report was prepared under Contract NOnr 225(71), NR.062-320
Office of Naval Research and
Grant GK 2506
National Science Foundation
This document has been approved for public release and sale;. its distribution is unlimited.
ABSTRACT
Two numerical techniques are utilized to study the motion of two-dimensional, finite-amplitude water waves by using an electronic digital computer. The nonlinear properties of water waves are of primary
interest.
The first part of the work introduces the Stanford-University-Modified MAC (STJNMAC) code which is proposed as a valid tool for analyzing
incompressible, viscous flows with a free surface under transient con-ditions. The method is applied to the study of the solitary wave run-up
on a vertical wall. The results are compared with the available experi-mental data and give a much better prediction of the wave run-up than
the existing analytic theory.
in the second part, Newton's process of successive corrections is applied to solve steady-state potential floWS with free surface and gravity. A specific application to the analysis of solitary waves is made and all the wave characteristics are in excellent agreement with experiments.
ACKNOWLEDGMENT S
We wish to express special gratitude to Mr. Jacob E.
Frotmn, IBMResearch Division, San Jose, California, for his many valuable ideas on
numerical methods and computer programming.
Appreciation is also
ex-tended to Mr. Shen-Fu Liang, Department of
Chemical Engineering, Stanford
University, for useful information on numerical
method's.
The senior
author ii especially'grateful to his wife,
Jing-Kuan Chän, for her
encour-agenient
nd help in preparing the draft of this report.
Finally the authcrs thank the School of
Engineering and' the
Asso-ciate Dean for Research, Dr. W. R. Rambo,
for the grants.that supported
the major part of the computations.'
-TkBLE OF CONTENTS
Page
INTRODUCTION . . 1
WAVE RUN:-UP ON A VERTICAL WALL 5
2.1 Review of Related Investigations 5
2.2 The Marker-And-Cell Method 9
2.2.1 Governing Equations and Auxiliary Conditions . 12
2.2.2 Finite-Difference Representations 17
2.2.3 Computational Procedure 21
2.3 Modifications of the MAC Method 22
2.3.1 The Use of Irregular Stars Near the Free Surface 23
2.3.2 Extrapolation of the Velocity Fields 25
2.3.3 Calculation of Particle Velocities 25
2.3.4 Acquisition of the Undefined Variables 26
2.3.5 Calculations with the Modified MAC -- S1ThAC . 27
2.4 Stability Considerations 28
2.5 Results and Conclus±ons 31
2.6 Future Applications 36
THE SOLITARY WAVE AS A STEADY FLOW . 38
3.1 Review of Related Investigations 38
3.2 The Successive Correction Method 41
3.2.1 Formulation of the Problem . . 42
3.2.2 Transformation of the Problem in the Complex
Potential Plane 43
3.2.3 The Correction Equations 50
3.2.4 Solution Techniques 54
3.3 Results and Conclusions . . . 60
3.4 Future Applications 65
LIST OF TABLES
Table Page
2.1 Solitary wave run-up by the SUMMAC method . 69
2.2 Solitary wave run-up by Peregrine's method .
...
693.1 Statistics of a typical computation by the successive
correction methOd 70
3.2 Convergence of different guesses of the solution y(,) 71
LIST OF ILLUSTRATIONS
Figure - Page
2.1 Definition sketch 72
2.2 Wave run-up on a vertical wall C R/H ) 72
2.3 Wave run-up on a vertical wall ( R/d ) 73
2.4 Position of field variables 73
2.5 The Marker-And-Cell setup 74
2.6 Linear interpolation for
uk
2.7 Irregular star near a curved boundary .
...75
2.8 The undefined variables u2 artv...
752.9 The effect of u on . 76
0 K
2.10 InterpolatIon near a maximum 76
2.11 Second-order interpolation for Uk . 77
2.12 Interpolation of u.. 77
13
2.13 Flowchart for SUNMAC 78
2.14 u/J distribution . . . 79
2.15 v/Jj distribUtion . .. 80
2.16 dIstribution 81
2.17 Wave deformations C
Hid
0.2 ) 822.18 Wave.deforrnations
( H Id = 0.5 ) ...82
0
2.19 Comparison of methods 83
2.20 Comparison of u distribution under the wave crest . . . 84
LIST OF ILLUSTRATIONS (continued)
Figure Page
2.22 A rectilinear dock 85
2.23 Generation of periodic waves 85
3.1 Definition sketch 86
3.2 Comparison of wave profiles 87
3.3 Comparison of u/C under crest 88
3.4 The physical plane 88
3.5 The complex potential plane . . . . 89
3.6 The finite difference mesh 89
3.7 The governing equations of c(q,) 90
3.8 Flow chart for the correction method 91
3.9 Example of the computed flow pattern 92
3.10 The effect of corrections on H . . . 93
p
3.11 Comparison of the wave speed 93
3.12 Comparison of u/C distribution . . 94
3.13 Comparison of v./C distribution 95
3.14 Comparison cf u/C under crest 96
3.15 Extrapolation to find the maximum H0/d 97 3.16 Extrapolation to find the maximum Q . . . 97
SYMBOLS ANI) NOTATION
A1, A2, A3, A1 Areas Used in weighted averages
Aw A constant In Laitone's solution of solitary waves
B BernoUlli constant
-o
C Wave speed
c
Correction of y()
D Divergence of the fluid velocity D/Dt Substantive derivative
d Still water depth
f( ) Function of
g Absolute value of
Gravitational acceleration inx-direction
5
Gravitational acceleration in y-directionH Depth of still water In Peregrine's equation of motion in Chapter 2; approximate value for H in Chapter 3
H Undisturbed wave height
0
H
Pressure heidäfëesur'face
p
h HOrizontal distance 9f a particle from the nearest point where u or
v iSefined
i Subscript denoting - or x-direction
J JacObian
Subscript denoting
4-
or y-direction k Subscript identifying a marker particle L Wave length; a characteristic length2. Vertical distance of a particle f.ro the nearest point
SYMBOLS AND NOTATION (continued)
m, m1, m3, m1 Coefficients in the correction equations
m The x-component of a unit tangential vector at free surface
m The y-component of a unit tangential vector at free surface
Normal direction
n A superscript
The x-component of a unit normal vector at free surface
n The y-component of a unit normal vector at free surface
O( ) Order of magnitude
p Pressure
Q Flow rate
q Fluid speed
ciN Velocity component normal to a fixed solid boundary
Velocity component tangential to a fixed solid boundary
Vertical rise of water surface at a wall in Chapter 2.; also source term of Poisson's equation in Chapter 2
Residual
t Time
U Speed of uniform flow at infinity u Velocity component in x-direction
Uk Velocity component in x-direction for the kth particle
u False Uk
v Velocity component in y-direction n
Vk Velocity component in y-direction for the kth particle x Abscissa of cartesian coordinates
Location of the kth particle in x-direction y Ordinate of cartesian coordinates
YBI;X) Boussinesq's solitary wave profile
Location of the kth particle in y-direction y Height of the free surface in Laitone's formula
Z The constant 2 (_i + in Chapter 2; the ratio
(/5)2
in Chaper 35t Time increment
Sx Horizontal width of a cell c5y Vertical height of a cell
Spacing of grid points in 4-direction Spacing of grid points in iL-direction
s Ratio of wave amplitude to still water depth n Height of the free surface, measured from still
water level
V
p
4)
V2
SYMBOLS AND NOTATION (continued)
Irregular-star lengths
Kinematic viscosity of fluid Density of fluid
Ratio of still water depth to wave length
Ratio of pressure to fluid density ( p/p )in Chapter 2; velocity potential in Chapter 3
p/p at free surface Stream function Relaxation parameter Laplacian operator
1 INTRODUCTION
An understanding of wave actions is of vital importance in the design of many coastline and harbor structures. A number of failures in the past stemmed from inadequate design which was the consequence of man's limited knowledge about the great ocean waves. The destructive waves are usually large in amplitude and their properties change as they
approach.the shore. For example, the underwater topography may cause waves travelling nearly parallel to a shoreline to change direction and
focus their destructive forces on a seemingly safe harbor nearby. In the deep ocean an observer may notice Only a foot or two of surge as a tsunami or seismic sea wave passes. Upon reaching the shoreline, this same wave grows to monstrous breakers that bring death and serious damage to man's works along the coast.
A review of the literature of coastal engineering (see, e.g., the Proceedings of the Tenth Conference on Coastal Engineering, 1966) indicates that the action of water waves in the near-shore shoaling region or the surf zone is known only in terms of the results of experiments and model tests on specific configurations. In addition, the difficulties of making experimental measurements under transient conditions prevent us from acquiring more than a superficial knowledge of the total wave action. It would appear that the consequences of the true or nonlinear properties of propagating and deforining water waves are neither fully understood nor accounted for in design practice.
Because of their finite amplitudes, the mathematical description of the large water waves in shallow water is necessarily nonlinear
(Stoker, 1957). The importance of nonlinearity excludes the possibility of empipying the linearized theory which is based on the assumption of small amplitudes (Lamb, 1945). Most nonlinear analytic studies of water waves use perturbation expansions about a base flow such as a uniform
flow with hydrostatic pressure distribution. Stokes (1880), for example, studied waves of small but finite steepness by expanding the velocity potential about the still water level. The solution was obtained by
(1891) used a different method to attack the problem of steady-state solitary waves in a two-dimensional, infinitely long channel. A trial function for the velocity potential was selected to obtain an equation for the free surface. Then the parameters in the trial function were determined by minimizing the pressure at the free surface (Wallace, 1963). In recent years, a higher-order approximation theory was developed by, Laitone (1960) in a study of solitary and cnoidal waves; he obtained good agreement with experimental ate in several of the wave characteristics. Wallace (1963) modified McCowan's approaci and studied the behavior, of a solitary wave on various slopes. Dean (1965) presented an excellent theoretical model for steady-state nonlinear ocean waves. Monkmeyer and Kutzback (1965) gave a higher-order theory for deep water waves. Thus the nonlinear behavior of finite-amplitude water waves has received con-siderable emphasis in recent investigations.
Despite the fact that an increasing number of nonlinear theories are available, applications to practical problems are seriously restricted by the limiting assumptions inherent in these theories. The most common limitations include admissible boundary geometries, restrictionto Infinite domain, nonuniform satisfaction of the essential free-surface boundary conditions, etc. It appears that comprehensive theoretical tools are needed to provide adequate descriptions of the steady and unsteady motion of finite-amplitude water waves in the near-shore zone.
In vIew of the continuing and extremely rapid development of high-speed digital computers, we turned to the development of two essentially numerical methods as tools for wave analysis. The purpose of this work
is to demonstrate the validity of these numerical methods by applying them to simple,. yet practical, problems whose solution properties are well
known and experimentally documented. The employment of our numerical methods in more complex problems is discussed at the close of the
presen-tation.
The general basis of numerical methods is the representation 9f the governing differential equations of a mathematical model by finite differ-ence quotients. Using these difference representatIons, we can carry. out calculation of wave motions on a computer in what amountS to a series of
-2-numerical experiments. The essential difference between numerical and physical experiment can be summarized as follows. In a physical experi-merit the number and type of measurements that can be made are limitedby instrumentation capabilities, interference by measurement dei'icés in the flow and the cost of sophIst±cated instrumentation. For example,it is difficult, if not impossible, to get a complete velocity map of the flow field as a function of time [cf., Miller and Zeigler (1964) in which only the horizontal component was measured]. On the other hand, in a (rimierical) simulation one can, if he wishes, obtain as either printed or graphical output the fluid velocity, pressure and acceleration fields, the free-surf ace configuration, and plots of the individual water particle motions at every time step. Thus, it is nOw necessary to take an entirely different view of an experiment. Rather than designing an experiment with a view toward demonstrating one point or achieving one type of data record, we must now design our numerical experiments so that maximum significant use
can be made of all the data obtainable. We are led to simulation then of generally useful configurations rather than specific design configurations, at least In the beginning. However, the numerical- approach has its own difficulties, too. The finite-difference representation of partial
-differential equations is itself an approximation technique. What's more, numerical instabilities may lead to totally misleading resUlts and reduce
the value of computation to nil. .-Therefore, the validity of -a numerical scheme must be verified by comparing the results with theoretical results or experimental data.
The present work deals with the motion of a solitary wave. - A
solitary wave is a wave consisting of a single elevation of fluid, of height not necessarily small compared with the total depth of the fluid. If properly started, this wave may travel for a considerable distance along a uniform canal, with little or no change of form (Lamb, 1945)-.. Following Jordaan's (1965) argument, Street and Canifield (1967) suggested that, to study the. Shoaling, breaking and run-up character-ist4cs of: the dispersive, catastrophic ocean waves, it is reasonable. to begin in the laboratory by generating and studying nondispers-ive waves such as the
results and experimental data on this wave, in addition to its practical. importance, we have also selected the solitary wave as the medium for the illustration an4 testing of our numericalmethods presented in this york.
The present study is divided into .two parts. The first part is concerued with the transient motion and deformation of a wave as it pro-pagates toward a vertical wall in a channel, of constant depth. in design-ing breakwaters and similar protective structures, the engineer must know
the maximum vertical rie of water, called the run-up, as his design wave reaches the wall to avoid overtopping and consequent disasters. Unf
or-tunately, as Street and Camfleld (1966) pointed out, the run-up criterion. widely used. in engineering practice tends to seriously underestimate the run-up of large ocean waves. Accordingly, the objectives of this part were two-fold. First, we wished to show that the numerical, technique
provides satisfactory predictions of the run-up on a vertical, wall. Secondly, we wished to implement and verify the Marker-And-Cell (MAC.)
numerical technique (Welch, et al., .1.966) as a valid tool for theanalysis of two-dimensional,, nonlinear wave motions.
in the second part of this work, a useful numerical technique, known as the Newton's method, is introduced with an application to the problem of a' steady-state solitary wave. The main concern is with large-amplitude solitary waves propagating in a channel of constant .depth and infinite length., such waves being useful as input data to the MAC cal-culations, for example. To date, no existing theory. is known to accur-ately predict the properties of a solitary wave with a very large ampli-tude. However, we shall show that the present approach can produce strikingly good results, provided some techniques are available to solve the mathematical model accurately enough. The technique presented here employs a method of successive corrections to transform a nonlinear problem into a linear problem. The use of this method to solve the steady-state .so'iitary wave problem was first suggested by Professor Theodor Strelkoff of the University of California, Davis, during his visit at Stanford University in the Spring of 1968.
The computations were performed on an IBM 360/67 system. Graphic outputs of the wave profiles were obtained through the use of a CALC0
2. WAVE RUN-UP ON A VERTICAL WALL
2.1 REVIEW OF RELAT INVESTIGATIONS
The investit-jons of Wallace (1963), Street and Camfield (1966), and Peregrine (1967) are directly related to the present study and are discussed below. The first two investigations were ultimately concerned with the primary design quantity R/H0 --the ratio of wave run-up to the undisturbed wave height. However, in reviewing their work we shall'
com-ment on other features as well because of their relevance to the present results.
Extending McCowan's (1891) treatment of the steady-state solitary wave, Wallace (1963) proposed an 'approach to the problem of describing
the continuous deformation of a solitary wave in a channel of cQnstant depth. This extension was made possible by constructing a fixed region--an infinite strip--whIch contains the fluid under the wave as a subregion and the mathematical sources which produce the flow field. The fluId motIon was assumed to be invicid and irrotational. Therefore., the
veloc-ity potential in the fluId is a harmonic function that satisfies the Laplace equation V2 = 0 . Wallace let the free surface be defined by
y = (x,t) + d (Fig. 2.1). Then at the free surface, 4 must satisfy
the Bernoulli equation
-+
(er]
± gy = constant (2.1)under the dynamic condition p/p = 0 , and the kinematic condition
D
Dt - n(x,t)] = 0
Here d is the depth of the undisturbed water far from the wave crest,
p and p are the pressure and the fluid density, and g is the accel-eration due to gravity.
V
The method. of Solution required fIrst the selection of a trial function for the velocity, potential . This trial functIofl contained
several parameters to be determined at later steps. Under the assumptions
of the problem, the complex potential field W = + i is analytic, where the stream function i4 is the conjugate harmontc function of
As a consequence of this analyticity, the differential equation (2.2) for the kinematic condition at the free surface can be written in a form which permits. an exact solution for the shape of the free surface in terms of the t-rial function for . Finally, the Bernoulli.equationwas
used to determine the parameters in the trial functipn for so that the presure variation at the free surface Is minimized.
For a solitary wave in a channel of infinite length, the trial function used by Wallace is the field of a single dipole moving in the direction of the wave. To describe the wav deformation near a vertical wall, two dipoles of equal strength an speed, but moving toward each
other, were used. The values of R/H0 computed y this method are shown in Figs. 2.2 and 2.3.
In the Wallace method the condition p/p 0 is not satisfied uniformly along the free surface. In fact, p/p = 0 only at infinity. At the crest of the wave, the pressure was made to vanish only to the
order ri3 . Consequently the approximation deviates farthe and farther
from the true solution as the amplitude of a wave increases.
Street and Camfield (1966) ran a series of experiments in a two-dimensional wave channel to study solitary wave deformation on plane slopes ranging from horizontal to vertca1. The test objectives included verification of deformation theories, delineation of the shoaling pro-cesses associated with various beach slopes, ad examination of the limit-height-wave concept. These experiments were run in the Stanford Wind, Water-Wave Research Facility (Hsu,1965). 7aves were generated by a vertical piston-type wave plate. The motion and position of the wave
plate were controlled through a iiydraulic-servo-electrbnic system by an externally applied voltage. It was assumed that the initial generated wave can he represented by
I
13H\½
= H sech2
L
<_.2.) (2..3)
applied to the wave generating system to prodUe such a wave. Wave heights and shapes were recorded or. strip charts by use of capacitance wave-height gages and Sanborn recording osciliograph systems. Defàrma-tion details near the slope were recorded on 16 mm movie films.
Of particular pertinence to the present analysis were the run-up data on a vertical wall. These data, plotted in Figs. 2.2 and 2.3, indicate that Wällacé's theory tends to underestimate the run-up for higher waves. It was found ecperimentally that for very low values of H0/d the value of R/H on a vertical wall is approximately equal to
2.0 . However, as H/d increases R/H0 also increases. Van Porn
(1966), making use of these data, suggested the equation
= 2.0 +
(2.4)
to describe the run-up on a vertical wall. This simple, empirical
equation fits the experimental data far better than the theoretical pre-dictions of Wallace.
Peregrine (1967) derived a set of equations of motion for long waves in water of varying depth. The esseice of his method is the ex pansion of variables such as the Surface shape ri , the pressure p
the horizontal velocity u , etc.-,- -in terms of the powers of a parameter --the ratio of amplitude to depth--under the assumption that
c =a
where a is the rat-io of water depth to wave length. Substituting these expansion expressions into Euler's equations of motion and the continuity equation in an integrated form yields equations (for long waves on a beach) that are accurate to the second order in
I-n two dimensions these equations are
T) 1 2 3u
/H a2u
. 1 2H u + u + -= H 232+
H 23x 23 t+
2 [(H + )u] = 0-1
(2.5)
The surface position is described by n ,' u is the average horizontal
velocity, and H = H(x) is the depth of still water. All nonlinear terms are included in Eqs. (2.5), consistent wlth.the usual requirement that E << 1 and second-order accuracy. Unlike the linearized equations
of motion (Lamb, 1945) or the finite-amplitude shallow water equations (Stoker, 1957) which assume that the surface slope of the waves Is very small, Peregrine's method provides a better approximation for steeper waves and Is computationally stable for the unsteady propagation of
finite-amplitude waves. Although the wave grows steeper when it runs up a vertical wall, this theory can be used to calculate R/H0 for a limited range of relatively small, but finite-amplitude, waves. At Stanford
University we have developed a computer program to solve Eqs. (2.5) numerically by the method outlined by Peregrine (1967). Some of the
results are presented in Figs. 2.2 and 2.3.
The computational methods of the present study require prescription of an undisturbed initial wave located far enough from a vertical wall so that subsequent deformations can be analyzed by a 'computer program. In terms of' expense of computation the most economic way to obtain such a starting wave is through some necessarily approximate, but satisfactory, analytic solution for the solitary wave in an infinite].y long channel of constant depth. Laitone (1960) used a perturbation method to obtain solutions to models of cnoidal and solitary waves approximately. He extended the first-order shallow water approximation to obtain the second and third approximations to shallow water theory in order to get' the solution to the second approximation to the cnoidal wave and the solitary wave. These results can be summarized as follows:
'For the, surface profile, he obtained
= 1
+ (_2-) sech2'(Awx - Ct)
-sech2(Aüx - Ct) [i - sech2(Awx - Ct)]
,The velocity and pressure fields were given by
11
/H\I
/H U IOil
ii
0-()2[
9()2]
IH\
/H \ 2IH\r
__i1_aj(\ +(2.1I_2
8\d/
2\d/\d/
\d/L
)(]sech2(Awx
- Ct) sech(Awx - Ct) IH \3/2- ..JT-#) (i-) sech (Ax - Ct) tanh(Awx - Ct)
(2.7) sech2(Awx -(2.8) and -
(H)2
1)+(_
)2] [2 sech2(Awx - Ct)- 3 sech(x - Ct)
(2.9)Here Awx .(x/d)J. (Hid) [i - .. (H/d) , g
Il
is theacceler-ation due to gravity in the y-direction., y9 is the free surf ace éleva-tion and Ct is a product of time and wave celerity that sets the initial location of the wave crest.. These equations are used herein because they agree more favorably with laboratqrymeasurements of several key wave
characteristics' than other theoretical results (Wiegel, 1964).
2.2 THE MARKER-AND-CELL METHOD
The numerical scheme employed in this part of the work is the Marker-And-Cell (MAC) method that was originally developed by Welch, et al., (1966) at the Los Alamos Scientific Laboratory, University of California. The Navier-Stokes equations are used to compute unsteady, incompressible viscous fluid flows with or without a free surface. In its presently implemented form the MAC method is suitable for analyzing two-dimensional flows. Welch, et a].., (1966) gave a detailed description of their original method, flow charts for computations and several ex-amples of numerical solution. Our discussion is concerned with a brief review of the method and Its implementation, after which we deal with
the specific application to low viscosity fluid motion in water waves and its attendent problems.
One feature that distinguishes the MAC method from others is that its finite-difference equations rigorously conserve momentim. Also, the calculationál procedure has been arranged in such a way that it is self-correcting. The errors made at any time step, regardless of their source, tend to be corrected in the next step. Furthermore, the use of hypothet-ical"marker" particles to mark the free-surface location makes it
possible to carry out the computations in a well-defined domain.
The solution technique and unique fluid representation method are described for a two-dimensional flow by Welch, et al., (1966) as follows:
"There are, in effect! two coordinate systems used in MAC-method calculations: The primary one covers the entire domain of interest with a rectangular grid of cells, each of dimensiois Sx by 6y . The cells are numbered by
indices i and j , with i counting the rows in the
x-direction and
j
counting the rows in the y-direction. The field-variable values describing the flow field are directly associated with these cells."Actually, the true fluid would, in general, have a different set of field-variable values for every infini-tesimal point in the fluid. The representation used for computing, however, must be restricted to a finite number of values, each apprOximating an average through the
inmiediately adjacent region. It follows that the accuracy of the representationdepends strongly upon the fineness of the mesh compared to the macroscopic structure of-the flow.
"In addition to the. primary-coordinate system of finite-difference cells, there is a coordinate system of particles whose motions describe the trajectories of fluid elements. These particles serve two purposes: First, they show which cells are surface cells, into which the surface boundary conditions should be applied. Second, they show the. motion of the fluid and all its distortions as it passes through the computing region. For the first purpose, it would be necessary to have particles only near the fluid surface. The second purpose, however, is also of. considerable value, particularly if facilities are available for easily
plot-ting particle positions at various times through the prog-ress of the calcülàtiôn. We have found that complete con-figuration plots carry in concise form much of the inf or-mation require4 for the analysis of results. .
"It should be emphasized that these particles serve only as asslèss markers of thecenters of mas of the
and enter into the calculations only insofar as they' designate which are the surface cells.
"The cell-and-particle system enables an instantaneous representation of the fluid for any particular time during the evolution of the dynamics. In addition, it i necessary. to have a means of actually calculating the changes with time of the fluid representation What is needed is a com-puting technique whereby the prescribed initial conditions can develop, within the limitations imposed by the boundary conditions, into that subsequent set of configurations that most nearly represent the behavior of a true fluid.
"As in most other fluid dynamics computing methods for transient problems, the MAC teclnique works with a time cycle, or 'movie, frame,' point of view. This means that the calculation proceeds through a sequence of cycles, each advancing the entire fluid configuration through a small, but finite, increment of time, 5t . The results of each
cycle, act as initial conditions for the next one, and .the calculation continues for as many cycles as the investigator wishes. Each cycle is itself subdivided into phases:
The pressure for each cell is obtained by solving a finite-difference Poisson's equation, whose source term is a function of the velocities. This equation was derived subject to the requirement that the resulting momentum equations Should pro-duce a new velocity field that satiSfies the In-compressibility condition.
The full finite-difference Navier-Stokes equations are used to find the new velocities throughout the mesh.
The marker particles are moved to their new posi-tions, using for their velocities simple inter-polated values from the nearby cells.
Bookkeeping processes are accomplished related to the creation or destruction of surface cells,, the input or output of particles, the advancement of a time counter, printing or plotting results, and. numerous similar matters.
"By the end of the cycle, the results have been arranged in the computer memory in sich .a way that the next cycle can inediately begin.
'There are. two types of off-line prints, one presenting information pertaining to each marker particle, and the other containing cell data, such as velocties and pressures.
The prints are made only occasionafly, according to .a pre-. determined print interval incorporated in the program or
upon demandby the operator through, the usepf sense switches." The original version of the MAC method was used in the early stages of the present investigation... The results were no,t satisfactory. The computed free-surface profiles were very 'ragged, iictiiig considerable
errors originating from the treatment of free-surface boundary conditions. Several modifióations have been made Sitice then and.the result Is a quite satisfying procedure. As these changes have a distinct effect on the accuracy of solutions, the changes are described in a separate section.
2.2.1 Governing Equations and Auxiliary Conditions
The well-known Navier-Stokes equations--the equations of motio which form the basis of the science of mOdern fluid mechanics--are used
in the MAC nethod. When the continuity equation is included and the suitable boundary and initial conditions are specified., the Navier-Stokes equations are sufficient to determine fluid flow patterns uniquely.
Historically, the Navier-Stokes equations were first. derived by M. Navier in 1827 and by S. D. Poisson in 1831 onthe basis of an argument
involving intermolecular forces. The same equations were derived later without using such 'a hypothesis by B. de Saint Venarit in 1843 and G. G.
Stokes in 1845. Their derivations were based on the assumption that the normal and shearing stresses are linear functions of the rate of def or-mation of a fluid element--an assumption consistent with Newton's law of fluid friction--and that, for an incompressible fluid, the thermodynamic pressure is equal to the negative of one-third of the sum of the normal stresses.
It is evident that the assumption of linearity is actually merely an extension of the. stress-strain relations, used in solid-body mechanics. Therefore the Navier-Stokes equations do not necessarily give a true description of the motion of a fluid. They must be verified by experi-ment, whtèh remaIns the only true test of the vaiidty of a theoretical model. 'Indeed, a greatnumbet of problem solutions, such as the exact
solution of laminar flow through a circular pipe and the approximate solutions of boundary-layer flows, agree. so well with experiment that the general validity of the Navier-StokeS equations can hardly be doubted
(Schlichting, 1968).
For an incompressible flow the Navier-StokeS equations in two dimensions are
-Dv Dv Dv
- + U - + V - =
Dt Dx Dy U V L' sjT ' gL
+
fD2v
D2V\ Dy yand the conservation of mass equation for incompressible fluids is
+ !Y.
=o
Dx
y
Here, x and y are the coordinates of a. fixed Cartesian system;. u and v are the velocity components in the x- and y-directions, respec-tively; is the ratio of pressure p to the constant density p of the fluid (i.e.,
E p/p );
and are the x and y components of the body acceleration (gravity); and v is the kinematic viscosity of the fluid.In handling a family of dynamically similar problems, it is con-verlient to nondimensionalize the variables in Eqs. (2.1Q)-(2.12). By
introducing dimensionless parameters Eqs. (2.10)-(2.1.2) can be written in the form y (2.11) (2.12) ._1 t
-0
(213)'g
'/L7j'LjL
where g g and L is a characteristic length. In our study of' the sol±tary wave problem, for example, the water depth in the undisturbed region can be used as the characteristic length., that is, L. = d . Notice
that the form of Eqs. (2.10)-(2.12) is not changed by the nondimensional-ization of the variables.
Using Eq. (2.12), we can write the following relations:
Du2 Duv Du Du I Du
-
Dx+ = U 'fl + V + U I
-Dy Dx Dy \DX Du= u - + v -
Dx Dy(32u
D2u (2.10)Biw By2 By
(u
av
Bx By Bx By
\BX
ByBy
=
u+V--By
Thus Eqs. '(2.10) and (2.11) can be written in the alternat±ve form
+ = By
---+g +v
Bx '+i+_ .it+g
+v
Bt Bx By' By y fB2u B2u + (B2v B2v (2.14) (2.15)if the field variables u , v and are defined at the locations
shown in Fig. 2.4, then it can be shown, by integrating aver a fixed con-trol volume, that the finite-difference representations of Eqs. (2.14) and (2.15) conserve momentum rigorously
while
the finite-difference repre-sentations of Eqs. (2.10) and (2.11)'do not. In addition, the illustratedarrangement of variables is particularly useful for setting up a test of the degree to which the finite-difference form of the continuity equation is satisf1e4.
In the MAC method, velocity components and pressure are used as the dependent variables. Equations (2.14) and (2.15) provide a means of computing Bu/Bt and Bv/Bt so that the data of u and v can be up-dated. 'But this updating procedure requires knowledge of
3/Bx
andThus, further derivations are necessary to obtain a governing equation for the pressure field. First, we define
D =
3x By
and the continuity equation becomes
D = 0 - (2.16)
a2/x3t
2/tx and a2/yt
E a2/t3y and then adding the resulting equations, we have + - R (2.17) where -22
v2 Dfa2D
2D' R = + 2 + + - v +Equations (2.14)-(2.17) const-itute the basic equations from which a finitedifference scheme can be readily developed. It is necessary to retain the terms containing D in Eq. (2.17) because the combination of truncation Errors and other inaccuracies lead to a nonvanishing as the flow field evolves. It is desirable then to devise a scheme that aims at making D vanish at every time step. This is very important in maintaining the accuracy and stability of numerical solutions.
To complete the formulation of the problem, proper initial and boundary conditions must be imposed. The initial distributions of u v and , plus the free-surface. position, are required to start the calculations. One way to obtain the initiàl condition s to start with a sJthpe floc fiEld. For example, one thight begin with a uniform flog or even with a fluid at rest. Then moving obstacles or other types of dis-- turbance may be applied to the flO so that thE desired flow patterns
wil-1 eventually de4ielop. Another method is to use some analytic solution as the starting point.
The boundary conditions can be grouped into four categories for most problem of intErest.. These are () a solid wall that is impermeable,
(2) an inflow boundary where the velocities are prescribed, (3) an out-flow boundary through which the flUjd leaves the computation region, and
(4) the free surface. Only (1) and. (4) are employed in our work and discussed here.
Boundary conditions at an impermeable solid wall are the easiest to apply. Let and be the velocity components, normal and tan-gential respectively, to a solid wall. The bouidEry conditions are
(1) = 0 and 0 if no slip is allowed between the 'lu4 and the wall, or
(2) 0 and
q/N = 0
if the expected boundary layer is sothin that the f]uid is practically slipping along the boundary. The conditions pertaining to the free surface present a special problem. The free surface as defined here is the interface between a liquid and a massless gas. Generally speaking, the stress tangential to the 'surface must vanish and the stress normal to the surface must exactly balance any externally applied normal stress. Following Welch, et al., '(1966S we let n and fly be the components of a unit outward vector
normal to the surface and the surface, elevation be described by the function n(x,t) ; then,
Also, if m and thy are the correspondjng components of a unit
tangen-tial vector,
m =n.
andm =n
x y y
In addition, let at) be an externally applied pressure. Then, with-in the flUid at the free Surface, the tangential-stress.condition is
1 . u / \13u 3v\ 3v I
'i2nm +lnm +nm.ii+---i'+2nm i0
L X
x x \ x yy xjy
x/ y y 3yjand the normal stress condition can be expressed as
+2vin2+nn(.+)+fl2
a
L x ax x y ay 3x 'y ay
For incompressible fluids with low viscos±ty, such as water, it is sufficiently accurate to use only the ing.le condition .
Undr
[i + (an/x)2]_
2.2.2 Finite-Difference Representations
In order to solve the above partial differential equations by a finite-difference scheme, the computation region is divided into a number of rectangular cells (Fig,. 2.5). Each cell is named or "flagged" as being of a. certain type, depending on whether the cell represents a solid wall, a free surface, empty space, a cell full of fluid, and/or a boundary
or corner cell. By using these flags, the pioper boundary conditions can be applied to the pertinent cells.
The positioning of the field variables u , v and . in the MAC
method is quite different from other numerical schemes. Instead of de-fining all these variables at the same location, namely, at the nodes of the rectangular mesh, U is defined at the center of the right-hand and left-hand sides of the cell, while v is defined at the center of: the upper and lower sides. Then is defined only at the cell centers
(Fig. 2.4). As mentioned in the preceding section, this Special arrange-ment is needed to rigorously conserve moarrange-mentum. and to provide an easy way
to test the continuity equation D = 0
Now suppose the term
u/t
is represented by the first-order finite=difference approximationn n+1 n
kit)
(2. 18)where refers to. the nth time step and t is the time increment. This approximation is obviously crude, but not uncommon in the numerical treatment of many parabolic or first-order initial-value problems. One way to efne the approxImation is to iterate between the nth and the n+lth steps before proceeding to.the step n+2. However, this would demand a considerable amount Of computer time and storage. The prograrnthing work would also be quite formidable. In the MAC method the simple formula
Eq. (2.18) 'is used. The permissible size of' 6t is discussed' in Sec. 2.4, in connection with the problem of numetical stability. "
usual circumstances
a = 0 but le can also specify = f(x,t) along the free surface to generate some types of waves.
Terms that involve only the space derivatives can be approximated
more accurately. For example,
-I
n+1 I iiVjj+½hjjf½+
t
+ g + UU.
3fu\
26x i+½j U+ Uj½j -
2u1 :i.+-j-
(6x)
Applying similar procedures to each
term
in Eqs. (2.14)-(2.17), we obtain the finite-difference representationsn+1
Uj - U+1j
(uv)j½ -
(uv)jj
ui+½jUj.j
+t
5y + +ii xi+1j
+ V Pij -i1±1 c5y(U
-
18-(V
Vi+1+ + Vj_jf
-+ cSx +Uj½j -
2ti1.
5xz(UV)j_½j
-
(uv)jj
'Sx 3 + Vjj - 2vjj.+½ yTL1
(ii1
+ i-1 ij+1 + ij-113
Z\
x 6y 13(2.19)
(2.20)(2.22)
(2. 23)
1+1+
- 5y + Ui.f½j+l + 5yL-
2u)I1
(2.21)
Here
and
- /1. z
=
In the above equat-ions, variables with the superscript n+1 are related to the n±lth time step. Variables lacking a superscript are evaluated at the nth step. One of the essential features of the MAC method is the presence of the D terms in Eq. (2.24). Because of the inherent error.s in finite differences and roundoff by the computer, these D terms do not vanIsh in general. Thus, an effort aimed at making
= 0 at the n+lth step must be incorporated into the calculations at the nth time step. Indeed, Eq. (2.23) is obtained by writing Eq. (2.17)
in finite differences and then setting D1 = 0
Only the boundary conditions relevant to the solitary wave problem are presented hEre. For other types of problems special boundary conditions may arise, but. the treatment is analogous (Welch, et al., 1966). In Fig.
2.4, cell (i,j) is next to a solid wall and cell (i-1,j) is a boundary cell. For the "free-slip" boundary we have
a.
Vjlj+½
b u2 = 2
j_ij
)
u
.-u.1.
V-v
- i+½3 ij+½
ij--ij -
c5y V. .=V
i-lj-%
(vi+)2
+ (v)2 -
2(v..)
½
(2.. 25)+
xy
[uv)jj
± (uv) -(uv).
- (uv)i1j..]D.1
± D1 - 2D
+ - 2D6t SXZ
(uv)j_._½ = 0 2v /
-(U.
X \ 1-t-3
-
Uj½j)
=4
+ Ui_
=4 (v..
+ Vjj...11
(uv). = uiI.½iI
(2.26) C. (uv)i_½j = 0 , (uv)j__½ = 0 d.c_11 =
- g.Sx
e..Ui½1
= 0f. D1 = D1
If the thickness of the boundary layer is greater than x , then the "no-slip" boundary conditions must be used.:
v11
=. -u2i-ij
=u.
13 C.(uv).
0 ; = g.5x -= 0D1. = D11
Some variables in Eqs. (2.21)-(2.24) are not directly defined in
the i.i , v fields. They must be obtained by an interpolation formula.
Simple averages were used in the MAC method. For example,
)
As we shall see later, this linear interpolation is not satIsfactory. It may introduce significant truncation errors in a region where the
u/3x and v/y terms are large (see Sec. 2.4).
NExt consider the boundary condition at a free surface. In the original MAC method
= a the prescribed surface pressure, is assigned to the surface cells. Since thE field is defined at thE center of cells, the condition
= a' is applied at the cEnter of surface cells
rather than at the exact location of the free surface. Serious errors may result from this arrangement. It may even creatE the impression that
the MAC method is highly unstable. We shall return to this point in Sec. 2.3.1, where modifications are suggested.
2.2.3 Computational Procedure
The initial condition provides the starting values of ü and v Since the initial position of the free surface is known, we can obtain the pressure field by using Eq. (2.23). But, the value of at a given cell cannot be calculated directly from this equation. because the
values at the neighboring cells arE also involved. The field of
values must be obtained simultaneously. A direct matrix inversion method is rather time-consuming, and the direct elimination method is difficult to use because the geometry of the fluid domain is complicated. In this respect the. Gauss-Seidel iterative (relaxation) method (McCracken and Dorn, 1964) proves to be very useful, The field can be roughly guessed at first--a hydrostatic distribution would be appropriate, for example. Then, proceeding from one cell
to
another in a systematic manner, Eq. (2.23) is used to replace the ol.d values of . The iterationistermi-nated when the new values do not. differ much from the old values. Because either or
/N ,
the normal derivative, is specified along the boundarIes, convergence of this iterative scheme is guaranteed(McCracken and Dorn, 1964).
Solving the pressure field consumes a great portion of the com-puting time. While Welch,, et al., (1966) make no recommendations on detailed computational procedures, it seems desirable to achieve some efficiency on the iterations. Various techniques are available to
n+1 ±1
Xk =Xk+Uk..cSt
n-I-i .n .
-=
+ Vkt
where Uk and Vk denote the velocity components fOr the kth particle.. The values Uk and Vk are obtained by interpolating from the u and v fields, respectively. For example, if particle k 'is situated within
the zone bounded by the points where u1 ., u2 , u3 andu are defined
(Figo 2.6),. uk can be computed by
A1u1 + A2u2 + A3u3 + Aku
Uk =
SXcS)7(2.29)
A similar formula is used to calculate
vk .
Equation (2.29) represents a bivariate linear interpolation which is inadequate under certaincircumstances.
NOW we have, completed a cycle of computation, i.e., we have update4 the u ,.. v and data to a 'new state. The, foregoing procedure can be repeated for as many time steps as necessary for the investigation.
2.3 MODIFICATIONS OF THE MAC METHOD
Several diffict4ties arose when we made. the first attempt to cal-culate the 'motion of: large-amplitude water waves b the MAC method. The well-shaped initial wave became irregular after a few time steps of
com-putation. Also, the. u and v 'fields. did not vary as' smoothly as they should. Because a small viscosity was used this experience seemed to
22
-= old
+ w(d 'told)
(2.27)'where the relaxation parameter w varies from 1.5 to 1.8 f Or best results. After solving the pressure'fie.ld, we can use Eqs. (2.21) and (2.22) to calculate the new distributions of u and v . At. the same time. a
new free surface geometry is produced by moving the marker partiëles to thejr new positions. The following formlas are
used
for this purpose:I
confirm the statement that the MAC method is numerically unstable. for computing low viscosjty fluid flows (Hirt, 1968).
Usually, a formal numerical StabiIty analysis is performed on the governing equations of the problem, suh as Eqs. (2.21) and (2.22). How-ever, improperly applied boundary conditions may be a more serious source of errors than the governing equatioiis.. We have found by experience that the original MAC method can be modified to yield very sat1sfactory results. Most of the changes have to do with the free-surface condition This
particular condition is inherently difficult to apply in the rectangular mesh system due to the complexity of the surface geometry.
The modifications presented here. were developed by heuristic rea.soning since a rigorous numerical analysis is diff:icult., if not im-possible. Nevertheless, the modified scheme has been teste4 and has produced good results.
2.3.1 The Use of Irregular Stars Near the Free Siurf ace It is not proper to impose the condition
a at the center of surface cells. Only at the free surface, which does not necessarily co-incIde with surface cell centers, is = . This difference affects.
the accuracy of the solution of the field wkich in turn influences the calculatiOns of the U and. v fields. Unfortunately, features near a curved boundary cannot be adequately resolved by a finIte-dIfference scheme unless "irregular stars" (Fig. 2.7) are used. Therefore, we derived a finite-difference form of Eq. (2. 17) and used it to apply
at the exact location of the free surface by employment of irregular stars. Let r1 , ii3 and be the lengths Of thE' four legs (Fig. 2.7) and
,
be the values of at the ends of the four legs. Using the Taylor series expansion (Allen, 1954), we can write
1
(P)
±
(..)
+ o
(nfl
iJ ii
(2.30)
Equation (2.17) now becomes
(}4)+
:Lj Finally we get - -2 3Lf I +13
L2
+ ii 2(nn1n3) ((-'- ±n)
+ + Ri.] (2.37)by substituting Eqs. (2.34) and (2.35) into Eq. (2.36).
It should be noted that Eq. '(2.37) reduces to Eq. (2.23) when
applied to aa interior cell, where =
= i+1 ' 3 2 = ij+1
and
= . Thus, Eq. (2.37) isa more general formula for solving the f±eld.
= 'ij + 2 1
Y)jj
= 2L(2 +
flLf) L / 'z)=-R.
y J 24 -2 + (2.32) + (2.35) (2.36) ij n±() - 0(n)
(2.33)We neglect terms of the order n , , etc., or higher, Then, by
eliminating the (3/ax) term between Eqs-. (2.30) arid (2.31)! we Obtain
n1n3(n1+ 3) [n
+ n1- (n1
:1) (2.34)2.3.2 Extrapolation of the Velocity Fields
The quantities u and v are not defined outside çhe fluid region, but they are needed to carry out the computations using Eqs. (2.21)-(2.24) near the free surface and to move the marker particles. For the case
shown in Fig. 2.8, u2 is set equal to. u1 and v1 is set equal to
v2
in the MAC method. This arrangement makes = 0 , but this choice of u2 and v1 is entirely arbitrary because D1 0 n general in a sur-f ace cell. Also, consider the case shown in Fig. 2.9. There, u is set equal to u1 in the unmodified MAC method. The result is an erroneous velocity profile in the extrapolated region. Since the particle is to be moved according to the Interpolated velocity from the u and v fields, the particle k will assume the erroneous velocity u instead of the more accurate value u.These problems can be resolved in a fairly simple way. We have only to calculate the undefined u and v values by a proper extrapo-lation formula. For this purpose the Gregory-Newton backward interpola-tion formula can be used (Jeppson, 1966) to give
fu\
1 /11 3 1= -s;,
Y
-3u2 + -
U1 +
11 2
-i- u1 - 5u2 + 3u3 u1
Similar considerations are applied to the extrapolation of v and to other Orientations of the free surface. In a given case we extrapolate the u and v fields far enough into the empty space so that the subsequent Eqs.
(2.40)-(2.42) can be applied near the surface.
2.3.3 Calculation of Particle Velocities .
The formula used in MAC to obtain particle velocities from the u and v fields is a bivariate linear interpolation forthUla [see Eq. (2.29)]. Generally speaking, this is adequate
it
a region where u and v change slowly and inonotonically. When the particle is situated near the maxima(2.39) Using (Du/y)1 (u0_u2)/2y in Eq. (2.38) and rearranging, we have
)
or minimg of the velocity fields (Fig. 2.10), however, the linearly inter-polated partiole velocity can be quite erroneous.
To derive a more satisfactory interpolation formula,. consider the particle k in Fig. 2.11. If the particle lies in the shaded area, we
can find the velocity component uk for the p.artjc by making 4 Taylo series expansion about the point 0 . Neglecting the third and higher
order terms, we have
= u +
h()
+ + r[h2()
+2h2() +
± I h\Iu1 - u3\ (ui.1-j
(u3 +
(2.41) u 0+
L\6x/\II 2 /J + I--H--
\y/\
2 + (u1 + u3 - 2u +()2
(u2 + u - 2u)+U7 - U6 -
us)] (2.40)In a similar manner, a formula can be derived to calculate vk, the y-component' o the particle velocity.
2.3,4 AcquIsition of the Undefined Variables
It was noted in Sec. 2.2.2 that undefifled variables are calculated by simple averages as given by Eqs. (2.26). Again, this is a linear
in-terpolation that introduces errors in a region where u or v changes rapidly.
Suppose we are to calculate from the known values at the neighboring points (Fig. 2.12). Because of the symetry of the data
points about the interpolated u1 , the Lagrangian interpolation form4a
If the (i+lj)th cell is a boundary cell, a Taylor series can be used to achieve 2
uij + u3 -
2u ) (x) -[ujij
+ 3u3
-.Ui½j]
(2.42)2.3.5 Calculations with the Modified MAC--SUAC
The flow chart in Fig. 2.13 gives an outline of the computational procedure of the Stanford-University-Modified-Marker-And-Cell Method
(SUMMAC). The general procedure is essentially the same as the one used. in the original MAC method (Welch, et al., 1966) and discussed in Sec. 2.2.3. However, some additional features are needed to implement the modifications presented n Secs. 2.3.1 - 2.3.4.
First, a set of input data is used to specify the key parameters of the problem, such as the mesh size, x , y , t , viscosity, etc.
Then the cell setup is established by assigning to each cell one or more numerical codes to specifyits type. The Ind.t-ial condition is Set by calculating the marker partic]e positions and the u and v distribu-tions at t = 0 . Now, we begin the cycle with a procedure that ref lags
the cells. This procedure has no effect on the first cycle. In the second and future cycles, cells near the free surface must be ref lagged from cycle to cycle because of the continuous change of the surface geometry.
TO calculate by Eq. (2.37), we must know which cells contain an irregular star and the leg lengths n . To establish
this, we sequentially number (from left to right) the particles that mark the free surface. The irregUlar-star cells are defined as those cells that contain surface particles and that have their centers lying
Uj.j
u3
on the fluid side of the free surf ace. Then,, if the free surface is rep-, resented by many short line segments that connect the surface partc1es,. the distance from the cell center to the free surface gives the length of a leg of the irregular star (Fig. 2.7). The actual computation is not as formidable as it appears. The terms D.. and R. are needed in
compu-- 13 13
ting . As usual we obtain D.. by Eq. (2.25) and R. by Eq.
(2.24). Then the values of are solved simultaneously by the iterative method which was outlined in Sec. 2.2.3, except that Eq. (2.37) is used
for an irregular-star cell and Eq,. (2.23) is used for an interior cell. Now data of u and v are updated by Eqs. (2.21) and (2.22). Remember that these velocity components are evaluated on the borders of cells' (Fig. 2..4). In the SU'NAC method we compute u and v only at locations' where they are defined, i.e.,,within the fluid,. For applica-tion of Eqs. (2.21)-(2.24.) near the free surface, the fields of u and v have to be' extended to the empty space above the surface. Equation
(2.39) is used for -this purpose.
Finally, the marker particles are moved to their new locations by the interpolated velocities obtained from Eq. ('2.40). Then the next time cycle begins with a reflagging of the cells and the process cited above is repeated.
2.4 STABILITY CONSIDERATIONS
The most crucial. consideration in the development of a numerical technique is whether or not the proposed' scheme is computationally stable. The finite-difference approximatiän,.nd roundoff in the computer neces-sar.ily introduce some errors in the course of' a èomputation. These
errors may be magnified and accumulated to such an extent that the numbers, obtained-are completely meaningless. However, it is possible to work out schemes that' keep the errors within reasonable bounds.
The.. nonlinearity of, the Návier-StOkes equations that are used in the MAC or SiJIAC methods makes it difficult to' perform rigorous stability analyses. -For linear' finite-difference equations with'constant coeff i-cients, stability can be determined by using a, Fourier' method proposed by, von Neumann (Hirt, 1968). Unfortunately, mpst equations of physical
simple heuristic method was proposed by Hirt (1968) for investigating the computational stability of such nonlinear finite-difference equations.
Hirt's method is based on a rather novel idea. A finite-difference equation is reduced to a differential equation by expanding each of the difference-function terms in a Taylor series. The lowest order terms in the expansion must represent the original differential equation being approxImated. 11 higher order terms constitute the truncation errors caused by the finite-difference approximation. Hirt (1968) showed that the stability of a finite-difference equation can often be determined from an examination of these truncation errors. Following the procedure just outlined above, the finite-difference equations (2.21) and (2.22) can be reduced to the differential equations, respectively, as follows:
Bu
T+;+-+X+X= \V---u
Bu2 Buy I (St 2 5x2Bu\
B2u j)-2-+ Itv_V
(St 2 (Sy2 Bv\ B2u-
4)z
By Buy By2 / (St ' 2(Sx2 Bu\ B2v
+--+---+
Bt Bx By By+g = tv--u
y ', 2--p---
4 Bx/ Bx I (St 2 Ôy2 Bv\ B2v 2 (2.43) (2.44)When comparing Eqs. (2.43) and (2.44) with the Navier-Stokes equations (2.14) and (2.15), we find additional terms n Eqs. (2.43) and (2.44). Those tetms involving (St result from using a first-order
approximation to represent Bu/Bt and Bv/Bt [see Eq. (2.18)], while terms containing (Sx or (Sy stem from evaluating undefined variables
by simple average formulas sucl as Eqs. (2.26). The effec.t of these additional terms is to contribute negative diffusion coefficients, i.e., negative viscosity. It can be seen from. Eqs. (2.43) and (2.44) that the finite-difference scheme will yield growingly unstable solutions if the physically meaningful. viscosity is smaller than the truncat&on error terms. Thus., according to Hirt, the stability conditions for, the MAC
method should be 2u 2 2v' c5t <2 and Sx < '¼X
-(2.45)where u 'is the average maximum fluid speed and eu/ax is the average maximuth velocity gradient in the direction of flow.
Equations (2.45) imply that very fine time steps and space incre-ments 'arerequired in computing low viscosity fluid motions. However, the actual relevance of these conditions to the MAC method Is open to
speculation and further discussion. For example, in our early attempts to apply the MAC method to the study of water waves we were plagued by the fast development of instabilities which were manifested by vary irregular surface profile and field variable distributions. Because of' the small viscosity that was used ( v/dJj = 1.8 x
io6
, it wassus-pected that the MAC method would not be suitable for computing flows with very small viscosity. Then, we tested the same problem by using a very 'large viScOity with v/d/j = 0.10 , which was more than that re-quired to meet the conditions of Eqs. (2.45). The resultIng fluid motion was found to be thuch'-. slower because of the damping effect of a large viscosity. The instabilties remained, however. We were therefore led
to examine the MAC method more carefully. The modifications which have been made were discussed earlier in Secs. 2.3.1 - 2.3.5. We shall
demon-straté the improvement in Sec. 2.5 by comparing the parallel results given by the MAC and the STJNAC methods.
The refinement made in Sec. 2.3.4 (in the evaluation of undefined variables) eliminates the truncation terms involving' x -or y in 'Eqs. (2.43) and (2.4'4). This occurs simply because the refinement
con-sists of using a second-order approximation whose truncatidn error does not contribute to the diffusion terms of Eqs. (2.43) 'and (2.44). Thus, with the SUMMAC method, the' only stability conditions that need consider-ation are these réIat'éd to the select-ion of '6t . This is particularly
true if th" computation 'gOes through a great number of time increments. Because' a large viscOsity was found unnecessary in our numerical experiment,
we presume that t is restricted only by the, stability criterion
CSt 2cSxSy 6x + cSy
where C is the speed of wave. This criterion, called the. Courant con dition, restricts the distance a wave travels in one time step to less than one space interval.. The application of the Courant conditidn to nonlinear systems Of equations is based on a heuristic extension of ]near concepts rather than on direct, rigorous theory.
By using the Courant condition to choose 5t , we have obtained very satisfactory results provided x and 6y arE of such sizes that
important features in the fliuid motiOn ate sufficiently resolved. It
appears that instability problems have been avoided in large part by our ref inethents. in the MAC method. However, we believe that the free-surface modifications, not the reduction of truncation errors, are the key changes leading to stable computations in this case. A self-correcting mechanism, naxely the presence of D terms in Eq. (2.24), operates in the MAC method,. This cotrective effect, that results from setting Dr' 0 in Eq.' (2.24)
n n+1
and from the interaction between the equations used to find , u
and v , tends to cancel out truncation errors, generated 4.n the flow
interior. Hirt's stability analysis that led to Eqs. (2.43) and (2.44) did not include the effect of' D terms. The formidable task of fully analyzing thi, particular feature of the MAC method remains to be done. One way to test the effect of the D terms n Eq. (2.24) is a direct numerical experiment on two parallel computer programs, one of whith con-tains the P terms while the other does not. In the present cases,. com-putations without D terms, i.e.., assuming 0 in Eq. (2.24), were completely unusable.
2.5 RESULTS AND CONCLUSIONS
As an application of the SU1AC method, thE run-up of a solitary wave on a vertical wall has been calculated for a range of. I/d
Laitone's (1960) formulas, Eqs (2 6)-(2 9), were used to compute the
initial conditions of the incident wave. The setup of the computation region is shown in Fig. 2.1.
Theoretically, the length of a solitary wave is infinite and Lal-tone's formulas hold for an infinitely long channel only. Because the computations must be done in a finite domain and the fluid at a distance from the wave crest is essentially still,, it is desirable to define a finite, practical length of the solitary wave. Street and Camf.ield (1967) suggested an effective wave length L by taking L/2 to be the
distance from the wave crest to the section where r = 0.001H .
Sub-stituting this criterion into Eq. (2.3) gives
10.1 (r
(2.67)EquatiOn (2.47) was found to be more stringent than necessary for the present study. The.main consideration is that the two vertical walls which constitute the boundaries of the computation region should be far enough from the initial wave crest so that the motion of a solitary wave into .the still water in frOnt of a vertical wall can be closely
approxi-mated. For this purpose, a less restrictive criterion, say ri = 0.0111
was used in Eq. (2.3) to obtain
=
6.90(-)
(2.48)It is seen that the practical length of'a solitary wave increases as the amplitude decreaSes. For the reason discussed above, the two ver-tical walls have to be located at least L/2 away from the Initial wave crest. When Hid = 0.1 , for example, the value L/2. as given by
Eq. (2.48) is about lid . As shown in Fig. 2.1, the walls were set lOd
from the initial wave crest'. Since L is even shorter for higher ampli-tudes, this setup is adequate for H0/d >0.1 . We have tried to set the left-hand wall at 5d behind the initial wave crest. Thi arrangement resulted, in a slight reduction of the run-up on the right-hand wall because
the moving water in the region beyona the lefthand wall that could have contributed to the run-up was eliminated by the left-hand wall. Thus, it is important nt to set the walls too close to the initialposition of the wave crest.
A mesh of 40 cells in the x-direction and 25 cells In the y
direction was used to represent the computation region. We used (Sx = 0.5
and (Sy = 0.1 for the case H0/d 0.2 . The reasOn for using much
smaller vertical than horizontal spacing i's that we know S priori, that the field variables change more rapidly in
tlle
vertical direction than in the horizontal direction. As a experiment,a
finer mesh, with ox = 0.25 and (Sy = 0.08 , was also used to compute the same case.. The difference in the results was, however, insignificant. The Courant condition, Eq.(2.46), was used to determine the tIme increment (St . The speed of
wave C s approximately (1 + 0.2Y = 1.095
for Hid = 0.2
. Then,from Eq. (2.46),
(St < 2(0.5) (0.1)
1.095 0.5
+ oi
0.152 (2.49.)We chose a moderate value of (St = 0.10 for our computations. We also trIed (St 0.05 and &t.= 0.20 . The use of (St = 0.20 , which
vio-lates Eq. (2.49), resulted in a significantly different solution with a sign of instability in the distributions of the field variables. When (St = 0.05 , we found a small amount of improvement on the solutIon
obtained by using (St = 0.10 . The difference is not significant and we
are led to believe that further reduction in (St is not necessary. For
the viscosity, we used the nondimensional term vidj = 1.8 x
1o6
We were able to obtain stable solutions for this value of viscosity.The quantities Rid rand R/H , as computed by the SITh4MAC method,
are shown n Tab]-e 2.1 for various values of H0/d , They are compared
with experimental data (Street and Camfield, 1966) and Wallace's (1963) theory. in Figs. 2.2 and 2.3. For comparison purposes, the values of RiH0 and Rid
have
also been calculated according to Peregrine's (1967) method using Eqs (2 5), and the results are listed in Table 2 2 These values are plotted on Figs. 2.2 and 2.3. It is seen that the quantityR/H as pedicted by this method is about 2.0 which is little different from the linear theory (i.e., R = 2H by superposition). Therefore, Peregrine's method does not actually give an adequate description of the nonlinear fluid motions for the present case. Thus, the SUfAC method yields.results that are in better agreement with experiment than the
other methods being compared. This is what we expect of numerical methods in which the governing equations an4 boundary conditions are more or less uniformly approximated by the finite-difference scheme.; Wallace's theory diverges from the true solution as the wave nears the wall as a result of the nonuniform approximation on the surface (see Wallace, 1963).
Figure 2.2, in which R/H is plotted against H0/d , is a very sensitive
plot. Better agreement among different methods being compared is found in Fig. 2.3, where R/d is plotted versus Hid
The contour lines of the u , v and fields are shown in
Figs. 2.14, 2.15 and 2.16, respectively. These plots describe the motion of a solitary wave with Hid = 0.5 as it approahes the wall at
t/..JT/
= .5.4 . These are examples of the information that can be obtainedfrom a numerical experiment. The plots presented here were done by hand. We could have used a graphic display device, such as an IBM 2250, to re-veal these internal flow details. However, as our primary phasis was on the development and testing of the SUIAC method, we postponed the
incorporation of this added feature.
in the present problem we used only three rows of marker particles near the free surface. Only the row that marks the free-surface posit-ioi has real significance in our examples and is shown. However, more
par-ticles can be used in the interior Of the fluid if details on the particle movement there are desired. Generally, we use as few particles as
pos-sible to conserve the storage of the computer. Successive stages of wave deformation are shown in Fig. 2.17 for the case H0/d = 0.2 and in Fig.
2.18 for the case H0/d = 0.5 . These free-surface configurations were
obtained by plotting the surface particle coordinates using a CACOMP
plotter. .
As will be discussed in Chapter 3, the velocity distributions described by Eqs. (2.7) and (2.8). are not valid for a solitary wave with
