This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.

flows, the inertial force arises due to the acceleration of the pore fluid and such a

force can be quite large and can even be dominant when the permeability is high.
To model the unsteady porous flow, a number of models have been suggested.

The first model found in the literature containing the inertial force term was
the one by Reid and Kajiura (1957). It is the direct extension of Darcy's model:

Vp = P+ (2.6)
K, nat

The inertial term in this model is in fact the local acceleration term in Navier-Stokes
equation with the fluid velocity expressed as a specific discharge velocity q. It is

clear that this term contains only the inertia of the fluid but not the inertia induced

by the fluid-solid interaction, or the added mass effect.

The extension of Dupuit-Forchheimer model leads to Polubarinova-Kochina's
model (Scheidegger, 1960; McCorquodale, 1972). It has the form

Vp = aq+ bq2 + c (2.7)
at

8
where a, b and c are empirical coefficients.
Murray (1965), when studying the viscous damping of gravity waves over a
permeable seabed, proposed

Vp = u+f pn (2.8)
k at

where k is the permeability of porous media (defined differently from Kp) and 6

is the spatially averaged microscopic velocity of the pore fluid, related to q by
q'= ni. Comparing to Eq.(2.7), the coefficient c is equal to p in Murray's model.

McCorquodale in 1972 further modified the expression of this coefficient as

P (2.9)
n

and his model (McCorquodale, 1972) becomes

-Vp = (a + + b) 4+ (2.10)
n at

At about the same time as McCorquodale, Sollitt and Cross (1972) also extended
Eq.(2.2) to include the effects of unsteadiness. In their model, the inertia resistance

consists of two components, one is the inertia of the pore fluid, and the other one
is the inertia induced by the virtual mass effect due to the fluid-solid interactions;
the model reads

+ pC+ 1- n ai
Vp = (n- + n U) + p(1 + Ci ) (2.11)
Kp TKp n at

where u is again the spatially averaged microscopic velocity vector and Ca is the
virtual mass coefficient.

The value of the virtual mass coefficient Ca in Sollitt and Cross's model was

assumed to be zero in their actual computation because it was unknown at the time.
It has been so assumed in almost all the later applications of this model (Sulisz,

1985, etc.). When comparing to the experimental data, Sollitt and Cross found that

the correlation improved by taking nonzero values for Ca and asserted that "one

9
cannot predict the magnitude of this coefficient a priori because the virtual mass

of densely packed fractured stone is not known" (Sollitt and Cross, 1972, p1842).
They also pointed out that "Evaluation of Ca, however, may serve as a calibrating
link between theory and experiment in future studies" (same page).
Hannoura and McCorquodale (1978) made an attempt to determine the value

of Ca by laboratory experiment. They measured the instantaneous velocity and

the pressure gradient in the experiment, and calibrated the data with their semi-

theoretical porous flow model:

Vp = (a + b q)q+ p(l + Co) (2.12)

The results so obtained for Ca scattered in a range of -7.5 +5.0.

Dagan in 1979, based on the microscopic approach, arrived at a generalized

Darcy's law for nonuniform but steady porous flows (Dagan, 1979):

Vp = q- V2q (2.13)
Kp Kp
where is a constant coefficient depending only on the geometry of the media and
it was defined as
dL (2.14)
80
In applying this model to the wave-porous-seabed interaction, Liu and Dalrym-
ple (1984) added an acceleration term to the Dagan's model and it becomes

Vp = -pq+p-q V' (2.15)
K,, nat K,
It was found, by performing a dimensional analysis (Liu, 1984), that the last
term in Dagan's model is a second order quantity in comparison to the other terms.
The inertial term in this model is the same as that of Reid and Kajiura and that of

Sollitt and Cross where Ca is set to be zero.

Based on the phenomenological approach, Barends (1986) added another porous
flow model to the list:

