A PARABOLIC EQUATION METHOD IN POLAR
COORDINATES FOR WAVES IN HARBORS
Haruhiko Kaku and James T. Kirby
Technical Report UFL/COELTR/075
Coastal and Oceanographic Engineering Department
University of Florida
Gainesville, FL 32611
May, 1988
Abstract
This report details the development of a parabolic approximation method for surface
water waves in polar coordinates, and its application to the prediction of wave conditions
inside a breakwater harbor. The parabolic approximation is developed by employing a
coordinate transformation to the linear mildslope equation, after which the approximation
is obtained using standard splitting methods. The linear model is developed in both small
angle and largeangle forms.
The basic linear model is modified to include the effect of wave nonlinearity (using a
Stokes thirdorder formulation) and the effect of bottom boundarylayer dissipation.
The model is employed to model the decay of waves incident on a breakwater harbor
formed by two vertical walls lying on constantangle radii of a circle. Model results are
compared to laboratory data, and reasonable agreement is found between data and the
largeangle nonlinear model.
Contents
1 INTRODUCTION 1
1.1 MildSlope Equation .................... ........... 2
1.2 Parabolic Approximation of the Wave Equation ................ 4
2 THE PARABOLIC APPROXIMATION IN ALTERNATE COORDINATE
SYSTEMS 9
2.1 General Coordinate Transformation of the Reduced Wave Equations .. 9
2.2 Parabolic Approximation in the Alternate Coordinate System ........ 11
3 APPLICATION OF THE PARABOLIC APPROXIMATION IN THE
POLAR COORDINATE SYSTEM 13
3.1 Parabolic Form in Polar Coordinate System ................ .... 13
3.2 MildSlope Equation in Polar Coordinate System ............... 16
3.3 Additional Physical Effects ................... ....... 17
3.4 Applicability of the Parabolic Approximations ................. 20
4 NUMERICAL SOLUTION OF THE PARABOLIC WAVE EQUATIONS 22
4.1 Numerical Scheme for the Parabolic Approximations ............. 22
4.2 Examples of Numerical Solution ................. ...... 27
4.2.1 Effects of Wave Nonlinearity and Laminar Bottom Damping . 27
4.2.2 Noise in the HigherOrder Approximation ............... 28
4.2.3 Effects of a Channel between the Breakwaters ............. 28
5 LABORATORY EXPERIMENT 45
5.1 Facility and Apparatus ..................... ......... 45
5.2 Experimental Procedure .................... ......... 48
5.3 Wave Record Analysis .................... .......... 49
6 COMPARISON BETWEEN NUMERICAL SOLUTION AND LABO
RATORY DATA 53
6.1 Isobe's Experiments .................... ........... 53
6.2 Comparison with Laboratory Data ................. ..... .. 54
7 SUMMARY AND CONCLUSIONS 66
List of Figures
1.1 Comparison of apparent wavenumber direction for each approximation with
exact value ...................................... 8
3.1 Geometry of physical domain ................... ....... 14
3.2 Comparisons between the exact forms (solid line) and the asymptotic forms
(dot line) for a) the Bessel function of first kind, b) the Bessel function of
second kind and c) the absolute value of the Hankel function. . ... 21
4.1 Geometry of computational domain. . . . ..... 23
4.2 Wave amplitude variation in r direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.1m, Ao = 0.01m ......................... 29
4.3 Wave amplitude variation in r direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.lm, Ao = 0.02m ......................... 30
4.4 Wave amplitude variation in r direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.01m. . . . ...... 31
4.5 Wave amplitude variation in r direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.02m. . . . ...... 32
4.6 Wave amplitude variation in 0 direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = O.lm, Ao = 0.01m...................... 33
4.7 Wave amplitude variation in 0 direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.1m, Ao = 0.02m ................... .. 34
4.8 Wave amplitude variation in 0 direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, A = 0.01m. . . . ...... 35
4.9 Wave amplitude variation in 0 direction. Solid line, nonlinear with damping
effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.02m. . . . ... 36
4.10 Wave amplitude contour of the nonlinear solution; Ao = 0.01m, o = 0. 37
4.11 Wave amplitude contour of the nonlinear solution; Ao = 0.01m, 00 = 450 38
4.12 Noise arising in higherorder model with filtering process (Method 1). Solid
line, e = 0; fine dot line, E = 0.05; dot line, e = 0.1; h = 0.15m, Ao = 0.01m. 39
4.13 Noise arising in higherorder model with filtering process (Method 2). Solid
line, e = 0; fine dot line, e = 0.05; dot line, e = 0.1; h = 0.15m, Ao = 0.01m. 40
4.14 Depth contour and breakwater configuration. . . . .. 41
4.15 Wave amplitude variations over the channel after a gap of breakwater. Solid
line, higherorder nonlinear solution; fine dot line, lowestorder nonlinear
solution; dot line, higherorder nonlinear for plane bottom (h = 0.15m);
Ao = 0.01m, ro = 0.615m. Arrows indicate the position of the bank of the
channel . . . . . . . . .. 42
4.16 Wave amplitude contour over a channel; lowestorder approximation. .... 43
4.17 Wave amplitude contour over a channel; higherorder approximation. ... 44
5.1 Geometry of the wave tank for constant depth case (in metric unit) ... 46
5.2 Geometry of the wave tank with channel (in metric unit). . ... 47
5.3 Wave gage response and calibration curve fit by L.S.M. . . .... 48
5.4 Flow chart of the experimental procedure. . . . . ... 50
5.5 Flow chart of analysis ................... ............ 51
5.6 Typical wave record and its energy spectrum. . . . .... 52
6.1 Geometry of the wave tank of Isobe's experiments (in metric unit). . 55
6.2 Comparison of higherorder nonlinear (solid line), higherorder linear (fine
dot line), lowestorder nonlinear (dot line), Isobe's model (dash line), and
laboratory data (x).................... ............. 56
6.3 Plane bottom: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear
approximation (dot line), and laboratory data (x). Test No.l; T = 0.49sec.,
Ao = 0.0085m ................... ................ 57
6.4 Plane bottom: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear
approximation (dot line), and laboratory data (x). Test No.2; T = 0.49sec.,
Ao = 0.017m .. .. .... ... ... ... .. ... .. .. .. ... .. 58
6.5 Plane bottom: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear
approximation (dot line), and laboratory data (x). Test No.3; T = 0.74sec.,
Ao = 0.009m ............ ... ... ........... ... ... 59
6.6 Plane bottom: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear
approximation (dot line), and laboratory data (x). Test No.4; T = 0.74sec.,
Ao = 0.016m ..................... ............. 60
6.7 Plane bottom: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear
approximation (dot line), and laboratory data (x). Test No.5; T = 0.74sec.,
Ao = 0.022m .................... ..... ....... 61
6.8 Channel: comparison of higherorder nonlinear approximation (solid line),
higherorder linear approximation (fine dot line), lowestorder nonlinear ap
proximation (dot line), and laboratory data (x). Arrows indicate the posi
tions of the bank of channel. Test No.6; T = 0.49sec., Ao = 0.0085m. . 63
6.9 Channel: comparison of higherorder nonlinear approximation (solid line),
higherorder linear approximation (fine dot line), lowestorder nonlinear ap
proximation (dot line), and laboratory data (x). Arrows indicate the posi
tions of the bank of channel. Test No.7; T = 0.49sec., Ao = 0.0128m. . 64
6.10 Channel: comparison of higherorder nonlinear approximation (solid line),
higherorder linear approximation (fine dot line), lowestorder nonlinear ap
proximation (dot line), and laboratory data (x). Arrows indicate the posi
tions of the bank of channel. Test No.8; T = 0.74sec., Ao = 0.0105m. . 65
Chapter 1
INTRODUCTION
One of the most important tasks for coastal engineers and scientists is to obtain mathe
matical models which provide accurate predictions of wave conditions in coastal or offshore
areas. One such model which has served as the basis for extensive recent investigations is
the mildslope equation for linear, harmonic surface gravity waves, first derived by Berkhoff
(1972). The mildslope equation describes the simultaneous refraction and diffraction of
arbitrary surface wave motions in a twodimensional domain in situations where the lo
cal water depth varies slowly with respect to the surface wavelength. In its elliptic form,
however, the equation is only suitable for use in regions of small spatial extent due to
practical limitations on computing times, and thus further approximate formulations have
been developed. Fortunately, waves in open water in the coastal zone may be character
ized by a dominant propagation direction, as utilized in ray methods and other refraction
schemes. This geometrical property of the waves can be used as an appropriate reference
frame for determining a coordinate system in which the wave field is described, and also
establishes a framework for developing approximate models. The parabolic approximation
method is based on the restriction that forwardscattered waves propagate principally in
the prespecified direction with all transverse phase variation absorbed into variation of the
complex wave amplitude.
Various types of parabolic approximations for surface wave motion have been developed
and verified in comparison with experimental data. Using a WKBexpansion approach, Mei
and Tuck (1980) derived a parabolic equation governing the leadingorder amplitude for a
forwardscattered wave in constant water depth. This equation was extended to include
nonlinear effects by Yue and Mei (1980), who studied wave diffraction by a thin wedge
and found that nonlinear diffraction of waves at grazing incidence was governed by a cubic
Schr6dinger equation.
For a nonconstant but slowly varying depth, Radder(1979) and Lozano and Liu (1980)
first derived a parabolic equation by using a multiplescale perturbation method. Kirby
and Dalrymple (1983) extended the variable depth approximation to include the cubic non
linearity for Stokes waves, and showed that the nonlinear model predicted the development
of wavejump conditions and reductions in amplitude in the vicinity of caustics in com
parison to linear results. Their model was verified by comparing with the laboratory data
of Berkhoff, Booij and Radder (1982) (Kirby and Dalrymple, 1984). A model of similar
form to that of Kirby and Dalrymple was obtained by Liu and Tsay (1984), who studied
wave focussing over a semicircular shoal. The numerical solutions were compared with the
laboratory data of Whalin (1971) and also showed the importance of nonlinear effects in
the focusing zone. Alternatively, a higherorder approximation was obtained by splitting
method approximation by Booij (1981), which allowed the mathematical model to treat a
wider range of wave directions with respect to the principle propagation direction. Further
analysis and verification of Booij's method (and extension of his model to include the cubic
nonlinearity) can be found in the work of Kirby (1986).
The primary purpose of this study is to outline the development of the parabolic ap
proximation in alternate coordinate systems which are either conformal or nonconformal
with respect to transformation from the Cartesian coordinates. Advantageous coordinate
systems are those in which long coastal structures lie along a coordinate line. In the case of a
conformal transformation, such a boundary then imposes only normalderivative boundary
conditions on the modelled domain, and the domain does not expand or contract with dis
tance along the boundary. Our attention is focused on the mildslope equation or its reduced
form for the constant depth case, and we further consider only applications employing a
transformation to polar coordinates. Additional results employing elliptic coordinates are
being prepared and will be reported separately (Kirby and Kaku, 1988).
1.1 MildSlope Equation
In order to lay the foundation for the derivations below, we provide a derivation of the
mildslope equation in this section. For irrotational, harmonic linear waves propagating
over the seabed, where the water depth is given by z = h(x, y), the velocity potential #
satisfies the Laplace equation,
92 4
V + = 0 ; h(x,y) < z < 0, (1.1)
where Vh = a+ 8 (x, y) are the horizontal coordinates and z is the vertical coordinate.
At the free surface, the combined linearized boundary conditions yield
g;z = ; z =0, (1.2)
where w is the angular frequency of the harmonic wave motion and g is the gravitational
acceleration. The bottom boundary condition is given by
= Vh' Vhh ; z= h(x,y), (1.3)
where Vh is the horizontal gradient operator. We assume that the velocity potential P
can be separated at leading order into a vertical part (the crossspace dependence) and a
horizontal part with a time dependence (the propagating mode):
4)= F(z,h)qt(x,y,t),
where q is the surface displacement which, for a harmonic wave, may be written as
O= _(Xy)eiwt
where 4 is the twodimensional surface fluctuation. The time dependence may be dropped
from the velocity potential, leading to a problem in spatial coordinates alone. We then
rewrite the velocity potential as
P = F(z,h) (x,y). (1.4)
The vertical part F is a function related to local depth h which is assumed to vary only
slowly in space in accordance with h(x,y). Assuming that the model is being used to
represent the propagating components of linear theory, we may take F to be given by
coshk(h + z)
F (1.5)
cosh kh '
where k(x, y) is wavenumber satisfying the linear dispersion relationship:
w2 = gk tanh kh. (1.6)
In order to obtain an equation governing the twodimensional function 4, the Laplace equa
tion (1.1) is integrated over the depth after multiplication with the function F, which gives
f (FV i + Fz) dz = 0. (1.7)
Applying Green's formula to the second component of the integral gives
(Fh #F2) dz = [Fz Frz] z=h. (1.8)
After manipulating (1.8) using (1.2) and (1.3) and substituting the result in (1.7), we obtain
Vh. F2dzVh + 4 F2F dz + k2 F2 dz = 0. (1.9)
h hh
We now assume that the bottom topography changes gradually so that
a2F I aF
V2F = (Vhh)2 + V2h << k2F. (1.10)
The second term in (1.9) may then be neglected. Using (1.5) together with the dispersion
relationship (1.6), the integral of F2 over the depth reduces to
F 2 dz = o cosh2 k(h + z) dz C2 1 + 2kh
F2 dz = dz = 1+ (1.11)
h h cosh2 kh 2g sinh 2kh g
where C = w/k and Cg = C(1 + 2kh/ sinh 2kh)/2. C is the wave phase speed and Cg is the
group velocity. Substituting these results in (1.9) finally gives
V (CCgV) + k2CCgo = 0, (1.12)
where subscript h has been omitted from the gradient operator for further convenience.
This derivation is similar to the one given by Smith and Sprinks (1975).
Let us consider some limiting cases where the equation is exact. For the constant depth
case, C and Cg are constant. Equation (1.12) then becomes the Helmholtz equation,
V20 + k2 = 0, (1.13)
which governs linear wave diffraction. For the shallow water or long wave case, C = Cg,
the mildslope equation reduces to the long wave equation for a single frequency,
gV (hVW ) + w2 = 0. (1.14)
1.2 Parabolic Approximation of the Wave Equation
When we treat the surface wave field in coastal and offshore areas, it is often possible to
restrict attention to waves propagating principally in a uniform direction. The parabolic
approximation of the reduced wave equation is based on the assumption that the variation
of wave amplitude in transverse direction of wave train is much more rapid than in the
propagation direction. As a numerical problem, the parabolic equation method has the
advantage of being an initial boundary value problem which may be solved by marching
over the domain of interest, and relatively rough gridsize may be used.
The twodimensional function (zx,y) in the mildslope equation (1.12) describes the
surface wave motion. 4 can be given by the product of wave amplitude and phase term,
= + ., = ds, (1.15)
where A is a wave amplitude, k is the wavenumber vector and g is the position vector.
We assume that the wave propagation direction may be identified as the x direction.
Recalling the fact that the wavenumber k = k is a function of the slowly varying water
depth h(x, y), the phase function may be defined to be only a function of x after defining
some modified wave number k, which is allowed to vary only in the x direction. (1.13) is
then rewritten as
S= 2Ae" + c.c., = kdx, (1.16)
where A must be a complex wave amplitude to absorb the difference between the modified
and actual phase, which gives
A = Aei('). (1.17)
The parabolic model of the mildslope equation may be derived by defining a small
parameter 6 related to the rate of amplitude modulation and also to the rate of change of
depth over the space of a wavelength. We assume that the derivatives in terms of space
(x, y) have the following explicit scale relations,
a ~ +62 (1.18)
ax ax ax'
9 a
~ 8 (1.19)
ay OY'
where X and Y relate to the slow spacial variations of wave field, which gives
X = 62x,
Y= y.
These scales are based on Yue and Mei's approach. The wave amplitude A thus may be
defined as a function of X and Y.
Kirby and Dalrymple (1983) chose a scale for Vh as 0(62) in order to restrict a bottom
variation to be locally flat up to the second order in 6; consequently, we adopt the following
relationship,
ak (k1.0
a_ 62 k' (1.20)
09X aX'
a(CCg) 2 a(ccg) (1.21)
6 (1.21)
ax ax
Substituting equation (1.16) together with (1.18)(1.21) into the mildslope equation gives
CCg(k2 E2)A+
62 [2icCCgAx + i(kCC,)xA + (CCgAy)y] +
64(CCAx)x = 0. (1.22)
A parabolic approximation may be obtained by neglecting the term of order 0(86). For the
case of constant depth, this equation reduces to
62[2ikAx + Ayy] + S4Axx = 0. (1.23)
The same form as (1.23) may be obtained by substituting equation (1.16) with (1.18) and
(1.19) into the Helmholtz equation (1.13). The O(62) parabolic form of (1.23) is the simplest
approximation of the wave equation as given by Mei and Tuck (1980).
On the other hand, a number of methods for deriving parabolic models have been
proposed to make the models more accurate. Booij (1981) gave a higherorder splitting
method to derive a parabolic approximation of the Helmholtz equation and the mildslope
equation with currents. We apply the same procedure here to the mildslope equation
without current. The differential equation,
a i a )+c~ =O, (1.24)
can be split exactly into two decoupled equations describing a transmitted wave field and
a reflected wave field:
=irf+, (1.25)
iIc, (1.26)
ax
where
S =g + +,
in which superscript + indicates the transmitted wave and indicates the reflected wave.
Considering the special case where the reflected waves may be neglected then gives
S = d. (1.27)
We rewrite the mildslope equation as
s + (CC Z + 2 = 0, (1.28)
where the operator c2 is given by
X' = k2 1 + kCC a(CCg) (1.29)
We define the function as
= co, (1.30)
where e is initially unknown. Substituting (1.30) into (1.24) gives
S+[2'x '. 4+ [ eKIeX 2] 0.\( (1.31)
In the third term of (1.31), the second x derivative of e and the product of ex and .,x may
be neglected on the assumption that they are second order in the slowly varying quantities.
Then, comparing (1.31) and (1.28), the second term of (1.31) gives
& ,z (CC,)x
2 K (CC, '
or
e2
= CC,. (1.32)
Substituting (1.29) into (1.32) yields a pseudo operator given by
S(kCCg)'1 1+ k2CC (CC+,' (1.33)
Then, equation (1.25) becomes
S{(kCC,) [1 + (CC,]y }
3
=ik(CC,) [1+ k2CCgay(CC, ) (1.34)
The pseudo operators ( and i are related to the variables k and (CC,) and to the y
derivatives. The second term inside bracket for both operators, (CCg)/(k'2CCg), is
assumed to be much smaller than unity through the parabolic approximation. Using the
power series expansion and retaining the first two terms then gives
S(kccg) [+ 4k2c (CCcco )j = ik1(CCg) [ + 4k2 (CC,), ,
or
( kCC,)+ (kCC,), (k'CC), 3i
S+ 2kCC, 8k3 (CC,)2 4k' (CC,)2 4kCCg ,
+4k2CC (CC ), = 0. (1.35)
In this equation, a crossderivative term (CCg4y)y, appears which expresses a diffraction
effect in the x direction, and thus partially recovers the effect of the Ox term. The equation
can be approximated for more severe restriction on the size of yderivative terms. In (1.34),
we may consider the approximation for & as
O= [(kCC)1]= 1 + 2CCg (CC Y)]'. (1.36)
This means the z derivative of [1 + (Z g(CCg )]1 is neglected. Then (1.35) reduces to
(k )1 t
S+ 2kCCg ik 1+ + (y(CC, ) (1.37)
2kCC, I k2g
Taking the first two terms of power series expansion for righthand side of (1.37) then gives
+ [k ik] 2kc (CCg'P") = 0. (1.38)
The wave amplitude equation can then be obtained by substituting (1.16) into (1.35) and
(1.38). The lowerorder approximation (1.38) gives
2ikA, + i(CCg +2k(k ) A+ g(CC,Ay)y = 0. (1.39)
This equation has a very similar form to (1.22) after neglecting the term of order O(64).
Since the modified phase function V was applied before the process of parabolic approxi
mation in (1.22), the local wavenumber k was replaced by k. For the case of plane bottom,
however, if the x coordinate is taken to coincide with the direction of bottom slope, then
k = k, and both equations become identical.
The higherorder approximation (1.28) yields a higherorder wave amplitude equation:
[ iCCg,]
2ikA + [i(kcC + 2k(k k) A
Si(kCCg)z ik. 1 ( CA)
4k2(CCg)2 2k2CC, 2kCC (3
1
+ 2 (CCAy), = 0. (1.40)
It is noted that, for constant depth, (1.35) reduces to
3i 1
ik 4 + 4 Y + 4k2 O, = 0, (1.41)
and (1.38) reduces to
Oz ikO 2kiyy = 0. (1.42)
Both equations are the same forms as were derived by Kirby (1986) or Booij (1981). Booij
also derived a criterion for the accuracy of each approximation. The wavenumber k is
expressed by its x and y components, I and m respectively as
k2 12 +m2,
or
1 m i
k 1 (k)]. ](1.43)
Assuming that the y component ofk is very small; i.e., m << k, (1.43) can be approximated
by a lowestorder binomial expansion
1 1m
Ik [1 2 k )2]. (1.44)
0
0
o
1.43
o 1.45
0
6.
L0 1.44
cos'() )
C
(n)
0
0
0.00 30.00 60.00 90.00
sin1 ()
Figure 1.1: Comparison of apparent wavenumber direction for each approximation with
exact value.
This is equivalent to (1.42) in the order of approximation. The form equivalent to (1.41) is
 4 (1.45)
k 1 71_1)2
This is the first Padd approximant of the right hand side of (1.43), as shown by Dingemans
(1983).
The difference between the three equations (1.43), (1.44) and (1.45) can be evaluated
for fixed directions of the wavenumber vector k. The larger the direction of k with respect
to the z axis, the larger the error in each approximation. Figure 1.1 shows a comparison
of each approximation. If we take the maximum error as 5%, 42.60 of wave direction are
allowable for (1.44) and 55.90 for (1.45) (Kirby, 1986).
Chapter 2
THE PARABOLIC
APPROXIMATION IN
ALTERNATE COORDINATE
SYSTEMS
Most existing mathematical models for the spatial evolution of a wave field are carried out
in Cartesian coordinates. In practical situations such as the study of wave fields inside a
harbor or estuary, some difficulties may arise in applying parabolic models in Cartesian
coordinates to the problem because of complicated boundary shapes and also because of
a loss of correspondence between the chosen propagation direction and the actual physical
propagation direction. There are thus advantages to developing the mathematical models in
alternate coordinate systems, which may allow a wide range of problems to be solved easily.
Since the coordinate system may be chosen such that boundaries may be set on constant
values of the coordinates, the problems can be solved with simple boundary conditions.
In addition, the diffraction effects by some obstacles in the wave field cause the diffracted
waves to propagate along the boundaries, that is, the propagation direction may also be set
on the alternate coordinate. A parabolic approximation thus may be derived reasonably.
In this chapter, the mildslope equation and the Helmholtz equation in alternate coor
dinate systems will be discussed.
2.1 General Coordinate Transformation of the Reduced Wave
Equations
We consider the transformation of the mildslope equation (1.12) into an alternate coordi
nate system (u, v), related to (x, y) by,
( x (xu, v) (2.1)
y y (u, v)
The transformation performed here is assumed to be onetoone in the domain of relevant
wave motion. The Jacobian matrix of the transformation is defined by
[dx x v x du (2.2)
dy Yu yv [dv (2
1
The corresponding Jacobian determinant J is given by
J = Xuyv yuXV. (2.3)
The first derivatives of u and v with respect to x and y are given by
Xv
Yu
vx = UT
J
V = (2.4)
= 7
The first and second derivatives in terms of x and y are given by
8 y, 8 yu 8
ax Jau Jav'
a a a (2.5)
ay J au J 8v' )
a2 1 [ a2 a2 a2
2 = J (~)j2 + (Jy~ ) ZJU ( )
22 auav
S = [(Jxz) + (Jz) (2Jzzu,)
ay2 J3 aU2 U 2 auaU
+(Jx,ZuU JuXu J5x,, + J,uZ)
Bu
+(Juzux, Jxvzx + JxuXuv Jvx ) (2.7)
Substituting (2.6) and (2.7) into (1.12) gives the mildslope equation in the alternate coor
dinate system (u, v):
c(CCgo)u 7 [(CCt0u)v + (CCgOv)u] + f(CCgo ),
+J2CC, [(V2u) + (V2v), + k2, ] = 0, (2.8)
where
= y + X,
i=y +xv,
7 = YuYv + XXv,
V2u = u, + utY = [xv (ayuu 2yuv + yvv) yu(cauu 2xuv + xvv)],
V2v = VI + v = [u(a1uu 2u + ) Xu(a 2u )
V2 zz + gy = u y uscu 2,yxuv +)6xvv) xu(ceyuu 2yyuv + Pyvv)] *
For the case of constant depth, (2.8) reduces to
auu 2,yu + fvv + J2 [(V2U). + (V2V)O, + k2 = 0. (2.9)
This is the Helmholtz equation in the alternate coordinate system.
2.2 Parabolic Approximation in the Alternate Coordinate
System
We now consider a parabolization of (2.8) using the WKB approach of chapter 1. Assuming
that waves are propagating in the +u direction so that the phase accumulates along lines
of constant v, the phase function 0 can be obtained according to
S= kds= kx, y du k du, (2.10)
in which s is physical distance along a wave ray in the propagation direction. The phase
function may alternately be defined to be a function of u alone by adopting an appropriate
average over v. Then, the function 4 in new coordinates (u, v) is given by
2= Ae' + c.c., = kdu. (2.11)
Following the procedure in chapter 1, we set the scales as
a a+ 2 a (2.12)
a a
~ (2.13)
Tv 8V*
The wave amplitude A is then given by
A = A(U,V). (2.14)
Substituting (2.11) into (2.8) using (2.12) and (2.13) then gives a wave equation with terms
ordered by small parameter 6:
aCCg(k )2A + J2CCg [(2)ikiA + kA+
[ [jikf3 {2CCgAv + (CCg)vA} + J'CC (V2v)Av +
62 [ai { (CCAkPA)u + CCgk62Au} + P(CCAv)v + J2CCg(V2u)Au +
6' [y {(CCAv)y + (CCAu)v}] +
64 [a(CCgAu)u] = 0. (2.15)
A parabolic equation may be obtained by neglecting the term of order O(64) in (2.15).
However, the resulting equation has a complicated form because of the existence of (and
ambiguous scale position of) crossderivative terms, and is thus much more difficult to
interpret than the usual Cartesian parabolic approximation. On the other hand, for the
case of a conformal coordinate transformation, the CauchyRiemann conditions are satisfied:
Mu = Yv,
Xv = yu.
Therefore,
a= =J,
y = V2u = V2v = 0.
Then, (2.8) reduces to
(CCgu,) + (CC,,), + k2JCCg~ = 0 (2.16)
while (2.9) reduces to
Ouu + ~vv + k2Jo = 0. (2.17)
This is a variable coefficient Helmholtz equation. Also, (2.15) reduces to
CCgA(Jk2 (kJ)2+ 2 i(CCskJA)u + iCCkJAAu + (CCgAv)v +
64 {(CCgAu)u} = 0. (2.18)
A conformal transformation thus leads to a model equation of simpler form. This result
leads to some computational efficiency, and alsoexplicitely maintains the scale relations used
to derive the original parabolic approximation in Cartesian coordinates.
Chapter 3
APPLICATION OF THE
PARABOLIC
APPROXIMATION IN THE
POLAR COORDINATE
SYSTEM
3.1 Parabolic Form in Polar Coordinate System
In the remainder of this study, we will consider the physical domain shown in Figure 3.1,
where two semiinfinite breakwaters are set on the lines of constant angle 0a and Ob. The
waves entering the harbor through the gap are assumed to propagate radially in the +r
direction so that the general concept of a parabolic approximation may be applied in a
polar coordinate system (r, 0).
In a polar coordinate system, r and 0 are related to Cartesian (x, y) with a coincident
origin by
r = (2 + y2),
0 = arctan(),
x
x = rcos0,
y = r sin 0.
(3.1)
The Jacobian determinant and coefficients become
J = XrvY yrXz = r,
a = r2
S= 1,
y = 0,
V2 = 1
r
V2v = 0.
From (2.9), the Helmholtz equation written in polar coordinates is
r2 rr + ro+ ose + (kr)2o = 0.
(3.2)
Breakwater segment
Incident wave crests
Figure 3.1: Geometry of physical domain.
It is noted that the transformation between Cartesian coordinates and polar coordinates
is not conformal. We consider another coordinate system, which represents a conformal
transformation, given by
u = In ( 2 ) = roeu cos 0,
0 = arctan ( y =roeu sin 0, (3.3)
where ro is a constant representing a finite minimum value of r. For this choice of coordi
nates, the Jacobian determinant becomes
J = re2". (3.4)
Using (2.17), the Helmholtz equation in (u, 0) plane becomes
Ouu + 400 + (kroe")20 = 0. (3.5)
The coordinate system (3.3) is equivalent to a polar coordinate system which has been
stretched in the r direction. Equation (3.5) is transformed into (3.2) by letting
r = roeU, (3.6)
then,
Ou = ru/r = rjr,
'uu = r(rrr. + Orb),
and (3.2) is recovered by substituting these values into (3.5).
We now consider the parabolic approximation of the Helmholtz equation in polar coor
dinates. For the assumption that the principal wave direction is in the +r or +u direction,
the splitting method shown previously can be applied to either (3.2) or (3.5). Here, to show
the resulting form briefly, we consider (3.2) rewritten as
Orr + Lr + ,20 = 0, (3.7)
r
where i is a pseudo operator given by
x= k 1+ (k)2 (3.8)
1 (kr)2 92 *
Following the method of section 1.3, the higherorder approximation is given by
Or + ik ] [8r3k + 3ir] 0 + 4( re = 0 (3.9)
The lowestorder approximation is also obtained by following the same procedure as was
done for Cartesian coordinates. The final form is given by
Or + [I ik] 0 ~ 00 = 0. (3.10)
1
3.2 MildSlope Equation in Polar Coordinate System
The transformation of the mildslope equation from Cartesian to polar coordinates may be
carried out in two steps in order to apply the splitting method to a simple equation. First,
the coordinates (3.3) are used to develop the transformation in conformal form. Then,
after being approximated into a parabolic equation, the final model is reached in polar
coordinates.
Following the procedure presented in previous chapter for a conformal transformation,
the first and second derivatives of 0 and (CCg) are obtained according to
z = 7(cos 04u sin q0e),
qy = (sin Ot, + cos Oe0),
(CC,), = [cosO(CC,), sin O(CC,)e],
(CCg)y = 1 [sin O(CCg)u + cos 0(CC,)e],
4XS = UU,
1
Oyy = j4o,
then, the mildslope equation in (u, 0) plane becomes
(CCgOu)u + (CCgoe) + k2JCCg = 0. (3.11)
We now utilize the splitting method to derive a parabolic approximation. Equation (3.11)
is rewritten as
u CC + (g + k2J 1+ kj (CC )] = 0. (3.12)
Under the assumption that the principal wave direction is in +u direction, the lowestorder
approximation is given by
[(kCCgJ 2 ),u
4 + 2kJCCgJ0 2C(CC9) = 0. (3.13)
S2kCCJ'2 2kCCJ2
The higherorder form is
S(kCCgJ)
4u + 2kccgJu ikJ 4
2kCC9JJ
+[(kCCJi) (k2GCJ) 3 ]
8k3(CC)2J 4k4(CCg)2J 4kCCJ (CCs,)
+ 4k2 (CCgo)e = 0. (3.14)
By letting r = roe", the equations are transformed into a polar coordinates space, giving
S+ [ + i kCCk) 2kCCgr+ (CCg4) = 0, (3.15)
1 2kCCg 2r 2kCC~r2
from (3.13), and
[(kCC,), 1 .
[2kCC, 2r 
( kCC), (k'CC,), 3 3i (cC
+ 8k3(CC,)2,r2 4k4(CC,)2r2 8k2CCgr3 4kCCgr2
1
+4k2CCgr2 (CCg O),r = 0, (3.16)
from (3.14).
We assume that the wave amplitude is proportional to the wave surface motion function
0 with an averaged phase function
= 2Aei fdr+c.c. (3.17)
2
We define the modified wave number k according to
( 1 f bk(r, O)dO. (3.18)
) b Of .
Substituting (3.17) into (3.15) and (3.16) gives the wave amplitude equations
2ikCCgAr + [i(kCC)r + ikCC +2kCC,(k ) A+ (CCgAo) = 0, (3.19)
which is the lowestorder approximation, and
2ikCCgAr + [i(kCC,)r + liC + 2kCC(k ) A
i(kCCg)r ikr 3i 1
k2(CCg) + 2k22 + 4k3 (3k k) (CCAe)o
4r2 2k22 4kr3 2kr2
+ 2k (CC, A$)$, = 0, (3.20)
as the higherorder approximation
3.3 Additional Physical Effects
The mathematical models discussed previously do not account for various effects which
modify a physical wave field, including energy dissipation at the bottom or sidewalls, and
wave nonlinearity. With regard to the first effect, Booij (1981) has suggested that the term
iwwo might be added to the mildslope equation, in which w is an unspecified damping
factor. The equation becomes
V (CCgV)) + (k2CC, + iww)= = 0. (3.21)
The effect of the extra term can be illustrated by considering a conservation law of wave
energy. Letting 0 have the form
S= ae (3.22)
2
where a is the real wave amplitude and 0 is the real phase function. Substituting (3.22)
into (3.21) gives
CCV2a + 2iCCgV Va + iCCgaV20 CCa(Vb)2 +
k2CCga + iwwa + V(CC,) Va + iaV(CCg) V V = 0. (3.23)
The real part of (3.23) yields
k2 (Vb)2 + V(CCVa)= 0. (3.24)
CCga
For the case of small diffraction effects, the third term of (3.24) is small, and the eikonal
equation may thus be obtained in the form
k2 = (Vt)2, (3.25)
The imaginary part of (3.23) is
2CCgVp Va + CCaV24 + wwa + aV(CCg) Vb = 0. (3.26)
Using (3.25) together with the vector manipulation:
Cgk = kCg k = kg,
where k = Ikj, then gives
2Cg Va + aV C, + wa = 0. (3.27)
Multiplying this equation by a yields the energy transport equation:
V (a2C) = wa2. (3.28)
The term w thus represents the energy dissipation rate. Dalrymple, Kirby and Hwang (1984)
have investigated the effect of various types of localized energy dissipation and have provided
corresponding expressions for the term w. For the case of energy dissipation resulting from
laminar bottom boundary layer damping, the lowestorder parabolic approximation for
water of constant depth becomes (following Dalrymple et al.)
4ikc3 iv
2ikA, + AY + A=0, (3.29)
sinh 2kh + 2kh 2w (3.29)
where v is the kinematic viscosity. Equation (3.29) is for constant depth case and is identical
to the parabolic form of (1.23) except for the last term, which causes damping of the
wave amplitude. However, it is well known that the bottom boundary layer also causes a
frequency shift of the waves. Liu (1986) used a perturbation method to derive the viscous
damping coefficient, which had both an imaginary part and a real part. The term obtained
by Liu's method may be represented in the mildslope equation by simply letting w be
complex. The corrected form of the parabolic approximation becomes
4ks(1+i) (
2ikA + A + sinh 2kh + V = 0. (3.30)
sinh 2kh + 2kh 2w
In the third term of (3.30), the real part of the damping coefficient causes a phase shift
while the imaginary part causes wave amplitude damping. The magnitude of the viscous
term is relatively small compared to the wave motion itself, and thus, for the case of variable
depth, the term may be added to the parabolic approximation of the mildslope equation
without producing any extra terms concerned with the depth variation.
The amplitude dispersion caused by increasing wave steepness is effective, especially for
regions where focused waves exist, or where the diffraction effects are strong, such as in
the vicinity of the gap in a breakwater. Linear wave models tend to overpredict the wave
amplitude in those areas. Kirby and Dalrymple (1983) derived a parabolic approximation of
wave equation for mildly varying topography and weakly nonlinear waves using a multiple
scale perturbation method developed by Yue and Mei (1980). The equation was based on
Stokes wave theory and had the form of a cubic Schr6dinger equation given by
2ikCCgAz + 2k(k !)CCgA + i(kCCg)zA + (CCgAy), kCCgK'IA12A = 0, (3.31)
where
K, 3C cosh 4kh + 8 2 tanh2 kh
Cg 8 sinh4 kh
(3.31) may be applied within the region where the Stokes theory applies. This region
may be delimited in terms of the Ursell number, Ur, given by
k AI
U (kh)3
U, increases as waves enter shallower water, and, when Ur becomes 0(1), the wave field
is described by cnoidal wave theory rather than Stokes theory. In the region Ur < 0.32, the
Stokes theory gives reasonable results (Isobe et al., 1982). The linearized form of (3.31) has
the same form as (1.39). In the process of coordinate transformation, since the nonlinear
term does not contain any derivative terms, the parabolic models derived previously may
be modified to include this term directly. The final parabolic amplitude equations with
bottom damping effect and amplitude dispersion terms are
2ikCCgA, + i(kCCg), + ikCC + 2kCCg(k ) A
+ 1CCgAe)o kCCgK'IAjIA + kCC+ A = 0, (3.32)
; sinh 2kh + 2kh 2w 0
as a lowestorder approximation, and
2ikCCgAr + [(kCCg)r + ikCC + 2kCC(k A
r
i(kCCg)r ir + 3ik (
4k2CCgr2 2k2r2 4kr3 2 (3k (CCA)
(CCgAe)e, CCK'A2A+ 4k3CCg (1 + i)
CCgA)r kCCgK A sinh 2kh + 2kh = 0, (3.33)
as a higherorder approximation. These equations will be solved numerically and verified
by means of comparisons with laboratory data in following chapters.
3.4 Applicability of the Parabolic Approximations
The Helmholtz equation in polar coordinates (3.2) reduces to a Bessel's equation when the
problem is radially symmetric (1 = 0.). The solution for outgoing waves is given by
S= cH )(kr), (3.34)
where c is constant and H(1) is the Hankel function of the first kind of order zero, which is
given by
H(l)(kr) = Jo(kr) + iYo(kr), (3.35)
where Jo and Yo are the Bessel functions of the first kind and the second kind of order zero
respectively. For large kr the Bessel functions may be approximated as
Jo ~ cos(kr ),
r2 a.
Yo ~ sin(kr ). (3.36)
The function 4 is then given by
= [c ei'r. (3.37)
The term inside bracket may be replaced by a wave amplitude, given by
A(r) = c e
Vnkr
I 1
= cr2. (3.38)
This is a solution of the reduced form of (3.19) for constant depth and for radial symmetry,
which is
ik
2ikAr + A = 0. (3.39)
r
It is thus apparent that the parabolic approximation may only be applied at large kr where
the asymptotic form of the Hankel function is valid. Figure 3.2 shows the comparisons
between the exact forms and the asymptotic forms of the Bessel functions and the Hankel
function. In the plots of the Hankel function, for kr less than 2, a considerable amount of
error can be recognized in the asymptotic form. This effect appears to interfere with the
modelling of some of the longer wavelength cases in Chapter 6, for which kro at the harbor
entrance may not satisfy the far field condition for all the angular components involved.
o
o
0O (a)
0
0
0.00 2.00 4.00 6.00 8.00 10.00
o
C
0
0
Yo (b)
o 2
0
o
0.o00 2.00 4.00 6.00 8.00 10.00
kr
Figure 3.2: Comparisons between the exact forms (solid line) and the asymptotic
forms (dot line) for a) the Bessel function of first kind, b) the Bessel function of
second kind and c) the absolute value of the Hankel function.
21
Chapter 4
NUMERICAL SOLUTION OF
THE PARABOLIC WAVE
EQUATIONS
4.1 Numerical Scheme for the Parabolic Approximations
The parabolic equations obtained previously may be solved numerically as an initial value
problem, with initial values of A specified at r = ro and appropriate boundary conditions
applied at 0 = Oa, Ob. Implicit schemes for the diffusion equation such as the CrankNicolson
scheme are thus applicable. The computation domain, r > ro and a > 0 > Ob, is partitioned
into a rectangular grid where each section has size Ar in the r direction and A0 in the 0
direction. The parabolic equation can then be discretized onto a 6 point segment with 2
points in the r direction and 3 points in the 0 direction (see Figure 4.1). The computation
proceeds in the +r direction by solving n 1st order simultaneous equations at each ith row.
The CrankNicolson scheme is fully spacecentered on the middle point marked by x
in Figure 4.1. Then the first derivative of A with respect to r and second derivative with
respect to 0 are given according to
Aj 1
r A Ar),
82A 2(A)_ [Ai+ 2A+ + A.i) + (A+ 2A + A._)]
a02 2(AG0)2 [ + +iAS1 2+i A1)]
and A is given by
Ai = (A+' + Aj).
2 "
The nonlinear term may be approximated in an iterative manner:
1
Al2Ai = (I +1I2A+ + lA !2).
For the first iteration, A. is used for Ai+1. Then, after the first iteration, the result of
the former iteration is used as A.'1. The number of iterations may be determined by
convergency of A+'1. In the higherorder approximation equation, a crossderivative term,
Aoer, appears. This may be expected to play a role of diffraction effect in the r direction
S. i+1.
rtAr
ro 0. '"O
AO
... ..... i + 1
S+l ......
r 0.A r ... .....
ro j _l
Oi+ "'"
Figure 4.1: Geometry of computational domain.
23
instead of Arr, which was neglected in the process of deriving the approximation. The
numerical scheme for this term becomes
Aoorj = (AO r (A+ 2A+' + Alj) (A 2A + Af ]
This term may introduce spurious behavior in the computational results due to the non
physical behavior of Fourier components with short transverse wave length; see Kirby
(1986). As a result, a filtering process may be necessary. We consider two different methods
of filtering in the numerical solution. The first one involves including the relationship,
A*j = eAj_, + (1 2)A+ + eA+.,, (4.1)
in the scheme. Here, the parameter e is selected as the filter weighting. Then, the first
derivative with respect to r is written as
8A= 1 [(cA+_1 + (1 2e)A.+1 + A )(eA + (1 2e)A + eA+,)]
9r Ar c"I
(Method 1).
For the second method, we also consider (4.1), but the numerical scheme does not
include this relationship. For each grid row, the amplitude variation is obtained using the
basic CrankNicolson scheme. Then, (4.1) is applied to smooth the result before the next
computational step (Method 2).
The initial condition at the breakwater entrance, r = ro, is given by a plane incident wave
which has uniform amplitude and phase variation corresponding to the physical distance
and propagation direction. The phase function is then given by
= k(cos Ooz + sin 0oy)
= kro cos(8 00),
where 00 is the propagation direction of the incident wave. For constant depth at the break
water entrance and under the assumption that phase accumulates along the +r direction
after the entrance, the initial amplitude A(ro, 0) may be given according to
A(ro, O)eikro = aeikro cos(00o)
then,
A(ro,0) = aeikro[cos(e0o)1]
where a is the real value of the wave amplitude. For rigid and impermeable breakwaters,
the boundary conditions are
OA
0 = 0 ; = 0a,eb.
In the finite difference scheme, these are approximated by
Ai+ Ai+l+ A.+
i=1 lj=2 j=n j
Grid size may be determined in accordance with the physical wavelength, and the step size
of a single gridspace should be significantly smaller than a wavelength. In the region far
from origin, because our mathematical models are carried out in a Polar coordinate system,
care is necessary with the size of AO, that is, the computation result may not be well enough
resolved at large r if A0 is taken to be too large at ro. Thus, AO should be chosen carefully
in accordance with the overall computational domain.
In the numerical computations, the NewtonRaphson technique was used to determine
the wave number k from the linear dispersion relationship.
The parabolic models obtained previously yield the difference equation,
S(COE1A*+ COE A*) + (COE2+A*+ + COE2jA*j)
( 2 {[COE3 +1G3'+~Aij+l COE3 +G2i+1A*+ + COE3t1'G1l+'1A+]
+ [COE3j+G3'A,+ COE3'G2' A* + COE3 G1'A ]}
+2() Ar { [COE4 G3i+1A COE4 +'G2s+1A*+ + COE4t iG1+1AA+^]
[cOE4+ zcyG3 COE4 G2IA* + COE4 _,G1 'A_]}
J1* 2 (COE5j A* + COE5A*)
+1 (COE6+1A*,+1 + COE6 A*j) = 0 ,
where
COE1I = 2i(kCC,);,
i(kCC,).
COE2 = i(kCC,)rj + ri + 2(kCC,) (k F),
i(kCC,)r, ik, 3 3ki 
+ +
4(k2CC,)i(r)2 2(k.)2(r,)2 4kf(ri)3 2k.(ri)2'
COE3 ; higherorder approximation
COE3. =
1
S(rt)2, ; lowestorder approximation
COE4 M = 2k)(r ; higherorder approximation
CE0, ; lowestorder approximation
< 0, ;lowestorder approximation
c cosh 4(kh)' + 8 2 tanh2(kh)
(kj) (C 2)  8 sinh4(kh) .'
COE58 = ; coefficient of nonlinear term
0, ; linear solution
S4(kc)3(1 + i)(CC,) / \
sinh2(kh)i + 2(kh)i \ '
COE63 = ; coefficient of bottom damping effect term
0, ; without bottom damping effect
G1' = (CC,)'_I + (CC,)',
G2' = (CC9,)_j + 2(CC,)} + (CC,)+,
G3' = (CC)) + (CC,)z+,,
( (kCC,)_+1 (kCC,))
(kCCg)rj Ar '
kj = (i )
( EAj1 + (1 2E)Aj + eAi+1, ; with filtering process of Method 1
SAj. ; without filtering process
or filtering process of Method 2
For the constant depth case, the corresponding coefficients become,
COE1' = 2ik,
COE2. = i,
/ i 1
4kc(ri)3 (ri)2
COE3 = ; higherorder approximation
(ri)2' ; lowestorder approximation
V V )2
S 2k(ri, ; higherorder approximation
COE42 = 2k)
0, ; lowestorder approximation
C2 8 sinh(kch)
4 C cosh4(kh) + 8 2 tanh2(kh)
COE5 = ; coefficient of nonlinear term
0, ; linear solution
S 4k3(1 + i) V\
sinh 2kh + 2kh '
COE6 = ; coefficient of bottom damping effect term
0, ; without bottom damping effect
4.2 Examples of Numerical Solution
In order to illustrate the effects of various terms in the mathematical models, some numerical
results are discussed here.
4.2.1 Effects of Wave Nonlinearity and Laminar Bottom Damping
Figures 4.24.9 show the comparisons of the solutions obtained from the linear model without
damping effect, the nonlinear model without damping effect and the nonlinear model with
damping effect. The lowestorder of approximation is used and the breakwater configuration
follows the experiments for constant depth as explained in chapter 5, with ro = 1.23m and
0a = 0b = 45. Water depths h = 0.15m and h = 0.lm are examined. The conditions for
the incident waves are T = 0.5sec, and Ao = 0.01m and 0.02m. The kinematic viscosity
v = 1.003 x 106m2/sec was used for all numerical calculations. Figures 4.24.5 show the
wave amplitude variations along lines of constant 0, 0 = 00, 0 = 22.50 and 0 = 450, and
Figures 4.64.9 show the variations along lines of constant r, r = 1.83m (r/ro = 1.49),
r = 2.43m (r/ro = 1.98), r = 3.03m (r/ro = 2.46). The Ursell numbers for each incident
wave amplitude and water depth are Ur = 0.0339 in h = 0.lm and U, = 0.0111 in h = 0.15m
for Ao = 0.01m and Ur = 0.0678 in h = 0.lm and Ur = 0.0222 in h = 0.15m for Ao = 0.02m.
The bottom damping effect is very small in any case. Because the effect is accumulated
with the distance of wave propagation, its effect is more apparent in areas far from the
breakwater entrance. As is evident in the mathematical model, the damping effect term in
the models is a function of depth for constant k and w; then, for the same wave motion, the
effect becomes more pronounced in shallower water. This may be explained by the fact that
the relative thickness of the bottom boundary layer with respect to water depth increases
with decreasing water depth, and thus the damping effect also increases.
The effect of wave nonlinearity is effective for both cases of water depth. This effect
causes amplitude reductions and phase shifts near the peak and increases amplitude near the
breakwaters, indicating an increase in diffraction effects. Especially for higher incident wave
amplitude and larger r area, there are considerable amounts of difference in the maximum
wave amplitude between linear and nonlinear solutions.
The wave amplitude contours of the nonlinear solutions for Ao = 0.01m are also given
in Figures 4.10 for the incident wave direction Bo = 0 and Figure 4.11 is for 8O = 45.
4.2.2 Noise in the HigherOrder Approximation
In the higherorder approximation, components with large transverse wave number appear.
These components propagate in the +0 direction more rapidly than realistic longer compo
nents. Figure 4.12 shows the wave amplitude variations of higherorder approximation along
the 0 direction with the solutions including filtering process Method 1. High frequency wave
components can be seen as unrealistic noise. These components are shifted in phase by the
filtering process, but the magnitudes are not effectively reduced. In the filtering process
Method 1, the relationship (4.1) is included in the difference equation. On the other hand,
with the filtering Method 2, the effect of relationship (4.1) is accumulated with the number
of rows calculated. The high frequency noise is thus filtered (see Figure 4.13). This method,
however, also has the effect of reducing the amplitude of low frequency components since
the filter is not conservative of wave energy. In the section r/ro = 2.46 in Figure 4.13, it
is apparent that in the vicinity of the center line, 0 = 00, the filtered results fall below
the unfiltered ones. The computational results, therefore, tend to underestimate maximum
wave amplitude, especially in areas far from the breakwater entrance, where damping due
to the filtering process has accumulated. In calculations of our models, the filter weight
parameter e is selected between 0 and 0.15 depending on incident wave characteristics.
4.2.3 Effects of a Channel between the Breakwaters
In order to examine the effects of varying bottom topography, we introduce a channel
along the harbor centerline 0 = 00. The configuration of the channel and the breakwaters
follows the experiments described in chapter 5. Figure 4.14 shows the depth contour and
the breakwater configuration (see Figure 5.3 for detail). The wave field after a gap of
breakwaters is affected by means of combined effects of refraction and diffraction. Since
the mildslope equation may be applied on slowly varying bottom, a bank of the channel
which forms the transition of two different depth regions, needs to be mild. In the parabolic
models, where the wave propagation direction is assumed to be principally in the r direction,
the above restriction is effectively more severe in the r direction than in the 0 direction.
In addition, the waves radially propagating toward the r direction after diffraction by the
breakwater edges are refracted toward the shallower region by the channel banks. This
causes the wave number component in the 0 direction to be changed somewhat. Figure
4.15 shows the normalized amplitude variations of the lowestorder and the higherorder
nonlinear models in the 0 direction. The solution without a channel, h = 0.15m, is also
plotted in Figure 4.15. The corresponding amplitude contours are given in Figures 4.16 and
4.17. The plots indicate that the higherorder approximation reduces the wave amplitude
near the line of 0 = 0 faster than the lowestorder approximation. This is due to the
higherorder model's increased ability to allow refraction of waves towards larger angle 0.
o
C0
.*
0
o
o
0
0 = 0
t
").00
o
o
*C
Ao
S
0.40
0.80
1.20
2.00
0 = 22.5
0.40
0.80
1.20
1.60
2.00
0 = 45
0'.0o 0o.4 0.80 1.20 1.60 2.00
r(m) (after breakwater gap)
Figure 4.2: Wave amplitude variation in r direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.lm, A0 = 0.01m.
o..............
.. .. .. .
0
f)
A
Ao
0
0
I o =I u
in
o
1 
Ao
0o
IoI
0
0 = 0
4
00
0.40
0.80
1.20
1.60
2.00
0 = 22.5
o 1
ob. 0
0.40
0.80
1.20
1.60
2.00
0 = 45
0
0
_
Ao lo
.oo0 o0.40 0.80 1'.20 1'.60 2.00
r(m) (after breakwater gap)
Figure 4.3: Wave amplitude variation in r direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = O.lm, Ao = 0.02m.
;I I I .
.......... .
...........
.Z
o = 00
0=0
:oo o:o oso 1.20 1.60 2.00
Ao
Ao
0
0"
.oo o.uo o'.o I'.20 1.60 2.00
0 450
0
T.oo 0.40 0'.80 1'.20 1.60 2.00
r(m) (after breakwater gap)
Figure 4.4: Wave amplitude variation in r direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.01m.
0
in
o
A
Ao a
o
0
*
0 = 0
't. oo
o
.*
IA
lAo
0'
0.410
0.80
1.20
1.60
2.00
0 = 22.5
0
0J
.1
b.oo
0
U)
o
.
Ao
o0.14
0.80
1.20
1.60
2.00
0 = 450
0 I
01
o.oo
0.40
o0.80
1.20
1.60
2.00
r(m) (after breakwater gap)
Figure 4.5: Wave amplitude variation in r direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.02m.
.........
.. .. .. .
.. .. .. .. .. .. ..
.,
ri
'
'


`
'~
o
o
................................
r = 1.49
0 ro
Ao
.o00 10.00 20.00 30.00 40.00 50.00
oro
in
Ao
cboo ib. 00 20. 3 4b.oo s.o00
= 2.46
o
Ao
%.oo Ib.00 2b0.0 3b.oo 0b.00 50.00
0 (degree)
Figure 4.6: Wave amplitude variation in 0 direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = O.1m, Ao = 0.01m.
Ln
S= 1.49
o ro
A
0
'. O0 10.00 20.00 3b.00 u40.00 50.00
oro
62
a00 I0.00 20.00 30.00 i0.00 50.00o
 2.46
o
m.0oo 10.00 20.00 30.00 40.00 50.00
0 (degree)
Figure 4.7: Wave amplitude variation in 0 direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.1m, A0 = 0.02m.
1b.oo 20.00 3b.00 40.00 50.00
).00 10.00
0 (degree)
Figure 4.8: Wave amplitude variation in 0 direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.01m.
 = 1.49
ro
).00
L = 1.98
ro
0
o
o
o
CD
'. oo
 = 2.46
ro
50.00
i
...
i
S= 1.49
oro
oo 10b.oo00 2b.00 3b.00 b0.00oo 5b0.00oo
0
*n
0
= 1.98
Ao
o
C.
.00 10.00 20.00 30.00 q0.00 50.00
r = 2.46
oro
0
0 .. ....
10 10.00 20.00 30.00 40.00 50.00
0 (degree)
Figure 4.9: Wave amplitude variation in 0 direction. Solid line, nonlinear with
damping effect; fine dot line, nonlinear without damping effect; dot line, linear with
damping effect; h = 0.15m, Ao = 0.02m.
Figure 4.10: Wave amplitude contour of the nonlinear solution; Ao = 0.01m, 00 = 0.
37
Figure 4.11: Wave amplitude contour of the nonlinear solution; Ao = 0.01m,
O0 = 45.
38
= 1.49
Ao, ,;"
.OO 10.00 20.00 3.0.00 1. 00 50.00
o
o
0
S.oo0 10.00 20.00 3b.00 b.00 50.00
.
o. r =1.98
Ao
C 7.
Tb.oo ib.oo 20.00 3b.00 .00oo sb.oo
0
o (der2.46
0
SI I I
%.00 Io.oo 20.00 30.00 Do.oo s5o.oo
0 (degree)
Figure 4.12: Noise arising in higherorder model with filtering process (Method 1).
Solid line, E = 0; fine dot line, e = 0.05; dot line, E = 0.1; h = 0.15m, A0 = 0.01m.
~" = 1.49
AA''
.00o 10.00 2b0.00 30.00 '0.00 50.00
S= 1.98
C ro
.A
Ao
C)
ID
.o00o 10.00 2. 30.00 0.0 .oo 50.00o
C 2.46
0 0
oo.oo 10.00 20.o0 3.0o ob.oo s5b.oo
0 (degree)
Figure 4.13: Noise arising in higherorder model with filtering process (Method 2).
Solid line, E = 0; fine dot line, e = 0.05; dot line, e = 0.1; h = 0.15m, Ao = 0.01m.
h = 0.075m
breakwater
channel
breakwater
nciaent waves
Figure 4.14: Depth contour and breakwater configuration.
h = 0.15m
h = 0.075m
" = 2.22
3.00
 = 3.44
ro
20.00
__ I_
... .. .. ... .. .... .. ..
10.00 20.00 30.00 40.00 50.00
= 4.66
ro
0.00 10.00 20.00 30.00 'b.00 5b.00
0 (degree)
Figure 4.15: Wave amplitude variations over the channel after a gap of breakwater.
Solid line, higherorder nonlinear solution; fine dot line, lowestorder nonlinear so
lution; dot line, higherorder nonlinear for plane bottom (h = 0.15m); Ao = 0.01m,
ro = 0.615m. Arrows indicate the position of the bank of the channel.
0
o
Ao
0
01
00oo
Z ,!
I I
Figure 4.16: Wave amplitude contour over a channel; lowestorder approximation.
43
Figure 4.11: Wave amplitude contour over a channel; higherorder approximation.
44
Chapter 5
LABORATORY EXPERIMENT
In order to verify the accuracy of the parabolic models described in the previous chapters,
a series of experiments were performed at the Coastal and Oceanographic Engineering
Laboratory, University of Florida. The experiments fell into two categories depending on the
type of topography used. The first category consisted of tests in constant water depth. The
wave field in the region between breakwaters is then described by the Helmholtz equation
with constant coefficient in Cartesian coordinates. For the second case, a channel was
included in the region between the breakwaters, and the mildslope equation was expected
to describe the wave field.
5.1 Facility and Apparatus
Figure 5.1 illustrates the wave tank configuration for the tests of the constant depth case.
The wave tank was 7.24m x 7.24m (23.8ft x 23.8ft). The tank topography consisted of
a horizontal part and a sloping part. The water depth at the wavemaker was 0.45m and
decreased with a slope of 1 : 10 up to the horizontal part, which had a water depth of
0.15m. The wooden breakwaters, which were lin thick, were placed on the horizontal part.
The breakwater enclosed a 90 sector with a gap of 1.74m width. The other ends of the
breakwaters were connected to the corners of the tank so that wave propagation around the
ends of the breakwaters could be prevented. Sand, rocks and glued fiber mats were placed
at the endwall and the corners of sidewall and breakwaters in order to damp reflected
waves inside and outside the breakwaters. A flaptype single paddle wavemaker generated
the incident waves, whose direction was normal to the gap. A polar coordinate grid was
drawn on the bottom. Each breakwater was placed on a line of constant 0, and the origin
of the coordinates was on the horizontal part.
For a second test of the parabolic model, a submerged channel was built in the wave
tank as shown in Figure 5.2. The geometry of the channel was symmetric about the center
line 0 = 0. A 1 : 3 sloping bank was used as the transition between the two constant depth
areas h = 0.15m and h = 0.075m. The channel was 0.87m wide at its deepest part, which
was also the gap width between the breakwaters.
Sample stations were then established at fixed locations on the (r, 0) grid. Wave records
were obtained at each predetermined sampling station using capacitancetype wave gages
mounted on steel tripods. The analog signal from the wave gage was converted into a digital
signal by means of an A/D converter driven by a MICRO PDP11 digital computer. The
computer also controlled the measuring system, including calibration of gages, computa
Wavemaker
2.69
I
1.74
2.69
2.69 0.87
( I 1.
Figure 5.1: Geometry of the wave tank for constant depth case (in metric unit).
Breakwater
Breakwater
r = 0.075 I h = 0.075
0.225 0.225
2.96 0.87 2.96
I I
Figure 5.2: Geometry of the wave tank with channel (in metric unit).
3.185
0.435
0
(cm)
0
o
Z
(cm) o
0
0
I0_ I I !____ I
150.00 175.00 200.00 225.00 250.00
A/D OUT x10
Figure 5.3: Wave gage response and calibration curve fit by L.S.M.
tion of calibration coefficients by leastsquare method, timing the sampling frequency and
period, and storing the data into a specific file. Since the capacitancetype wave gages used
did not give a linear response with variation of water surface elevation, a quadratic approx
imation was utilized for gage calibration. Figure 5.3 shows plots of values of gage readings
versus water surface elevation during static calibration, and the quadratic approximation
determined by leastsquares method. The calibration curve agreed with gage values with
errors less than 0.05mm.
5.2 Experimental Procedure
Because of the limited size of the tank, some physical effects caused by the tank boundaries
such as crosswise rereflected waves and currents on the down wave side of the breakwaters
were generated by the incident waves and often grew with time. A modulation of the
incident wave amplitude possibly due to the reflected waves or to complex motion of water
behind the wave maker was also observed several minutes after turning on the wave maker
and, therefore, every sampling was taken within the first minute of each run. Up to 5
minutes with the wavemaker off was then allowed between two consecutive samples in order
to let the water calm down. Figure 5.4 shows a flow chart of the experimental procedure.
A 20Hz sampling frequency was chosen for each sample and the sampling period was
Table 5.1: Summary of experimental conditions
Test No. Tank Bottom ro(m) Wave Period(sec) Ao(m) kAo Ur
1 0.49 0.0085 0.144 0.0089
2 0.017 0.288 0.0176
3 Plane 1.23 0.74 0.009 0.077 0.0364
4 0.016 0.137 0.0647
5 0.022 0.188 0.0888
6 0.49 0.0085 0.144 0.0089
7 Channel 0.615 0.0128 0.217 0.0132
8 0.74 0.0105 0.090 0.0425
selected between 15sec and 30sec depending on the period of the incident waves and on the
location of the sampling station. The longer the wave period was, the longer the sampling
period should be in order to recognize unexpected disturbances in enough number of waves.
The stations far from the gap and near the breakwater sometimes required longer sampling
period due to the presence of relatively larger amplitude modulations relative to the incident
wave amplitude. Table 5.1 shows a summary of conditions for each experiment.
5.3 Wave Record Analysis
Because of the existence of long waves at the natural frequency of the tank and higherorder
harmonics, the wave record analysis must resolve the extra frequency regions not considered
in the theory, although the waves beyond the gap showed quite stable trend in time and
space. In analysis of the records, a Fourier transform was applied to isolate the fundamental
frequency band, and then wave amplitudes were calculated for each station. Figure 5.5
shows a flow chart of wave record analysis. A typical wave record and its energy spectrum
are shown in Figure 5.6. It is apparent that the record includes higher harmonics as well
as some broadbanded noise. Since the theory here is for the amplitude of the fundamental
frequency component, it is appropriate to filter out the range of frequencies corresponding
to bound harmonics of the fundamental component generated by the wavemaker. Selecting
the frequency band 1 H. ~ 3 H, in Figure 5.6(b) and transforming into timedomain gives
the modified wave record shown in Figure 5.6(c). This procedure is seen to reduce the
apparent modulation of the wave amplitude, leading to more stable estimates of uniform
wave height.
We note that several of the experiments were run with incident wave steepnesses large
enough to promote the onset of BenjaminFeir sideband instabilities. The form of the energy
spectrum displayed in Figure 5.6(b) is consistent with the presence of growing side bands;
this mechanism could account for most of the remaining amplitude modulation apparent in
the filtered wave record shown in Figure 5.6(c).
read gage signal
and
surface elevation
convert digital signal
into surface elevation
Figure 5.4: Flow chart of the experimental procedure.
Figure 5.5: Flow chart of analysis.
51
6.00 9.00
w (H.)
12.00 15.00
(c)
12.00 15.
12.00 15.00
T(sec)
Figure 5.6: Typical wave record and its energy spectrum.
( )
(cm)
17
(cm)
3.00
i,
n I
Chapter 6
COMPARISON BETWEEN
NUMERICAL SOLUTION AND
LABORATORY DATA
In this chapter, we discuss the wave phenomena on the down wave side of the breakwaters
by means of comparisons of the mathematical models and two sets of laboratory data. The
first set of data is taken from the paper of Isobe (1987). The second set of data was obtained
in this study, as described in the previous chapter.
6.1 Isobe's Experiments
The first set of data considered here was obtained by Isobe (1987), who studied the wave
refraction and diffraction problem by employing a rayfront coordinate system which is
determined by considering not only refraction effects by a bottom topography but also
diffraction effects by structures in the relevant domain. The geometry of the wave tank
is shown in Figure 6.1. Isobe derived a parabolic approximation in rayfront coordinates
following the method of Tsay and Liu (1982). The resulting equation is equivalent to
equation 3.19 in the order of approximation, which is a lowestorder linear approximation.
Since we use a polar coordinate system and the boundaries are placed along the lines of
constant 0, our computational domain is restricted to have straight boundaries at each side
of the breakwater harbor. We thus omit the second parts of the breakwaters, which are
perpendicular to the shoreline, from Isobe's experimental configuration. The computational
domain has the r distance the same as the first parts of the breakwaters, which are inclined
to the shoreline. The domain is shown in Figure 6.1 by dash lines. The conditions of
the incident waves are 9.1cm wave height, 0.83sec wave period and 180 wave propagation
direction in the deep part of the tank. At the entrance of the breakwaters, the refracted
incident wave angle is 16.030 using Snell's law. The incident wave directions in our domain
at the first row of the domain are not constant, because the entrance line of the domain is
not parallel to the depth contour. The computed results, however, show that assuming a
constant value 16.030 as the incident wave direction of the first row makes little difference
from the result using values at each point from Snell's law. The results shown in this
chapter are obtained by using constant incident wave direction 16.030. Considering the
validity of the nonlinear models, at the entrance of the breakwaters where the water depth
is approximately 0.2m, the incident wave has the wave steepness and the Ursell number as
kAo = 0.3 and U, = 0.13 respectively using the incident wave amplitude in the deep part of
the tank, which is used in the numerical computations. The Ursell number indicates that
the nonlinear models may be applied to the problem.
Figure 6.2 shows the comparison of the parabolic approximations and the laboratory
data at the line BB'. The numerical results of our models were obtained by taking corre
sponding points near the line BB'. The lowestorder nonlinear approximation gives higher
wave amplitude than the other approximations and laboratory data for the most part, es
pecially on the left side of the peak, where the assumed wave propagation direction differs
from the real one the most causing the worst error in the lowestorder model. The higher
order approximations give reasonable predictions in this area, and the nonlinear solution
shows an appropriate amplitude reduction in the vicinity of the peak, where Isobe's model
and the present linear higherorder model overpredict the results. These results point out
the importance of nonlinear effects in determining the pattern of wave height in a diffracted
wave field.
6.2 Comparison with Laboratory Data
The second set of experiments was described in chapter 5. The laboratory data and corre
sponding numerical solutions of the higherorder nonlinear, the higherorder linear and the
lowestorder nonlinear models are given in Figures 6.36.10. In the higherorder numerical
solutions of longer wave period, the large wavenumber components still exist, although the
filter weight parameter E = 0.15 for filtering process was used. However, the qualitative
features of the solutions are apparent enough to compare to the laboratory data for most
cases. Thus, we did not use larger c values in order to avoid losing physical meaning in the
solutions.
In general, for the constant depth case, results of the higherorder approximation agree
well with laboratory data. In the section r/ro = 1.38 (Fig. 6.36.7), the laboratory data
indicate that the diffraction effects occur at smaller 0 than the values predicted by the
lowestorder model and are well predicted by the higherorder approximation except the
values at 0 = 00 for T = 0.49sec, where the laboratory data give much smaller amplitude
than any numerical solutions. This may be due to the reflected waves from the end wall of
the tank and the incident waves forming a modulated wave envelope along the line 0 = 0.
The laboratory data along the line 0 = 00, thus, may include larger error than that of the
other position, even though a long period of time was used for wave record sampling. In
the section r/ro = 1.87 for T = 0.49sec (Fig. 6.3 and 6.4), the higherorder approximation
predicts the wave amplitude very well within 0 < 0 < 200 and slightly underpredicts the
laboratory data at 0 = 30 and 0 = 40. In the section r/ro = 2.2 for T = 0.49sec, the lower
incident wave amplitude case (Fig. 6.3) shows that the laboratory data has milder amplitude
variation in the 0 direction in the vicinity of the center than any numerical solutions while
the higher incident wave amplitude case (Fig. 6.4) shows reasonable agreement with the
higherorder nonlinear or lowestorder nonlinear approximation.
In Figures 6.56.7 (T = 0.74sec), it is shown that the difference between the lowestorder
approximation and the higherorder approximation is more obvious at the section of small
r than at the large r. The large modulation of amplitude with 0 at small values of r is
thought to be a result of problems with small kro and invalidity of the Hankel function
asymptote. For this case (T = 0.74sec), kro = 10.4. The wave propagation direction at the
large r is supposed to be closer to the direction of the r coordinate than at the small r so
0.60.50.8
WHtt
Figure 6.1: Geometry of the wave tank of Isobe's experiments (in metric unit).
2.3
0.9
1.4 9.0
0.9
3.5
3.6
h=0.4
Ih=0.24
_I
x ....I. ".X
o
0
1i.5o 1.00 b.so '.00oo 0'.so '.oo '.so
B (m) B'
Figure 6.2: Comparison of higherorder nonlinear (solid line), higherorder linear (fine dot
line), lowestorder nonlinear (dot line), Isobe's model (dash line), and laboratory data (x).
that the error involved in the lowestorder approximation becomes small. The laboratory
data are well predicted by the higherorder approximation at r/ro = 1.38 and r/ro = 1.87.
The effects of wave nonlinearity is not readily apparent for the T = 0.74sec cases.
When a channel exists behind the breakwaters, the wave propagation direction is affected
by refraction. Figures 6.8 and 6.9 show the comparison of the numerical solutions and the
laboratory data for T = 0.49sec. In the section r/ro = 2.22, there are obvious differences
between the lowestorder approximation and the higherorder approximation between 0 =
100 and 0 = 250, where the bottom profile changes from the flat part to the slope part
of the channel. The lowestorder approximation shows higher wave amplitude in this area
than the laboratory data while the higherorder approximations agree well with the data.
The difference of the real wave propagation direction from the assumed one, which is the
r direction, is not expressed properly by the lowestorder approximation. In the section
r/ro = 3.44, the numerical solutions agree well with the laboratory data for Ao = 0.0085m
(Fig. 6.8). Here wave nonlinearity is not effective enough to be examined in comparison with
the laboratory data, and neither is the difference between the lowestorder and the higher
order approximation. For Ao = 0.0128m (Fig. 6.9), however, the higherorder nonlinear
approximation seems to give the best fit. In the section r/ro = 4.66 for both incident
wave amplitude cases (Fig. 6.8 and 6.9), the laboratory data shows slightly different trend
from any approximations. Data points at 0 = 150, 200 and 30* for Ao = 0.0085m, and
at 0 = 150 and 20* for Ao = 0.0128m have much larger values than those in the section
r/ro = 3.44. This effect is not explained in terms of the forward propagating waves covered
by the parabolic approximation.
 = 1.38
ro

Ao I
0
0 ...... ro
I AoI I
o
0
o oo 1.o00 20.00 30.00 40.00 5.o00
*2
C =2.2
CD ro
3 X
o
.0OO 10.00 20.00 3b.00 40.00 50.00
0 (degree)
Figure 6.3: Plane bottom:. comparison of higherorder nonlinear approximation
(solid line), higherorder linear approximation (fine dot line), lowestorder nonlin
ear approximation (dot line), and laboratory data (x). Test No.1; T = 0.49sec.,
Ao = 0.0085m.
Co  = 1.38
A
Ao
X
.c00 10.00 20.00 30.00 40.00 50.00
S= 1.87
AO
A o I .... "'.4.
0
o.00 1b.00 2b.00 3b.00 4b.00 50.00
AOr = 0.017m.
0 .
%.oo 10.00 20.00 30.00 40.00 50.00
0 (degree)
Figure 6.4: Plane bottom: comparison of higherorder nonlinear approximation
(solid line), higherorder linear approximation (fine dot line), lowestorder nonlin
ear approximation (dot line), and laboratory data (x). Test No.2; T = 0.49sec.,
Ao = 0.017m.
.:... .. = 1.38
o ro
'.0oo 10.00 2b.00 3b.00 4b.00 5b.00
C.
o  =1.87
or
o 
0
T.oo 1.o00 20.00 3b.00 40.00 sb.00
Ao
0 = 2.2
.........
o C',
0
b.00 b10.00 2b.00 3b.00 40.oo 50.00
0 (degree)
Figure 6.5: Plane bottom: comparison of higherorder nonlinear approximation
(solid line), higherorder linear approximation (fine dot line), lowestorder nonlin
ear approximation (dot line), and laboratory data (x). Test No.3; T = 0.74sec.,
Ao = 0.009m.
C
Un
o x = =1.38
.o 1 20. 3 ro
Ao
o00 1i0.00 20.00 30.00 40.00 50.00
CD
0 = 1.87
o ro
S .. .............. ........
.2.2
oro
.oo0 10.00 20.00 30.00 40.00 50.00
0 (degree)
Figure 6.6: Plane bottom: comparison of higherorder nonlinear approximation
(solid line), higherorder linear approximation (fine dot line), lowestorder nonlin
ear approximation (dot line), and laboratory data (x). Test No.4; T = 0.74sec.,
Ao = 0.016m.
0
S.= 1.38
.0 10.00 20.00 30.00 40.00 50.00
C3 ..2.....2.
to
.. X
o
AO
0
A = 0.022m.2.2
o ro
1.00 10.00 20.00 30.00 40.00 50.00
ear approximation (dot line), and laboratory data (x). Test No.5; T = 0.74sec.,
A0 = 0.022m.
Results for the case of wave period T = 0.74sec is given in Figure 6.10. The higherorder
approximations give large amplitude modulation at the section r/ro = 2.22 and r/ro = 3.44.
For this case kro = 5.27. The parabolic approximation for this condition is thus no longer
reliable. In the section r/ro = 4.66, the numerical solution give smooth line, however, the
filtering process is accumulated through the computational domain, the result may not be
reliable since the earlier results of the domain do not have reasonable predictions.
0
o'
o
ao
o
.0oo
IA I
Ao
S= 2.22
ro
XX
Xx
1b.oo 2b.00 3b.00 40.00 50.00
" = 3.44
ro
.00
 = 4.66
ro
%t.oo b10.00 20.00 30.00 40.00 5C.00
0 (degree)
Figure 6.8: Channel: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear ap
proximation (dot line), and laboratory data (x). Arrows indicate the positions of
the bank of channel. Test No.6; T = 0.49sec., Ao = 0.0085m.
''
CI
Ao
0
0
0
C,
*
o
o_
CD
0D
'.00
Ao
00 20.00 30.00
 = 2.22
ro
).00
S= 3.44
ro
b1.00 2b.00 3b.00 40.00 50.00
S= 4.66
ro
'U.00 10.00 20.00 30.00 40.00 50.00
0 (degree)
Figure 6.9: Channel: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear ap
proximation (dot line), and laboratory data (x). Arrows indicate the positions of
the bank of channel. Test No.7; T = 0.49sec., Ao = 0.0128m.
I I
Ao
0X
.
.o00 1b.oo 2b.00 3b.00 0b.oo00 5.00
o
U)
S = 3.44
0o r
g ;" = 4.66
Ao
......... .
.... . ....7....... .. ....
X
.o00 10.00 20.00 30.00 40.00 50.00
0 (degree)
Figure 6.10: Channel: comparison of higherorder nonlinear approximation (solid
line), higherorder linear approximation (fine dot line), lowestorder nonlinear ap
proximation (dot line), and laboratory data (x). Arrows indicate the positions of
the bank of channel. Test No.8; T = 0.74sec., Ao = 0.0105m.
Chapter 7
SUMMARY AND
CONCLUSIONS
The objective of this study was to develop the parabolic approximation of the mildslope
equation and the Helmholtz equation in a polar coordinate system for application to diffrac
tion or combined refractiondiffraction problems for the wave field inside a harbor. A general
coordinate transformation of the reduced wave equation was introduced, and the parabolic
approximation of the mildslope equation in the alternate coordinate system was obtained
using WKB approach. Due to the existence of crossderivative terms resulting from the
transformation, the resulting parabolic equation assumes an inconvenient form. A restric
tion of the transformation to the conformal case reduces the complexity of the approximation
to the simple form found in the Cartesian case.
The parabolic approximations of the mildslope equation and the Helmholtz equation
for the transmitted wave field in a polar coordinate system were derived using a splitting
method. A lowestorder and a higherorder approximation were obtained. Wave nonlinear
ity and bottom damping effects were added in the models as a heuristic correction drawn
from previous results. The solution to the reduced form of the parabolic approximation for
the radially symmetric case was related to the large kr asymptote of the Hankel function,
and the applicable domain was restricted to that in which the asymptotic form of the Han
kel function was valid. The CrankNicolson technique was used to solve the approximations
numerically.
Numerical results indicated that the bottom damping effects were negligible for the
experimental conditions we examined, and the wave nonlinearity caused reduction of am
plitude near the centerline of the breakwater entrance and increased amplitude near the
breakwaters, due to enhanced diffrraction effects. The higherorder model was found to
increase the rate of amplitude reduction due to diffraction at small r in comparison to the
lowestorder model.
Laboratory experiments were conducted to verify the parabolic approximations. The
comparison between the parabolic models and the laboratory data indicated that the higher
order model predicted the wave amplitude variation well, especially when the effect of wave
nonlinearity was included in the model. The lowestorder model was least accurate in the
region where the wave propagation direction had large difference from the assumed one.
This trend was apparent in the vicinity of the harbor entrance and over the bank of the
submerged channel.
The results obtained by using the higherorder model were seriously contaminated by
large amplitude modulations which were too long to be sufficiently damped by the filtering
process. The numerical results for very small kro, where the asymptotic form of the Hankel
function is not obtained, lose physical meaning; this effect is accentuated in the higherorder
model, where the discrepancy between the near field form of the governing equation and the
parabolic model seem to enhance the generation of the spurious lateral Fourier components
in the higherorder model.
The numerical solutions were also compared with a model and laboratory data from
Isobe (1987). The comparison indicated the difficulty in applying the lowestorder model in
a polar coordinate system to the problem because of the difference between the real wave
propagation direction and the assumed one, while Isobe's model in a rayfront coordinate
system was in more reasonable agreement with the laboratory data. The higherorder
models, however, compared well with Isobe's data and showed improved prediction when
wave nonlinearity was included in the model.
Problems with the higherorder approximation resulting from the generation of spurious
lateral modes (Kirby, 1986), coupled with the fact that the higherorder approximation
provides only an incremental increase in angular validity of the parabolic method (see
Chapter 1) have led recently to the development of methods based on Fourier decomposition
which have no explicit angular restriction. This approach has been detailed for the case of
Cartesian coordinates and a flat bottom by Dalrymple and Kirby (1988) and is presently
being investigated in the context of the conformal coordinate transformation for application
to harbor problems.
Bibliography
[1] Berkhoff, J. C. W., Computation of combined refractiondiffraction, Proc. 13th
Coastal Eng. Conf., 471490, 1972.
[2] Berkhoff, J. C. W., N. Booij, and A. C. Radder, Verification of numerical wave propa
gation models for simple harmonic linear waves, Coastal Eng., 6, 255279, 1982.
[3] Booij, N., Gravity waves on water with nonuniform depth and current, Rep. 811, Dep.
of Civ. Eng., Delft Univ. of Technol., Delft, 1981.
[4] Dalrymple, R. A., J. T. Kirby, and P. A. Hwang, Wave diffraction due to areas of energy
dissipation, J. Waterway Port Coastal Ocean Div. Am. Soc. Civ. Eng., 110, 6779, 1984.
[5] Dalrymple, R. A. and J. T. Kirby, Models for very wideangle water waves and wave
diffraction, J. Fluid Mech., 192, 1988.
[6] Dingemans, M. W.,Verification of numerical wave propagation models with field mea
surements: CREDIZ verification Haringvliet, Rep. W488, part 1, Delft Hydraul. Lab.,
Delft, 1983.
[7] Isobe, M., A parabolic refractiondiffraction equation in the rayfront coordinate system,
Proc. 20th Coastal Eng. Conf., 306317, 1987.
[8] Isobe, M., H. Nishimura, and K. Horikawa, Theoretical considerations on perturbation
solutions for waves on permanent type, Bull. Fac. Eng., Yokohama Nat. Univ. 31, 2957,
1982.
[9] Kirby, J. T., Higherorder approximations in the parabolic equation method for water
waves, J. Geophys. Res. 91, 933952, 1986.
[10] Kirby, J. T., and R. A. Dalrymple, A parabolic equation for the combined refraction
diffraction of Stokes waves by mildly varying topography, J. Fluid Mech., 136, 453466,
1983.
[11] Kirby, J. T., and R. A. Dalrymple, Verification of a parabolic equation for propagation
of weaklynonlinear waves, Coastal Eng., 8, 219232, 1984.
[12] Kirby, J.T. and H. Kaku, Parabolic approximations for water waves in conformal co
ordinate systems, in preparation.
[13] Liu, P. L.F., Viscous effects on evolution of Stokes waves, J. Waterway
Port Coastal Ocean Engrg., 112, 5563, 1986.
[14] Liu, P. L.F., and C. C. Mei, Water motion on a beach in the presence of a breakwater,
1, Waves, J. Geophys. Res. 81, 30793084, 1976.
[15] Liu, P. L.F., and T. K. Tsay, Refractiondiffraction model for weakly nonlinear water
waves, J. Fluid Mech. 141, 265274, 1984.
[16] Lozano, C., and P. L.F. Liu, Refractiondiffraction model for linear surface water
waves, J. Fluid Mech., 101, 705720, 1980.
[17] Mei, C. C., and E. O. Tuck, Forward scattering by long thin bodies, SIAM
J. Applied Math., 39, 178191, 1980.
[18] Radder, A. C., On the parabolic equation method for water wave propagation,
J. Fluid Mech., 95, 159176, 1979.
[19] Smith, R., and R. Sprinks, Scattering of surface waves by a conical island,
J. Fluid Mach., 72, 373384, 1975.
[20] Tsay, T.K., and P. L.F. Liu, Numerical solution of waterwave refraction and diffrac
tion problems in the parabolic approximation, J. Geophys. Res., 87, 79327940, 1982.
[21] Whalin, R. W., The limit of applicability of linear wave refraction theory in a conver
gence zone, R.R.H713, U.S. Army Corps of Engrs., WES, Vicksburg, 1971.
[22] Yue, D. KP., and C. C. Mei, Forward diffraction of Stokes waves by a thin wedge,
J. Fluid. Mech., 99, 3352, 1980.
