|
REPORT DOCUMENTATION PAGE
1. Report No. 2. 3. Recipient's Accession No.
4. Title and Subtitle 5. Report Date
A Program for Calculating the Reflectivity
of Beach Profiles 6.
7. Author(s) 8. Performing Organization Report No.
James T. Kirby UFL/COEL-87/004
9. Performing Organization Name and Address 10. Project/Task/Work Unit No.
Coastal and Oceanographic Engineering Department
University of Florida 11. contract or Grant so.
336 Weil Hall N00014-86-K-0790
Gainesville, FL 32611 13. Type of Report
12. Sponsoring Organization Name and Address
Office of Naval Research Miscellaneous
14.
15. Supplementary Notes
16. Abstract
This report documents a program for computing the reflection of linear
water waves propagating over a one-dimensional, arbitrary topography. The
restrictions of linear theory are assumed to apply. Dissipative effects and
wave breaking are also neglected. The program is used in several examples
pertaining to the design and construction of artificial, reflective bar fields
on mildly-sloping beaches.
17. Originator's Key Words 18. Availability Statement
Nearshore bars
Wave reflection
19. U. S. Security Classif. of the Report 20. U. S. Security Classif. of This Page 21. No. of Pages 22. Price
Unclassified Unclassified 47 pp
UFL/COEL-87/004
A PROGRAM
FOR CALCULATING THE
REFLECTIVITY OF BEACH
PROFILES
BY
JAMES T. KIRBY
March 1987
Sponsor:
Office of Naval Research
Abstract
This report documents a program for computing the reflection of linear water waves prop-
agating over a one-dimensional, arbitrary topography. The restrictions of linear theory are
assumed to apply. Dissipative effects and wave breaking are also neglected. The program
is used in several examples pertaining to the design and construction of artificial, reflective
bar fields on mildly-sloping beaches.
1
Contents
1 Introduction 1
2 Theory 2
3 Numerical Approximations 6
4 Model Tests and Examples 10
4.1 Reflection from sinusoidal bars ................. ....... .. 10
4.2 Reflection from arrays of semicircular bars .................. 15
4.3 Artificial bar field on a sloping beach . . . ..... 24
4.4 Discussion of computational limitations . . . ..... 30
5 Program Listing 33
6 References 43
List of Figures
1 Definition of depth components ......................... 2
2 Sinusoidal bottom, m = 4, D/h = 0.32. Data from Davies and Heathershaw
(1984) . . . . . . . . ... 11
3 Sinusoidal bars as in Figure 2. 450 angle of incidence . . ... 12
4 Reflection and transmission for geometry for Figure 2; varying wave period
and angle of incidence........... .......... ......... 14
5 Reflection from row of semicircular bars. m = 4, regular spacing . 16
6 Bottom variation spectra for regularly spaced semicircular bars. ...... .. 17
7 Effect of irregular bar spacing on reflection by semicircular bars. . 20
8 Bar geometries corresponding to Figure 7. . . . .... 21
9 Spectra for irregular-spaced bar configurations. . . . ... 22
10 Reflection from regular array of semicircular bars, 2m spacing . ... 23
11 Configuration of artificial bars on a sloping beach . . ..... 24
12 Reflection and transmission for regularly spaced bars on a sloping beach... 26
13 Reflection and transmission for irregular bar spacing on a sloping beach. .. 27
14 Reflection from a sloping beach. 80m bar spacing. . . ... 28
15 Reflection from a sloping beach. 20m bar spacing. . . ... 29
16 Reflection and transmission for regular and irregular bar spacing on a sloping
beach, LLW+1.61m. ................... .......... 32
1 Introduction
The amount of wave energy arriving at a shoreline after propagation over offshore bathymetry
is affected by a number of mechanisms leading to loss of energy flux relative to values mea-
sured at an offshore station. The mechanisms of bottom damping have long been under
investigation, and may significantly reduce the amount of wave energy arriving at the surf
zone. Recently, the possible influence of the reflectivity of the beach profile has come under
closer scrutiny. In particular, several recent studies (Davies and Heathershaw (1984); Mei
(1985)) have shown that undular bottom features mimicking shore-parallel bars can lead
to strong reflection of incident surface waves. The existence of this mechanism provides a
means for approaching the design of artificial, bar-like structures intended to significantly
reduce wave energy arriving at the shore.
This report documents a program which may be used to calculate the reflection of
linear water waves by a specified offshore topography. By building up results for individual
frequency and directional components, the program can compute the transfer function for
the reflection of a spectrum with directional and frequency spread. The theoretical and
numerical approaches are summarized in sections 2 and 3. Section 4 then provides several
examples of the use of the program. A program listing and documentation are provided in
section 5.
2 Theory
The theory which serves as the basis for development of the computer program described
here is taken directly from the work of Kirby(1986). Kirby developed an extension to the
mild-slope equation of Berkhoff (1972) which allows the equation to be applied to bottoms
having a slowly-varying mean depth on which is superimposed a rapidly-varying but small -
amplitude depth change. The smallness of the rapid variation allows the bottom boundary
conditions to be expanded about the slowly varying mean depth. Following the notation
in Kirby, we take the total depth h'(x, y) to be expressed by
h'(x, y)= h(x,y) (x,y)
where h satisfies a mild-slope condition and S satisfies a small-amplitude condition. The
split of the total local depth into components is illustrated in Figure 1. The resulting
Figure 1: Definition of depth components
extension to the mild-slope equation is then given by
itt V (CCgV) + (w2 k2CCg)O + -9 V (6VO) = 0. (2)
cosh kh
Here, 4 represents the velocity potential in the plane of the free surface, and is related to a
three-dimensional velocity potential 0 according to
) = cosh k(h + z) ,y,t) (3)
cosh kh
Equation (2) is valid for a single frequency component of a linear wave train, with frequency
w and wavenumber k related by
w2 = gk tanh kh (4)
The remaining coefficients are determined from (4) and are given by
c = k (5)
C9 1 + sinh2kh (6)
8k 2k sinh 2kh
In the absence of currents, 4 is also simply related to the surface dispacement r through a
simple constant of proportionality, and hence t7 may be substituted in (2) in place of i. We
make this substitution. We then consider the case of purely harmonic motion and introduce
a spatial surface displacement q according to
tr(x, y, t)= @(x, y)e-iwt (7)
This substitution reduces (2) to the form
V (CCgV) + k2CCg cosh V (V) = 0. (8)
cosh2 kh
We now restrict the model to the case of one dimensional topography;
h'(x) = h(x) 8(x) (9)
The full elliptic problem (8) is further reduced to
CC, (V2 + k2) + (CCg,)zz -cosh2 ,skhS, = 0. (10)
cosh2 kh cosh kh
We now treat the two-dimensional surface i as a spectrum of directional wave components
which individually refract over the slowly-varying topography h(x) according to Snell's law.
For a given component which is incident at angle Oo in deep water, we have the relation
k sin 0 = ko sin 00 (11)
where
02
ko = (12)
g
We denote local wavenumber vector components according to
1 = kcos0 (13)
m = ksin0 = ksin00 (14)
where m is evidently a constant for all x. We also obtain the relation
k2 = 12 + 2 (15)
Noting that m < ko from (14), if decaying modes are to be neglected, we may write
formally as
Iko
f_= ( m)ei dm (16)
J-ko
Substituting (16) in (10) then yields a second-order ODE for the 7(m);
2 [, +c gm2
(CC,,), cosh2 h + csh2 Ih = 0. (17)
cosh2 kh I cosh2 kh
Equation (17) gives a well-posed problem for the reflection of each individual direction
component after the specification of boundary conditions. Formally, the problem may be
posed in the interval -oo < x < 0, where x = 0 represents the shoreline and the wave is
incident from the deep ocean at -oo. This problem is intractable. In order to simplify the
problem, we suppose that the incident wave condition is known at some finite distance xa
from the shoreline. Further, we neglect the region of the surfzone and take a second station
x2 to represent a position between the region of topography of interest and the breaker
zone. We then assume that wave energy transmitted past station x2 is lost to breaking.
This assumption implies the neglect of direct reflection at the shoreline, which may become
an important effect on steep beaches or on shores fronted by seawalls or other reflective
structures. We thus consider only the problem of (17) posed in a finite domain xz < x < x2
together with boundary conditions given at x1 and x2; the resulting reflection coefficients
will not contain the effect of shoreline reflection.
Boundary conditions are developed in the form of radiation conditions. At x2, we
assume that the wave is propagating out of the finite domain towards the shoreline, in the
+x direction. We thus take
tz = il(x)J; = x2 (18)
At the offshore station x2, the wavefield is composed of the incident component i, which
is assumed to be known, and a reflected component jr, which must satisfy a radiation
condition for propagation out of the domain:
,rx = -il(x))r; x = X1 (19)
Noting that
r = O (20)
we substitute (20) in (19) and obtain
= il(24i 4); x = z1 (21)
(17) together with (18) and (21) fully specifies the problem to be solved.
1
3 Numerical Approximations
We now develop a finite-difference scheme for solving the problem posed by (17), (18) and
(21). We drop the^superscript notation for ri(x, m) and further define the notations
p = CCg (22)
7 (23)
[= -cosh2 kh
We then discretize the domain xa < x < X2 into n points according to
x' = zx + (i- 1)Az; 1
where
n2 1
A- = (25)
n-1
All other coefficients and the dependent variable are defined in discrete form at the grid
points z'. We then develop a centered finite-difference scheme for (17) which is given by
2Ax2
2Az2
+ [(i)2p' + m2"i8i] i = 0. (26)
(26) may be written in the simple form
Aii-1 + Bii+ Ciri+1= 0; i= 2,...,n 1 (27)
where
Ai = pi+p-1 (S + '-1) (28)
B' = -(p'+l + 2pi + pi-1) + it(+1 + 25 + i-1)
+2Azx [(li)2pi + 2,i'] (29)
C' = p'+1 +p' i(+ + 5i) (30)
We note that, in the special case where waves are started in shallow water and then prop-
agated to deeper water, I given by
I = V/k2 m2 (31)
may become imaginary, which represents a turning point in the mathematical problem and
a caustic in the physical refraction problem. The problem is generally formulated to handle
this case by allowing the appropriate variables to have complex form.
We now turn to the boundary conditions. In order to simplify the application of the
boundary condition and the determination of reflection and transmission coefficients, we
assume that the input topography h(xz) has been modified to give a flat bottom at the
edges of the domain;
h1 = h2
h"-1 = h" (32)-
In keeping with this restriction, we require 6(x) to have compact support in the domain
zx < x < x2, so that the waves at the radiating boundary are not interacting with rapid
bed undulations. We thus require
81 = 82 = Sn-1 = Sn = 0. (33)
We then finite difference the transmitting boundary condition (18) in centered form accord-
ing to
n Mn-1 iln n in
Xl- t= + 7n-1 (34)
Ax 2 2
where In-1 = In due to (32) above. Let
a" = AIn. (35)
2
(34) may then be written as
ln(1 an) = n-l(1 + an) (36)
or, in the notation of (27),
B" = 1 a
A" = -(1+a") (37)
The reflective boundary condition is handled in a similar way. Defining
a = 2A1, (38)
we obtain
(1 + al')2 (1 a)1 = 2a(m)a1 [e2a' + 1] (39)
where we have assumed that the incident wave is specified by
li(x) = a(m)eilt(- 1) (40)
We then have
B1 = -(1-a1) (41)
C1 = l +a (42)
D1 = 2a(m)a1 [e2a + 1] (43)
The entire problem may then be written in the form of a linear matrix equation
Arq =D (44)
where D is a column vector with D2 D" = 0, rt is a column vector with elements r17 ,n,
and A is a tridiagonal matrix with diagonal vectors A', Bi and C'. The solution is obtained
using the double sweep algorithm as given by Carnahan, Luther and Wilkes (1969).
After determining the solution for rt, the reflection and transmission coefficients must
be determined. At xl, we obtain two estimates of the reflection coefficient R as follows.
Using (40), we may write rtr at xz as
f = rl -a(m)
r = ^-a(m~e2 (45)
We then take
a(m)
R2 1 (46)
a(m)
R is then estimated as the simple average of R1 and R2. Similarly, at x2 we obtain estimates
of transmission coefficent T according to
17"n-1
T, I -I-rI
a(m)
T2- (47)
2 (m) (47)
and then average Ti and T2 to obtain T. A test of the accuracy of the solution is obtained
by checking the conservation of energy requirement
R +T2 2c_ 0 = 1. (48)
The results for an entire frequency and directional spectrum may be built up by re-
peated application of the numerical scheme developed here. A program corresponding to
the present numerical algorithm is documented in section 5 below.
4 Model Tests and Examples
In this section, we provide an adequate set of tests of the present program in comparison
to previous results, and then test several cases of interest in the problem of constructing an
artificial nearshore bar field. The first set of tests pertain to the familiar problem of wave
reflection by a sinusoidal bar field, as studied recently by Davies and Heathershaw (1984),
Mei (1985) and Kirby (1986).
4.1 Reflection from sinusoidal bars
Davies and Heathershaw have obtained a detailed set of measurements of wave reflection
from a sinusoidal topography given by
6(x) = Dsin(Ax); 0 < x < L (49)
where A = 2r/lb, lb is the bar length, and L = mli. The bar field thus contains m
complete sinusoidal bars. Tests were run with 1b = lm, D = 0.05m and m = 2,4 or
10. Mean depth h(x) was uniform so that the bars were superimposed on a domain with
otherwise constant depth. Wave period T and water depth h were varied during the tests to
obtain the range of parameters desired. Figure 2 shows computed reflection (solid line) and
transmission (dashed line) coefficients for a case with m = 4 and h = D/0.32 = 0.15625m
for normally-incident waves. Also included in the figure is the laboratory data of Davies and
Heathershaw. Reflection associated with the presence of a beach downwave of the ripple
patch is estimated to be on the order of 0.2, accounting for much of the discrepancy between
theory and data away from resonance. Agreement between the approximate theory of Davies
and Heathershaw and the present numerical results is quite close, deviating the most at the
resonant peak. (This fact is probably due to Davies and Heathershaw's treatment of the
resonance in their theory; see the discussion in the introduction of Kirby (1986).) Agreement
between present numerical results and other sets of laboratory data presented by Davies
and Heathershaw was also good and additional results are thus not shown.
Existing theories for resonant (Mei, 1985) and non-resonant (Miles, 1981) reflection of
1.00
0.80
0.60
0.40 -
0.20 -
0.00 -
0.00
0.60 1.20 1.80 2.40
2k/ h
Figure 2: Sinusoidal bottom, m = 4, D/h =
(1984)
0.32. Data from Davies and Heathershaw
('-, - -
I
I
I 0 =0
I m =4
I
D=0.32
h
.I *
1 *
I *s *
AS
T/V 0** *
0MA^.
0-
II.
*e
gO
3.00
obliquely incident waves over topography S(x) all predict a zero in the incident-reflected
wave coupling at an angle of incidence of 45. This result indicates that the bar field is
transparent to waves at this angle, at least at leading order. Figure 3 shows an example of
results for a 450 angle of incidence and the geometry of Figure 2. The reflection coefficient
is non-zero but is not sufficiently large to reduce transmission by a significant amount.
1.00 r
0.80 1-
0.60
9 =450
m =4
S=0.32
h
0.40 I-
0.20 H-
0.00 L
0.c
r%~ -~I ~- ~ V I
0.60
1.20
2k/
1.80
2.40
3.00
Figure 3: Sinusoidal bars as in Figure 2. 45 angle of incidence
Figure 4 shows a contour plot of reflection and transmission coefficients for a range
of wave periods and angles of incidence for the geometry of Figure 2. The region of effec-
)0
CD
LD
LY)
Cf)
CD
Lr)
0.2
I i o
CD
Ln
I I I I c
O00800"0900 "T 00 "O 00 0
I3ONU 43AUM
Figure 4. a) Reflection coefficient
CO
LD
oo
SC 3
CZ)
]HOZ
angle of inci
r'i- ^ ) ) I- ;
0009008090fTO0O7 O0 0f
b) Transmission coefficient
Figure 4: Reflection and transmission for geometry for Figure 2; varying wave period and
angle of incidence.
tive reduction of transmission is localized to the patch 0 < 0 < 450 around the resonant
wavenumber, except for regions of strong reflection at large angles of incidence.
4.2 Reflection from arrays of semicircular bars
Construction of an artificial bar field which mimics the reflective behavior of a patch of sinu-
soidal bars will likely involve the placement of several discrete structures whose longitudinal
axes run parallel to shore and whose on-offshore spacing is controlled. In this section, we
look at several examples of bar fields consisting of rows of semicircular structures placed on
an otherwise flat bottom. The case of bars placed on a sloping bottom is left to the next
section.
For a given run, all bars are assumed to have the same radius r; 8(x) for x in the interval
of a bar location is then given by
(x) = [r2 I Xbl']1/2 (50)
where Xb is the bar center location. We first test a regularly spaced field of 4 bars with a
spacing of 1 m and a water depth of h = 0.15625m. We take r = 0.05m. This case thus
mimics the sinusoidal geometry for Figure 2, with bar crests having the same spacing and
the same vertical projection above the horizontal bottom. Results are shown in Figure 5,
where we retain the same definition for A as in the previous section.
The results for the field of semicircular bars exhibit several interesting features. The
reflection peak at 2k/A = 1 is still present and represents the Bragg-interaction of the
wave train with the fundamental spacing of the bar field. A second, stronger peak is
present at 2k/A = 2 and represents the Bragg interaction of the wave train with the first
superharmonic spacing of the bar field. These results need to be interpreted in light of
the Fourier transform of the constructed depth profile. The basic periodic interval for the
transform analysis is indicated in Figure 6a, and the FFT amplitude spectrum corresponding
to the interval is shown in Figure 6b. The spectrum of the bottom drops off quite slowly
above the fundamental component, which would be the only component present in an ideal
1.00
-- -.~ ~ N -
0 =00
0.80 4 semicircular bars
r=0.05m
1 m spacing
0.60 -
0.40 -
0.20 -
0.00
0.00 0.60 1.20 1.80 2.40 3.00
2k/
Figure 5: Reflection from row of semicircular bars. m = 4, regular spacing
, '
I-)1 mp ei 1 m bIo1 1 m--i-
a) periodic bottom interval
CD
1)
O
LD
CD
CD0
0
O
(CD
LO
0S -
0
00 0.0060.0080.00100.00
.0O 20.0 0O. 0 60. O0 80.0 100.0
N
b) amplitude spectra of bottom variation
Figure 6: Bottom variation spectra for regularly spaced semicircular bars.
1 1 H il l il li l
/^\
~1.. ....l.i I,~.
sinusoidal bottom. This indicates that the reflection response should be spread over a broad
band of wavenumbers, with intensified response occurring near the spectral peaks of the
bottom variation. (Compare to the case of a single discrete rectangular bump, for which
the spectrum is nearly white and the reflection coefficient varies nearly sinusoidally with
wavelength). The drawback to this particular type of broadening of the reflection response
is that much of the bottom variance is dedicated to reflecting wave components which are
outside the range of wave periods where reflection is to be concentrated. It would be more
efficient from a design point of view to devise a bottom configuration which concentrated
the bottom variance in spectral components corresponding to the band of surface waves
of interest. The practicality of constructing such arrangements, which would tend back
towards the sinusoidal configuration, is doubtful.
The amplitude variance associated with the first spectral component of the regularly-
spaced bar configuration corresponds to a sinusoidal bar with amplitude 0.008m. Reflection
from such a bar field was computed and is shown in Figure 5 as the dash-dot line. It is appar-
ent that most of the reflection near peaks of the response is associated with Bragg-reflection
from the individual bottom spectral components, but that non-resonant interaction with
other components of the bar field can overlap the region of resonant response, causing a
shift in the apparent peak amplitude. Kirby (1986) has shown that the reflection in the
vicinity of a resonant peak for one particular bottom component is almost equal to the
simple sum of the reflection contributions from all spectral components. (See, in particular,
Figures 5 and 6 in Kirby (1986).
The results above raise the question of whether bar fields consisting of relatively abrupt,
isolated structures should be constructed with a dominant spacing equal to the dominant
surface wavelength (or some other spacing), rather than to half the dominant wavelength
as would be the case for a simple sinusoidal bottom. This hypothesis will be tested further
below.
Another feature needing investigation is the effect of spacing the bars irregularly in
order to introduce an increased variety of spectral components in the bottom topography,
thus broadening the reflection response. Figures 7a-d show 4 examples obtained by altering
the spacings of the 4 bars corresponding to the Figure 5 example in such a way that the
total bar field width of 3m is maintained. The result for uniform spacing is included in
each figure as the dashed line. The bar field geometries corresponding to the 4 cases are
indicated in Figure 8.
Spectra of bottom variations for the 4 irregularly spaced bar configurations are shown
in Figure 9. The effect of irregular spacing is to produce a more densely populated bottom
spectrum, thus increasing the number of spectral components which individually take part
in Bragg-interactions with a broad surface wave spectrum. It is not immediately apparent
from the plotted spectra that an optimum choice of spacing distortion can be made, and
thus this problem remains to be treated from an analytic point of view.
We now return to the question of constructing a bar field with a dominant spacing other
than the 1/2- wavelength spacing indicated by the direct Bragg-resonance mechanism. We
first test a spacing of one wavelength (relative to the Bragg condition in Figures 2 and 5),
which, for the regularly spaced bars of Figure 5, corresponds to a 2m uniform spacing. The
resulting reflection response is plotted in Figure 10. Comparing these results to Figure 5,
we see that the main effect of lengthening the bar field by a factor of two is to sharpen
the resonant tuning associated with harmonics of the bottom variation. The resonant peak
at 2k/A = 1 (with A = 27r as in Figure 5) now corresponds to the first harmonic of the
basic periodicity, the fundamental now being located at 2k/A = 0.5. No shift in maximum
amplitude of the resonant peaks is achieved. The response may be characterized as being
less effective than the original Figure 5 response if a narrow band around 2k/A = 1 is the
principal concern, or more effective due to the extra peak at 2k/A = 1.5 if broad-spectrum
response is the major consideration.
a) 1.1-1.0-0.9m
0.2 I
0 I l I
0 1 2 3
b) 1.2-1.0-0.8m
0.2-
IRI 0 1 2 3
c) 1.3-1.0-0.7m
0.2 -
0 1 2 3
d) 1.4-1.0-0.6m
0.2
II
0 L v \ v
0 1 2 3
2k/ X
Figure 7: Effect of irregular bar spacing on reflection by semicircular bars.
----1 m --------m M ---- m
Figure 7a
--0.9m -- l--- m -- 1.1m
Figure 7b
4- 0.8m -- l----1m
1*-- 0.7m -- ----1
*- O.Xm- -1-----0
C>
m I1
Im 4- 1-4
1.2m---
Figure 7c
3m---
Figure 7d
m---
Figure 8: Bar geometries corresponding to Figure 7.
->* 1
C
C'
C
C
cr
c
CD
c
Ln
CD
cn
CD
CD
Ci
Lf
CD
J
3
1
^
)
)
)
)
>
0.00 20.00 40.00 60.00 80.00 100.00
(c) (d)
CD
CD
r'j
a:
o
.0 I 0 00 0 00 0 0 0 0 1 0 0 0 0 ,, 0.00 ,,,. 00 60. 00 0. 0. 0L
I 0 M
0.00 20. 00 40.00 GO. 00 80. 00 100.O00 0.00 20.00u40.0060.0080.00100.00
Figure 9: Spectra for irregular-spaced bar configurations.
o
,)
.-I
cr o
CD
in
m
ICI D
0.00 20.00 40.00 60.00 80.00 100.00
N
o
_1
d0 ----- --\. .. ?
0.80
0.60
0.40 -
0.20-
0.00
0.00 0.60 1.20 1.80 2.40 3.00
2k/X
Figure 10: Reflection from regular array of semicircular bars, 2m spacing
1.00
4.3 Artificial bar field on a sloping beach
We now turn to a particular application of the present theory to a proposed bar field
installation at Point Mugu beach. The profile geometry and bar field configuration are
taken initially from the report by DeVries (1987). The only initial alteration to the proposed
configuration is to impose a semicircular cross-section on the bars. The bar radius is taken
to be 2m (N 6 ft), corresponding to the 6 ft diameter proposed for circular obstacles. The
initial configuration is indicated in Figure 11. The computational scheme was initialized
hg =2.315m
2m
4.265m
f---- 40m ---
<4 130m t
Figure 11: Configuration of artificial bars on a sloping beach
by taking a section of beach extending 5m shoreward of the shoreward bar center and 5m
offshore of the offshore bar center. For an initial configuration of 4 bars at a nominal uniform
spacing of 40m (~ 120 ft), the horizontal extent of the computational domain is 130m. The
LLWL depths at the onshore and offshore ends are 2.315 m and 4.265 m, respectively, and
the beach slope is 0.015.
An initial run was made with no bars in place to determine the reflectivity of the finite
extent of the tested profile. Reflection was calculated for a range of wave periods from 5
to 20 seconds. The maximum reflection in this range was on the order of 0.02 and is not
dynamically significant.
The first set of runs with bars was made with a water level of LLW + 0.61m, which
corresponds to the level of the first profile break 100m seaward of the high-tide beach face.
Figure 12 shows the reflection coefficient (solid line) and transmission coefficient referenced
to the shoreward depth (short dashes) for the initial configuration of 4 bars at an equal
spacing of 40m, over a range of wave periods from 5sec to 20sec. The main Bragg resonance
peak is located near T = 14sec, and first and second superharmonic peaks are apparent
at T = 7.5sec and T ; 4.5sec. The reflection response has a sharp drop off in the range
T = 10 llsec, and the initial configuration is thus not well suited to the proposed storm
condition with T = 12sec. (Response centered on the annual mean of T = 14.69sec is
significantly better). 1
In order to shift the Bragg-resonance peak more towards the 12sec storm period and
broaden the overall response, a staggered array of 4 bars was tested next, with the shoreward
and seaward bars at the same location and an internal spacing of 30m 40m 50m. The
response resulting from this spacing is shown in Figure 13. Staggering the spacing had the
desired result of broadening the response characteristics with no concurrent reduction of the
magnitude of reflection. This approach thus appears to be quite successful in the sloping
beach application.
Two additional runs were performed in order to test the shifting of bar spacing as in the
previous section. Figure 14 shows the reflection response for 4 bars at a uniform spacing
of 80m (a rough one-wavelength/bar space in the 14sec band). The increased number of
spectral peaks in the response is apparent as in the previous section, but the overall response
drops off. This is likely to be due to the increase in depth over the seaward bars, which
renders these bars less effective.
In Figure 15, the results for a regular spacing of 20m are shown (again using 4 bars).
The main Bragg resonant peak is now shifted to ~ 8sec, and there is no spectral peak in the
10 15sec band. This case is not applicable to the West Coast site under investigation, but
'(Wave data values from DeVries, 1987)
1 II
r
r
/)
0.00
5.00 10. 00 15.00 20.00 25.00
WRVE PERIOD I(SEC)
Figure 12: Reflection and transmission for regularly spaced bars on a sloping beach. Four
bars, 40m spacing, LLW+O.61m
CO
0
0c
0
C0
a
CD
0
a
0
I /
---------------~
% ./
1- ----
II I~
Ii
0.00 5. 00 10.00 15.00 20.00 25.00
WHVE PERIOD(SEC)
Figure 13: Reflection and transmission for irregular bar spacing on a sloping beach. Four
bars, 30 40 50m spacing, LLW+0.61m.
CD
CD
CD
CD
CD
CD-)
CD
CD
a
0
0
CD
CD
0D
CD
1 11
^-I 1^--4-- - --- -
C'O
0 -
CD
CD
0-
CD
O
a
CD 'V
CD
ED -\-
0. 00 5. 00 10.O 15.00 20. 25.00
HRVE PERIOD (SEC)
Figure 14: Reflection from a sloping beach. 80m bar spacing, LLW+O.61m.
.rl
. -/- - L-- -
'I (
'I
t
*
''
5. 00 10. 00 15.00 20.00 25. 00
NRVE PERIOD (SEC)
Figure 15: Reflection from a sloping beach. 20m bar spacing, LLW+0.61m.
CD
Co
C
E D
0C
CD
O
c-I
CD
C1)
CD
CD
C_'
CD
O
ED
ED
0
0.00
I I -----------
gives an indication of what might be achieved on an East Coast site, where wave periods
are shorter. In this case, quite a strong reflection can be achieved due to the fact that water
depth over the bar field does not significantly change over the short spatial interval of the
bar field.
Finally, the effect of tidal elevation on the response of a given configuration must be
considered. Increasing tidal elevation causes greater submergence of the bar field and thus
reduces the effectiveness of the reflection process. The change in depth/bar spacing also
shifts the Bragg-resonant peak relative to a fixed range of wave periods. Figure 16 shows
the response of the 4 bar array with 40 40 40m spacing and 30 40 50m spacing for a
tidal elevation of LLW+1.61m. For the regular spacing, the peak is shifted down to ~ 13sec
from 14sec, and the overall magnitude of the reflection is dropped to ~ 0.2 from 0.4. For
the irregular spacing, the response is still broad in comparison to the regularly-spaced bars,
but the overall magnitude of the reflection is similar to the regular-spacing case, and the
downshifts in wave period are also similar.
4.4 Discussion of computational limitations
For the bar fields considered here, two factors are likely to conflict with the goal of compu-
tational accuracy of the results. First, for tidal stages corresponding to LLW, submergence
of the shoreward bar is slight and waves are likely to be breaking over the bar crest. This
effect could possibly be incorporated in the model by adapting a breaking wave scheme
from an empirical viewpoint, but strong nonlinearities would not be properly accounted
for. Secondly, for the lower water levels tested, the bar height/water depth is as large as
0.7, and the perturbation scheme adopted here is invalidated. Previous results obtained by
Dalrymple and Kirby (1986), using the boundary integral method for linear wave theory,
have indicated that the present method starts to overpredict reflection as the limits of it's
validity are exceeded. It may be necessary to develop the boundary integral method for
the cases studied here in order to obtain estimates of reflection for configurations having
shallow submergence.
CD
CD
CD
OD
--t
CO
CO
CD
C0
CD
OD
CD
O
CD1
0c
-- - - - -
-I.~ __ _
0.00 5.00 10.00 15.00 20.00 25.00
NRVE PERIOD (SEC)
Figure 16. a) 40 40 40m spacing.
Q
GD
S-- ------
jj., -- -- -- -- -- -- - -
/ \
00
CO
OD
CD
CD
WRVE PERIODISEC)
ZJ1
beach. 00Four bars, LLW+00 .61m.0.00 15.00 20.00 25.00
b) 30 40 50m spacing.
GDi
GDi
beach. Four bars, LLW+1.61m.
5 Program Listing
The program corresponding to the theory and numerical algorithm of sections 2 and 3 is
given here. The program is written using FORTRAN 77 and will compile under VAX VMS
FORTRAN or MICRO-SOFT FORTRAN with no alterations. The program is designed to
be run interactively, and prompts the user for all input data values in a self-explanatory
way.
The program is organized as a main program which calls a subroutine DEPTH to estab-
lish the one-dimensional topography h(x), and subroutine DELDET to establish the per-
turbation 6(x). A separate version of DEPTH and DELDET is supplied for each particular
topographic form, and the appropriate version thus needs to be linked to the main program
during compilation. Additional subroutines called by the main program are WVNUM,
which computes the wavenumber k at each grid point, and CTRIDA, which performs the
double-sweep solution for the complex-valued tridiagonal system of equations. The main
program and WVNUM and CTRIDA are included in the file PROFREF.FOR. Different
versions of DEPTH and DELDET are supplied in separate numbered files. Listings of all
programs are given below.
The program creates two data files. The file DEPTH.DAT contains the x-coordinate
and total depth h'(x) for each grid point. This file may be read by the statements
DO I=1,N
READ(unit,*) X(I),D(I)
END DO
where N is the number of grid points as specified in the input to the program.
The second data file REFLEC.DAT contains values of wavenumber k1, wave period T,
wave angle 0, transmission coefficient ITI, reflection coefficient IRI and the energy conserva-
tion test value for each computed reflection case. The values of transmission and reflection
are based on assuming a white incident wave spectrum, and a true reflected spectrum may
thus be computed by multiplying the incident spectrum at the offshore end of the bar field
by the values of IRI(T, 0). The data file REFLEC.DAT may be read by the statements
DO I=1,NP
DO J=1,NDIR
READ(unit,*) K,PERIOD,ANGLE,T,R,TEST
END DO
where NP is the number of wave periods and NDIR is the number of wave directions,
as specified in program input.
Main program.
C*--------------------------------------
C* REFLEC1.FOR
C*
C* PROGRAM COMPUTES THE REFLECTION COEFFICIENT FOR WAVES OVER AN
C* ARBITRARY BOTTOM WHICH VARIES ONLY IN THE X-DIRECTION. THE
C* APPROXIMATION USED IS GIVEN BY KIRBY(1986), JFM, AND ASSUMES
C* A SLOWLY VARYING TOPOGRAPHY WITH SUPERIMPOSED, FAST SMALL-
C* AMPLITUDE VARIATIONS.
C*
C* THE PROGRAM SOLVES THE FULL ELLIPTIC PROBLEM USING CENTERED
C* FINITE-DIFFERENCES AND A TRIDIAGONAL ELIMINATION SCHEME. ONE
C* DIRECTIONAL AND FREQUENCY COMPONENT IS HANDLED AT A TIME. THE
C* PROVISION FOR BUILDING UP SPECTRA BY MEANS OF MULTIPLE RUNS IS
C* INCLUDED. WAVES ARE ASSUMED TO BE UNBROKEN, NONLINEARITY IS
C* NEGLECTED, AND BOTTOM DAMPING IS IGNORED.
C*
C* JAMES KIRBY, JANUARY 1987
C*-------------------------------------------------------------------
PARAMETER(IP=1000)
COMMON/BLOCK1/H(IP),DEL(IP),K(IP),GAM(IP),P(IP)
COMMON/BLOCK2/A(IP),B(IP),C(IP),D(IP),V(IP),PHI(IP),ONE
REAL L,K,M
DIMENSION ANGLE(20)
COMPLEX ALF1,ALFN,A,B,C,D,PHI,ONE,V,L1,L2
OPEN(8,FILE='REFLEC.DAT',STATUS='NEW')
C*-------------------------------------------------------------------
C* ENTRY OF MILDLY SLOPING BOTTOM
C*-------------------------------------------------------------------
CALL DEPTH(L,N,DX)
C*------------------------------------------------------------------
C* ENTRY OF BOTTOM PERTURBATION INFORMATION
C*------------------------------------------------------.............
CALL DELDET(L,N,DX)
C*--------------------------------------...............--------------
C* DEFINE NONVARIABLE CONSTANTS
C*-------------------------------------------------------------------
ONE=CMPLX(1.,0.)
G=9.80621
PI=3.1415927
DO 1 I=2,N
D(I)=CMPLX(0.,0.)
1 CONTINUE
C*-------------------------------------------------------------------
C* SPECIFY WHETHER MULTIPLE RUNS ARE TO BE DONE
C*.-------------------------------------------------------------------
WRITE(*,90)
90 FORMAT(' ENTER NUMBER OF PERIODS, TMIN,TMAX')
READ(*,*) NP,TMIN,TMAX
DT=(TMAX-TMIN)/(NP-1)
PERIOD=TMIN
WRITE(*,91)
91 FORMAT(' ENTER NUMBER OF DIRECTIONAL COMPONENTS')
READ(*,*) NDIR
DO 92 IDIR=1,NDIR
WRITE(*,10)IDIR
10 FORMAT(' INPUT ANGLE FOR ', 12,'TH COMPONENT')
READ(*,*)ANGLE(IDIR)
92 CONTINUE
DO 300 INP=1,NP
OMEG=2.*PI/PERIOD
C*-------------------------------------------------------------------
C* SET UP PARAMETERS FOR A GIVEN FREQUENCY RUN
C*------------------------------------------------------------------
DO 2 I=1,N
CALL WVNUM(H(I),K(I),OMEG)
TKH=2.*K(I)*H(I)
P(I)=OMEG*OMEG*(1.+TKH/SINH(TKH))/(2.*K(I)*K(I))
GAM(I)=G/(COSH(K(I)*H(I))**2.)
2 CONTINUE
DO 3 I=2,N-1
C(I)=P(I+1)+P(I)-GAM(I)*(DEL(I+1)+DEL(I))
A(I)=P(I)+P(I-1)-GAM(I)*(DEL(I)+DEL(I-1))
3 CONTINUE
C*-----------------------------------------------------------------
C* SET UP PARAMETERS FOR A GIVEN DIRECTIONAL RUN
C*-------------------------------------------------------------------
DO 290 IDIR=1,NDIR
THET=PI*ANGLE(IDIR)/180.
M=K(1)*SIN(THET)
L1=CSQRT(CMPLX(K(1)*K(1)-M*M,0.))
L2=CSQRT(CMPLX(K(N)*K(N)-M*M,0.))
ALF1=L1*DX*CMPLX(0.,1.)/2.
ALFN=L2*DX*CMPLX(O.,1.)/2.
B(1)=-ONE+ALF1
C(1)=ONE+ALF1
D(1)=2.*ALF1*(CEXP(2.*ALF1)+ONE)
A(N)=-ONE-ALFN
B(N)=ONE-ALFN
DO 4 I=2,N-1
B(I)=-(A(I)+C(I))+CMPLX(2.*DX*DX*((K(I)*K(I)-M*M)*P(I)+M*M*
1GAM(I)*DEL(I)),0.)
4 CONTINUE
C+*-------------------------------------------------------------------
C* OBTAIN WAVE FIELD
*--------------------------------------------------------
CALL CTRIDA(1,N)
DO 5 I=1,N
PHI(I)=V(I)
5 CONTINUE
C+-------------------------------------------------------------------
C* REFLECTION AND TRANSMISSION COEFFICIENTS
C*-------------------------------------------------------
T1=CABS(PHI(N))
T2=CABS(PHI(N-1))
R1=CABS(PHI(1)-ONE)
R2=CABS(PHI(2)-CEXP(2.*ALF1))
WRITE(*,*)R1,R2,R1-R2
WRITE(*,*)T1,T2,T1-T2
T=(T1+T2)/2.
R=(Ri+R2)/2.
TEST=((K(N)*P(N))/(K(1)*P(1)))*T*T+R*R
WRITE(*,*)T,R,TEST
C*-------------------------------------------------------------------
C* STORE RESULTS AND MOVE ON TO NEXT COMPONENT
C*-------------------------------------------------------------------
WRITE(8,*) K(1),PERIOD,ANGLE(IDIR),T,R,TEST
290 CONTINUE
PERIOD=PERIOD+DT
300 CONTINUE
CLOSE(8)
STOP
END
Subroutine WVNUM
C*------------------------------------------------------------------
C*
SUBROUTINE WVNUM(D.K.S)
Ca
C* CALCULATE WAVENUMBER K ACCORDING TO THE FORM
C*
C* S*S-G*K*TANH(K*D)=0.
C*
C* WHERE
C*
C* D = LOCAL WATER DEPTH
C* S = ABSOLUTE FREQUENCY
C* G = GRAVITATIONAL ACCELERATION CONSTANT
C*
C* SOLUTION BY NEWTON-RAPHSON ITERATION USING ECKART'S
C* APPROXIMATION AS A SEED VALUE.
C* ----------------------------------------------------------------
REAL K,KN
G=9.80621
PI=3.1415927
K=S*S/(G*SQRT(TANH(S*S*D/G)))
DO 1 II=1,20
F=S*S-G*K*TANH(K*D)
FP=-G*TANH(K*D)-G*K*D/(COSH(K*D)**2.)
KN=K-F/FP
IF((ABS(KN-K)/KN).LT.EPS)GO TO 2
K=KN
1 CONTINUE
T=2.*PI/(SQRT(G*K*TANH(K*D)))
RETURN
2 K=KN
RETURN
END
Subroutine CTRIDA
C*-------------------------------------------------------------------
C*
SUBROUTINE CTRIDA(IF,L)
C*
C* TRIDIAGONAL MATRIX SOLUTION BY DOUBLE SWEEP ALGORITHM. PRESENT
C* SUBROUTINE ADOPTED FROM THE SUBROUTINE DESCRIBED IN:
C*
C* CARNAHAN, LUTHER AND WILKES, APPLIED NUMERICAL
C* METHODS, WILEY, 1969
C*
C* MODIFIED TO HANDLE COMPLEX ARRAY COEFFICIENTS AND SOLUTION
C* VALUES. INPUT AND OUTPUT ARE
C*
C* A,B,C = COEFFICIENTS OF ROW IN TRIDIAGONAL MATRIX
C* D = RIGHT HAND SIDE VECTOR OF MATRIX EQUATION
C* V = SOLUTION VECTOR
C* IF,L = BEGINNING AND END INDICES OF POSITIONS IN THE
C* DIMENSIONED RANGE OF THE COLUMN VECTOR TO BE
C* CONSIDERED.
C*--------------------------------------------------------------------
PARAMETER(IP=1000)
COMMON/BLOCK2/A(IP),B(IP),C(IP),D(IP),V(IP),PHI(IP).ONE
COMPLEX A,B.C,D,PHI,ONE,V
COMPLEX BETA(IP),GAMMA(IP)
C*--------------------------------------------------------------------
C* COMPUTE INTERMEDIATE VECTORS BETA AND GAMMA
C*--------------------------------------------------------------------
BETA(IF)=B(IF)
GAMMA(IF)=D(IF)/BETA(IF)
IFP1=IF+1
DO 1 I=IFP1,L
BETA(I)=B(I)-A(I)*C(I-1)/BETA(I-1)
GAMMA(I)=(D(I)-A(I)*GAMMA(I-1))/BETA(I)
1 CONTINUE
C*--------------------------------------------------------------------
C* COMPUTE SOLUTION VECTOR V
C*--------------------------------------------------------------------
V(L)=GAMMA(L)
LAST=L-IF
DO 2 K=1,LAST
I=L-K
V(I)=GAMMA(I)-C(I)*V(I+1)/BETA(I)
2 CONTINUE
RETURN
END
Subroutine DEPTH
a) Version DEPTH1.FOR; flat bottom
C* DEPTH1.FOR
C*
SUBROUTINE DEPTH(L,N,DX)
C*
C* DEPTH MODULE 1: SIMPLE UNIFORM DEPTH HO
PARAMETER (IP=1000)
COMMON/BLOCK1/H(IP),DEL(IP),K(IP),GAM(IP),P(IP)
REAL L,K,M
WRITE(*,1)
1 FORMAT(' ENTER UNIFORM DEPTH HO')
READ(*,*) HO
WRITE(*,3)
3 FORMAT(' ENTER DOMAIN LENGTH L, NUMBER OF POINTS N')
READ(*,*)L,N
DX=L/FLOAT(N-1)
DO 2 I=1,N
H(I)=HO
2 CONTINUE
RETURN
END
b) Version DEPTH2.FOR; sloping bottom
C+--------------------------------------
C* DEPTH2.FOR
C*
SUBROUTINE DEPTH(L,N,DX)
C*
C* DEPTH MODULE 2: SIMPLE PLANE SLOPE
PARAMETER (IP=1000)
COMMON/BLOCK1/H(IP),K(IP),GAM(IP),P(IP)
REAL L,K,M
WRITE(*, 1)
1 FORMAT(' ENTER DOMAIN LENGTH L, DEEP DEPTH H1,SHALLOW DEPTH H2')
READ(*,*)L,H1,H2
SLOPE=(H1-H2)/L
WRITE(*,2)
2 FORMAT(' ENTER NUMBER OF POINTS N')
READ(*,*) N
DX=L/FLOAT(N-1)
DO 3 I=1,N
H(I)=Hi-SLOPE*FLOAT(I-1)*DX
3 CONTINUE
H(1)=H(2)
H(N)=H(N-1)
RETURN
END
Subroutine DELDET
a) Version DELDET1.FOR; semicircular bars
C* DELDET1.FOR
C*
SUBROUTINE DELDET(L,N,DX)
C*
C* DEL MODULE 1: SEMICIRCULAR BARS
C*-------------------------------------------------------------------
PARAMETER(IP=1000)
COMMON/BLOCK1/H(IP),DEL(IP),K(IP),GAM(IP),P(IP)
REAL L,K,M
DIMENSION X(IP)
DO 1 I=1,N
DEL(I)=0.
X(I)=FLOAT(I-1)*DX
1 CONTINUE
WRITE(*,2)
2 FORMAT(' ENTER NUMBER OF BARS')
READ(*,*)NBARS
IF(NBARS.Eq.O)GO TO 7
WRITE(*,3)
3 FORMAT(' INPUT XO,RO FOR EACH BAR [ONE PAIR PER LINE]')
DO 5 J=1,NBARS
READ(*,*)XO,RO
DO 4 I=1,N
IF(ABS(X(I)-XO).LE.RO)THEN
DEL(I)=DEL(I)+SQRT(RO*RO-(X(I)-XO)**2)
ENDIF
4 CONTINUE
5 CONTINUE
C* STORE TOTAL DEPTH
7 OPEN(9,FILE='DEPTH.DAT',STATUS='NEW')
DO 6 I=1,N
D=H(I)-DEL(I)
WRITE(9,*) X(I),D
6 CONTINUE
CLOSE(9)
RETURN
END
b) Version DELDET2.FOR; sinusoidal bars
C*-------------------------------------------------------------------
C*
C* DELDET2.FOR
C*
SUBROUTINE DELDET(L,N,DX)
C*
C* DEL MODULE 2: SINUSOIDAL BARS
C*-------------------------------------------------------------------
PARAMETER(IP=1000)
COMMON/BLOCK1/H(IP),DEL(IP),K(IP),GAM(IP),P(IP)
REAL L,K,M
DIMENSION X(IP)
PI=3.1415927
DO 1 I=1,N
DEL(I)=0.
X(I)=FLOAT(I-1)*DX
1 CONTINUE
WRITE(*,2)
2 FORMAT(' ENTER NUMBER OF BARS')
READ(*,*)NBARS
IF(NBARS.EQ.O)GO TO 7
WRITE(*,3)
3 FORMAT(' INPUT XO, BAR LENGTH, AND BAR AMPLITUDE')
READ(*,*)XO,BARL,BARD
DO 4 I=1,N
IF(((X(I)-XO).GT.O.).AND.((X(I)-XO).LE.(NBARS*BARL)))THEN
DEL(I)=BARD*SIN(2.*PI*(X(I)-XO)/BARL)
ENDIF
4 CONTINUE
C*-------------------------------------------------------------------
C* STORE TOTAL DEPTH
C*-------------------------------------------------------------------
7 OPEN(9,FILE='DEPTH.DAT',STATUS='NEW')
DO 6 I=1,N
D=H(I)-DEL(I)
WRITE(9,*) X(I),D
6 CONTINUE
RETURN
END
6 References
1. Berkhoff, J.C.W., 1972, "Computation of combined refraction diffraction", Proc.
18th Int. Conf. Coastal Eng., Vancouver, 471-490.
2. Carnahan, B., Luther, H.A. and Wilkes, J.O., 1969, Applied Numerical Methods, John
Wiley
3. Dalrymple, R.A. and Kirby, J.T., "Water waves over ripples" Journal of Waterway,
Port, Coastal and Ocean Engineering 112, 309-319.
4. Davies, A.G. and Heathershaw, A.D., 1984, "Surface wave propagation over sinu-
soidally varying topography" Journal of Fluid Mechanics 144, 419-443.
5. DeVries, J., 1987, "Bragg reflection notes on field study preparation", Naval Civil
Engineering Laboratory, unpublished planning document.
6. Kirby, J.T., 1986, "A general wave equation for waves over rippled beds" Journal of
Fluid Mechanics 162, 171-186.
7. Mei, C.C., 1985, "Resonant reflection of surface water waves by periodic sand bars"
Journal of Fluid Mechanics 152, 315-335.
8. Miles, J.W., 1981 "Oblique surface-wave diffraction by a cylindrical obstacle", Dy-
namics of Atmospheres and Oceans 6, 121-123.
|
RSS
HIDE
