UFL/COEL95/013
NUMERICAL MODELING OF WAVES IN THE NEARSHORE ZONE WITH PERMEABLE STRUCTURES by
Santiago Alfageme Thesis
1995
NUMERICAL MODELING OF WAVES IN THE NEARSHORE ZONE
WITH PERMEABLE STRUCTURES
By
SANTIAGO ALFAGEME
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF ENGINEERING
UNIVERSITY OF FLORIDA
1995
ACKNOWLEDGMENTS
I would like to express my sincere appreciation and gratitude to Dr. Hsiang Wang, my advisor, who gave me the opportunity to come to the University of Florida and enroll in the master's program. He, along with the rest of the professors in this department, specially the other two members of my committee Dr. Robert G. Dean and Dr. Robert J. Thieke, has guided me through numerous research and class work problems. Appreciation is also extended to Dr. Masahiko Isobe for his comments and help during the past few months.
The author would also like to thank the other members of the research group he has worked for during the past two years, Dr. LiHwa Lin, Dr. Taerim Kim, and "future Dr." Xu Wang, their friendship, help and suggestions made those two years infinitely easier.
Gratitude is extended to the staff of the Coastal and Oceanographic Engineering Laboratory at the University of Florida and the staff of the Department of Coastal and Oceanographic Engineering who helped the author in many ways. Special acknowledgment is given to Becky Hudson, Helen T. Twedell and John M. Davis for their patience with me and the rest of the students.
Thanks also go to my fellow students Mike Krecic, Paul Devine, Rob Sloop, Jie Zheng, Yingong Li, and all the rest; their company and friendship have made my stay in Gainesville an memorable experience.
Finally I want to thank the most important people in my life, my parents and two siblings, who have supported my goal of getting this degree even though they did not like the idea of me being so far away from home. And to the person that has made that distance bearable, my "novia", Leah Gambal, my deepest love and gratitude.
TABLE OF CONTENTS
ACKNOW LEDGM ENTS .............................................. ii
LIST OF FIGU RES .................................................... v
A B STR A CT ......................................................... viii
CHAPTERS
I INTRODU CTION ................................................. 1
1. 1 Literature Review ............................................. 3
1.2 Sum m ary of contents .......................................... 6
2 WAVE ENERGY PROPAGATION THROUGH VARYING CURRENT FIELD 8
2.1 Governing Equations and Boundary Conditions 2.2 Wave Energy Equation ..................
2.2.1 TimeAveraged Wave Energy .......
2.2.2 TimeAveraged Wave Energy Flux .
2.2.3 TimeAveraged Wave Energy Equation 2.3 W ave Action Equation ..................
......................
. . . . I . . . . . .
.....................
.....................
......................
......................
3 MATHEMATICAL WAVE MODEL ................................. 20
3.1 M ild Slope Equation ........................................... 21
3.2 Derivation of Governing Equation ................................ 25
4 NUMERICAL SCHEME AND TESTING OF WAVE MODEL ............. 28
4.1 Numerical Scheme of Wave Model ............................... 28
4.2 Stability A nalysis ............................................. 32
4.3 Boundary Conditions .......................................... 34
4.3.1 Lateral Upwave and Downwave Boundary Conditions .......... 34
4.3.2 Boundary Conditions for Groins ............................ 37
4.4 Description of the W ave M odel .................................. 42
4.5 Testing of the W ave M odel ...................................... 45
4.5.1 W ave Shoaling and Refraction ............................. 45
4.5.2 Wave Diffraction and Reflection ........................... 49
4.5.3 W aveCurrent Interaction ................................. 59
5 SURF ZONE M ODEL .............................................. 62
5.1 TimeAveraged Wave Energy Equation in the Surf Zone .............. 62
5.2 Wave Action Equation in the Surf Zone ............................ 65
5.3 Wave Height Transformation in the Surf Zone ...................... 66
6 M ODEL APPLICATIONS .......................................... 72
6.1 W ave Field over a Paraboloidal Shoal ............................. 72
6.2 Wave Field around a Detached Breakwater ......................... 76
6.3 W ave Field around Groins and Jetties ............................. 81
6.4 W ave Field around Permeable Groins ............................. 83
6.4.1 Single Groin in Constant Water Depth ....................... 83
6.4.2 Single Groin on a Plane Beach ............................. 88
6.4.3 W aves around a Field of Groins ............................ 90
6.5 Preliminary Results from a Circulation Model ....................... 91
6.5.1 Circulation Results for a Plane Beach ......................... 95
6.5.2 Circulation Results for a Groin Field ......................... 98
7 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY .. 101
7.1 Summary and Conclusions ..................................... 101
7.2 Recommendations for Further Study ............................. 103
7.2.1 W ave Propagation M odel ................................ 103
7.2.2 Boundary Conditions ................................... 104
7.2.3 Breaking Criteria ....................................... 105
7.2.4 Radiation Stresses ...................................... 105
APPENDIX
A REFLECTION AND TRANSMISSION FROM POROUS STRUCTURES
UNDER OBLIQUE WAVE ATTACK ................................ 107
BIBLIO GRAPH Y .................................................. 113
BIOGRAPHICAL SKETCH ........................................... 117
LIST OF FIGURES
4.1 SubGrid system for wave model ................................. 30
4.2 Location of grid points where the wave potential
will be solved using a boundary condition .......................... 35
4.3 Description of wave field around a groin ........................... 37
4.4 Schem atic diagram .............................................. 41
4.5 Computational domain of the model ............................... 43
4.6 Comparison of wave shoaling .................................... 47
4.7 Comparison of wave refraction, t9=20 .. ............................. 47
4.8 Comparison of wave shoaling/refraction, Oi = 20 .. ..................... 48
4.9 Comparison of wave refraction, Oi = 30 .. ............................. 48
4.10 Coordinate system for analytical diffraction/reflection solution ........... 50
4.11 Grid and structure layout used for first diffraction test.
The values of x and y refer to model coordinates ...................... 51
4.12 Comparison of wave diffraction/reflection coefficient.
Top 0= 5', bottom 0= 15 .. ....................................... 52
4.13 Comparison of wave diffraction/reflection coefficient
along five different sections, 0= 5 .. ................................. 53
4.14 Comparison of wave diffraction/reflection coefficient
along five different sections, 0= 15 .. ................................ 54
4.15 Comparison of wave diffraction/reflection coefficient
along five different sections, 0= 30 .. ................................ 56
4.16 Comparison of wave diffraction coefficient, 0= 90 .. ................... 57
4.17 Comparison of wave diffraction coefficient, O= 700......................... 58
4.18 Conditions of collinear and shearing current .......................... 59
4.19 Comparison of collinear wavecurrent interaction ...................... 60
4.20 Comparison of shearing wavecurrent interaction,
amplification of wave height ..................................... 60
4.21 Comparison of shearing wavecurrent interaction,
change in wave angle .......................................... 61
5.1 Comparison between breaking analytical and numerical results.
Test 1: T= 3 s, Test 2: T= 5 s, Test 3: T= 7 s......................... 71
6.1 Bathymetry for paraboloidal shoal configuration,
after Ito and Tanimoto (1972).................................... 73
6.2 Wave height comparisons between computed and
experimental results by Ito and Tanimoto (1972) ....................... 74
6.3 Numerical model results for shoal configuration
used by Ito and Tanimoto (1972).................................. 75
6.4 Numerical model results for detached breakwater test ................... 78
6.5 Comparison of location of breaker line behind a detached breakwater,
experimental data after Watanabe and Maruyama (1986) ................ 79
6.6 Comparison of distributions of wave height around a detached
breakwater, experimental data after Watanabe and Maruyama (1986) ........79 6.7 Numerical model results for jetty test............................... 81
6.8 Comparison of location of breaker line around a jetty,
experimental data after Watanabe and Maruyama (1986) ................ 82
6.9 Comparison of crossshore distributions of wave height around a jetty,
experimental data after Watanabe and Maruyama (1986) ................ 82
6.10 Free surface and amplification factor around one groin, Case 1 ............ 85
6.11 Free surface and amplification factor around one groin, Case 2 ............ 85
6.12 Free surface and amplification factor around one groin, Case 3 ............ 86
6.13 Free surface and amplification factor around one groin, Case 4 ............ 86
6.14 Amplification factor along four sections parallel to groin ................ 87
6.15 Comparison of model results for the cases of fixed and variable
reflection/transmission coefficients................................ 89
6.16 Offshore incident wave angle influence on
reflection and transmission coefficients............................. 89
6.17 Model results for groin field #1................................... 92
6.18 Model results for groin field #2................................... 92
6.19 Wave height along the upwave side of the groins ....................... 93
6.20 Wave height along the downwave side the groins ...................... 93
6.21 Circulation in a plane beach: Case 1 ............................... 96
6.22 Circulation in a plane beach: Case 2................................ 97
6.23 Baltic coast bathymetry and groin location........................... 98
6.24 WaveCurrent field ........................................... 100
6.25 Mean free surface and Currents.................................. 100
A. 1 Schematic diagram of porous structure............................. 107
Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Engineering
NUMERICAL MODELING OF WAVES IN THE NEARSHORE ZONE WITH PERMEABLE STRUCTURES
By
Santiago Alfageme
August 1995
Chairman: Dr. Hsiang Wang
Major Department: Coastal and Oceanographic Engineering
This study presents a numerical model that predicts the transformation of monochromatic waves over complex bathymetry and includes refractive, diffractive and reflective effects induced by structures placed along the wave main direction of propagation. Particular emphasis is placed upon permeable structures which are common in the nearshore zone. Since the model has been developed with the idea of describing wave propagation in the surf zone, and in anticipation of the possibility of combining it with a circulation model, wavecurrent interaction effects are included in the derivation of the governing equations, namely the 'wave energy conservation equation', 'wave action equation' and 'kinematic conservation of intrinsic angular frequency'. Finite difference approximations are used to solve the governing equations, and the solution is obtained for a finite number of rectilinear grid cells that comprise the domain of interest. The solution technique employed, based on a marching Gragg' s method, is direct, and eliminates the computer storage problems
associated with large matrices inherent in other methods. This efficiency is obtained at the cost of resolution of the wave components normal to the incident wave direction.
Boundary conditions for all sides of the numerical domain are derived, as well as the necessary conditions to reach a solution on grid points next to permeable structures. The effects induced by different permeability values of the structure on the wave field are well simulated through the use of the appropriate boundary conditions that include calculated reflection and transmission coefficients, which in turn are based on the incident wave field and the physical characteristics of the structure.
A surf zone model based on the work of Lee and Wang (1993) is derived. The model takes a different approach to surf zone dynamics by showing that wave action is conserved in the surf zone. This, together with the wave energy equation, allows for the derivation of a surf zone model that requires fewer empirical coefficients and is better suited for a general three dimensional bathymetry.
Model results are compared first with data from analytical solutions. Then, the capability and utility of the model for real coastal applications are illustrated by application to three different experimental test layouts, where ineasured data were available. Finally, to study how the model performs when solving problems that include permeable structures, several cases are presented, from a single groin in constant water depth, to a shoreline protected by a groin field. The preliminary results from a circulation model coupled with the wave model are shown as well.
CHAPTER 1
INTRODUCTION
Coastal engineers and scientists face an increasing demand for more accurate solutions to a wide range of coastal problems, many of them involving the transformation of wind generated waves as they propagate from deep water, over complex bathymetry, into shallow water where they may be further affected by man made structures. The nearshore wave field will generally need quantitative evaluation in order to design those coastal and harbor structures, to study and comprehend coastal processes, and to predict nearshore currents and the resulting sediment motion. Thus, there is a need for an efficient and accurate method to estimate nearshore wave conditions with arbitrary shoreline shape and nearshore bathymetry, as well as the presence of coastal structures.
Several mathematical models have been proposed over the years, as well as their corresponding numerical solutions. Most of them are capable of producing reliable results under certain conditions, but with a series of drawbacks as well. That is, most of the models are restricted to a certain range of imposed conditions.
The main objective of this study was to develop a model that would effectively simulate the nearshore wave field in the presence of structures aligned approximately in the main wave direction such as a field of groins with different permeability values. Over the years, different kinds of groins have been used for coastal protection. There are a great variety of shapes and materials used for groin construction; traditionally these structures consist of rubble mound. The permeability of this kind of groin could be very different
2
depending of the size of the stone used, the dimensions of the groin, the use of different layers with different stone sizes across the groin, etc. Completely impermeable groins have several disadvantages, including the following:
* They are very expensive.
0 They tend to induce deep erosion pits seaward of the groin head because of locally
increased currents.
* They may stimulate ripcurrents and seaward loss of sand.
* Extensive leeside erosion may occur on the adjacent unprotected shoreline.
Other types of groins include the ones composed of metal sheets that are driven into the beach, or even made out of precast concrete units. An interesting solution often used in places like the Netherlands and Germany consists of single or double permeable rows of wooden piles perpendicular to the beach. This system costs only 10 to 25% of the cost of rubble stone groins used in those countries. In contrast to the impermeable groins, which seek to form a complete obstruction to the longshore current and sediment transport, highly permeable groins like pile screens are meant as an artificial hydraulic resistance in order to gradually reduce the longshore rate of sediment transport. The advantages of this kind of groin can be summarized as follows:
* Low cost structures.
* Less pronounced current concentration seaward of the head.
0 Reduced tendency to form ripcurrents.
0 The aforementioned ability to reduce longshore transport in a gradual way, thus
reducing the leeside erosion behind massive structures. An adequate screen field
3
layout will distribute the total sediment deficit over a larger area, thus decreasing the
rate of recession of the shoreline.
Because of the different effects induced by completely impermeable groins, and groins with different degrees of permeability, it is important to develop a model that includes the correct boundary conditions at the grid points next to these structures, so that their effects can be properly accounted for. Wave propagation models reviewed by the author usually do not allow for reflective boundary conditions, while a few do by imposing a complete reflection condition at the upwave side and zero transmission on the other side of the structure. In this thesis, adequate boundary conditions that allow for different values of reflection and transmission coefficients are derived and implemented in the model.
To solve the wave propagation problem with the inclusion of wavecurrent interaction effects, such that the model can eventually couple with a nearshore circulation model, two equations governing the dynamics are needed, the wave action equation and the wave energy equation. Another basic equation to be satisfied is the conservation of intrinsic wave frequency. All of them are formally derived first in this study.
Since the model is to be applied in the nearshore zone, a surf zone model is necessary. In this work, the model first developed by Lee and Wang (1993) is used. The model is based on wave action and wave energy conservation and fully includes wavecurrent interaction.
1. 1 Literature Review
Much of the early work done in the area of wave propagation was based on wave refraction theory and entailed the construction of wave rays using linear potential wave theory. This approach is an extension of Snell's law from optical wave refraction to the
4
analogous problem of water waves propagating over straight and parallel contours. Thus, an important assumption made was that the contours of the given bathymetry were locally straight. These ray methods do not yield solutions at prescribed points, but rather along the wave rays only. During the 1960s and early 1970s the linear wave refraction problem was solved in a more efficient way through the use of computers. Noda et al. (1974) was the first to devise a numerical scheme that solved the wave energy equation and used the irrotationality of the wave number to obtain wave height and direction at points on a rectangular grid. Background currents were also included in the formulation.
A serious limitation of the refraction models is the omission of diffraction effects, which creates unrealistically large wave heights when the wave rays converge to a small area. The problem of water waves diffraction was first studied analytically by Penney and Price (1944), who showed that Sommerfeld's solution of the optical diffraction problem is also a solution of the constant depth water wave diffraction problem. Again, Penney and Price (1952) solved the constant depth diffraction problem for breakwaters assuming a semiinfinitely long, infinitesimally thin barrier. In order to include diffraction by coastal structures in existing wave propagation models, first a solution to the diffraction problem was given, using tabulated analytical solutions for constant water depth, and the results from the ray method were superimposed after a few wave lengths.
Berkhoff (1972) derived an elliptic equation describing the complete wave transformation process for linear waves, including both refraction and diffraction, from deep to shallow water in terms of velocity potential function with the assumption of mild slope. The elliptic equation was first solved using finite element methods with good success but at high computational effort. Also, the treatment of boundary conditions is generally difficult.
5
Since then many others have solved the wave transformation problem for complex but idealized bottom configurations. Ito and Tanimoto (1972) proposed a numerical method for harbor wave analysis that allowed easy specification of boundary conditions, but requires modification to describe wave shoaling. Both Radder (1979) and Booij (1981) derived parabolic approximations to the original mild slope equation, the latter extended the scope of the equation to include the effects of currents. Kirby (1984) made corrections to the equation derived by Booij. Ebersole (1985) suggested a finitedifference solution to the elliptic equation in a manner suitable for large coastal areas. Watanabe and Maruyama (1986) proposed a time dependent version of the mild slope equation consisting of two first order equations, which are obtained by separating the original mild slope equation in terms of the water surface elevation and the flow rate. Madsen and Larsen (1987) also developed a model based on a system of hyperbolic equations similar to those governing horizontal flow. The model used an ADI algorithm to iteratively converge to a stationary solution. The difference with the set of equations used by Watanabe and Maruyama, as well as Copeland (1985), was that the harmonic time variation was extracted from the equations, so that the method converges to a final solution more rapidly.
Panchang et al. (1988) presented a method to solve the boundary value problem of the elliptic mild slope equation in large domains based on a marching 'Error Vector Propagation' method that solved backscattering and propagation in the x direction, thus overcoming the major limitation of the parabolic approximation.
Lee and Wang (1993) evaluated the performance of several numerical wave models, including two models based on a set of hyperbolic equations, two based on an elliptic equation and finally a fifth using a parabolic approximation. They concluded that no single
6
model outperformed the other and thus the selection of a model for application depends on the intended purpose.
More recently, in an attempt to model irregular nonlinear waves in very shallow water, models based on Boussinesq equations have been used (Abbot et al. (1978)). Those equations include nonlinearity as well as frequency dispersion. Furthermore, they operate in the time domain, so that irregular waves can be simulated.
Theoretical solutions for reflection and transmission coefficients and for the wave field have been derived for porous structures using different methods of approach (e.g. Sollit and Cross, 1972; Madsen 1974). Recently Darlymple, Losada and Martin (1991) have obtained results for transmission and reflection coefficients for obliquely incident waves. In their studies the equivalent work condition by Lorentz is established yielding a potential flow problem. Matching pressure and horizontal mass flow at the region interfaces solutions of reflection and transmission coefficients are then obtained. They also showed how neglecting the nonpropagating evanescent wave modes the solution technique can be simplified for all practical applications.
1.2 Summary of Contents
It has been mentioned already that the objective of this study is to develop a numerical model that simulates the nearshore wave field in the presence of porous structures such as offshore breakwaters, jetties and groins.
In Chapter 2 the equations governing the wavecurrent interaction problem are derived. Basic assumptions and fundamental equations are first presented, from them dynamic equations describing wave energy transport and wave action conservation as well as the kinematic intrinsic frequency conservation equation are obtained. Chapter 3 derives
7
both the hyperbolic and elliptic mild slope equations that govern the refractiondiffraction of waves. Based on the equation of hyperbolic type, and eliminating the harmonic time dependence, an elliptic equation more suitable for a numerical solution is presented. Chapter 4 describes in detail the finite difference form of the governing equation and the numerical solution by a forward marching Gragg's method. Boundary conditions are established on all sides of the numerical domain as well as on grid points immediately next to possible structures. A section in this chapter gives a general description of the wave model, detailing its components and how they interact until a final solution is reached. The resulting model is tested against several analytical solutions to basic wave propagation and wavecurrent interaction problems. Chapter 5 introduces the theory for wave transformation in the surf zone based on the work of Lee and Wang (1993). The surf zone model is implemented in the numerical model and computed values are compared to analytical solutions for a beach with straight and parallel contours. Chapter 6 demonstrates the applicability of the model based on the results obtained for several situations in which measured data are available from the literature. The model is also used to obtain the wave field around permeable groins in constant water depth and on a planar beach. The case of a groin field is also simulated with the model. The preliminary work and results leading to the complete incorporation of a nearshore circulation model are summarized too in this chapter. Chapter 7 discusses the results, the advantages and drawbacks of the model, and makes several recommendations for further study. One Appendix, where the theory for a linear wave impinging obliquely on a vertical sided porous structure is presented, is included.
CHAPTER 2
WAVE ENERGY PROPAGATION THROUGH VARYING CURRENT FIELD
When waves propagate through a region with currents, wave energy propagation speed is no longer equal to the wave group celerity, but it is modified by the current effect. Other characteristics of the wave train will also be altered, including the wave height, length and period. This situation is commonly observed in regions with strong currents, like tidal inlets. In this chapter, two fundamental dynamic equations governing the behavior of wavecurrent interactions are derived. They are the wave energy equation and the wave action equation.
In Section 2. 1, the basic problem formulation and boundary conditions are given. In Section 2.2, the depthintegrated wave energy equation is derived and the exact forms of wave energy and wave energy flux are presented in waveaveraged quantities. Section 2.3 introduces the wave action equation and the conservation of intrinsic frequency equation. It will be shown how the wave energy equation originates from the conservation of energy, while the wave action equation is derived from the free surface boundary conditions.
2. 1. Governing Equations and Boundary Conditions
The first assumption made is that the fluid motion is irrotational and thus a velocity potential 0 exists and water particle velocities are given by Vq5. The kinematic and dynamic boundary conditions to be satisfied at the free surface, z= q, are, respectively,
27t + VhO'Vhr 0,z =0 (2.1)
1
+ l(VO)2 + gz = C(t) (2.2)
2
where C(t) may depend on t, but not on the space variables. We may take C(t) 0 without any essential loss of generality. The subscripts t and z indicate the differentiations with respect to time and zaxis, respectively. The symbol Vh representents the horizontal gradient operator.
The cartesian coordinate system is used with origin at the still water level, x(x,y) in the horizontal plane and z directed vertically upwards. The velocity vector, U(u,v,w), is related to qbby
u, v and wax ayaz
The velocity potential and free surface displacement are assumed to be composed of current and wave components.
O( x,z,t) = 0b(x;z, t) + eow(X,z,t) (2.3)
q(x, t) = r7(x; t) + erlw(x, t) (2.4)
where e is an undefined factor used to separate the current from the wave part of the velocity potential. The symbol ';' is used to express that the current component of the velocity potential as well as the free surface vary slowly over time when the time scale is much larger that the wave period. Also, qb, can account for small vertical variations of the current. Equations (2.1) and (2.1) are expanded in a Taylor series around the mean water level z= q,
['it+ VhO'Vh'J VZ]z:.c + V2
17t+Vh h zJzc+Cw[ (t +VhOVh)+ h z=Q+ =.. 0
lot + l(VOq) +zl t Z]z +... 0
[ 32 z +8 b + )2 1gz= 0
Substituting equations (2.3) and (2.4) into the above equations gives [(47,)t+Vh c'Vh7c( c)z z=, t (w)+Vh c'Vh7w +VhgwVhzc(w)2+ wh 'c= c=C (2.5)
( Vc) c)2 + g[c] z= + e(0w)t + c w VwV + z=.=O (2.6)
226
The above equations are now separated for wave and current parts and truncated to retain only the first order terms, 0(e):
at~
( = 0 +c Vhc (2.7)
S 2 (2.8)
g 2
(2.8)
(O~w)z = ] + (V2" I + Vh95. Vh 17
w D t c whh) w q (2.9)
1 = D (2.10)
g Dt
where D/Dt /at+Vh ocVh.
Note that the last two terms in Eq. (2.9) are, in general, of higher order than the first one, only when diffraction occurs they become important. For slowly varying water depth, the wave part of the velocity potential may be written as
(Pw(X,z,t) =f(z) (x,t) (2.11)
where f(z) = cosh k(h+z)/cosh k(h+ q) is a slowly varying function of x, k is a real wave number. For progressive waves, the velocity potential and free surface displacement can be written in terms of wave averaged, slowly varying quantities as cw(x,z,t) = f(z)A(x;t) ie if (2.12)
i/w(x,t) = a(x;t) e i (2.13)
where a is commonly defined as the wave amplitude. The phase function is defined as Vf= (Kxwt), where K is the wave number vector including the diffraction effects, and aW is the absolute frequency. The relation between a and A can be established by the dynamic free surface boundary condition specified in Eq. (2.10), which, after substituting Eqs. (2.12) and (2.13) into it, yields
Dq5f
Dt
gaeip= { +V} {Aiei'}
8t
= oAeip+ I + VA}iei, (2.14)
where od is the intrinsic wave frequency including the diffraction effects. Its value is defined as o= UK,and U as Vh 0.
Equation (2.14) states that a and A should have a phase difference unless we impose the condition
o +UV.A =0 (2.15)
at
Then, the relation between A and a can be given by the following familiar expression
a
A= ga (2.16)
Similarly, substituting Eqs. (2.12) and (2.13) into the kinematic free surface boundary condition given by Eq. (2.9) yields o = gktanhk(h + io7) gVA. V (2.17)
CA
8a +V"(Ua) +AK.Vi=0 (2.18)
at
Again, the last term in both the above equations reflects the wave diffraction effect and, under normal circumstances, is of a higher order.
2.2. Wave Energy Equation In this section, the Eulerian expression of energy equation, which governs the local balance between the rate of change of energy and the divergence of energy flux at a point, will be presented. The Euler equation for incompressible and inviscid flow is pDU =V(p + pgz) (2.19)
Dt
taking the scalar product of U(u,v,w) with the respective terms in Eq. (2.19) and then summing the three components yields
p q2 = U V(p +pgz) (2.20)
Dt 2
where q2=(uZ+v2+w2). With the use of the continuity equation, the mechanical energy conservation equation becomes
[1] + V [U( +p +pgz)] = 0 (2.21)
at 2 2
Integrating Eq. (2.21) over the water depth,
1[j2 ]+V[U(P 2+p+pgz)]Idz=0 (2.22)
fht2 ]2
The first term inside the integral represents the volumetric rate of energy change and the second term gives the energy flux through the enclosing surface. The depthintegrated wave energy equation has been derived from equation Eq. (2.22) by many including LonguetHiggins and Stewart (1961) and Witham (1962). A brief account is given here.
Using Leibnitz's rule, Eq. (2.22) can be written as
f [q]dz U( +p+pgz)]dz
f at 2 2
h h
_[U(p 2 +p +pgz)]7. Vh q [U( +p + pigz)]h hh (2.23)
22
+[w( +p+pgz)]h =0
2 h
Substituting the kinematic boundary conditions at the free surface and at z=h,
w  U"Vh r=0 (2.24)
at
W h +UVhh =O (2.25)
and letting p be zero at the water surface, Eq. (2.23) becomes [ 2]dz + pg a +v. fU(pq 2 +p + pgz)dz 0 (2.26)
t f t at 2
h h
Defining the total energy, E, and the total energy flux, F, as
E = [ dz + P 2 (2.27)
2 2
h
77 2
F= f U(2 +p+pgz)dz (2.28)
h
Equation (2.26) yields the well known energy equation given by LonguetHiggins and Stewart (1961) and Witham (1962), 8E
at +VhF =0 (2.29)
at
Here, E is the energy density in the water column per unit surface area and F is the energy flux through the vertical surface enclosing the water column.
2.2.1 TimeAveraged Wave Energy
The energy per unit surface area given in Eq. (2.27) is expanded in a Taylor series with respect to the mean water level z= c
17c
E=p ( 2 2 '2 h 2 1 h272
E p J lI[(Vh ')2 +GP()')2>+pg[(q~c +eg.])n 2] +pe.7[(Vhq)2+ (q5)2]r
h
=P f [(Vh+ h 2 +(q5z + eq5Z)2 +lpg[(17,e7)2 h 2]
h
2
+ pCeli(Vh c h c w21O (2.30)
Taking time average over the wave period and collecting terms associated with current (O(1)) and wave (O(82)) separately, we obtain the following pair of equations,
Ec =P f 1[(V 2 c5) 2 + (,)z + 1 pg(7c h)2 (2.31)
h
E f I(Vh )2 + (w)2Z+ pgq 17+ Pw [Vh c.Vh w] (2.32)
h
where Ec and E are, respectively, the mean values of current and wave energy. The mean wave energy density can be expressed in terms of the slowly varying quantities by substituting Eqs. (2.12), (2.13) and (2.17) in Eq. (2.31): E 'H 2 (2.33)
8
where H is the wave height defined as twice the wave amplitude a.
A few general remarks regarding the definition of wave energy density are made here: i. E is the energy density directly associated with the O(e) fluctuation motions only.
It does not take into account contributions associated with mean water level change. ii. The last term in Eq. (2.32) represents the contribution due to wavecurrent
16
interaction. LonguetHiggins and Stewart (1961) and Phillips (1977) all included this term in the mean flow energy rather than in the mean wave energy. However, this
is truly a O(e) term from the fluctuation motion.
iii. In the absence of current, E reduces to the conventional definition of energy density
in a wave field.
2.2.2 TimeAveraged Wave Energy Flux
The energy flux expression given in Eq. (2.28) can be given, using the Bernoulli equation, as
77
F =p pfU8'dz (2.34)
Introducing the current and wave components defined in Eqs. (2.3) and (2.4) and expanding in a Taylor series with respect to the mean water level z= qk, we have
1Ca a
F p f V(qc + + egbw)dz M]. [Vh(O, +00)4c +ebw)], (2.35)
h
Taking time averaging over the wave period and collecting terms of 0('E), we have the mean wave energy flux,
P = P W adz pw[vhocW+vhowato, (2.36)
h
Substituting Eqs. (2.12) and (2.13) into Eq. (2.36) it can be obtained that,
'P 'CiOt k]E (2.37)
The above equation contains an additional term, q5,ck/a, that does not appear in the conventional form of this equation. This term will be of importance only in the case of unsteady currents, and will not be included in the derivation of the governing equations for this model. As explained above, the definition of E is different than the conventional term too.
2.2.3 TimeAveraged Wave Energy Equation
Now Eqs. (2.33) and (2.37) can be substituted in Eq. (2.29), obtaining,
aE + Vh (Cg+Ui Oltk)E = (2.38)
at (A
2.3 Wave Action Equation
Subtracting Eq. (2.10)xpgir/w from Eq. (2.9)xpg4, and ignoring the higher order diffraction effect, we obtain
(p~a ) A 2pqw;b ,~~ (2.39)
 ( P 7w 06w) + V7 (Uf, Pr w)7, P 117w lc g 17, =0 (.9
where U. is the current velocity of the mean flow at the water surface level. Substituting Eqs. (2.12) and (2.13) into Eq. (2.39), the following equation is obtained:
(Bie 2i f) + (UJsBie 2i) + aBe i2 V gktanhk(h + 0) where B is defined as B _pg H 2
8 a
Expanding and separating the harmonic motions,
(2.40)
ie2i [ aB + V'(jsB) + aBei2 2 ktanhk(h + = + 1] 0
at S 21
which yields the dispersion relation, o2 = gk tanh k(h + q,), and the following wave action equation:
aB
+Vh [U sB] =0 8t
(2.41)
If we substitute Eq. (2.16) into Eq. (2.18), the following equation is obtained,
aA +Vh [Us aA] = 0 8t
(2.42)
Eliminating A from Eqs. (2.42) and (2.15), we arrive at the final equation that governs the change of the intrinsic wave frequency in a current field:
Ou
+Vh [Us] =0
at
(2.43)
which is termed as the 'the kinematic conservation equation of intrinsic frequency', or simply the 'conservation equation of intrinsic frequency'. It is possible to derive two other
19
kinematic conservation equations. The phase function was defined as t= (k.xw6t), where k is the wave number after dropping the diffraction effects. Thus we can obtain k and wo as Ot
Then, if we take the curl on k, we find that Vxk =0 That is, the irrotationality condition on k. It also follows that At+V6=O
at
known as the 'kinematic conservation equation of the wave number'.
CHAPTER 3
MATHEMATICAL WAVE MODEL
A significant advance in the field of wave modeling was made by Berkhoff (1972, 1976), who derived the combined refractiondiffraction mild slope equation. This equation eliminated the usual problem encountered in refraction studies caused by caustics. It can be used for a wide range of ocean wave frequencies, since it converts, in the limit, to the deep and shallow water equations. The usefulness of Berkhoff's equation in yielding good simulations of wave behavior in a wide variety of situations has been proven by many investigators. Berkhoff, Booy and Radder (1982) have used it to study short wave propagation over arbitrary variations in bottom topography. Berkhoff (1976) and RottmannSode and Zielke (1984) employed it to study wave propagation into harbors. However, when attempting to simulate wave conditions in large coastal areas, the numerical solution becomes a formidable task in terms of computer time and memory.
In view of this problem, it became necessary to use the parabolic approximation (Radder, 1979) of the mildslope equation. With this method a solution is rendered feasible in fairly large open coastal areas, but there are two limitations to it, the waves must have a dominant direction of propagation, and reflections are neglected.
The extension of the mild slopeequation to cover wavecurrent interaction has been developed by Kirby (1984), and the corresponding equation eventually yielded the energy equation suggested by Bretherton and Garrett (1969):
DVCg~ 2(o22CCg) bn2 + (v.D)O) V( ccQV + (o kcc) o Dt2 Dt
where q is the complex potential at the water surface.
In Section 3.1, the mild slope equation revised by Kirby (1984) will be obtained from the energy equation, assuming that the steady state condition of wave and current field could be retained for the waveinduced currents which are usually accompanied by the waves. In Section 3.2, the governing wave equation of the wavecurrent interaction model will be derived from the mild slope equation. The solution obtained does not rigorously satisfy the boundary condition for a sloping bottom, but even the firstorder effect of the slope is negligibly small. Booij (1983) and Smith and Sprinks (1975) numerically examined the validity of the mildslope equation and concluded that it is applicable to waves propagating over slopes as steep as 1/3, and even to waves propagating across as step.
3.1 Mild Slope Equation
The mild slope equation has been derived by the variational principle and Green's second identity. In this section, however, the mild slope equation will be derived directly from the energy equation as described below. Introducing the velocity potentials and free surface displacements composed of current and wave components as given in Eqs. (3.33.4), Eq. (3.26) becomes
f[ (V e + e)2]dz + g(c + erlw)(rc+ e77w) ft 2 at
h (3.1)
Vh f[(Vc + eVkc)(,Oc + O)Odz 0
h
Taking Taylor expansion about z= qc,
(VqC + eW)2 a (VOq + 0w)2
t f 2 te 2 ]at+g(+e)(+e)
h
vh. f[(Voc + eV)(4 + eqw)t]dz + eq[(Vc, + eVw)(c + e~,),] = 0
Collecting the terms of O ) and ignoring the longterm fluctuation of current field:
Collecting the terms of O(i) and ignoring the longterm fluctuation of current field:
a (V,)2 a a
f[ ] ]dz + [ RwVc.Vw]+ g Vh w[V 4wt]dz +V c ww.
h h
0 (3.2)
Substituting Eq. (2.11) into Eq. (3.2),
Sff2dz[
h
+ 7a g
+7 at
+ f 2dz [ )
2 h 2
]+ [ 7,V ; wj
+at
Vh f2dzVv w, + V 4 j= 0
h2
where, substituting the expression forf
f2dz
h
Sdz =
CCg
g
o k2CCg
Therefore, Eq. (3.3) becomes
ccg a (Vh 2
[W]
g at 2
+ ag t
+ o2k2CCg a )2
+ []~ ] M~,+ V V < ]o
g at 2 at
Vh wt cV w wt =0
(3.3)
(3.4) (3.5)
23
Hereafter, we omit the subscript denoting the wave mode, w.
CCgVh Vh t + t(02 k2CCg) b + g(cV'.V)t + g 2t
 tVh (CCgVO) CCgVq5 Vh qt gVh (VOc 7 t) = 0
The first and sixth terms offset each other and expanding the third and last terms,
t(02 k2CCg) t + gctV'V + g+V'Vt + g2t
 Vh(CCgVO) g Oh c(V 7) VcVh = 0 and also offsetting the third and seventh terms, o,(d2 k 2CCg) + gtv Vc.V + g2,
Oth.(CCgVO) gt hV(V c ) 0 Substituting the dynamic free surface boundary condition of O(e) given in Eq. (2.10),
t(d k 2CCg) q + g,7Vq .Vo g( qt + Vc.Vo) 1t
tVh'(CCgV) g0 tVh'(V' q) =0 t( k 2CCg) q g t h ttV(CCgV) gqtV(V5c74) =0
Combining the second and fourth terms by use of the kinematic free surface boundary condition of O(e) given in Eq. (2.9), (02 k 2CCg) Vh(CCgV) gq = 0 (3.6)
Now we consider two kinds of expressions of z. The first is obtained by substituting the dynamic free surface boundary condition, Eq. (2.10), into the kinematic boundary condition,
(3.7)
g Dt2 Dt J
and the second expression is simply obtained from Eq. (2.11) as
o2~
qS~q5
(3.8)
which eliminates the slowly varying motion through the mean surface. With the first relation, we finally get the mild slope equation as derived by Kirby (1984);
D 2 + (Vh.J)D Vh.(CCgVh ) +(o2 k2CCg)q0=O
Dt2 Dt
and with the second expression, we get the mild slope equation of elliptic type: Vh'(CCgVq5) +k2CCg =O
which has the same form as the case of no current.
Equation (3.9) can be fully expanded, becoming,
2 2 2+(u2 CCg) +(v2 CCg)2 =
t 2 (Vhxyt + 2u O
8t2 aa X2 ay2
where
= (VhU) D +V h h(CCg) (o2 k 2CCg)
Dt
J=ui+vj
(3.9)
(3.10)
(3.11)
25
It can be easily seen that Eq. (3.11) is a three dimensional second order partial differential equation of hyperbolic type. This type of equation is also known as the KleinGordon equation or telegraph equation.
It can be shown that Eq. (3.10) reduces to the wellknown Helmholtz equation in the case of deepwater,
Vq5+k2=0 (3.12)
and in shallow water Eq. (3.10) becomes the long wave equation, k 2CCg(_) +V(r )Vh(CCg) 0 (3.13)
a a7
3.2 Derivation of Governing Equation In this section, the governing wave equation of the numerical model selected is derived from the linearized mild slope equation Eq. (3.9) for waves interacting with currents. It is rewritten below with the symbol descriptions
D2q5 + (V'.J)Dq V'(CCgV) + (o2 k 2CCg) q0 Dt2 Dt
where is the twodimensional complex velocity potential,
C is the relative phase velocity (a/k),
Cg is the relative group velocity (eoaOk),
ais the intrinsic frequency (oi=gktanhkh),
w is the absolute frequency,
k is the wave number,
h is the water depth,
Uis the steady current velocity vector (u,v), and V = a/Ax + a/ay, omitting the subscript h.
The values of a and k are determined by the Doppler relation, 6 = a+ k.
If we express the twodimensional velocity potential in the linear stationary field as: (x,t) = O(x) e iot (3.14)
where q is the surface potential in steady state. Substituting Eq. (3.14) into Eq. (3.9) gives,
 2iaUV V + (U+V)(Vq) + (U.V)(iw + .V )
V(CCgV) +(o2 k 2CCg)= 0
Rearranging and factoring out e iwt,
io[2UVO + (UV) b] + (JV)(IVq) + (.V)(..V#6) (3.15)
V.(CCgVq) + (d2 k 2CCg)q5=O
The above equation is valid as long as the wave motion is time harmonic. The threedimensional (x, y and t) problem has been reduced to a twodimensional (x and y) one with an equation of elliptic type.
If we expand Eq. (3.15), the first term will result in,
iw(24xu + 2 yv+u x + Vy= io(ux + Vy) # i62u x io2vy
The second and third terms will give,
(UV)(u/x +VOy) + (U +Vy)(Uvx +vy) ~Y (2ux+vu +uv )x +u2 xx+2uv +(uvx+2vVy+u ) +V2y
And the fourth and fifth terms,
(CCg)x x CCg xx (CCg)yy (cCg) yy + 2 k C) 2 5 Combining all terms we get,
(u 2 CCg)x + {2iou +2uux +uyv + UV (CCg)x }x
+ 2uvx + { 2i(v + 2vvy + uvx+ uxv(CCg)y} y
+ (v 2 _CCg)yy + {i(ux + 2 k 2CCg } = 0
To study the nature of the above equation we can rewrite it as, aqxx +2bxy + CO, F( qx, qy,,x,y)
where,
a = u2CCg b = uv c = v2Ccg
in general, CCg is >> U2=u2+v2, and, ac b 2 = (CCg)2 CCgU2>O
Thus, Eq. (3.16) is an elliptic partial differential equation. In Chapter 4, a numerical method to solve this governing wavecurrent interaction equation is presented.
(3.16)
CHAPTER 4
NUMERICAL SCHEME AND TESTING OF WAVE MODEL
In the last chapter a governing equation of elliptic type was derived from the original mild slope equation of hyperbolic type. A general solution to this kind of governing equation is not available for a wave field in the vicinity of structures in the nearshore zone; therefore, numerical computations are needed. In this chapter, a numerical method is developed based on the work of Lee (1993). The method falls under the general category of finite difference schemes. Also, the boundary conditions used for the edges of the computational domain and the grid points next to a groin are obtained. In the next section a general description of the wave model is given. The last part of this chapter includes comparisons between the results obtained from the wave model and a series of analytical solutions to basic wave propagation and wavecurrent interaction problems.
4.1 Numerical Scheme of Wave Model
In considering the different possible numerical approaches to Equation (3.16), the author tried to use one that would combine a minimum computational effort and with results that would include refraction, diffraction as well as reflection due to structures. Several mathematical models have been proposed for combined refraction and diffraction of simple harmonic linear waves. All of them have advantages and drawbacks depending on the kind of problem that one is trying to solve. The solution to the mild slope equation derived by Berkhoff (1972) is very accurate and complete, but the treatment of boundary conditions is generally difficult. The numerical method proposed by Ito and Tanimoto (1972) allows easy
29
specification of boundary conditions but requires modification to describe wave shoaling. A parabolic approximation to the mild slope equation is a popular approach to alleviate some of the problems addressed above, but it does not include reflections from structures such as jetties and groins. The numerical solution adopted by the author will include such reflections without making the computational effort excessive.
Equation (3.16) can be treated as an ordinary differential equation for 0, as given below, so that the ydirectional difference operator, D, is explicitly approximated by using finite differences,
(u2 CCg) + { 2iwu +2uux +uv +u v(CCg)x}x
+2uvDy(A)+ { 2iv+2vvy +uvx +uxv(CCg)y}Dy( ) (4.1)
+(v 2 CCg)D(q)+{ iW(ux + Vy) + o2 d k2CCg} =0
where, the x axis is selected for the main direction of wave propagation. We convert the above equation to a pair of firstorder equations by simply defining the x derivative as a second function.
lx 21 [{2io)u+2uux+u v+y (CCg)x} 1 + 2uvDy(,)+ ()]
CCg u 2
where
7 = {2iev+2vvy +uvx +uxv(CCg)y }D( )+ (v2CCg)Dyy( ) +{ i(ux + y) + o2 k 2CCg}
In this study, the ordinary differential equations are numerically solved by Gragg's method, whose main algorithm for a differential equation q5x(x) =f(x, O(x)) is given as
30
Yl fi1 +hf(xi1,qi 1) yj+ = yi1 +2hf(xilI +jh,yj) ,j= 1,2,..,n1 i (Yn + Yn1 + 2hf(xi, Yn))/2
where h is a subgrid space defined as h = Ax/n as shown in Figure 4.1.
WAVES
V y axis
i1
i +1
, Main Grid
  SubGrid Figure 4.1: SubGrid system for wave model.
If we obtain the truncation error of the above equations as a power series in h: Error= a h 2i
i=1
That is, it only contains even powers of h. If n is even, the result for n steps would be q, and for n/2 steps 9,/2. Combining both solutions:
4 3q5 2
3
31
where q5 is now fourth order accurate, but requiring about one and a half derivative evaluations per step h instead of four.
The governing wave equation used, as well as the numerical solution used are capable of solving the reflected wave field generated by structures posed along the x direction, that is, the reflected potential propagates in the direction of positive x like the incident does. However, the numerical solution is limited in the sense that it does not include reflected waves propagating in the offshore direction, like the ones created by a detached breakwater, which would be posed along the y direction. The model is capable though, of computing the reflected potential backwards. This reflected potential would be obtained as the conjugate of the incident potential at the grid points immediately next to the structure. A solution of this kind would work well in the case of seawalls or 'infinitely' long detached breakwaters, that is, as long as the only potential propagating offshore is due to reflections and thus easy to obtain from the incident potential. In the case of a finite detached breakwater though, the diffraction effects produced by the structure and propagating offshore from its position are impossible to simulate with the method used. Nevertheless, since the objective of this model is to simulate the wave field in the nearshore area around a field of groins, the solution adopted fits the goal. However, this does not mean that detached breakwaters and Tgroins cannot be included in the model, since the diffraction effects induced by such structures downwave of their location will be accurately reproduced.
Once the wave potential is obtained the remaining characteristics of the wave field are obtained as follows. The wave angle is calculated by
6tan '(K/K)(4)
(4.2)
where
Kx = Im { / Ky = DyS (4.3)
where, S= K x =tan [Im( q)/Re()].
The wave height is calculated easily by
H = 2 Re { _ }2 + Im{Z}2 (4.4)
H=2e 1 II0 (4.4)
g g
4.2 Stability Analysis It is possible to roughly obtain a stability condition for the method used by performing avon Neumann stability analysis on the Helmholtz equation (Panchang et al., 1988),
V2q5+k 20=0
which has a solution for constant k ~ e imx e iny
such that m2 +n2= k2. The solution for b can be written as: O(p,q) = e im(p1)AXe in(q 1)Ay = p e in(q 1)Ay Expressing the Helmholtz equation in finite difference form and using q(p,q), P20p +p_ 2cosnAy2 2p
(Ax)2 (Ay)2
33
Defining the amplification factor as G= Op+I/op, the above equation gives G2+BG+I =0
where B=k2Ai2 (2(Ax/Ay) sin a)22, and a=nAy/2. When the amplification factor satisfies a quadratic equation such as the one above, it is convenient to use the formulas derived in CushmanRoisin (1984) for investigating stability. In this case those formulas yield the following condition for stability (Panchang et al., 1988): [B]2
The right inequality is always satisfied if kAx :< 2, or Ax< L/t (where L=27t/k)
The left inequality is always satisfied if Ay L/ht
If the marching Gragg's method described in the preceding section is used with these stability criteria, a stable scheme results. While the condition obtained for Ax is similar to commonly obtained stability criteria for many partial differential equations, the condition for Ay is unusual, and it has both advantages and disadvantages. Most important it explains why a marching method can be used for a wave equation, but is inherently unstable for most other commonly occurring elliptic problems. For instance, for the Poisson or Laplace equations, where k2=0, or L, it is obvious that any finite grid size Ay will result in an
34
unstable scheme. In addition to making the scheme stable, this type of condition reduces the number of computational points required, which is important when trying to model the wave field in the nearshore area, which the area of interest will be usually larger in the y direction than in the x direction.
There are two disadvantages as well: the grid points being spaced L/rT length units apart, information about depth variations between grid columns is not well represented in the model, and the wave components in the y direction are not well resolved. The first disadvantage is not a major drawback, since in coastal areas we do not expect rapid variations of bathymetry along the y axis. The second one makes the numerical scheme approximate, in that, as in the parabolic approximation, the dominant wave direction is still x. This problem will be described latter when testing the model against analytical solutions.
4.3 Boundary Conditions
Radiation boundary conditions are defined at both upwave and downwave sides of the computational domain. In addition, each grid row is divided in ng+1 subdomains, where ng is the number of groins crossing that row (see Figure 4.2). Each one of these subdomains will need lateral boundary conditions. The upwave side corresponds to either the lee side of a groin or simply a virtual boundary, and the downwave side to either the front side of a groin or again a virtual boundary.
4.3.1 Lateral Upwave and Downwave Boundary Conditions
Side boundaries of the computation area are generally virtual open boundaries established for performing the numerical analysis. The presence of such virtual boundaries should not affect the solution in the study area. In other numerical models of this kind it is assumed that there is only one wave train crossing the side of the domain. A simple radiation
35
Wave Direction 0 o0o o o olnternal boundary grid points
M W M X aVirtual boundary grid points 10 a
2 0 i . . . .: . e . .... . . . . i.... . .. .... . . 0 . . . . . .. . . .
20 0 0
o o o
0 0 0 0a
o 0 00 a
o 00 a
0 STRUCTURE
0 0 o 0
30 0 9
30 . . . .. ... 0 o . . .. . . .
z b sle o
=,o ik# hr sn=kso (45
a a
ag a aa
0 0 0
0 0 0 0 a
0 0 0 30 4 0 0 6
Fiur 9Loaio o ri pitswhr th wav poeta
 k#k oR (4.6
wil b sove uigboay codtin
boung con n is ued n wit csebo a
= i ( o
000 0 0 0
in cid ena l................ ...............
Assmig... wecnrert th goudr conito as:
0 0 0 0
o 0 0 0
o 0 0. o o0
60 0 0o 0 0 ......
0 1 0 0 30 4 0 06
0y 0 o 0o
0 0 to 20b 30 40b 50 60~ti' 46
will be sleusnabouny codtin
Assunr codiio is k uwed an thwit case budr odto s
yY Y 0, 45
'30 iky + ORi=kyAqq_ 1 ik'YO R (4.6)
ay
36
This expression is only valid when the absolute value of the wave number in the y direction is the same for both incident and reflected wave fields, which is the case for problems involving reflecting structures posed in the x direction and with straight parallel contours along the y direction. The incident potential will be obtained at the upwave lateral boundary using Snell's law. The boundary condition can then be expressed in finite difference form to relate the wave potential at the first and second columns of the grid. The value of q5 on j=1 will be known and thus we can obtain qb for j=0.
Using finite central differences on row i yields
i
+ O _iky(20, 2 J)
Ay 2
The lateral boundaries are assumed to be located between the grid columns 0,1 and NY, NY+ 1. Therefore, q5 values at the upper lateral boundary can be given as
Siky2Ayqj 1 + ikyAy/2
1 ikyAy/2 1 kyy/2 (4.8)
on row i. The same procedure is used for lower boundary, the only difference is that there is no reflected wave field escaping the domain through this boundary, thus we can substitute q5, for q5 in Eq. 4.7 and using finite differences yields
1 +ikyAy/2 (49
1 ik Ay/2 95NY (4.9)
37
The two boundary conditions can also be further approximated using higher order terms from the Taylor expansion of the y derivatives.
4.3.2 Boundary Conditions for Groins
We will derive the boundary conditions for the general case when waves are incident on both sides of the groin and from different directions. The groin itself will separate two regions, (I) and (II), see Figure 4.3.
B
Groin /
N /yyaxis Region (1) /NRegion (11)
7 /
7 /
03
Fiur 4.:Dsrpino/ ae ilrudagon
Region ~ Fiur 4.3: Decrpto of wave fiin ih rnmtel d difaround av in. d
The total wave potentials for those regions are, respectively,
Region (1): 0I) + ,+q' Region (II): =/~~ 0' ',+
where the velocity potential components can be expressed as:
a= e ik(ycosa+xsina)
S= Rae i[k(ycosr+xsina)+,r] r ik(ycosa+xsina) r= Rae' Prae
i Ta e [k(ycosa+xsina)+6t] a e ik(ycosa+xsina)
Assuming the depth in regions (I) and (II) is the same, thus, the wave numbers are the same too, then:
Se ik(ycosd'+xsina')
'Pr = R'a'ei[k(ycosa'+xsina')+6'] ra e ik(ycosa'+xsina')
 T _,a i[k(ycosa'+xsina')+6't] ta eik(ycosa'+xsina')
where:
,r = Re t = e r = R' e ,
At = T' e i 't
Now we can obtain the wave potential and its derivative in the y direction at y=O and y=B.
0(i) y= = (1 + r)aeikxsina + Ita 'e ikxsinad a0l ly=o = ikcosa(1 r)ae ikxsina ikcosaj'ta 'e ikxsina' ay lr.ae
0fH)Iy=B = (1/y' + p'ry')a'e ikxsina' + ytae ikxsina
0 ) y= i kcosa'(1/y '+P'ry')a 'e ikxsina' + i kcos ay tae ikxsina
Wy
Where: y e ikBcosa r e ikBcosa'
We can rewrite the above expressions as:
= (1 +6r)Y +ftY' y
= ikcosa(1 fir)Yikcosa'f,6'tY'
() = (l1/y'+'ry')Y'Y'+ tY OpB= ikcosa'(lly'+/'ry')Y' +ikcosayY Where: Y = ae ikxsina, a Y'e ikxsina' Now we can combine equations (4.10) and (4.11):
(5.10)xik' (5.11): (5.10)xik'y + (5.11):
ik' = ik'y(l+fr)Y + ik'y'tY' S= ik(1f)Yik' B'Y'
ik'yO + = [ik' (1 +fir) +ik(1 fir)]Y
From equations (4.12) and (4.13):
(5.12)xik :
(5.13):
(5.12)xi ky (5.13):
ik = ikyf/,yY'+iky(1/y' + f'ry')Y'
of,' = ik ftyY+ik (1/y'+I'ry')Y' ik y = [ik(1ly' +f'rY') +ik'(1/y'fl'rY'
Combining (4.14) and (4.15) to obtain:
01_ (1+fr) [ik' + Mo]+ 6[ikt . ']
=ik'y(1 +,Or) +iky(1 fr) ty 001 iky(lly'+pg'ry ')+ik'y(l/y'f8'ry') y
Rearranging terms in the equation above:
(5.10) (5.11)
(5.12) (5.13)
(5.14)
(5.15)
ik (1 fr) (1+fir) ) +
ik( +f) +iky(1fr) i( + k'(1 )+ik) +(1+lr) +6YO ikyft f, 61
(II) ,.
ik(1/y'+p'ry') +iky(1/r''rr) yB U
Combining (4.14) and (4.15) to obtain:
( ) t [ir on l1/ '+P r [ (I ) (TI.1]
ik'y (1+fr) +ik(1 fr) [lk'yO('O +q)Yo] + ik(1/r'+p'rr')+ik'y'(1/y'p6'rr') [1/yOB OyB
Rearranging terms:
ik'yf, r +t t? /
ik'y(1 +fr)+ik(1frr ) + i k'y(1 +9r) +ik(1 +r) yO
ik (1y'f'rr') ( _1) (lly'+ 6'ry') m = 0(5.17)
i k y ( l ly ,'+ p 'r y ') + i k y ( 1 / _' 'r y ) i k ( / y ,'+ r y ) + i k y ( 1 / y ,' f 'r r ) y B
If we assume a'= 900 for pure diffraction waves in downwave side, then, k' =0 and f'r=f't =0 and substituting in Eqs. (4.16) and (4.17):
(O (1 fr) o = i k ( l + .r) I)
(1O fl
(M) tY (I)
(1 fr)
(5.18)
(5.19)
However, equations (4.16) and (4.17) are to be used in the general case when a' 900. For a single groin, we can assume that the incident potential in Region (II) is only
(5.16)
ik (1/Y'+p'ry') + ik'(1/y ' 'ry) B
41
due to the diffracted wave field. The wave fronts corresponding to it are perpendicular to the groin, which means that we can in fact assume a'= 90o. In the case of a groin field, this assumption may create problems, since the reflected waves from a groin 'downwave' may be part of the incident potential in Region (II) for the previous groin in the upwave direction. In that case a' will be different than 900 .
Equations (4.18) and (4.19) will now be expressed in finite difference form. The upwave groin boundary condition will be located between grid columns m and m+ 1, while Eq. (4.19) will be located between columns n and n+l, the groin width, B, is then equal to Ayx(nm), (see Figure 4.4).
y axis
Upwave Downwave
7rB
y axisGr n
0 + 40, I I
4Ia I
] .Groin
Upwave Downwave
I I B I Il d,"
I I I I
SI 4
m m+1 n n+1
in in+1 n +1
Column Index
Figure 4.4: Schematic diagram.
(1 /r) Ay
1 +iky
(1 +,0) 2
/1 12 ( (5.20)
(1 /r) Ay
(1 + r) 2
To obtain q5,, we will use /b+,1 and b, which are already known since the upwave boundary of the groin has already been solved.
Ay (1/i,) Ay
95" y A i~~ (5.21)
'kl(1,6r) Ay
if the complete boundary equations are to be used, it will be necessary to solve a system of equations involving unknowns 40, and 0, as well as k,' which has to be assumed to be equal to icr. This assumption is necessary when solving the wave field numerically. This is because the numerical solution gives the total wave potential, which in the lee of a groin could be the combination of transmitted, diffracted, incident and reflected potentials, consequently, in general we do not know in what direction the wave field represented by Oi is propagating,
Appendix A includes the theoretical derivation of transmission and reflection coefficients for a permeable structure. The case solved represents the simple situation in which there is no other potential in region (I) than the transmitted one. But the coefficients obtained in that simple situation are also valid when solving the general case that includes incident potential on both regions. This is true when assuming linear wave theory, the general case is just the superposition of two simple cases like the one solved on Appendix A. The coefficients can be calculated independently for each incident wave field.
4.4 Description of Wave Model
This section describes the numerical model that has been constructed to calculate the wave field within a given domain. The domain will be usually bordered by a shoreline along one boundary and an offshore region along the opposite boundary, although a completely
43
'open' region is also possible. The set of governing equations is solved using a finite difference method which requires that a computational grid system be constructed. Definition of the grid scheme is shown in Figure 4.5. A grid mesh comprised of (NX+l)x(NY+I) rectangular grid cells of constant size (Ax by Ay) will be considered.
i = 1 .. . . .. . . . !. . . . . . . i . . ... . . .. . . . . . . . . ..
i Wave Direction
i=2
13"
i= 3 . . . .: . . . . . . . . . . . ... . . . . .. . . . . . . . ...
. . . . . . .... . . . .. ...... . . . .. . . . . . . . ..
.. . . .. . . . . .. . . . . .... . . . . . .. . ... . . .
i x . . . . .. .. . . . . . . .. . . . . . . .. . .
i=NX1
i=i xX ...... ......... .. ......... i ....... ....!.....
i=NX
j=O j=1 j=2 j=3 j=NYI j=NY j=NY+I
y axis
Figure 4.5: Computational domain of the model. Specification of the wave parameters at the offshore grid row, height (H), angle (6) and period (T) is the only wave input required for the model when no structures are present. In the case when one or more structures are situated in the area of interest, their dimensions, location as well as physical characteristics, porosity (e) inertial coefficient (S) and friction or damping coefficient (f) will have to be given. The model will calculate the corresponding reflection (R) and transmission (T) coefficients depending on the characteristics of the impinging wave train. It is also possible to directly assign the values of the reflection and transmission coefficients if known in advance.
Based on the given wave period and the specified depth matrix describing the bathymetry in the area of interest, the wave number, k, for each grid point is computed using the dispersion relationship. The dispersion relationship already includes the effects of a possible current, which can be also given in a matrix form for every grid point. With the
44
computed wave number and period we can calculate the phase speed (C) and wave group speed (Cg). Once all these values are known, the subroutine containing the numerical solution to the finite difference equation governing wave propagation can be called. The subroutine contains an algorithm that ultimately calculates the complex steady free surface potential (q5) at each grid point and row by row, starting on row one and 'marching' towards row NX+ 1. Each row will have two boundary grid points at j = 0 and j= AIX+ 1, the appropriate boundary conditions are then solved relating the wave potential at those grid points to the already solved at j = 1 and j = NX.
When the domain to be solved includes structures the method of solution varies slightly, new boundary conditions in the numerical 'subdomains' are necessary, the expressions governing those boundary conditions as well as the external ones were presented earlier in this Chapter. The so called subdomains are simply the two divisions that a structure crossing a row will create in that row. The problem then becomes a matter of bookkeeping and adequate indexing. When the row being solved is limited by a structure on one or both of its sides, a call to a subroutine that computes the reflection and transmission coefficients used in the boundary conditions is made. The subroutine not only calculates R and T but also gives the wave number inside the porous structure, the real part of that wave number corresponds to the propagating wave mode inside the structure, from it we can obtain a phase speed (Cpt,r) and wave group celerity (Cgpor,) inside the structure which will be also necessary when solving the governing wave propagation equation on the boundary grid points.
During the first call to the wave propagation subroutine the model is run over every groin in the domain but without considering breaking, once the values of the wave height (H) and angle (0) are obtained at every grid point from the complex free surface potential, the
45
breaking and surf zone subroutine is called. This subroutine numerically applies the surfzone model derived in Chapter 5. The subroutine first checks every grid point for a simple breaking condition relating water depth and wave height, when the condition is fulfilled the model keeps track of the grid point were it occurred, thus obtaining a 'breaking line' that connects all those grid points. After the model has located the grid points were the waves break, for every column the model obtains a modified group celerity at every grid point with an index greater than the index where the waves break for that column. The modified group celerity (Cg*) includes the dissipative effects present in the surf zone.
Once the breaker line and the modified wave group celerity are determined the wave propagation subroutine is called once again, the waves are propagated in exactly the same fashion as in the first call, with the only difference that when the grid point being solved is beyond the breaker line, the modified group celerity is used instead of the one obtained from the dispersion relationship.
After the second run of the wave propagation algorithm is over, wave angles and wave heights are obtained at every grid point.
4.5 Testing of Wave Model
The model will be tested against theoretical and experimental results for problems that include wave shoaling, refraction, diffraction, reflection, wavecurrentinteraction, etc.
4.5.1 Wave shoaling and refraction
In shallow water waves transform over a sloping bottom, and if incident waves are normal to a beach with straight and parallel contours, the change in wave height is solely due to the change in water depth. This transformation is called wave shoaling. The wave celerity also changes due to the change in water depth or due to a current field. The gradient of wave
46
celerity along the wave crest results in a modification in wave direction. This process is known as wave refraction. To solve the wave shoaling process we make use of the energy conservation equation, while to calculate the refraction we assume a gradually sloping bottom, constant period and thus irrotationality of the wave number. In the particular case of a simple bottom topography in which the contours are parallel to the y axis the irrotationality of the wave number yields Snell' s Law.
Numerical results from the wave model were compared to the analytical solution based on energy conservation and Snell's law. The input data are given as follows:
NX NY Ax (m) Ay (m) Hi (m) T (sec) Slope 60 20 .1524 .3048 .04 0.8 1:20
In order to test how the model reproduces shoaling without refraction, the wave direction is given as perpendicular to the beach. Next, the waves are run at an angle with respect to the depth contours, and from the results we can obtain, dividing by the previously calculated shoaling coefficient, a refraction coefficient.
As shown on Figures 4.6 to 4.8, the model produces results that agree well with the theory for both the refraction and shoaling coefficients. However, if the incident angle is increased too much, due to the nature of the numerical method applied, the refraction coefficient differs considerably from the exact solution: this problem is shown on Figure 4.8. This might induce errors at the boundaries as well. In the numerical scheme, the boundary condition is based on Snell' s law, which will yield the correct incident wave angle along the boundary. However, due to numerical error, the wave potential escaping the domain might propagates in a slightly different direction when the angle becomes large; this will induce unrealistic values of the wave height and direction at the boundary.
47
SHOALING TEST
Theory
0 0000 Numerical Model
0
u.8 s
0.6 0.8 1 1.2 1.4
kh
Figure 4.6: Comparison of wave shoaling.
1.6 1.8 2
REFRACTION TEST
1.15 Theory
1.1 0 0 0 0 0 Numerical 1
1.05
1. ..O"
0.95 0.9 0.85
0.8 I I I
0.6 0.8 1 1.2 1.4
kh
Figure 4.7: Comparison of wave refraction, Oi= 20'.
1.15
1.1 1.05
0.95
0.9 0.85
SHOALING plus REFRACTION
1.15 Theory
1.1 0 0 0 0 0 Numerical Model
1.05
0.95
0.9 0.85 0.8
0.6 0.8 1 1.2 1.4 1.6
kh
Figure 4.8: Comparison of wave shoaling/refraction, 0/= 20'.
1
0.95 0.9 0.85 0.X
REFRACTION TEST
Theory
o0 0 00 Numerical Model
0000 00000 ooo o
0.6 0.8 1 1.2 1.4
kh
Figure 4.9: Comparison of wave refraction, O= 30'.
4.5.2 Wave Diffraction and Reflection
Diffraction is a phenomenon by which wave energy diffuses or flows transversely to the main direction of wave propagation. This allows the waves to propagate into sheltered regions like the lee of a breakwater, jetty or groin.
Penney and Price (1952) derived a solution for constant water depth based on the Helmholtz equation and similar to the Sommerfeld theory of the diffraction of light. The expression they obtained governs diffraction of small amplitude waves around a semiinfinite thin breakwater and propagating through water of uniform depth. It has been already confirmed that this analytical solution agrees well with hydraulic model experimental results. The diffraction coefficient, H/Hi, that is, the ratio of diffracted to incident wave amplitudes, is given as follows in terms of the polar coordinates r and a (Figure 4.10):
KHH ( 4kr sin 2) jos(Oa 14sin e ikrcos(6+a)1 (5.22)
K =H/i1I sin esi
d T 2 7tU
where
2 f 2
Equation (4.22) includes both the effects of diffraction and reflection from the structure posed along the x direction.
Equation. (4.23) can also be expressed as
I(A) = 1 +CQ) +S(2) + i C(A) S(1) (4.24)
2 2
50
in which the Fresnel integrals, C(I) ans S(A) are defined as A ir2A iA
C(A)=fcos TEAd,,A, S(A) fsin EAd (4.25)
J 2 J 2
0 0
X0
x/
(XI Structure
r
S/
(Xy)
i/Waves
Figure 4.10: Coordinate system for analytical solution.
The model capability of handling wave diffraction was evaluated by comparing the numerical results with the above analytical solution. Since the model is limited in terms of the wave incident angle, two different structure layouts where used. In the first one (see Figure 4.11) the given conditions were Hj= 1 m, depth=5 mn and0L =30.0 m, which corresponds to a wave period, T, of 4.96 seconds. The grid used had dimensions NX= 60, NY= 60, Ax= 10.0 mand Ay= 15.0 m. The structure was arranged along the x axis, that is the main axis of wave propagation. Therefore, the numerical solution included not only the diffracted wave field behind the structure, but the combination of reflected, diffracted and incident wave fields in front of it.
3 .. . . . . .. . . .
Sec. 3
Sec. 4
1 2 . . . . . . . . . . . . . . . . . .
STRUCTURE
1 5 . . . .
See. 5
0 3 6 9 12 15 18 21 24 27 y/Lo
Figure 4.11: Grid and structure layout used for first diffraction
test. The values of x and y refer to model coordinates.
Figure 4.12 shows the results for two different incident angles, 50 and 150. Figures 4.13 and 4.14 include comparisons between theory and the numerical solution along five different cross sections. In all the aforementioned figures x and y as well as angle 0 refer to the coordinate system for the analytical solution. The results seem to compare well for these small angles of incidence, of special significance is the agreement obtained in the region upwave of the structure were reflected and diffracted wave fields interact, since the final objective of this model is to simulate the wave field around a field of groins. Such groins would be posed along the x axis and the reflected wave field that they generate would greatly affect the resulting nearshore circulation.
In the upwave region (Figure 4.14, Sec. 1) the amplification factor (H/H,) actually becomes larger than two, this is due to the diffraction effects that are superimposed over the reflected and incident waves. Such effects will tend to disappear away from the tip of the breakwater.
8
8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 y/Lo
12"
030.
00.
16" \
18
8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 y/Lo
Figure 4.12: Comparison of wave diffraction/reflection coefficient. Top O = 50, bottom O= 15'.
2.5
2
'1.5
0.5 0 000 Numerical Model
18 16 14 12 10 8 6 4 2 0 1
x/Lo
18 16 14 12 10 8 6 4 2 0 1 x/Lo
F
eeeGeee~e2~
1 1 I I I I I I I I
Sec. 3
8 7 6 5 4 3 2 1 012345678 y/Lo
2.5 . .
2 Sec. 4
i1.5
0.5
8 7 6 5 4 3 2 1 012345678 y/Lo
8 7 6 5 4 3 2 1 012345678 y/Lo
Figure 4.13: Comparison of wave diffraction/reflection coefficient along five different sections, 0= 50.
2.5
2
1.5
0.5
n
2.5
2
0.5
1
0.5
2.5
2
0.5
18 16 14 12 10 8 6 4 2 0 1 x/Lo
2
0
2.
0.
.5
2 .5
1
.5
01.
18 16 14 12 10 8 6 4 2 0 1 x/Lo
5
1 I I I I
8 7 6 5 4 3 2 1 012345678 y/Lo
2.5
2 1.5
1
0.5
0
8 7 6 5 4 3 2 1 012345678 yILo
2 Sec. 5
1.5
I L
0.5 0
0 I,
8 7 6 5 4 3 2 1 012345678 yILo
Figure 4.14: Comparison of wave diffraction/reflection coefficient along five different sections, O= 15'.
0 0 00 Numerical Model
Sec. 2
55
For larger angles of incidence the model starts to give inaccurate results, this problem is illustrated in Figure 4.15. The extension of the diffraction zone seems to fall short of the theoretical solution, the same can be said about the reflection region on the other side of the groin. The variations of wave height in the latter region along a section perpendicular to the groin increase with increasing angles of incidence. The model is not capable of simulating these variations since there is a limit to the size of the grid the y direction, Ay.
In order to simulate cases with large angles a different layout is used in which the breakwater is posed along the y axis. The input conditions are the same as for the first layout with the exception of Ax which is equal to 5 meters now. Figure 4.16 shows the comparison for waves approaching a breakwater at 90', the model seems to agree reasonably well with the theory for this case. Figure 4.17 shows the results when the waves approach the breakwater at a 70 degree angle and directed away from the breakwater. The results for this case are not as good, most likely due to the limitations of the numerical model for higher angles of incidence. The results can be improved by reducing Ay, but as explained earlier in this chapter, there is a numerical limit to that reduction. It is noticeable how the diffraction bar that starts at the tip of the breakwater is higher than what the theory predicts. This effect is not observed when the structure is posed in the y direction as in the first layout used. As explained in the first part of this chapter, the numerical solution is based on a marching method along the x axis, this means that diffraction and reflection effects that propagate in the x direction from the body and tip of the breakwater are not obtained.
2.5
2
0.5
2.
1.
0.
18 16 14 12 10 8 6 4 2 0 1
x/Lo
5
2
5 Sec. 2
1 5
18 16 14 12 10 8 6 4 2 0 1 x/Lo
2.5
2
'1.5
0.5
8 7 6 5 4 3 2 1 012345678 y/Lo
2.5
2 Sec. 4
'1.5
0.5
8 7 6 5 4 3 2 1 012345678 y/Lo
2 .5 . . . .
2 Sec. 5
1.5 C
1 oOO
0.5 I0
0.
8 7 6 5 4 3 2 1 012345678 y/Lo
Figure 4.15: Comparison of wave diffraction/reflection coefficient along five different sections, 0=30o.
. 000 Numerical Model Sec. 1
Sec. 3
0 00 0900
] I I I I I I I I
57
DIFFRACTION COEFFICIENT
5 4 3 2 1 0 1 2 3 4 5
y/L
WAVE FRONTS
01 1 1 1
5 4 3 2 1 0 1 2 3 4 5
y/L
Figure 4.16: Comparison of wave diffraction coefficient, 0= 90'.
58
DIFFRACTION COEFFICIENT
5 4 3 2 1 0 1 2 3 4 5
yAL
WAVE FRONTS
5 4 3 2 1 0 1 2 3 4 5
y/L
Figure 4.17: Comparison of wave diffraction coefficient, 0= 700.
4.5.3 WaveCurrent interaction
When surface waves of any kind propagate over the surface of a medium in steady but nonuniform motion, they tend to undergo changes in length, direction and amplitude. In the case of water waves propagating over a current field the common assumption was, at first, that the wave energy is simply propagated with a velocity equal to (U + Cg), where Cg is the vector group velocity and U is the total stream velocity, and that no coupling between the waves and the current took place. However, as described in Chapter 2, waves are modified to a much greater extent than would be predicted using that assumption when riding through a nonuniform wave field.
Collinear Current Shearing Current
y axis y axis
Waves Waves
Lon .4
Current XCurrent
Figure 4.18: Conditions of collinear and shearing current.
The performance of the model in terms of wavecurrent interaction is tested through two different cases, one with collinear wave and current directions, and the other with a shearing current that produces refraction of the wave field. Figure 4.18 describes the two
60
situations tested. The analytic solution to both cases is given by LonguetHiggins and Stewart (1961). The comparisons with their solutions are shown on Figures 4.1921. The given conditions are H =0.1 m, depth=3 m and T= 1 sec. In the case of the shearing current the initial angle between the wave fronts and the current is 30 degrees.
The grid used had dimensions NX= 101, NY=21, Ax=0.1 m and Ay=0.6 m. The model solution for the problem of collinear wavecurrent compares very well with the theory.
COLLINEAR CURRENT(Amplification)
3
2.5 'T=1 sec
depth=3 m
2
_ Theory
1.5 0000 Numerical Model
0.5
0
0.2 0 0.2 0.4 0.6 0.8
U/Co
Figure 4.19: Comparison of collinear wavecurrent
interaction.
SHEARING CURRENT (Amplification) 2i 1 1
1.
Theory
0.5 0 0 0 0 Numerical Model
0 1
0.2 0.1 0 0.1 0.2 0.3
U/Co
Figure 4.20: Comparison of shearing wavecurrent
interaction, amplification of wave height.
T=1 sec
5 depth=3 m
Initial Angle=30 degrees
61
SHEARING CURRENT (Angle change) 50
45 T=1 sec depth=3 mn
40 Initial Angle=30 degrees
a: 35 Theory
0o o o o Num. Model
25
20
0.2 0.1 0 0.1 0.2 0.3
U/Co
Figure 21: Comparison of shearing wavecurrent interaction,
change in wave angle.
In the case of the shearing current, the model predicts adequately the amplification factor, but the change in wave angle deviates considerably from the theory, especially for
high positive currents.
CHAPTER 5
SURF ZONE MODEL
The flow properties in the surf zone are extremely complex due to the strong interactions among motions induced by waves, currents, and turbulence. The present knowledge on surf zone dynamics is still inadequate and most of the models are rather rudimentary. Most of the available models that describe wave breaking are based on an approximation to wave energy conservation. These models can be classified into two groups: one based on the similarity of the wave breaking process with other hydraulic phenomena such as a hydraulic jump (Daily et al., 1984), a tidal bore (Battjes and Janssen, 1978), etc., and the other is based on estimation of the eddy viscosity (Mizuguchi, 1980) or turbulence (Izumiya and Horikawa, 1984).
In this chapter a simple surf zone model proposed by Lee and Wang (1993) is presented. The model is based on the consideration of wave energy balance and wave action conservation so that the wavecurrent interaction is fully taken into account. The model is presented in analytical form for the case of two dimensional gradually sloped bottoms.
5.1 Time Averaged Wave Energy Equation in the Surf Zone
It is assumed here that the surf zone maintains a wavelike periodic motion that is quasistationary when timeaveraged over a wave period. The turbulent motions occur on a much smaller time scale and thus can be treated as dissipative forces represented by turbulent friction terms which include the eddy viscosity. Therefore, we can use the Navier Stokes' equation to describe the flow in the surf zone:
U +V +p+gz (VxU) xU +UV2xU (5.1)
at 2 p
where p is the kinematic viscosity coefficient. Let the surface displacement and the velocity vector, U(u,v,w), be decomposed into mean value, wave and turbulent fluctuations, which are distinguished by subscripts c and w, and prime respectively, thus,
7 = + 77' 17 +7,v + 17' (5.2)
U = U+U' UC+Uw +U' (5.3)
where the superscript is used to denote turbulent averaging. After turbulent averaging, Eq. (5.1) becomes
a' +V +P+g (VxU)xU + vV2xU (5.4)
at 2 p
where v is the total viscosity including the eddy viscosity due to turbulence. The subscript is omitted thereafter.
Taking the scalar product of U(u,v,w) and the respective terms in Eq. (5.4), and summing the products gives the energy equation with a dissipation term:
2 [ p
h I h
by applying the kinematic boundary conditions at the free surface and the bottom,
w 7 7 UVh7 =0 at
w h + UVhh = 0
and setting p equal to zero at the surface, we obtain,
r7 2 7 22
(q )dz+ gh + Vhf[U( + +gz)]dzf vUV2Udz=O (5.6)
8t 2 8t 2 p
h h h
using the Leibnitz' rule of integration. This equation is similar to the wave energy equation given by Eq. (2.26) with the additional dissipative term. This volumetric dissipative term is replaced here by equivalent energy flux term introducing a head loss term in the context of Bernoulli sum, i.e., 17 17
D = f vUV2Udz =Vh fglUdz (5.7)
h h
where 1 is defined as the head loss due to turbulence. Eq. (5.6) can then be expressed as 17 2 r/ 2
S(q )dz + gha + Vhf [U(L + +gz +gl)]dz = 0 (5.8)
t 2 at h 2 p
h h
Following the same procedure used in Chapter 3 by taking time average over the wave period, an energy equation similar to that of Eq. (3.38) can be obtained, BE
 +Vh(UhE) =0 (5.9)
at
65
in which the transport velocity, Uh, is the counterpart of (Cg+UJ) in the non dissipative case and can be represented by
Uh = Cg +U + CgD (5.10)
The first two terms in (5.10) constitute the transport velocity due to nondissipative forces, whereas the last term manifests the effect due to the dissipative forces. The latter, in general, should be negative indicating a reduced energy flux due to dissipation. In theory, it can be estimated from the timeaveraged energy dissipation term as follows: Vh.(CgDE) =D (5.11)
where D is equal to the time average of D.
5.2 Wave Action Equation in the Surf Zone In this section, the wave action equation given in Eq. (2.41) is shown to be also valid even in the presence of strong turbulence. It is assumed here that the dynamic free surface boundary condition given in Eq. (2.10) is still valid with the inclusion of a head loss term. Following the approach by Kirby (1983), a virtual work term proportional to W(Dqw/Dt) is introduced to represent the head loss, where W is a positive undefined coefficient indicating the strength of the dissipation. The dynamic free surface boundary condition then becomes
7 (1 + W) 0 +Vh95'Vh'Ow (5.12)
g
which, after substituting Eqs. (2.12) and (2.13) into it yields aei_ (1 + [adAei+ iei{ A +js.VA =0}]
g at
66
where ad is again the intrinsic frequency including the diffraction effects. Following the same reasoning as in Chapter 2 we obtain
a
A g9 where aD=(1 +W)(aUS'k). (5.13)
a D
a +_s.VA =0 (5.14)
at S
The subscript, s, denotes the mean water surface level. The kinematic free surface boundary condition also yields
2
UoD = (1 + W)gk tanhk(h + 7c) (5.15)
a +V' (Ua)=0 (5.16)
a~t
Combining Eq.(5.14) and Eq. (5.16), we obtain the following wave action equation which is valid in the surf zone: g H2)=0 (5.17)
at 8OaD 8Sgy
This equation enables us to estimate the surface velocity in the surf zone once the wave height decay rate can be established.
5.3 Wave Height Transformation in the surf zone In this section an approach developed by Lee and Wang (1993) to wave breaking and wavecurrent interaction in the surf zone is presented.
The wave energy equation given in Eq. (5.9), when expressed in terms of wave height, can be written as,
a o H2 o H2
(pg )+Vh1(Cg+U+CgD)pg ] =0 (5.18)
at UD 8 OD 8
Assuming that the surf zone retains a quasisteady state when integrated over a wave period, then the slowly varying flow properties become time independent, and the absolute frequency becomes constant. Therefore, Eq. (5.18) becomes pg H2]
Vh [(Cg +U+CgD) p =0 (5.19)
8 OD
Now we apply the wave action equation in the steady state, and the wave energy equation in the surf zone is reduced to gH2]
Vh' [(Cg +Cg D) ] =0 (5.20)
8 OD
We will assume that the dissipative term in equation (5.21) is proportional to the group velocity at the breaking point, Cgb, CgD = 6Cgb (5.21)
where f is a positive coefficient. Substituting in Eq. (5.20)
 pg H2]0 (5.22)
8 OaD
where the Cg* is the group velocity in the surf zone and is calculated by
Cg* =Cg /JCgb (5.23)
This energy transformation model given in its final form by Eq. (5.22) has only one unknown coefficient, 8.. The model is applicable to three dimensional bathymetry and any incident wave angles.
In the case of twodimensional bathymetry of uniform slope, it is possible to obtain an analytical expression that describes the wave height transformation in the surf zone. Considering that x is normal to the beach and directed onshore, Eq. (5.22) becomes a Cg1 Xgb) =0 (5.24)
The x component of the above equation results in, H2 ~
cosH(,oCgb Cg) (5.25)
Using the dynamic free surface boundary condition given in Eq. (5.13),which states that UD is proportional to the wave height in the surf zone, equation (5.25) can be written as
H & (5.26)
cosO(/JCgb Cg)
with ,H determined later. The wave angle 0 can be obtained using Snell's law as
0 sin1 (C sin 0/Co)
(5.27)
69
where Ca,=&
Hb HbCgb csO(,oCg)
where Cg'=Cg/Cgb. After applying boundary conditions at breaking point, H/fl =1 and d/db= 1, and H/Hb=H', where Cg becomes zero, we obtain H H'
H oHs (5.29)
Hb cosl[1Cg'(1H's/cos~b)]
Now, using a shallow water approximation to Cg',Eq. (5.29) yields H H's
Hb cosO[1 Vd(l H's/cosOb)] (5.30)
Where d'=d/db. Only one of the three parameters H's, PH and /3 is independent. If the value of H', is obtained from an experiment, the other two are solved by, cos9b
16= COS 0b
cos0b H's (5.31)
,PH COS ObH's HCgb (5.32)
cosOb H's
Lee and Wang (1993) showed good agreement between the present theory and the laboratory experiments by Horikawa and Kuo (1966). The values of H', are evaluated from the data, they depend on the slope of the beach face and can be closely approximated by (tan a)'/2
70
where a is the slope of the beach. This idea coincides with the fact that wave height decay across the surf zone is strongly influenced by the beach slope.
In the numerical wave propagation model in which this theory is to be included Eq. (5.22) will be used since full three dimensional cases will be solved. The model treats the value of C9C9b as the new group velocity inside the surf zone and solves the wave propagation equation using it. Figures 5.2 to 5.3 compare the numerical results with the theoretical curves obtained with Eq. (5.30) for a I to 20 slope with H,'=0.22, that is, '8= 1.28.
The figures show how the numerical results oscillate around the analytical solution, and the amplitude of the oscillations increases with increasing wave period. The author has compared the solutions from other numerical wave propagation models that include this surf zone theory and the deviations between the analytical and numerical solutions were smaller. Therefore it is reasonable to conclude that the fluctuation of the wave height after breaking is entirely due to the numerical scheme used.
BREAKING TEST 1
d/db
BREAKING TEST 2
d/db
BREAKING TEST 3
0 0.2 0.4 0.6 0.8 1
d/db
Figure 5.1: Comparison between breaking analytical and numerical results. Test 1: T=3 s, Test 2: T=5 s, Test 3: T=7 s.
CHAPTER 6
MODEL APPLICATIONS
In the previous chapters a numerical model describing wave propagation in the nearshore zone has been developed and tested against several analytical solutions. The model will be now used to simulate the wave field in various cases where experimental data are available in the literature.
The first section of this chapter evaluates the performance of the model for the case when waves propagate over a paraboloidal shoal surrounded by a region of constant water depth. In the second section, the model is compared with measurements of wave height obtained in a laboratory experiment simulating the nearshore wave field around a detached breakwater on a sloping beach. The third section concerns the wave field in the vicinity of ajetty placed perpendicular to the shoreline, again against data from a laboratory experiment.
In the last section, the model application is illustrated by a field of groins in constant and sloping bottoms. It will be shown how the model can simulate the wave field around such structures allowing for different dimensions and permeabilities.
6.1 Wave Field over a Paraboloidal Shoal
Since the conventional wave refraction theory is based on the geometrical optic approximation, it fails to predict wave height at and near caustics where ray intersections occur. The effect of diffraction should be included in the analysis of waves in the vicinity of ray convergence. The basic equation for the present model includes this effect. As a demostration, the model is applied to wave propagation over a submerged shoal with
73
concentric circular contours where the conventional refraction theory indicates the formation of caustics. Hydraulic model experiments for this situation were conducted by Ito and Tanimoto (1972). All tests were conducted for non breaking waves. The arrangement of the shoal in the numerical computation is shown in Figure 6.1.
Cross Section
y/Lo Plan View
81
0 1 2 3 4 5 6
yILo
Figure 6.1: Bathymetry for paraboloidal shoal
configuration, after Ito and Tanimoto (1972).
The water depth around the shoal is constant with a value of 15 cm, the top of the shoal is
at a depth of 5 cm. The initial wave length on the upwave row of the numerical is 40 cm, the
rest of the input conditions are given as follows
NX NY Ax (m) Ay (m) Hi (m) T (sec) 0 (o) 63 17 .05 .15 .02 .5107 0
The results of the numerical calculation and the hydraulic model tests are presented in Figure 6.2. Figure 6.3 shows three dimensional plots of both free surface and wave height. It clearly shows how the wave fronts bend and wrap around the shoal. The maximum wave height is generated right behind the shoal, where the wave rays would be crossing according to conventional wave refraction theory.
3
2.5
. 2
1.5
0.5
0
3
2.5
..2
1.5
0.5
1 2 3 4 5 6 7
x/Lo
0 2 3 5 6
y/Lo
3
2.5
2 Sec. 3 0/'
~1.50
1 0
0.5 1 1 1
0 1 2 3 4 5 6
y/Lo
Figure 6.2: Wave height comparisons between computed and experimental results by Ito and Tanimoto (1972).
S Numerical Model
O000Measured 0 0 %
See. 1
Sec. 2 0
0o
75
FREE SURFACE (ETA/ETAi)
x/Lo
WAVE HEIGHT (HIHi)
6 0
x/Lo
Figure 6.3: Numerical (1972).
model results for shoal configuration used by Ito and Tanimoto
y/Lo
y/Lo
76
6.2 Wave Field around a Detached Breakwater A detached breakwater is also a common coastal structure. The numerical model is tested here to predict the wave field, breaking and surfzone decay around one of these structures.
Watanabe and Maruyama (1986) performed a hydraulic model test that included a detached breakwater made of a steel plate 2.67 m long and placed at a water depth of 6.0 cm, and parallel to the shoreline on a plane beach made of mortar with a slope of 1150. The breakwater was designed to almost perfectly reflect the incident waves. The waves were incident normal to the breakwater and the contour lines, and the deepwater wave height and period were 2.0 cm and 1.2 sec., respectively. The depth at the offshore grid row of the numerical domain was 12 cm, the rest of the input data were given as
NX NY Ax (m) Ay (m) Hi (m) T (sec) 0) 60 60 0.10 3814 .0202 1.2 0
The results of the numerical calculation for both free surface and wave height are presented in Figure 6.4. The location of the computed breaker line is compared with that obtained in the experiment in Figure 6.5. The overall agreement seems to be fairly good, with the exception of the area located in the middle, behind the breakwater, where the predicted breaker line is closer to shore than the measurement. Two explanations can be given for this difference. One related to the fact that the computed results do not include nearshore circulation, which in this area would be directed offshore against the incoming waves, thus, the waves break earlier. A second reason could be related to the model's partial inability to fully simulate waves that do not propagate mainly in the x direction. Comparisons of crossshore distributions of wave height are shown in Figure 6.6. The
77
results compare reasonably well, except the area affected by the reflected waves from the breakwater. As explained in the last chapter, the model is not capable of simulating diffracted waves propagating in the negative x direction. Although it would be possible to include the reflected wave field from the breakwater, in the present model is not included since usually our main interest for a detached breakwater is the wave condition on the downwave side.
6.3 Wave Field around a Groins and Jetties
Groins and jetties, in general, are shore perpendicular structures used on the coast for multiple purposes. Sometimes jetties are constructed to keep a tidal inlet in place so as to stabilize the inlet entrance. An important use of groins is as terminal structures for beach nourishment projects, with the idea of maintaining a wide beach for a longer time by holding the sand in place. They might also be used to simply slow down erosion along a natural shoreline. In any of these cases, to accurately predict the effects of these structures on the shoreline, the first step is to predict the wave field around them.
Watanabe and Maruyama (1986) constructed a hydraulic model to test this situation, the basic layout of the test was the same as the one described in section 6.2, but instead of a detached breakwater, the experiment set up included a jetty made of steel plate extended to an offshore distance of 4 in, where the water depth was 8 cm. The deep water incident wave angle is also different, 60' for this case. After applying Snell's law and calculating the refraction and shoaling coefficients between deep water and the offshore grid row at a water depth of 12 cm, the rest of the input conditions were given as
AW NY & (m) &y (m) Hi (m) T (sec) 0
60 60 0.1 0.4 .0143 1.2 28.3
78
FREE SURFACE ETA
3 2 1 0 1 2 3 4
Distance from center of breakwater y (m)
4 3 2 1 0 2 3 4
Distance from center of breakwater y (m)
Figure 6.4: Numerical model results for detached breakwater test.
WAVE HEIGHT H (cm)
4 3.5 3 2.5 2 1.5 1 0.5 0 0.5 1 1.5 2 2.5 3 3.5 4
Distance from center of breakwater y (m) Figure 6.5: Comparison of location of breaker line behind a detached breakwater, experimental data after Watanabe and
Maruyama (1986).
5 x m Nueric M o d e l 1 i
4 x=4 m Numerical Model
3 0 0 0 0 0 Measured
20 00P
 0Po 0o 00
. . .CO o .. .
4 3.5 3 2.5 2 1.5 1 0.5 0 0.5 1 1.5 2
Distance from center of breakwater y (m)
2.5 3 3.5 4
3 3.5
Distance onshore x (m)
5 I
4 y=2 m
3 000 00..0 _.0 0
"I:2 O00000000000 00 0
0
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
Distance onshore x (m)
Figure 6.6: Comparison of distributions of wave height around a detached breakwater, experimental data after Watanabe and Maruyama (1986).
80
The computed results for both free surface and wave height are presented in Figure 6.7. The location of the computed breaker line is compared with that obtained in the experiment in Figure 6.8. The overall agreement is acceptable, although there is some discrepancy, in particular, near the jetty on both sides. Part of this deviation is attributable to the effect of the nearshore current which has not been incorporated here. Figure 6.9 shows comparisons of crossshore distributions of wave height. They reveal that the model predicts with sufficient accuracy the wave field around the jetty, especially taking into account that the model is based on linear theory.
One of the limitations of the model can be identified studying Figure 6.9 in more detail; the wave height along y = 4.8 m does not seem to decrease as much as the data when approaching the shoreline at x = 6 m. This is thought to be due to two different reasons: (1) the breaking criteria is simple and probably inaccurate when more than one wave field is superimposed. Watanabe, Hara, and Horikawa (1984) presented a new breaking condition and diagrams in which the main governing parameter is the ratio of orbital velocity at the crest to the phase velocity, in place of the conventional depth to height ratio. Under this kind of breaking criteria, waves would break in cases where the conventional depth limiting breaking criteria has not been reached yet; (2) the surf zone model used, which affects the wave height reduction rate after breaking might not be adequate for the case when incident and reflected waves are combined.
81
FREE SURFACE (ETA)
0
1
4
5
6
0
8 9 10
WAVE HEIGHT H (cm)
6 1 I I I I 3 I I
0 1 2 3 4 5 6 7
Distance alongshore y (m)
Figure 6.7: Numerical model results for jetty test.
8 9 10
1 2 3 4 5 6 7
Distance alongshore y (m)
82
1 i I I I
Wave Direction
1.5
1.5 Numnerical Model
2 O O O O Measured
2.5
s 3
3.5 .
4 O
0 0 0 0 0 00
5 0 0 0 0o
40000 0 0 0 (0
5.5
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 Distance alongshore y (m)
Figure 6.8: Comparison of location of breaker line around a jetty, experimental data after Watanabe and Maruyama (1986).
5 i I 1 i I I
4 y=5.2 m Numerical Model
0 0 0 0 Measured
3
S2 000000 00no
i nO0000000000000000 0
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 Distance onshore x (m)
5
4 y=4.8 m
S
0 II I I
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 Distance onshore x (m)
5 i i i i i I
4 y=7.0 m
3 0000
2 ~ ~~ 0 0 0nannpoonoo
2 00 00000000000 0
1
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 Distance onshore x (m) Figure 6.9: Comparison of crossshore distributions of wave height around a jetty, experimental data after Watanabe and Maruyama (1986).
83
6.4 Wave Field around Permeable Groins The main objective of this thesis was to obtain a wave propagation model that could simulate how a field of permeable groins alters the wave propagation in the nearshore region. The basic numerical solution and appropriate boundary conditions were described in Chapter 5. In this section it will be shown how it is possible to assign different physical characteristics to a groin and obtain the corresponding wave field
6.4.1 Single Groin in Constant Water Depth
To first study how the characteristics of the groin affect the amount of energy that is reflected and transmitted, and thus the wave field around the structure, a simple case is presented. The model was run over a domain with dimensions 600x900 meters and a 400 mn long groin located at y=450 m. The rest of the necessary input conditions are given as,
NrX NY Ax (m) Ay (m) Depth (in) Hi (m) T (sec) 0 (') 60 60 10.0 15.0 5.0 1.0 7.0 20.0
Several cases were run for this configuration with different values assigned to the physical properties of the groin. The following table includes the values used as well as the reflection and transmission coefficients calculated by the model and employed in the boundary conditions around the structure, Case 1 Case 2 Case 3 Case 4 Inertial Coefficient, S 1.0 1.0 1.0 1.0
Friction Coefficient, f 0.0 1.0 0.5 0.01
Permeability, e 0.0 0.1 0.5 10.95
Reflection Coefficient, R 1.00 0.66 0.14 0.05Transmission Coefficient, T 0.00 0.12 0.36 0.96
84
The model results are presented in Figures 6.10 to 6.13, the contour lines represent the instant water free surface and show the bending of the wave crests behind the groin as well as the formation of short crested waves perpendicular to the groin in the upwave area. The relationship between the calculated and the incident wave height (H/Ho) is shown in shades of gray with a scale on the right hand side of each figure. Figure 6.14 presents the wave height change along four sections parallel to the groin, the first two at 6.5 and 6.5 meters and the others at 67.5 and 67.5 meters. It can be seen that H/Ho always increases for increasing R values in the upwave side of the groin (negative y). However in the lee of the groin H/Ho does not seem to always decrease with decreasing transmission coefficient. For Case 3, which has a larger transmission coefficient than Case 1, H/Ho values at y = 6.5 are greater initially than the values from Case I near the head of the groin, but become smaller towards the root of the groin. This effect is thought to be due to the superposition of transmitted and diffracted wave fields on the lee side of the groin. At y = 67.5 m, further downwave where the diffraction effect becomes less important, H/Ho for Case 3 is clearly greater than that for Case 1. Case 4 shows the results when the groin is close to being almost 'transparent' to wave propagation.
The special case of groins composed of pile screens can also be simulated approximately by choosing the right physical characteristics for the groin. Here, an inertial coefficient, S, of around 2 would seem to be more suitable. The permeability will also be considerably higher than that of a rubble mound groin, with values usually around 0.8 to 0.9. The friction factor on the other hand will be reduced. Ideally f should be a value dependent upon the wave field also, thus, successive iterations would lead to the rightf value based on the Lorentz' s equivalent work concept..' In this model, for the sake of computational time, such dependency is not included, andf is a given input.
30( 40(
200 150 100 50 0 50 100 150 200
Distance from center of jetty y (in) Figure 6. 10: Free surface and amplification factor around one groin, Case 1.
HO
300 ~350
400
200 150 100 50 0 50 100 150 200
Distance from center of jetty y (mn)
Figure 6.11: Free surface and amplification factor around one groin, Case 2.
H/Ho
300 35C
400 450 500
600 0' _........ _0
200 150 100 50 0 50 100 150 200
Distance from center of jetty y (m)
Figure 6.12: Free surface and amplification factor around one groin, Case 3.
H/Ho
2
200
1.8
1.6
300 1.4
5 350 1.2
4001
4500.
0.6
0.2
550 0.2
600 0
200 150 100 50 0 50 100 150 200
Distance from center of jetty y (m)
Figure 6.13: Free surface and amplification factor around one groin, Case 4.
2.5 2.25
2
1.75
o 1.5 1.25
1
0.75 0.5
0.25
i i i Ir
100 150 200 250 300 350 400 450 500 550 600
x (M)
00 150 200 250 300 350 400 450 500 550 600 x (M)
2.5 2.25
2
1.75 o 1.5 1.25
1
0.75 0.5
100 150 200 250 300 350 400 450 500 550 600
x (M)
2.5 2.25
2
1.75 1.5
1.25
0.75
0.5 y=67.5 m
0.25
100 150 200 250 300 350 400 450 500 550 600
x (M)
Figure 6.14: Amplification factor along four sections parallel to groin.
 Case 1 ........ Case 2
. .Case 4
y=67.5 m
 Case 1 Case 2
 . Case 3 .. Case 4
y7.5m
. .. .... . . .. 
6.4.2 Single Groin on a Plane Beach
For the case of a varying bottom topography the computation of reflection and transmission coefficients for a porous groin becomes more complicated, as transmission and reflection coefficients are functions of wave properties which, in this case, change along the length of the structure. At first it was decided not to include these changes in the model. However, after further consideration it was decided to include the complete solution to compute reflection and transmission coefficients to examine the effect. A case was run with constant reflection and transmission coefficients based on the groin characteristics and angle of incidence at the offshore tip of the groin. The groin was located in plane beach with a 1 to 50 slope. Then, the full wavedependent solutions of transmission and reflection coefficients were used. The results are compared in Figure 6.15. The table below gives the computed TIR values at the tip and the base of the groin for the case of 6i = 20'.
S f~e 6 Oi Obas Rtp Tti, Rbase Tbs
2.0 j1.0 0. 0 16.55 2.68 0.24 0.07 0.80 00
It can be seen that the reflection coefficient increases from 0.26 at the offshore end of the groin to 0.80 at the its base; this effect will always occur when groins are placed across a surfzone, thus, it seems important to compute R and T coefficients based on a porous flow model rather than just using constant values.
Using the full solution to compute R and T coefficients not only accounts for changes in their value along the groin, but also simplifies the input given to the model. Figure 6.16 compares the cases when same R and T coefficients are used for two different offshore wave angles and when the model computes those coefficients depending on the value of Oi. The
89
fixed coefficients correspond to an offshore angle, 0, of 20 degrees, and the variable values
to 0i=5.
1.75 Fixed R,T
1.5 .... Variabe R,T
1.25
0.75
0.5 Front side of groin
0.25
00 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo
2
1.75 Fixed R,T
1.5 ....... Variabe R,T
1.25
0 .75 ... . . . . . .
0.5 Lee side of groin
0.25'0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo
Figure 6.15: Comparison of model results for the cases of fixed and
variable reflection/transmission coefficients.
2
1.75 Fixed R,T
1.5 ....... Variabe R,T
1.25
0.75
0.5 Front side of groin
0.25
00 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo
2
1.75 Fixed R,T
1.5 ... Variabe R,T
1.25
0
0.75
0.5 Lee side of groin
0.25
00 065 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo
Figure 6.16: Offshore incident wave angle influence on reflection and
transmission coefficients.
6.4.3 Waves around a Groin Field
Groin fields are used frequently as coastal defense systems. Their performance will depend very much on their ability to disrupt both wave energy propagation and longshore currents. The applicability of the model to this kind of situation will be shown in the following examples. Here, three groins with arbitrary dimensions are introduced. An unlimited number of groins can be introduced in the computational domain as long as physically permissible. The input wave conditions, as well as the computational domain dimensions are given as follows,
NX NY Ax (m) Ay (m) Slope Hi (m) T (sec) 0
60 70 5.0 15.0 1/60 1.0 7.0 20.0
The three groins are located at y 225, 525 and 825 meters, with lengths 100, 150 and 125 meters respectively. The dimensions are definitely larger than what is usually constructed on real coastlines, such that the effects of groins are exaggerated to aid in visual appreciation in the figures. Anyway, structures up to 300 meters long have been used in places like the Mediterranean Spanish coast. Figure 6.17 shows the results when the groins are completely impermeable, referred as Field 1. On Figure 6.18 the case when the three groins have different physical characteristics, Field 2, is presented. The characteristics of the groins for this case are given as:
S f e
Groin at y=225 m 1.0 1.0 0.4
Groin at y=525 m 1.0 1.0 0.2
Groin at y=825 m 1 1.0 2.0 0.1