V = ( + g (2.16)
n at K(q)

10
and K(q) was defined by Hannoura and Barends (1981, see Barends, 1986) as

/ pd (2.17)
K(q) = g (2.17)

with
Ca(1 n)
~ n5
where C is the drag coefficient related to the porous Reynolds number (R =q q

d,/v), a and 3 are two constants and d, is the relevant grain size.
With so many models, it is difficult to decide which one is most appropriate for

the porous flow problem in this study without further analysis.

2.2 Wave-Porous Seabed Interactions

The computation of gravity wave attenuation over a rigid porous bottom has
been performed by a number of investigators. The classic approaches by Putnam

(1949) and Reid and Kajiura (1957) was to assume that the Laplace equation is

satisfied in the overlying fluid medium and that the bottom layer can be treated
as a continuum following Darcy's law of permeability. In Reid and Kajiura's pa-
per, although the inertial term was included in the porous flow model, they later
neglected it and concluded that Darcy's law is adequate for sandy seabeds. Under

the assumption of low permeability, Reid and Kajiura found, for infinitely thick

wave number and the damping rate, h is the water depth and v is the kinematic

viscosity of the fluid.
Hunt (1959) examined the damping of gravity waves propagating over a perme-

able surface using Reid and Kajiura's porous flow model and retained the inertial

term to the end. The thickness of the permeable beds was infinite and water above

11

the permeable surface was assumed viscous as opposed to the inviscid approach by

the previous investigators. The tangential velocities at the surface of the permeable
bed were set to be zero and the continuity of the normal velocity and the pressure
were observed on the same boundary. The resulting dispersion equation was quite
lengthy. Under the consideration of sand beds of low permeability, the results for
kr and k, are

kr = k + 2k (2.20)
k, = ko 2koh + sinh 2koh

2ko aK, K2
ki = 2k o (- + ko (2.21)
2koh + sinh 2koh v 2(

with
a2 = gko tanh koh (2.22)

Murray (1965) investigated the same problem by the use of stream function and
the porous flow model given in Eq.(2.7). Instead of letting the tangential velocities

at the interface be zero as done by Hunt, Murray assumed that the tangential
velocities were finite and continuous. The first order solution of the dispersion
equation for the spatial damping was

2k a + 2koN
2koh + sinh 2koh

where k = ko + ki with ko being the initial wave number for the impervious bottom

as defined in Eq.(2.22) and
N= naKp
N =

For fine sand, N -+ 0 and
i/-v
k, = __k_ V a (2.24)
2kh + sinh 2krh
and the temporal damping for known k in such case is

a = -i-= kh- (2.25)
v2 = vsinh 2kh

12

Liu (1973) added an interface laminar boundary layer in his solution for the

same problem. In his paper, the water in the fluid domain was regarded as inviscid
and the porous flow was assumed to follow Darcy's law. The imaginary part of k,
after the boundary layer correction, was given by

2k, )
ki- =2k (aKp/v + k, ) (2.26)
2k,h + sinh 2k,h FT2a

with k, remaining the same as that in Eq.(2.18). It is noted that this is the same ex-

pression as that given by Hunt (1959). The correction term of Eq.(2.26) to Eq.(2.19)

is found to be insignificant unless Kp is extremely small.

Dean and Dalrymple (1984) later obtained the expressions for shallow water

waves over sand beds:

k, [1 1- h/g g 2 (2.27)
f T 4v h
(1 u2h/g) aK,
ki = ( ) (2.28)
2h V

All the solutions above are for infinitely thick seabeds. Liu (1977) solved the

problem for a stratified porous bed where each layer has a finite thickness and

different permeability. It was found that the wave damping rate was insignificant

to the permeability stratification while the wave induced pressure and its gradients

are affected significantly by stratification. The porous flow model employed in the

analysis also obeys Darcy's law only.

Liu and Dalrymple (1984) further modified the solution, for a bed of finite

thickness, by replacing Darcy's law with Dagan's unsteady porous flow model, which

partially included the effects of unsteady flows. The dispersion relationship from
the homogeneous problem-without laminar boundary layer corrections-was found

16
this method to porous breakwaters with constant elements, which assumes that

the physical quantities remain constant over each individual element. The whole
computation domain was divided into three regions, two fluid regions and one porous

region. In the porous domain, Darcy's law was used and in the two fluid regions,
two artificially defined vertical boundaries were placed at the offshore and inshore
ends. The patterns of the vertical distribution of the velocity potential at these two

boundaries was assumed undistorted by the presence of a structure. The reflected

and the transmitted potential functions outside the computational domain were

given by

,(, z) = [e(-')+ Aoe-k(-)] (z +h) (2.35)
cosh kh

et(x,z) = Boe-i*(+I,)coshk'z + h) (2.36)
cosh k'h'

where I and 1' are the distances from the left and the right vertical boundaries to

the origin, k and k' are the wave numbers in the reflection and the transmission
regions, respectively. The reflection and transmission coefficients are then

K, =1 Ao (2.37)

Kt =| Bo 1 (2.38)

The agreement of the transmission coefficient for impermeable floating struc-

tures with the experiment data was fairly good, but the comparison of Kt and Kr
for a sloped-face permeable breakwater was not as satisfactory.
Finnigan and Yamamoto (1979) further added the modulated standing wave

modes to 4, and Ot, while keeping the constant element and Darcy's law unchanged

a few empirical formulae for transmission coefficient have been established from ex-

perimental data by different authors (see Baba, 1986). For example, the equation

given by Goda (1969) for the transmission coefficient due to overtopping is

Kt = 0.5[1 sin ( F) (2.43)
2 a H,

18
with a = 2.2 and 3 = 0.0 ~ 0.8 for vertical face breakwaters, where F is the depth

of submergence of the breakwater crest and HI is the incident wave height. Another
equation for impermeable submerged breakwaters was given by Seelig (1980):

F F
K,= C(1 ) -(1- 2C)F (2.44)

where R is the wave run up given by Franzius (1965, cited in Baba, 1986)

R = HIC,(0.123 I)(c2v d+c3) (2.45)

with Cx = 1.997, C2 = 0.498 and C3 = -0.185. d is the total water depth and

0.11B
C = 0.51- 1B (2.46)

where B is the crest width and h is the height of the breakwater. Several other
empirical equations and methods are also available for designs.
The latest theoretical approach for impermeable submerged breakwaters was
given by Kobayashi et al. (1989) by using finite- amplitude shallow-water equations.

The numerical results were found to agree well with the experimental data by Seelig

(1980) for such structures.
However, for permeable structures, no parallel empirical formula could be found

in the literature besides the nomograph by Averin and Sidorchuk (1967, cited in
Baba, 1986). The research for such structures is still in the stage of physical exper-

iments.

Dick and Brebner (1968) carried out an experiment on solid and permeable
breakwaters of vertical faces. In the experiment, it was discovered that, over a

certain wave length range, the permeable breakwater was much better than the

solid one of the same dimension in terms of wave damping, and that the permeable

breakwater had a well defined minimum value for the coefficient of transmission. It

was also found that a substantial portion, 30% ~ 60%, of the wave energy of the

transmitted waves was transferred to higher frequency components. The equivalent

19
wave height for the transmitted waves in his paper was defined as

Heq = 2fT () dt (2.47)

Due to the smallness of the submergence of the breakwater crest, it is felt that

the majority of the waves in the experiment were breaking waves.

Dattatri et al. (1978) tested a number of submerged breakwaters of different

types, permeable, impermeable, rectangular and trapezoidal. One of the conclusions
was that the important parameters affecting the performance of a submerged break-
water are the crest width and the depth of submergence. The transmitted waves

were found "to be irregular though periodic". They also reported that porosity and

wave steepness did not have significant influence on the transmission coefficient,

which contradictory with the observations by Dick et al. (1968).

Seelig (1980) conducted a large number of tests on the cross sections of 17

different breakwaters for both regular and irregular waves. Most of the breakwaters

tested were rubble-mound porous structures with multi-layer designs. Beyond the

experiments for impermeable breakwaters by which Eq.(2.44) was obtained, he also
tested these permeable structures as submerged breakwaters. Since "no generalized

model was available" at the time for such breakwaters, the experimental data with

regular waves were compared to Eq.(2.44). It was found that the formula is quite

conservative in estimating the transmission coefficient for permeable submerged

breakwaters. The phenomenon of wave energy shifting to higher order harmonics

in the transmitted waves was again reported without further explanation.

Baba (1986) conducted an experiment on concrete submerged breakwaters and

compared the data with the four widely used computational methods for wave

This velocity also known as the discharge velocity, is related to the actual dis-
charge over a unit surface area, or, Q/A. These apparent properties are of final

interest in engineering. Since we are now dealing with a continuum, we have

Dq 8q aq
+ = Vq n" *tf (3.6)
Dt at at

here it has been assumed that the convective acceleration is negligible.

The various force terms are established as follows:

(A). Drag Force

The drag force consists of two terms: that due to laminar skin friction and that

due to turbulent form drag; the former is proportional to the velocity and the latter

is proportional to the velocity squared. Their functional form, when expressed in

terms of the final field variables, may be written as

FDz = [Ax(qz nu,) + B. q nu, [ (q= nu.)]dx nA dz (3.7)

The coefficients A, and B. are to be determined later.

(B). Inertia Force

The inertial force can be treated as an added mass effect and expressed in the

following form:
FI, = Cz( ni,,)dx nA dz (3.8)

The coefficient C, will be determined later.
(C). Body Force
In the pore fluid, there is no horizontal body force, and the vertical body force
is balanced by the static pressure gradient and, thus, vanishes. In the solid skeleton,
the vertical body force is the net weight. If the solid skeleton is subjected to an
unsteady vertical motion, this term should be included; otherwise, it can be ignored.
Substitution of Eqs.(3.5) through (3.8) into Eq.(3.3) results in the following
differential equation:

aP
=- A(qf nu,,) + Bze q nu, I (q, nu,,)
ax
+ (+ Cz)(4 -n i.) (3.9)

This is the basic equation of pore fluid motion. In a general case, it is coupled

with Eq.(3.2) and must be solved simultaneously. Only a special case with no

movement of the solid skeleton will be analyzed here. For such a case, Eq.(3.9)
reduces to
BP p
S= Aqz + B qz I q9 I + ( + C.)4. (3.10)
ax n
The force balance in the z-direction leads to

aP p
Arq + B, I q I q + (- + C,), (3.11)
8z n

3.2 Force Coefficients and Simplifying Assumptions

The evaluation of force coefficients involves a great deal of empiricism. However,

successful engineering application relies heavily upon successful estimation of these

coefficients. The part of the flow resistance which is linear in q clearly will lead to

24
a Darcy-type resistance law. The coefficient A is, at least, related to four founda-

tional properties, one of the fluid-the dynamic viscosity-and three of the porous

nificantly from a sphere shape. A new coefficient, Cm, can be defined such that

Sp n +C(1 n) 2
Cm, = p + Cz = p[ + C -- (3.23)
nAz nAz '

which has the same meaning as the mass coefficient in hydrodynamics.

Clearly the force coefficients in other directions can be defined in a similar

fashion. Equations.(3.10) and (3.11) can now be expressed in the following general

form:

A4+ | q 1 q+ Cm = -Vp (3.24)

Here, the carat indicates a second order tensor. This equation can be generalized

to include the motion of the solid simply by replacing q with q u,.

We now proceed to assume that the porous material is isotropic and that all

the coefficients reduce to scalar quantities and

nAz = Az = n

Recall that the convective acceleration in q has been assumed, in Eq.(3.6), to be

negligible, i.e.

Sq=
q= -
Bt

27
Substituting these conditions into Eq.(3.24), we obtain

(A + B q I)q+ Cm = -Vp (3.25)

This equation is similar to the well-known equation of motion in a permeable
medium such as given by Reid and Kajiura (1957) and others, with the exception
that the added mass effect is now formally introduced.

To seek a solution to this equation, the common approach is to ignore the

inertial term and linearize or ignore the nonlinear term. Under such conditions,

Eq.(3.25) can be reduced to the Laplace equation by virtue of mass conservation.

To retain the inertia term, the most convenient approach is to assume the motion

to be oscillatory, as is the case under wave excitation. Now we let

q= q-x, z)e-it (3.26)

where q'is a spatial variable only, and a is a generalized wave frequency, which could

Obviously both Reynolds numbers signify. the relative importance of the inertial

force to the viscous force. The origins of the inertial forces are, however, different.

In R/, the inertia is of a convective nature and the resistance arises due to change

of velocity in space (fore and aft the body) whereas in Ri, the inertia is of a local

29
nature and the resistance arises due to the rate of change of velocity at a specific
location. The ratio of R, and R! is the Strouhal number, which clearly identifies

the different origins of the two inertial forces.
In the range of common engineering applications, the magnitudes of various
coefficients can be estimated as follows:

n s 0.3 ~ 0.6

Co ~ 0(1)

ao ~ 0(103)

bo ~ 0(1)

Therefore, Eqs.(3.33) through (3.35) reduce to

t 10-'R, (3.38)

-j I ~ 10-2Rf (3.39)

LfI R- (3.40)
fi RI
When Eqs.(3.38) to (3.40) are plotted on a Rf and Ri plane as given in Fig. 3.2,
we identify seven regions where the three resistance forces are of varying degrees of
importance. There are three regions where only one resistance force dominates, that
is, the dominant force is at least one order of magnitude larger than the other two.

There is one region where the three forces are of equal importance. Then, there are

three intermediate regions where two out of three forces could be important.
Since both Reynolds numbers are flow-related parameters, accurate position of
a situation within the graph cannot be determined a priori. However, a general
guideline can be provided with the aid of the graph or with Eqs.(3.38) through

(3.40). We give here an example of practical interest. In coastal waters, we com-

monly encounter wind waves with frequency in the order of O(1)(rad/sec) and the
nondimensional bottom pressure gradient defined as V(p/y) in the order of O(10-1).

106 N Kegion

fn > 0lfi t, f.
10o I > 10f

102 ,/.n f, A

1
L Region
102 / I Region

fA > lofn A f > loft
104 f, > 10f,

10-6
10-6 10-4 10-2 1 10 104 106

Figure 3.2: Regions with different dominant resistant forces

Under such conditions, the importance of the three resistance components for bot-

tom materials of various sizes can be assessed. Table 3.1 illustrates the results.

As mentioned in the previous section, Eq.(4.15) is an equation involving both

of the complex unknown variables of either a or k, and the complex coefficient

fo. The other necessary equation can be obtained from the linearization process
(This step is not necessary if C1 = 0, as in Liu and Dalrymple's solution). The
common method of evaluating the linearized resistance coefficient, fo, is to apply
the principle of equivalent work. This principle states that the energy dissipation

within a volume of porous medium during a time period should be the same when
evaluated from the true system or from its equivalent linearized system

(ED)I = (ED)ni (4.22)

where ED is the energy dissipation in a controlled volume during one wave period

and the subscripts I and nl refer to linearized and nonlinear systems, respectively.

It can be readily shown (Appendix A) that such energy dissipation (considered
as a positive value) can be expressed in the form of a boundary integral

ED= = ds d (4.23)

where ED is a complex energy dissipation function and .7 is the complex energy flux
normal to S, with the real parts of them being the corresponding physical quantities.
Here S is the closed boundary of the computation domain. The complex energy
flux function 7, can be expressed as (Appendix A)

1 e-2it (4.24
= (unp e-2i+ UnP) (4.24)
2

36

where u, is the normal velocity at the boundary S and p* is the conjugate of pore

pressure p, both of them are complex quantities; a, is the wave frequency, areal

value, and the subscript r is used to distinguish the complex a.

Physically, there are two classes of problems: standing waves of a specified wave

number, and progressive waves of a specified wave frequency. In the former case, k

is real, a is complex, and the solution has the following form

r (x, t) = a e'' cos kze-'"'' o, < 0 (4.25)

where ai is the imaginary part of the complex a and a, is the wave frequency. In

the latter case, a is real, k is complex and r becomes

ti(x,t) = ae-kizei(k'-at) ki > 0 (4.26)

where k, and ki are, respectively, the real and imaginary parts of k with k, being
the wave number.

For standing waves, the pore pressure function is

P, = p e-io" = Dcosh k(h + h, + z) cos kxe'ie-it'' (4.27)

Since P, is periodic in x, the boundary curve S for the contour integral in

f = -- iP a = C (4.32)
R a
In Eq.(4.31), the magnitude of wo is approximately taken as

o = Dfo) sinhkh, (4.33)
afo

with

D(fo) = (4.34)
cosh kh cosh kh,(1 tanh kh tanh kh,)

where the decaying wave amplitude aec'i has been approximately replaced by u,
which is the averaged wave amplitude in the period of [0,T].
With such approximations, f for the nonlinear system is no longer a function
of x and t. The approximations made above are not expected to cause significant
errors for the final results.
By introducing Eq.(4.29) into Eq.(4.28) and carrying out the contour integration
regarding f as a constant, the energy dissipation within that one-wave-long portion
of a seabed is obtained as

kLT oD'(f) t + D(f) I
ED = sinh2kh, T(t) + D( 2(t)] (4.35)
8p h f f

where Ti(t) and T (t) are the nondimensional functions of t generated by the inte-
gration with respect to time.
Substituting the linearized and the nonlinear resistance coefficients given in
Eqs.(4.30) and (4.31) for f in Eq.(4.35) respectively, equating the energy dissipations

38
by the two systems according to Eq.(4.22) and assuming that a and k are the same
for both systems, we have

Substituting the expressions for fl, f2 and wo in the above equation, it becomes
1 + Cd I kD(f) sinh kh (4.40)
R a |fol
It is clear, from the definitions of the parameters, that all the terms in Eq.(4.40)
are complex quantities in the standing wave case.
For progressive waves, the pore pressure function P, is given by Eq.(4.10) and
the contour S for the integration is the same as that for standing waves. The energy
dissipation is
t+T
ED = f ds dt
Jt s

= -Jt (- dz+ Fn dz + Lf 7nd)dt (4.41)
St ==0 =L =-h
By substituting the expression of 7, derived from Eq.(4.10) into the above equa-
tion, the summation of the first two integrals is found to be O(I ki 1) while the last

39
one is 0(1). If we assume that I k, j< 1, which is generally true in coastal waters

(see the results for the progressive wave case), then the energy dissipation becomes
= t+T ft
ED = t 7 dx dt (4.42)
t O
Eq.(4.42) will clearly lead to the same relationship as that given by Eq.(4.40)
except that the averaged wave amplitude d in D, in this case, is the spatially
averaged value of ae-k'" over [0, L].
In the progressive wave case, we assume that the wave period is given and that
there is no time dependent damping, i.e. a is real. Therefore, the first and the last
term on the R.H.S. of Eq.(4.40) are all real numbers. It is then obvious that

(fo)i = -/ (4.43)

and
1 Cd I kD(k, fo) I(4.44)
(fo)r = + (4.44)
R ao | fo
where the subscripts i and r are for imaginary and real parts, respectively.
To actually compute fo requires specification of two fundamental quantities of
n and d, (see Eq.(4.40) and the definitions for R, Cd and P). It is often more
convenient to specify the permeability parameter R = aK,/v (Rr = a,Kp/v for
standing waves), as opposed to specifying the actual granular size, d,. Under these
circumstances, we could replace Cd with the following expression (Ward, 1964):

Cd = a- = (4.45)

with
C bo(1 n) (4.46)
bao(l n)

Now Cf is a nondimensional turbulent coefficient. To be consistent with the
suggested values of ao and bo given in Eqs.(3.14) and (3.20), Cf should be in the
range of 0.3 ~ 1.1 for a porosity of n = 0.4. Equation (4.40) now becomes:
1 + C I kD(k, fo) sinh kh I
fo = + f(4.47)
R |Ra afo

40
This is the final form of the equation used for computing fo for both cases.
The dispersion equation, Eq.(4.15), which is coupled with Eq.(4.47), is solved
iteratively. The procedures to obtain the solution for the cases of standing waves
and progressive waves are different.
a) For standing waves a = a, + iai is complex, k is real and known. Equation
(4.15) can be rewritten as

where the superscript n denotes the level of iteration, and the criterion of conver-
gence is
fn) O-) and o'(n) I
I fo(n) I_ and a ((n)

with e being a pre-specified arbitrarily small number. It is set to be 1.0% in this
model.
b) For progressive waves k = k, + iki is complex and a is real and known. In
this case, the dispersion equation can be written as

where F, and F, are the real and imaginary parts of F and n indicates the iteration
level.
The criterion of convergence for such case is
fo") fc"-1) k ,(") f("-)
fI n)0 ) f ( and I k(-) -| (4.59)

with e being a pre-specified arbitrarily small number. It is again set to be 1.0% in
this model.
4.3 Results

The predicted damping rates kip from the solution of Eqs.(4.15) and (4.47) are
first compared to the laboratory data ki, by Savage (1953) for progressive waves

propagating over a sandy seabed. The experiment was conducted in a wave tank of
29.3 meters (96 ft) long, 0.46 meters (1.5 ft) wide and 0.61 (2 ft) meters deep. The
porous seabed was composed of 0.3 meters thick of sand with the medium diameter
of 3.82 mm. The water depths are h = 0.229 m. and h = 0.152 m. The data
for the water depth of h = 0.102 m was ignored for the reasons given by Liu and
Dalrymple (1984). The wave conditions and the comparison of the damping rates
are listed in Table 4.1. The parameters used in the solution of Eqs.(4.15) and (4.47)
are: n = 0.3, ao = 570, bo = 2.0 and Ca = 0.46 and the average wave height H in
the table was calculated according to
Ho HLg
H = (1 i
HoIn Ho
HLg
where Ho and HLg are the wave heights measured at two points of Lg apart with
Lg = 18.3 m (60 ft). The values computed by Liu and Dalrymple (1984) are also

42

listed. The relative error in the table is defined as
k Iko,
A% = i m | x100% (4.60)
kim

with ki being either kip or kiLD.

In this case, the errors are of similar magnitude between the present model and

the model by Liu and Dalrymple (1984). This is because the experimental values

fall in the region where the inertial resistance due to the virtual mass effect and the

turbulent resistance, neglected in Liu and Dalrymple's model, are unimportant.

In Fig. 4.2 through Fig. 4.5, the solutions of the complex dispersion equation are

illustrated graphically. The case of progressive waves are demonstrated first. The

The experiment was carried out in a wave flume in the Laboratory of Coastal
and Oceanographic Engineering Department of University of Florida. The flume
is about 15.5 meters long, 0.6 meters wide, 0.9 meters high and equipped with a

mechanically driven piston-type wave maker. All the tests were conducted with

standing waves.
5.1 Experiment Layout and Test Conditions

Figure 5.1 shows the experiment arrangement. A porous gravel seabed was
constructed at the end of the wave flume in the opposite side of the wave maker. A

sliding gate was positioned at one wavelength from this end of the tank to trap the
standing wave after a sinusoidal wave system was established by the wave maker.

The decay of the freely oscillating standing wave was then measured by a capacitance

wave gage mounted at the center of the compartment. The damping rate of the
porous seabed for each particular test condition was determined by applying a least
squares fit to the data according to the following equation:

1 H-
(oi) = I-n( ) (5.1)
Tj Hj-i

where j refers to j-th wave. The corresponding wave frequency was obtained by

averaging individual waves.

The contribution to the wave damping due to the side walls and the bottom
was subtracted from the data according to the following equation

(36 cases) is as high as 0.327. It can be seen from the figure that considerable

discrepancies occur for the damping rate, and the maximum relative error is about
87% of the experimental values. The prediction of wave frequency by their solution

is about the same as that by this study.

A test over a sand bed (d50 = 0.16 mm) was also conducted. Active sand move-
ment was observed along with the rapidly formed ripples. As a result, theoretical

predictions are several orders of magnitude smaller than those actually measured
(Table 5.8). The basic assumption of no particle movement was violated. How-

ever, the predictions of wave frequencies were as good as the coarse-grained cases.

The same test, with dye injection into the sand bed, disclosed that the flow inside

the sand bed was very slow, and indicated that the assumption of impermeability
invoked in the linear wave theory is justified in such situations.

When the thickness of a gravel bed is reduced to one grain diameter, the model
also failed to yield good results as shown by the values in Table 5.9. In this case, the
effect of the boundary layer at the interface, which was not included in the model,
obviously cannot be ignored.

The nature of a turbulent boundary layer over a thicker seabed was qualitatively
investigated for the case of h,= 20 cm, dso = 1.48 cm, h =20 cm and L =200 cm,

by dye studies. Turbulent diffusion was spotted over virtually the entire porous

domain and up to approximately 1.0 cm above the interface. The effects of turbulent

boundary layers to the wave damping are important in wave interaction with porous

seabeds and more theoretical and quantitative experimental research is needed.

CHAPTER 6
BOUNDARY INTEGRAL ELEMENT METHOD

The Boundary Integral Element Method (BIEM) is an efficient numerical method
for conservative systems such as those governed by Laplace equation. It is efficient
because the computation is carried out only on the boundaries rather than over the
entire domain as it would be in finite difference and finite element methods. Consid-
erable amount of computation time and storage space can be saved without losing
accuracy. In addition to the efficiency, the data preparation becomes much simpler
as compared to the other two methods. For non-Laplace equations, the application
of this method may become difficult and sometimes impossible. By the efforts of
many scientists, the scheme has been extended to more and more problems governed
by non-Laplace equations (Brebbia, 1987). In this chapter, we restrict ourselves to
the two dimensional Laplace equation only.

6.1 Basic Formulation

Let U and V be any two continuous two dimensional functions, twice differen-
tiable in domain D which has a close boundary C. According to the divergence
theorem (Loss, 1950; Franklin, 1944, cited in Ligget and Liu, 1983) for continuity
in a volume

f D(V f u)dA= 6. ds (6.1)

where V is the vector operator, i is any differentiable vector and n' is the unit
outward vector normal to C. If we define that v = UVV and then v7 = VVU in
Eq.(6.1), two integral equations can be generated in terms of U and V,

To apply this equation, we choose U as the velocity potential 4 and V as a free
space Green's function, G. Both of them satisfy Lapalce equation. Equation (6.7)
can then be rewritten, in terms of 1 and G, as

G(P a)(Q)
C[(Q) -aG(PQ G(P,Q) ]ds = 0 (6.8)

where P is a point in the domain DnC and Q is a point on C.
One of the free space Green's functions for two dimensional problems is

G(P,Q) = nr (6.9)

where r is the distance between point P and point Q and it can be expressed by

r = V(zp Xq)2 + (zp zq)2

on the x z plane.
Substituting Eq.(6.9) into Eq.(6.8), it becomes

#(Q) Br 8(Q)-
[ In r n( ds = 0 (6.10)
r 5n an

66
Despite the fact that the Green's function has a singularity at r = 0, the contour
integral in Eq.(6.10) exists and can be worked out by removing the singular point
from the domain.

In general, Eq.(6.10) becomes

(P) = In r ]ds (6.11)
cz i(P) T(Q a
r an an

after removing the singularity at point P. Here a is 27r if P is an interior point
and is equal to the inner angle between the two boundary segments joining at P if
it is on the boundary. Generally, point P can be anywhere in the domain DnC,

while Q is always on C due to the contour integration. In BIEM, since only the
a
boundary values of D and are solved, P has to be kept on the boundary C in
an
the process of computation. When 4 and on C are solved, the interior values
can be derived with Eq.(6.11) by placing P at the point of interest inside D.

The closed boundary C in Eq.(6.11) is the same as that in Eq.(6.10) except that
it does not pass point P as it does in Eq.(6.10). The singular point of the Green's
function r = 0 has been excluded by a circular arc of infinitesimal radius from the
domain DnC. Thus there is no singularity on the new contour C. Eq.(6.11) is the
basic equation for BIEM formulations.

67
quadratic element, special element and so on so forth. The classification of elements
is based on the type of the function used to interpolate the values over an element.
For example, on a constant element, the physical quantities and their normal deriva-

tives are assumed to have no change over the element; while on a linear element,

the quantities and their normal derivatives are interpolated by a linear function
between the values at the ends of the element, i.e.

= -N(e)D- +Ng(C)4y+1 (6.13)
-- (C) = N~ ()( )i + N ()() )j+ (6.14)

Cj C ej+i

where N1 and N2 are linear functions with the feature of

N (C) = 1 N1(e+1) =

Ng2() = 0 N2(e,+1) = 1

In the same manner, higher order elements can be defined by replacing N1 and
N2 with higher order functions. When the variation of the quantities is known,
the variation function can then be chosen as the interpolation function to obtain

a more accurate element. Such elements are usually classified as special elements.
Generally, it is difficult to judge which element is superior over the others with-

out analyzing the particular problem. As a rule of thumb, the higher the order of
an element, the more accurate the approximation would be for the same element

size. But for a higher order element, as a trade-off for precision, the formulation

would be much more complicated and tedious, and the computation might be much

more time consuming. The constant element works fairly well for problems with
smooth boundaries and continuously changing boundary conditions, but could gen-

erate considerable errors at 'corner points' where the boundary conditions are not

continuous. For boundary with such 'corners', the linear element is usually a better

choice.

/
/
-'ii
/
/

Pi-

Figure 6.1: Auxiliary coordinate system

6.2 Local Coordinate System for A Linear Element

To formulate a problem with linear elements, an auxiliary local coordinate sys-

tem, as shown in Fig. 6.1, has to be established to facilitate line integration over

the element. The two axis C r7 are perpendicular to each other with one lying on

the element and pointing in the direction of integration, from Pj to Pjy+, and the
other one being in the same direction as the outward normal vector in. Based on this

auxiliary coordinate system, the integration over each segment can be completed

analytically and be expressed in terms of 4., 4i+I and t7i, the distances from Pj,

P,+i and Pi to the local origin defined in Fig. 6.1. The values of j, (e+i and i7, are
determined by the global coordinate information of three points.

Since

+ 17+ 2 = rI-

(i+ Ai)X +171 = ~tij+i

where

r,+, = vi x+) + (z. z+)

^ l= \{i- Za+l)2 + (z.- *'~)

then
'- Yj (6.15)
2Aj

+l = ,i + A (6.16)

The value of ri is therefore

mI= ,/'- (6.17)

However, the numerical test showed that Eq.(6.17) could lead to significant error

of I ri due to runoff errors. A more accurate expression for I ri j can be obtained by
computing the distance from Pi to the straight line PyP,+i directly from the global
coordinates, i.e.
I A(x- x) + zi (6.18)

where
A = jz+1 z)
zy+l zy
The sign of ri depends on the relative position of Pi to the boundary line element

PjPj+1. If Pi is in the same side of the element with ii, r7i is positive, otherwise, it
is negative. Mathematically, it can be expressed by the following two equations for
vertical and nonvertical elements:

sign (7) = (S zi) Az =x (6.19)
(xj x,) Az

sign (ri7) = A Azi X,) A x+ (6.20)
Ax Azio

where

Ax = xY+ Xi

Az = zy+ z
Az
Azio = (X- r) + z z

Therefore

(6.21)

rni = sign (r1i) I 7, I

70
6.3 Linear Element and Related Integrations

By the definition stated in the last section, on a linear element, a quantity and
its normal derivative are interpolated by a linear function between the values at
the two nodes, P, and P,+l. For the velocity potential function 4 and its normal
derivative we have, in terms of q r] coordinates,
an
(+) = ( )+ ( f-iCi+) ,i i+1 (6.22)
Cj+1-

a = ([)i] + [C.a+l)+ e
n a a a (6.23)
j e
It is obvious that
a(D) = ; (-)(a) = )(-

) = +1; ) )= (an)

Substituting Eqs.(6.22) and (6.23) into Eq.(6.12), and noting that

r, = 7, + (6.24)
ri a= r (6.25)
an a 7i r,
the integration over segment j (between P, and Pj+1) can be carried out analytically
in the local coordinate system (Ligget and Liu, 1983):
f fi+ 9ri a j +r
l (r In r, )dan

= Il 4+ H +1,-K i( ), +- K (i+ (6.26)

where the superscripts and + refer to the positions immediately before and after
the nodal point denoted by the corresponding subscript. And

(6.27)

H -i = -ia1 + j+iy12j
cj 'ij ti

71
KI. .= I 1 2lI

Kl = -,j + 1j+4

K, = I2j- lSi2

1 1f i+1 =r
i+1 fi fi ri Kn

- ij7 +1I
2AE 171?+

= ~1
j1

fJ+

1ar dEr = r fi+l
ri 8n 2Ai

1 [t an-'( 1
A ji

(6.32)

- tan-'( )]
ri

= 1 fi+nrid.
fi+1 j 1i

= -{(, + J+1)[ln(o77 + -+1)- 1]

-(7,? + )[1ln(,? + ?) 1]}

1=In ri d(
*+ n( + -+

1 2
= {+I i(2 + 77 ) ln(7? + 2(J) 2+ )

+2rjj[tan-'l(e'+l) -tan--l()]}
th 77-

with

(6.28)

(6.29)

(6.30)

rii
2A,

I CJ+ i

;dC
r.

(6.31)

d xi
rT

(6.33)

(6.34)

72

for all i j and i j + 1 and

Ai = Ei+1 ti

If i = j, the above integrals have to be re-evaluated to avoid the singularity of

the Green's function,

1
I j = lim

1 -
Ij 1= lim
j+1 1i -
:i+: C- 'o

1
-= lim
ji+1 i e-o

Si+i1 =
Ji+1 1 ar 0

ei+, ri an

i+, In ri d
fi+
*(,+<

= (2nA 1)
4

= 1+ im t In r, de
$y~i-,

= lnA- 1

Similarly, when i = j + 1,

1 /f++i- 1 ari
= lim 5l dC = 0
S+1i C-0tI ri 8n

1 ,i+l-* 1 ar,
-= li~+ -od = 0;
j T (e--#40 ri an

= lim Inri d4
Tj+1 -i T ei

= (1-2 In A,)
4

(6.35)

(6.36)

(6.37)

I1+1.1

j13+11J

Ij111

(6.38)

(6.39)

(6.40)

(6.41)

= nA 1

n ri d

(6.42)

Summarizing all the Ij's according to Eq.(6.12) for node points i = 1 N, and

assuming that

(6.43)

(6.44)

a0
at all the nodes, a system of linear equations with unknowns of P and n can be

obtained:

N
Si= E Qin
j=1

P4,1 = H~,_N + Hli 6,1aj
?j., = Hijl + H ,i
^**2 HLi b ijci

1 ifi=j
.-- 0 if i # j

Sij = K, N + Kj,1

LiJ = Ki_ i- + Kaa

j = 2, 3,..., N

Equation (6.45) can also be expressed by a matrix equation

where

(6.45)

with

(6.46)

(6.47)

(6.48)

(an

74

where 4 and 4, are the vectors of order N x 1 containing the unknown T and L,

respectively. The coefficients in R and L are determined solely by the boundary
geometries.
6.4 Boundary Conditions

In Eq.(6.45), there are N linear equations and 2N unknowns. The additional

N equations necessary for obtaining a unique solution are usually introduced by

the application of the boundary conditions. There are essentially three types of

boundary conditions for boundary value problems.

I J = (6.49)

II = (6.50)
an an
III a- k= (6.51)
on

where D, f and k are known functions on the boundary.
In
In water wave problems, the first relation usually means that the pressure distri-

bution is known, and the second describes a given flux distribution on the boundary.
They are also referred as 'essential' and 'natural' conditions (Brebbia,1989), respec-
tively. The third one some times takes place when none of the quantities in the first

two equations are specified but only the relationship between them can be found.
A typical example of this is the free surface boundary condition in linear wave the-

ory. For nonlinear wave problems, the right hand side of the above equations may
include other known functions.

As might be noticed in Eq.(6.45) or Eq.(6.48) there are two unknowns for each

nodal point. Since we have obtained one equation for every node on the boundary,

only one of the three relations needs to be specified at each node, if the assumptions
in Eqs.(6.43) and (6.44) are strictly satisfied. However, for some nodal points where

the boundary changes directions, ('-)~ may not be equal to ( )+. For such cases,

more than one boundary conditions may be necessary for one node. The treatment

75

of such points will be discussed in the next chapter.

By introducing the boundary conditions at all N boundary nodal points, another

system of N equations with the same unknown in Eq.(6.48) is established. The

number of unknowns is now equal to the number of the equations. It is usually

convenient to eliminate N unknowns with the 2N equations, and the resulted matrix

equation can be expressed as
AX = b (6.52)

in which A is a known matrix of order N x N, X is the unknown vector containing

4 or or D for some part of the boundary and for the rest of it, depending
an an
on which one is not specified by the boundary condition; and b is a known vector

resulted from boundary conditions.

CHAPTER 7
NUMERICAL MODEL FOR SUBMERGED POROUS BREAKWATERS

In this chapter, a numerical model of BIEM for submerged porous breakwaters

is developed by using the unsteady porous flow model given in Chapter 3. The

basic function of the numerical model is to compute the wave flow field and related

quantities such as the wave form, the dynamic pressure and the normal velocity

along all the boundaries and so on. With these quantities, the wave transmission

and reflection coefficients and wave forces can then be obtained.
7.1 Governing Equations

The computation domain of the problem, as shown in Fig. 7.1, consists of two

sub-domains, the fluid domain and the domain of the porous medium. In the fluid

domain, the water is considered inviscid and incompressible. The flow induced by

gravity waves is assumed irrotational. Thus, the governing equation in this domain,

for the velocity function 4, is the Laplace equation,

V2 = 0 (7.1)

with fluid velocities being defined as

S= (7.2)
ax

w = (7.3)
az

While in the porous medium domain, the viscosity of the fluid cannot be ig-

nored since the flow is largely within the low Reynolds number region. The flow is

described by the linearized porous flow model,

pafof= -VP (7.4)

where P is the pore pressure, fo is the linearized resistance coefficient and q*is the

There are two vertical lateral boundaries, one is on the offshore side and one

is on the lee side of the porous domain, as shown in Fig. 7.1. In this study, the

radiation boundary condition for these two boundaries is adopted. Such a boundary

80
condition assumes that the waves at these two boundaries are purely progressive,
and that the decaying standing waves generated by the object inside the domain are
negligible at these boundaries. Comparing with the matching boundary condition
with the wave maker theory employed by Sulisz (1985) and some other authors,

the radiation boundary condition offers significant simplicity in programming and
provides sufficient accuracy with much less CPU time. In applying this boundary
condition, the two lateral boundaries have to be placed far enough from the struc-
ture. Numerical tests show that a distance of about two wave lengths from the toe

of the structure is more than adequate for this purpose.

At these two lateral boundaries, the potential functions for the transmitted and
reflected waves are assumed to have the forms

r (xz) = e'k(z+')R(z) (7.18)

t(x, z) = e-i'(L-')T(z) (7.19)

where T(z) and R(z) are two unknown functions of z; k and k' are the wave numbers
at the two boundaries respectively, and the subscripts t and r here refer to the
transmitted and reflected waves, respectively. On the lateral boundary of up-wave

side, the potential function for the incident wave is known:

,;(z, z) = e-'k(+')A(z) (7.20)

with
Hg cosh k(z + h) (7.21)
2a cosh kh
Therefore,

S= + 4' at x= -1 (7.22)

S= 4t at x= (7.23)

where I and 1' are, respectively, the distances from the origin of the coordinates to

the offshore and onshore lateral boundaries and 0 is the unknown potential function.

where h' is the water depth at x = 1'.
While on the lateral boundary of reflection (offshore side), x = -1,

r, = 4 4I (7.26)

a a a, a,._ ik, ik,
an az ax ax

= ik ik(4O 4i) = 2ik, ik4'

or
a 2ikd; ikO (7.27)
an
in which k is the incident wave number determined by

gk tanh kh = a2 (7.28)

where h is the water depth at x = -1.

7.2.2 Boundary Conditions for The Porous Medium Domain

In this domain, there are two types of boundaries, one is the common boundary

with the fluid domain and the other one is the impervious bottom. On the interface,

the boundary condition is the same as the one specified by Eq.(7.16) and (7.17).
On the impermeable bottom, the no-flux condition, in terms of the pore pres-

sure, is
ap 0 (7.29)
an
Again, this boundary can have any shape, it can be the surface of the domain

for core materials if it is considered as impermeable, etc.

82
7.3 BIEM Formulations

The formulations in the two domains are slightly different because of the dif-
ferent boundary conditions. In the fluid domain, due to the radiation boundary
condition on the lateral boundary of offshore side, the terms containing the incident
wave potential are introduced, which will form the RHS vector of the matrix equa-
tion. On the free surface, the normal derivative are expressed in terms of 4 owing
to the CFSBC. In the porous sub-domain, the matrix equation, after applying the
no-flux condition, has to be manipulated to match with the fluid domain.
7.3.1 Fluid Domain

According to the formulae given in chapter 6, Eq.(6.12) can be expanded, in
terms of node values, assuming that 4- = 0+ for all the nodes, as:

where a, is the inner angle of node i and the subscripts 1, Ni, N2, N3, Ne and N,,
refer to the nodes of 'corner points' as shown in Fig. 7.1. The superscripts and +
refer to the positions immediately before and after the nodal point in the direction

of contour integration, and the subscript n refers to the normal derivative. Here 0-
is not necessarily equal to 0+ for all the nodes.