Interactions of breaking waves with a current over cut cells

Tingqiu Li
*
, Peter Troch, Julien De Rouck
Department of Civil Engineering, Ghent University, Technologiepark 904, B-9052 Zwijnaarde, Belgium
Received 4 November 2005; received in revised form 4 October 2006; accepted 5 October 2006
Available online 22 November 2006
Abstract
By design of the external and internal wavecurrent generators, the objective of this paper is to extend our ecient
NavierStokes solver [T. Li, P. Troch, J. De Rouck, Wave overtopping over a sea dike, J. Comput. Phys. 198 (2004)
686726] for modelling of interactions between breaking waves and a current over a cut-cell grid, based on a dynamic sub-
grid-scale (SGS) model. This solver is constructed by a novel VOF nite volume approach, coupled with surface tension.
When studying waves following a positive current, our external generator creates the combined inow motions of waves
and a current, which is viewed as one type of wavy inow conditions. For cases of waves against strong currents, our inter-
nal generator describes the opposing current, by incorporating the source function to the continuity and momentum equa-
tions as a net driving force, acting on the uid elements lying within the nite thickness source region. The outgoing waves
downstream are dissipated with a breaking-type wave absorber placed in the tank extremity. Five test cases recommended
are of distinctly dierent applications of interest, characterized by overtopping of following waves over sloping and vertical
structures.
Under the grid renement eects, the results in 2D and 3D are in close agreement with the experimental data available
in terms of the surface wave. Additionally, the performance of convergence in computations is also investigated, including
full discussion for waves on beaches between 2D and 3D. By visualization of the motions that describe the physics of tur-
bulence, it has been shown that our solver can capture most of the signicant features in wavecurrent interactions varying
with three dierent current speeds (positive, zero, negative).
2006 Elsevier Inc. All rights reserved.
Keywords: The NavierStokes solver; A numerical wavecurrent generator; Following or opposing waves; A breaking-type wave absorber
1. Introduction
Surface wave propagation in water of varying depth is of particular interest in the study of nonlinear wave-
wave interactions and breaking wave-induced motions, especially coupling of steep waves with a structure.
For the latter, the major physical processes involved are extremely complex, characterized by reection, shoal-
ing, diraction and breaking, until waves overtop. Based on uid mechanics of turbulence in waves, the
dynamics of these transformations can be considered as small-scale turbulence, larger-scale coherent vortical
0021-9991/$ - see front matter 2006 Elsevier Inc. All rights reserved.
doi:10.1016/j.jcp.2006.10.003
*
Corresponding author. Tel.: +32 92645489; fax: +32 92645837.
E-mail address: tingqiu.li@UGent.be (T. Li).
Journal of Computational Physics 223 (2007) 865897
www.elsevier.com/locate/jcp
motions, low-frequency waves and steady ows [2]. This subject itself, therefore, is one of the challenging top-
ics in research and engineering. Importantly, the additional eects of a current substantially aggravate the sit-
uation that may change the coastal environments. Typically, the wave-induced currents alter the nearshore
transport of sand and contaminants at the seabed, as the mechanism of sediment movement is closely related
with pickup by wave action and transport by currents. Moreover, the combination of extreme wave and cur-
rent conditions is also used to establish the design loads for large oshore structures.
The relevant study mainly involves computational modelling of moving boundary problems associated with
an airwater interface (referred to as the free surface), separating the two segregated uids. There are many
ways available in theory for solution of this particular problem, which basically may be classied into two
types of categories: one- and two-phase uid model. In the case of the former, the free surface is assigned
as one of boundaries, and a necessity is to specify boundary conditions at this surface before driving impul-
sively from rest; but conversely for the latter the free surface is viewed as a material interface of a contact dis-
continuity in variable density elds, captured naturally without special provision as part of the solution.
Denitely, both methods are of their essential features. Given the exact location of interface at transient state,
for example, a high order of accuracy is easily preserved with gravity-dominated one-phase solver [32], as in
this case a non-smeared interface can be materialized. Additionally, computations are less expensive than a
two-phase model, which needs the considerable CPU time, because of the extra air domain involved. Never-
theless, it indeed oers the potential to situations where an air bubble is trapped in the ambient water [25,38]
and can well describe the behavior of the waterair mixture. Owing to the strong interface deformation and air
entrainment phenomenon involved, the question raised here is relevant to numerical diusion, primarily
caused by handling a sharp jump of the density over a narrow band relying on the fact: the density and vis-
cosity at a transitional layer between the two uids vary appreciably across the interface. If such dissipation
accumulates with time, this behavior could induce an impact on the accuracy and stability. Consequently, for
high density ratio cases convergence may be slow, and the modelling of surface tension may lead to accuracy
problems [36].
An earlier approach is based on the linear boundary layer approximation [20], when studying energy trans-
fer between waves and a current. Until recently, the most popular one is to introduce the radiation stress con-
cept [37], for instance, as the basis of the famous SWAN solver [4]. Within the framework of this technique,
the wave motions are distinguished and are in fact separated from the current motions. An alternative is to
account for surface-piercing numerical wavemakers for forcing surface waves. The emphasis is on modelling
of external or internal wavemakers plus a current, by prescribing inow conditions or adding a net driving
force. This provides a convenient framework point for investigation of wavecurrent interactions in a consis-
tent manner. More specically, solution of the NavierStokes equations can directly give the current state and
describe the wave pattern simultaneously. Importantly, distinction between waves and currents is unnecessary
in this theory, as done in Park et al. [43]. It is one of the very few NavierStokes simulations even in the case of
non-breaking waves. For situations where the positive (or negative) currents vary gradually in space and time,
a key point in computations is that mass conservation has to be held over a long duration simulation, as waves
develop rapidly in both scales. Accordingly, the numerical schemes proposed would identify the basic mechan-
ical design feature. To avoid waves reected in the near-inlet regions, one suggestion is to discretize a kind of
Orlanski condition on a grid for the velocity and total wave proles [30] or to apply for the internal wave-
maker [35]. In view of the computational issues mentioned above, therefore, a numerical wavemaker has ben-
eted greatly, especially it substantially diers from the radiation stress approaches.
In this paper, we mainly emphasize the design of an external and internal generator for wavecurrent cou-
pling, by using a dynamic Smagorinsky model. Theoretically, an external generator can be viewed as one type
of inow conditions into the major computational domain [1,48], which perturbs ows through the boundary
rather than the interior. It is a major dierence with an internal generator, in which a net driving force is added
to the source term of the governing equations for creation of the desirable adverse current. Both numerical
generators provide an ecient tool for the study of coupling of breaking waves with a current over a sloping
and vertical structure both in 2D and 3D, as addressed in waves following or opposing a current (mostly for
the former). This considerably extends our previous study [30], in which only a pure wave in 2D is investigated
with our developed NavierStokes solver, coupled to a static Smagorinsky model. The outgoing waves down-
stream are dissipated with the breaking-type wave absorber placed in the tank extremity.
866 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
Five test cases without the exact solution are applied in this study, some of which are related to dierent
applications in coastal and harbour engineering. Three cases rst recommended involve the modelling of fol-
lowing (or opposing) waves with various wave characteristics (see Table 1) in the open channel ow over con-
stant water depth, due to measurements available in the literature. In all test cases, particular interest is paid to
the near-boundary free-surface turbulent region, especially when validated our design wavecurrent genera-
tors. The last two cases are characteristic of waves breaking in water of varying and constant depth, respec-
tively, where steep waves overtop over the sea dike and vertical wall. Our results demonstrate that most of
typical features in the wave-induced motions are captured with our solver in the presence of a current or
not. Comparisons with the experimental data are also rather satisfactory in terms of the mean surface and
horizontal velocity proles, including prediction of the long-term propagation behavior of waves.
The outline of the present study is as follows. We rst describe our numerical methods, e.g. by introducing a
dynamic SGS model for the study of breaking water-wave problems, with particular emphasis on design of the
wavecurrent generator plus the wave absorber. A comparison of the calculated results in 2D and 3D with the
experimental data is represented next, followed by the concluding remarks.
2. Methodology
2.1. A MILES model
On the assumption that waves are travelling on a current in a numerical wave tank (NWT) that is of a rel-
atively long ume (see Fig. 1), the plane of z = 0 corresponds to the mean water level in the Cartesian coor-
dinate (x, y, z) system, where x is in the direction of the wave propagation, while (y, z) are toward the lateral
direction and positive upwards, respectively. Based on MILES, the mathematical model in this study is gov-
erned by the unsteady incompressible NavierStokes (NS) equations, as the concept of MILES (monotoni-
cally integrated large-eddy simulation) attempts to resolve accurately the NS equations. Namely
ou
ot

o F
i
F
v
F
a

ox

o G
i
G
v
G
a

oy

o H
i
H
v
H
a

oz
Q; 1
where the variables u = (0, u, v, w, a)
T
, in which (u, v, w) are the components of the velocity in the x-, y- and
z-directions, respectively, and a is the volume fractions in VOF (volume-of-uid). By dening the inviscid
uxes (F
i
, G
i
, H
i
), the viscous uxes (F
v
, G
v
, H
v
) and the other uxes (F
a
, G
a
, H
a
), we have
Table 1
The wave and current characteristics for Case 1 to Case 5: the wave height H, the period T
a
, the water depth d and a current speed U
c
Test cases U
c
(ms
1
) d (m) H (m) T
a
(s)
Case 1 (WCA4): Kemp & Simons (UK) 0.11, 0.0, 0.11 0.2 0.0394 1.0
Case 2: Klopman (Netherlands) 0.16, 0.0, 0.16 0.5 0.12 1.44
Case 3 (WC1): Fredse et al. (Denmark) 0.222, 0.0, 0.222 0.45 0.13 2.5
Case 4: a vertical wall problem 0.0, 0.10 0.7 0.16 2.0
Case 5: a seadike problem 0.0, 0.17 0.7 0.16 2.0
x
y
z
0
wave alone
direction
(inflow)
U
c
opposing current
following current
U
c
source function
absorber
outflow
a
x
z
0
following current
U
c
absorber
source
opposing current
U
c
b
Fig. 1. Denition sketch of the NWT for waves following or opposing a current. (a) 3D; (b) 2D (on the symmetry plane, y = 0), where the
origin is upstream.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 867
F
i

0
u
2
vu
wu
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; G
i

0
uv
v
2
wv
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; H
i

0
uw
vw
w
2
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; 2
F
v

0
s
xx
s
yx
s
zx
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; G
v

0
s
xy
s
yy
s
zy
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; H
v

0
s
xz
s
yz
s
zz
0
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; 3
F
a

u
1
q
p
0
0
ua
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; G
a

v
0
1
q
p
0
va
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; H
a

w
0
0
1
q
p
wa
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; 4
plus the source terms Q
Q
f
1
q
F
x
b
F
x
d
fu
1
q
F
y
b
F
y
d
fv
1
q
F
z
b
F
z
d
fw g
a
ou
ox

ov
oy

ow
oz
_ _
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
; 5
where g is the gravitational acceleration and f is the internal source function (see Eq. (50)) for creation of an
adverse current. (F
x
b
; F
y
b
; F
z
b
and (F
x
d
; F
y
d
; F
z
d
are three components of a body force
~
F
b
[30] and an extra dis-
sipative term
~
F
d
(see Eq. (52)) in the (x, y, z)-directions, respectively. For the latter, it is active on a limited
portion of regions (i.e. the so-called numerical damping zone, see Fig. 7) adjacent to the articial boundary.
For restricted domain simulations, this provides an alternative for dissipating energy of waves in numerical
models. p is the modied pressure that considers the eect of the subgrid kinetic energy k (from the trace
of s
ij
) by setting p P
2
3
k, in which P is the total pressure.
With the Einstein summation convention (such as s
xx
= s
11
), the viscous stress tensors s
ij
(i, j = 1, 2, 3) may
be evaluated by
s
ij
m
ou
i
ox
j

ou
j
ox
i
_ _
6
with the molecular viscosity m (a constant of 1.3 10
4
m
2
s
1
). The local density q and the viscosity m are lin-
early connected with the scalar value a
q aq
w
1 aq
a
; m am
w
1 am
a
; 7
where the subscripts (w, a) denote the water and air, respectively.
2.2. Solution procedures
Applying Gausss divergence theorem to Eq. (1), this yields the integral form as
o
ot
_
V
udV
_
S
F dS
_
V
QdV 8
868 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
for the surface S surrounding the domain of interest V. Accordingly, the discretization at the (n + 1)th time
level is driven by
ou
n1
ot

1
V

faces

F
n1
S Q
n1
9
over each cell of a hexahedron.

represents the summation of all cell faces on a grid, and

F includes the
inviscid, viscous and other uxes via

F F n F
i
F
v
F
a
n
x
G
i
G
v
G
a
n
y
H
i
H
v
H
a
n
z
10
with unit normal surface vector n, where (n
x
, n
y
, n
z
) are its three normal components outwards in the x-, y- and
z-directions, respectively. (F
i
, F
v
, F
a
) to (H
i
, H
v
, H
a
) are the uxes dened by Eqs. (2)(4).
By introducing the delta form du = u
n+1
u
n
for the velocity, the uxes

F
n1
in Eq. (9) can be linearized
locally with respect to time t, when dropping the high-order terms. Namely

F
n1

F
n

F
ou
du 11
only for the convective and diusion terms, resulting in Eq. (9) to be factored into the three one-dimensional
equations as follows:
A
i1
du

i1
A
1
p
du

A
i1
du

i1
DtR;
A
j1
du

j1
A
2
p
du

A
j1
du

j1
du

;
A
k1
du
k1
A
3
p
du A
k1
du
k1
du

;
12
based on the alternating direction implicit (ADI) approximate factorization with a local time step Dt. The sub-
scripts (i 1, j 1, k 1) represent two adjacent neighbours of cell (i), respectively, and the coecients
(A
i1
, A
j1
, A
k1
) typically involve a blend of the inviscid and viscous volumetric uxes. R in Eq. (12) is the
residual of the momentum equations dened as
R
1
V

faces
F
n
i
F
n
v
_ _
Q
n
_ _
13
with F
i
F
i
n
x
G
i
n
y
H
i
n
z
and F
v
F
v
n
x
G
v
n
y
H
v
n
z
, where the source terms Q
n+1
= Q
n
= (0, 0, 0, g)
T
,
indicating that the body force
~
F
b
, the dissipative term
~
F
d
, the source function f(x, y, z, t) and the other uxes
are retained in the explicit formulation during the implicit process.
When the contribution of the dissipative component F
x
d
and the source function f is incorporated into the
temporal velocity u (for example, for the momentum equations in the x-direction)
~u u
n
du Dt
1
q
F
x
b
F
x
d
fu
_ _
14
this results in the Poisson equation for the pressure as
r
rp
n1
q
n1

1
Dt
fr
~
~u f g 15
according to the projection algorithm. In our implicit stage, the convective uxes F
n
i
(see Eq. (13)) at the nth
time step are calculated by means of the ux-dierence splitting approach, yielding the uxes F
n
i
at the face
i
1
2
between cell (i) and its neighbour (i + 1)
F
n
i
_ _
i
1
2
u
n
Su
L
_ _
i
1
2
when u
n

i
1
2
> 0: 16
Under the nite volume approximation, two typical variables (u
n
, u
L
) are separately treated for maintaining
the second-order spatial accuracy, including enhancement of the numerical stability. First, the normal face
velocity u
n
(dened by u
n
= un
x
+ vn
y
+ wn
z
) is directly modied by our design combination of the second-
and fourth-order articial damping terms. Namely
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 869
u
n

i
1
2

1
2
u
n

i
u
n

i1
_ _

1
2qw
2
c
2
p
R
i
1
2
p
L
i
1
2
_ _
1

1
2qw
4
c
4
p
R
i
1
2
p
L
i
1
2
_ _
3
_ _
: 17
As expected, this correction helps to damp high frequency oscillations, since for problems of the convective
domination the linear interpolation of the face value results in an unbounded solution. Secondly, the relevant
face values (such as the left state u
L
) are estimated with the characteristic variables at adjacent cells, by incor-
porating the second-order essentially non-oscillatory (ENO) scheme for u
L
i
1
2
and a
i

1
2
L
, and the MUSCL
(monotone upstream-centered scheme for conservation laws) approach for p
i

1
2
R

3
and (p
i

1
2
L

3
, while
an one-order upwind discretization is used for p
i

1
2
R

1
and p
i

1
2
L

1
. Finally, Eq. (16) in the ENO context
is derived by
F
n
i

i
1
2

u
n
S
1
2
u
i
u
i1

1
8
o
2
u
ox
2
Dx
2
if ju
i
1 u
i
j 6 u
i
u
i1
j;
u
n
S
1
2
3u
i
u
i1

3
8
o
2
u
ox
2
Dx
2
otherwise;
_
18
in smooth regions. Based on Taylor series, the dierent leading-order truncation errors (
1
8
o
2
u
ox
2
Dx
2
,
3
8
o
2
u
ox
2
Dx
2
) are
obtained, dependent on the upwinding of the stencils. For Eq. (17), our damping terms with regard to the pres-
sure can be expressed by

1
2qw
2
c
2
p
i1
p
i

1
12qw
4
c
4
p
i2
3p
i1
3p
i
p
i1

_ _
; 19
where the scaling factor c
4
(c
2
= 1 in this case) is determined by the three mainly diagonal coecients
(A
1
p
; A
2
p
; A
3
p
) of Eq. (12), which is formulated according to [30], in case the second- and fourth-order coef-
cients (w
2
, w
4
) are given adaptively with the local diusion velocity. Both schemes can be considered as
one type of shock-capturing schemes [16], which provide a particular form of the functional reconstruction
for the convective uxes: by using TVD (total variation diminishing)-based schemes as ux limiters that are
implicitly implemented SGS models. Denitely, such features essentially construct one type of MILES for
the spatial discretization of the convective terms. Hopefully, the scheme rendered prevents unphysical oscilla-
tions in computations.
More importantly, a hybrid approach constructed by the weighted upwind scheme and blending algorithm
is devised for high resolution of a step prole of a in VOF, as opposed to geometrical strategies. It is a simple
but eective algorithm in multidimensional directions, characterized by the interfaces that are arbitrarily ori-
entated with respect to the computational grids. To achieve higher-order accuracy, for example, the volume
fractions a
i

1
2
n
in VOF are updated by adding the second-order correction that accounts for the ux limiting
a
n
i
1
2
Ca
n
i
1 Ca
L
i
1
2
if a
n
i1
or a
n
i1
6 0 20
implying that a bounded solution can be guaranteed during tracking. In particular, the high-resolution mono-
tone scheme for a
i

1
2
L
may adjust amount of the numerical diusion induced by the low-order scheme (this
feature is also suitable to Eq. (17), as done in the ux-corrected transport algorithms, by introduction of a
ux-limiter C. For the fully implicit cell-staggered nite volume (FV) method on non-uniform stationary
cut-cell grids, we refer the reader to our study [30,33] for considerable detail of this procedure. The following
is to describe a dynamic subgrid-scale (SGS) model, because MILES coupled to such model can give a better
accuracy for particular periodic breaking water-wave propagation problems, as guided by the results of
Appendix A (e.g. Figs. 25 and 26).
2.3. A dynamic Smagorinskys model
In the conventional LES (large-eddy simulation), the SGS stress b
ij
caused by nonlinear terms u
i
u
j
is
expressed as
b
ij
u
i
u
j
u
i
u
j
21
in case a low-pass lter is applied to the NS equations. The overbar represents the spatially ltered quantities
at a characteristic length scale of the small eddies D, based on the predened lter kernel (e.g. a top-hat shaped
870 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
kernel). In the present study, we address the issue of modeling the total stress term (i.e. a purely dissipative
term). On the assumption that the SGS stress is restricted to eddy-viscosity branch of a family, theoretically,
it may be cast into the constitutive relation as follows:
b
ij

1
3
d
ij
b
kk
2m
T
S
ij
with S
ij

1
2
ou
i
ox
j

ou
j
ox
i
_ _
; 22
where d
ij
and S
ij
are the Kronecker delta function and the resolved-scale strain rate tensor, respectively. Phys-
ically, the major role of the SGS model is to remove energy from the resolved scales. In this sense, Eq. (6) can
be rewritten as
s
ij
m
eff
ou
i
ox
j

ou
j
ox
i
_ _
23
as the isotropic part of the SGS stresses (
1
3
d
ij
b
kk
) has been absorbed to the pressure term. An essential point in
our MILES, now, is to evaluate the eective viscous coecient m
e
, assigned with m
e
= m + m
T
, by introducing
a dynamic Smagorinskys model (DSM) that captures the averaged motions without elaborate numerical
computation.
Owing to the simplicity, a Smagorinsky-type isotropic eddy-viscosity model remains in wide use [28], as
small scales tend to be more isotropic than the large-scale uctuations, which are simulated directly. Theoret-
ically, the DSM model eliminates the uncertainty in a static Smagorinsky model (SSM) [50], especially helps to
correlate with varying in dierent regions of the ow (see Fig. 17 for the eddy viscosity contours). An attrac-
tive benet is to avoid the free-surface boundary conditions specied, as the satisfactory conditions is usually
lack for the turbulent variables involved, for example, by using a k turbulence model. Consequently, a sim-
ple DSM is coupled to our shock-capturing scheme as the blending algorithm: MILES + SGS. Hopefully, it
becomes a promising LES approach for modelling of breaking water-wave propagation problems. One way is
capable of capturing the inherently anisotropic small-scale ow features, as in MILES implicit SGS models
provide an additional damping, caused by nonlinear high-frequency lters built-in the convection
discretization.
Within the framework of the eddy-viscosity technique, the introduction of the DSM is adaptively
to estimate the adjustable dimensionless parameter C
s
, related with the unknown eddy viscosity m
T
.
Namely
m
T
C
s
D
2
jSj; 24
where jSj is the Frobenius norm of the strain tensor rate S
ij
, setting by jSj

2S
ij
S
ij
_
. In this case, it is
evaluated with the magnitude of the vorticity vector by denition of jSj jr
~
V j, where r
o
ox
;
o
oy
;
o
oz
_ _
and
~
V u
~
i v
~
j w
~
k. Additionally, the constraint C
s
= max(C
s
, 0) is articially enforced whenever C
s
returns
negative values, implying that no mechanism for backscatter is allowed, as done in Smagorinskys constant
model [30]. It is absolutely dissipative in computations, probably guaranteeing numerical stability.
Based on the multiple ltering operators [17], the subtest-scale stress tensor similar to Eq. (21) also holds by
T
ij


u
i
u
j

u
i

u
j
25
provided that the test lter

D (denoted by a tilde) is introduced and used in fact as the second lter to the NS
equations. For situations involving the same shape in the lter kernel as the grid lter D, the lter width of

D is
typically a function of D. Owing to the scale-similarity in the algebraic identity, the resolved turbulent stresses
D
ij
are obtained by
D
ij
T
ij

b
ij

u
i
u
j

u
i

u
j

u
i
u
j


u
i
u
j


u
i
u
j

u
i

u
j
; 26
which removes the problem in estimate of the quantity

u
i
u
j
. Under a rened discrete lter given (see Eq. (34),
D
ij
can be exactly evaluated with the resolved scales, as it physically represents the contribution of the smallest
resolved length scales to the Reynolds stresses.
On the other hand, the anisotropic subgrid-scale stresses b
ij
and subtest-scale ones T
ij
are formulated by,
respectively,
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 871
b
ij

1
3
d
ij
b
kk
2m
t
S
ij
2C
s
D
_ _
2
jSjS
ij
;
T
ij

1
3
d
ij
T
kk
2

m
t
S
ij
2C
s

D
_ _
2
j

Sj

S
ij
:
27
Substituting Eq. (27) into Eq. (26), this yields
D
ij
2C
s
D
2

D=D
2
j

Sj

S
ij


jSjS
ij
_ _
: 28
To achieve more stable values of C
s
, Eq. (28) is replaced with a least-square manner that accounts for the ef-
fects of all the components [29]. This ensures the error E to be minimal
E D
ij
2C
s
D
2
P
ij

2
29
by forcing
oE
oC
s
D
2

0. Thus, we have
C
s
D
2

hD
ij
P
ij
i

2hP
ij
P
ij
i

30
with P
ij
dened by
P
ij

D=D
2
j

Sj

S
ij


jSjS
ij
; 31
where the variation of C
s
on the scale of the test lter width

D is ignored. Because of the lter-grid ratio
c

D=D 2 (usually and more smooth than c = 4, as shown in Fig. 2), the estimate of P
ij
is not necessary
to dene the values

D and D, respectively. This feature is useful in some cases, as the ow structure near a
wall seems very anisotropic, leading to the denition of the length scale for anisotropic grids or lters unclear
[11].
The sign
*
stands for the averaging operation over the homogeneous directions, which is introduced addi-
tionally, for example, by using the following trapezoidal lter:
h/wi
i;j;k

9
16
/w
i;j;k

3
32
f/w
i1;j;k
/w
i;j;k1
/w
i1;j;k
/w
i;j;k1
g

1
64
f/w
i1;j;k1
/w
i1;j;k1
/w
i1;j;k1
/w
i1;j;k1
g: 32
Eq. (32) implies that averaging in the wall-normal direction (i.e. along the y-axis) is performed over eight
neighboring grid points surrounding its center (i, j, k), which avoids very small denominators or negative
numerator values. One problem is in the near-wall regions, where turbulent ows are characterized by much
less universal properties. Owing to some mathematical inconsistencies with the formulation of ltering oper-
ation, it appears that Eq. (30) is not applicable to general geometries, as C
s
is actually a function of time and
position [18]. For the present water-wave propagation problems involved, however, we focus attention on the
-0.12
-0.06
0
0.06
0.12
0 5 10 15 20 25 30 35 40 45

(
m
)
t (s)
WG: x=1.0 (m)
= 2
4
Exp.
0
0.002
0.004
0.006
0.008
0.01
20 22 24 26 28 30 32 34
C
s
t (s)
WG: x = 1.0 (m)
= 2
4
Fig. 2. Computations with two dierent lter-grid ratios c over grids (251 40) for Case 5 (pure regular waves): the time evolution of the
surface elevation g and Smagorinskys parameter C
s
(on the free surface) with c = 2 and 4.
872 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
capture of high turbulence that dominates over the wave-induced motions, generated in waves breaking. In
this sense, such problem associated with the breaking waves may be thought of as a kind of isotropic decaying
turbulence, due to without the uctuation of walls (if assuming that far from walls the SGS ow is ideally
considered homogeneous isotropic). For this reason, the dynamic model can remove the mathematical incon-
sistencies, i.e. C
s
can depend only on time (see Fig. 2) during averaging, as conducted with Eq. (32). With pos-
sibilities of simultaneously handling ow and grid anisotropies, on the other hand, our MILES may give a
better approach for inhomogeneous turbulent ows near the free-surface (or solid walls), because of the aniso-
tropic features of the SGS modeling in MILES [15].
To ensure smooth C
s
, three times in cycle seem enough, based on our experience. When evaluating P
ij
(see
Eq. (31) with inclusion of the vorticity components, a favourable choice is with high-order numerical methods,
for example, by using a fourth-order central-dierence approximation in smooth regions. In this way, the rst-
order derivative,
ou
ox

i
in the ith direction (where u
i
denotes a variable on the grid, and Dx is the grid size in the
relevant direction) is given by
ou
ox
_ _
i

u
i2
8u
i1
8u
i1
u
i2
12Dx
; 33
while a central-dierence scheme is employed in the near-wall regions and adjacent to the free surface. To re-
move the high spatial frequencies, moreover, we utilize a seven-point discrete test lter with trapezoidal rule
integration [51] for estimate of a ltered variable on the grid

u (see Eq. (26), i.e.

u
i

1
64
fu
i3
6u
i2
u
i2
15u
i1
u
i1
u
i3
44u
i
g 34
in smooth regions but weights for symmetric lter with three- or ve-point vanishing moments elsewhere.
The ltering is performed over regions of a sequence of one-dimensional lter for 2D or 3D problems. It
can eliminate the highest wavenumber that could be sustained by the grid [40], as the function becomes
zero at the grid cuto in spectral space, despite the absence of the inconsistency in the eective lter width
for Eq. (34). Note that we do not introduce a weighting function that accounts for the non-uniform grid
spacing eects, as done in our ship ow studies for ltering the wave height over a stretched grid
[32].
Based on the DSM, a smooth function of C
s
varying with time is given, as illustrated in Fig. 2, where
the shape of C
s
forms a more regular step function, when c = 2. By comparison, the DSM produces the
better and less physical dissipation behavior for the free-surface proles, as in our study Smagorinskys
coecient is constant (C
s
= 0.01), according to the SSM. Interestingly, the computations are not much sen-
sitive to the lter-grid ratio c (although c = 2 looks better), especially a large variation on C
s
does not
have too much of a bearing on the results. One reason for that is there exists no net eect of the SGS
model due to m
T
= 0, in case the derivative of the velocity (
ou
oz
) is small, where C
s
may be huge. The relevant
studies in the literature are given with [39,45,52], including two new SGS models, a dynamic free-surface
function model (DFFM) and a dynamic anisotropic selective model (DASM), proposed by Shen and Yue
[49].
2.4. Initial and boundary conditions
2.4.1. Initial conditions
At t = 0, all simulations start from a state of rest in quiescent water via (u, v, w) = 0 plus the hydrostatic
pressure distribution in the entire computational domain. After that, the accelerator in uid xed at the inow
boundary will vertically oscillate the static water, and as a result of forcing surface waves propagating towards
the positive x-direction.
According to a at free-surface at t = 0, the volume fractions a in VOF are initially assigned by
a
1 in the full fluid;
gdz
k1
Dz
k
over mixed cell;
0 in the air;
_

_
35
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 873
where g is the surface elevation (setup and setdown) with respect to the still water line, Dz is the vertical grid
size and the index k represents the mixed cell pointing to the free surface.
2.4.2. Boundary conditions
On the assumption that the free surface only supports the normal stress of a constant surface tension, the
more accurate normal free-surface boundary condition for the pressure (despite the absence of wind) is
expressed by
p qm
eff
2
ow
oz

og
ox
ou
oz

ow
ox
_ _

og
oy
ov
oz

ow
oy
_ _ _ _
36
as in our solver surface tension has been handled with a body force. Of course, Eq. (36) has been used for the
study of turbulent free-surface ows in 3D modern ships [32], which is referred to as the Reynolds stress free-
surface boundary condition. It essentially represents the normal force only acting on the free surface under the
pressure surrounding air. A major benet is that such boundary condition oers the potential for the capture
of the (thin) free-surface wave boundary layer, especially in waves breaking.
At the seabed, the free slip conditions are enforced and the no-ux condition is used on a cut cell (dened
by a cell partially dry) that is utilized for handling an arbitrary geometry, while the no-slip conditions are
implemented in the lateral direction. Additionally, a kind of Orlanski linear open boundary condition is
employed at the outow boundary
ou
ot
C U
c

ou
ox
0; 37
where u = (u, v, w, a, p). C is the phase velocity of the long wave and its upper limit value is set by C

gd
p
in
water of nite depth d. Theoretically, Eq. (37) simulates the propagation of velocity waves out of the compu-
tational domain. Unfortunately, the energy of waves is reected there (see Fig. 9), implying that the use of
such condition is not expected with high absorption eciency (probably, due to approximation of the phase
velocity). By comparison, an alternative is to introduce the physical wave absorber device for a tank, which
would be given, following the explanation of our design wavecurrent generator as inow conditions.
2.5. An external wavecurrent generator
By superimposing a steady current onto waves, the combined motions on the inow boundary are easily
obtained with the linear wave theory, which turns out to give a substantial solution (i.e. inow conditions)
in regions of irrotational ow, as the nonlinear features of the free surface can be captured in the interior
of the uid. In the context of this theory, the reected wave elevation g
f
at the inlet can be estimated with
the incident wave elevation g
i
specied via
g
f
g
t
g
i
38
provided that the total wave elevation g
t
is generated. The subscripts (t, f, i) represent the total, reected and
incident waves, respectively. Since we concern the steady vertically-uniform current speed U
c
[44], the reected
waves g
f
can be satised the modied Orlanskis theory [42]. This yields
og
f
ot
C
x
U
c

og
f
ox
0: 39
Substituting Eq. (38) into Eq. (39), the resulting wave surface is given by
og
t
ot
C
x
U
c

og
t
ox

og
i
ot
C
x
U
c

og
i
ox
40
indicating that the net waves relative to the weak disturbance of the reected waves are created at the inow
boundary. Substantially, the nal result can be considered as the addition of driving force similar to the con-
tribution under the gravity acceleration. Mathematically, it is constructed as the 1D hyperbolic wave equation
for the total wave g
t
. C
x
in Eq. (40) is the component of the group celerity in the x direction, and for any wave
propagation its derivation is always as follows:
874 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
C
x

ox
r
ok

1
2k
1
2kd
sinh 2kd
_ _
x
r
; 41
subject to the intrinsic angular frequency x
r
(the zeroth order approximation) and the wavenumber k dened
by
x
2
r
gk tanhfkdg and k
2p
k
; 42
where the wavelength k is precomputed from Eq. (49) below.
When solving this type of initial value problems, applying the second-order explicit AdamsBashforth for-
mulation yields
g
n1
t
g
n
t
DtR
n

Dt
2
R
n
R
n1
43
with an adaptive time step Dt, that accounts for the variation of the subgrid spacing plus the ow in the inte-
rior of domain, and the residual R dened by
R
og
i
ot
C
x
U
c

og
i
g
t

ox
; 44
where the desirable incident waves are specied a priori at the location of the inow boundary.
According to our previous study [30], the use of Eq. (40) helps to guarantee a consistent manner in solving a
due to a well-suited boundary value at the mixed cell even in the presence of a current. In particular, it only
creates the very weak reected waves.
Based on the velocity potential in uid, the wave elevation g
i
at ghost cells is given by
g
i

m
j1
H
2
cos k x
k
4
_ _
x
a
t
_ _ _ _
j
45
along a vertical wavemaking boundary in 2D. With the help of Eq. (44), the principal source term eventually
becomes
og
i
ot
C
x
U
c

og
i
ox

m
j1
H
j
2
x
a
C
x
U
c
k
j
sin k x
k
4
_ _
x
a
t
_ _
j
_ _
46
for the combined motions of irregular waves and a current, as it evolves in time t (where the variable (x) is
xed). m denotes the total number of waves (obviously, the simplest is m = 1 for a target of regular waves),
and the subscript j means the jth component. Given the standard wave energy spectrum of irregular wave
trains, the desirable irregular waves can be created by linearly superimposing a series of regular waves to a
specied number of terms (m = 10, enough for our computations), as illustrated in Fig. 3, when the charac-
teristics of one basic single regular wave (the wave height H and period T
a
) is specied.
Owing to a Doppler shifted frequency (kU
c
cos b), the absolute frequency x
a
can be evaluated by
x
a
x
r
kU
c
cos b or x
a

2p
T
a
; 47
-0.25
-0.12
0
0.12
0.25
0 5 10 15 20 25 30 35 40

t

(
m
)
t (s)
Fig. 3. Generation of irregular waves at the inlet, based on the signicant wave height (H = 0.16 m) and the wave peak period (T
a
= 2.0 s)
in the presence of a current (U
c
= 0.14 ms
1
).
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 875
where b is an angle of the wave orthogonal with the current streamline. In our present study, it is set to zero, as
waves propagate into a current without obliqueness.
Given the knowledge of the discrete quantity for g
t
(see Fig. 3), it is easy to obtain the corresponding veloc-
ity distributions at the inow boundary. One way is that three components (u
t
, v
t
, w
t
) of the velocity for the
total waves can be represented by superposing linearly the steady current speed U
c
onto the velocity induced
by incident waves that absorb the Doppler shifted frequency. This yields
u
t
U
c
cos b

m
j1
H
2
x
a
kU
c
cos b
cosh kg
t
d f g
sinhkd
cos k x
k
4
_ _
x
a
t
_ _ _ _
j
;
v
t
0;
w
t

m
j1
H
2
x
a
kU
c
cos b
sinh kg
t
d f g
sinhkd
sin k x
k
4
_ _
x
a
t
_ _ _ _
j
:
48
Accordingly, the wavelength k that absorbs a current is derived by
k
gk
2p
tanh
2pd
k
_ _1
2
U
c
cos b
_ _
T
a
; 49
which is solved iteratively, based on NewtonRaphsons technique that can give rapid convergence, under the
period T
a
and depth d specied. As a consequence, Eqs. (40) and (48) construct the basis of our external gen-
erator as wavy inow boundary conditions for the wavecurrent coupling. The major mechanism is that the
eects of the steady current have been absorbed to the wavelength (see Eq. (49), related with the Doppler-
shifted linear dispersion relation) together with the velocity proles (see Eq. (48), corresponding to full poten-
tial theory), while the combined ows are realized by superimposing a current onto waves.
We rst apply our external generator for modelling of waves travelling on three dierent current speeds
(positive, zero, negative), which can be recognized in Fig. 4 (refer to Fig. 5 for its geometry). By observation,
the external generator did well in the case of waves following a positive current, as the positive net ux formed
by the boundary can maintain the energy balance in the NWT over a long duration simulation. But it fails to
work for waves opposing the adverse current, in case the current is assigned as a net mass transport. Some
reasons may be interpreted as the fact: the uid leaving through the inow boundary could drain water out
of the tank, leading to a continuous drop in the mean water level. To solve particular counter current prob-
lems, the water has to be fed into the ow domain by the outow boundary, as done in our study [34]. In this
way, the computational domain is usually coupled with a numerical damping zone, which costs the substantial
CPU time plus the heavy memory, resulting in convergence to be extremely slow in general. Hence, a hopeful
alternative is to apply our internal generator, as described below, which creates an adverse current (see Figs.
1416), based on the source function. A major feature is that our NWT can be made as much small as possible
(at least in length), because of the use of a breaking-type wave absorber. This further maintains the minimal
expense for computations without instability problem even for lengthy computations.
Fig. 5 demonstrates a typical example regarding wave-structure coupling problems with and without a posi-
tive current, by using the external generator. When the synthetic inow motions fully develop and propagate
along the tank, the following waves overtop over a vertical barrier in the front of a pier, and then the waves are
strongly reected at the seawall, as one case of the diraction of regular and irregular waves by a barrier sitting
-0.3
-0.15
0
0.075
0.15
0 2 4 6 8

(
m
)
t (s)
WG: x = 1.0 (m)
Fig. 4. Computations with three dierent currents (positive, zero, negative) over grids (251 40) for Case 4 (regular waves): the wave train
at a certain position, by using the external generator. U
c
= 0.1 ms
1
; U
c
= 0.0; U
c
= 0.1 ms
1
.
876 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
above the bottom (see [31] for its 3D computations). From observation on this gure, the presence of the cur-
rent causes substantial changes for the water-wave dynamics, especially in the near-barrier region, character-
ized by a horizontal large-scale eddy shed from the barrier (also see Fig. 26 for the irregular waves). Closed to
the barrier, a violent overtopping jet occurs and surface waves naturally break against this structure, leading
to high turbulence induced by the breaking waves, and as a result of increase of the wave-induced loading act-
ing on the pier (due to the sea level rising). Additionally, the overall shapes in the wave trains exhibit the typ-
ical feature of nonlinearity after perturbing from its rest state.
2.6. An internal generator
In our present study, the source function f in Eq. (5) is incorporated into the continuity and momentum
equations as one type of a net driving force that creates an adverse current, as opposed to our previous study
[34]. In this sense, it can be thought of as an internal generator. Given the adverse current speed U
c
, the func-
tion f is formulated as follows, based on amount of mass ux over the nite volume region specied. This
yields
f
d gbjU

c
j
l
x
l
y
l
z
cos b
d gjU

c
j
l
x
l
z
cos b; 50
where (d, b) are the constant water depth and a nite width of the tank, respectively, (l
x
, l
y
, l
z
) represent the
size of the rectangular source region in the corresponding coordinates (x, y, z) and b is the same as that of
Eq. (47).
Note that l
y
= b in Eq. (50), which is one case equivalent to the uniform shape specied in the lateral direc-
tion. It is clear that this approximation is straightforward in 2D. In addition, the source strength involves the
eects of the surface elevation g at a certain location. By comparison, we found that computations look much
more ecient, when the source region covers several cells rather than one cell on a Cartesian grid (see Fig. 13).
By rule of thumb, its length (l
x
) should be thin (two to three cells) or less than 6% of the wavelength, while the
height (l
z
) is restricted to distance: from the wave trough below to water depth (d/3) above. To suppress noise
from suddenly starting at t = 0, the modied speed U

c
is given by
U

c
0:5U
c
f1 tanh10pt pg; 51
in which tanh(t) is smoothly varied from 1 to 1 over a specied number of wave periods.
-0.2
-0.1
0
0.1
0.2
0 5 10 15 20 25 30 35

(
m
)
t (s)
WG: x=3.81 (m)
-0.7
-0.35
0
0.15
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 33.5 (s)
U
c
= 0.0 m/s
1.0 m/s
-0.7
-0.35
0
0.15
0.3

0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 33.5 (s)
U
c
= 0.1 m/s
1.0 m/s
Fig. 5. Following-current-wave eects over grids (251 40) for Case 4 (regular waves): the wave train at a certain location and the
instantaneous velocity elds. U
c
= 0.1 ms
1
; U
c
= 0.0.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 877
Now, let us observe the eects of the prescribed net driving force: the water surface above the source region
(located at x = 7.89 to 8.04 m, covering three cells in this case) responds rapidly (ranging from t = 0.1 to 0.4 s)
and arises to a certain height (see Fig. 6). Closed to the source region, an adverse current generated ows
towards two the ends of the tank, forming spatially-varying current proles in regions of interest, because
the current is strongly coupled with incident waves. Interestingly, the mechanism of current generation is
through the rapid growth of sea level, similar to the water wave generation by underwater explosion, as gravity
is only the restoring force.
2.7. A wave absorber
In the course of modelling of wave propagation, the numerical solution is probably highly sensitive to
waves reected from articial boundaries, as waves propagate over distance much longer than several wave-
lengths. To minimize spurious reection in computations, a less ecient method is to extend the computa-
tional domain until a large distance, in which some results are achieved, before the signicant reection
error propagates upstream. Additionally, well-suited boundary conditions have to be specied at open bound-
aries in order to minimize the computer resources required for modelling of an open system. Because no exact
absorbing conditions are available at open boundaries, two basic categories recommended are numerical and
physical techniques, respectively. In the following, we emphasize how to specify suitable absorbing conditions
at the outow boundary for the conned tank (see Fig. 1), as in this case a well-designed open boundary con-
dition is essential for the semi-open NWT (for example, to eectively reduce the size of the computational
domain in length).
2.7.1. Numerical techniques
Two classes of the absorbing boundary conditions are usually applied for handling articial boundaries.
One is based on the wave equation [10] or its zeroth approximation like the open condition [42], where only
the waves (asymptotically ideal for low frequencies) leaving the outlet boundary are described. The other relies
on the sponge layer or articial damping zone technique (applicable to absorption of short waves), by means
of the addition of a forcing term, and as a result of reduction of oscillations in the solution near the boundary
[24]. Mathematically, it is substantially served as a function of the numerical dissipation, by incorporating the
appropriate articial damping terms
~
F
d
F
x
d
; F
y
d
; F
z
d
f0; 0; Sx; tgfu; v; wg
T
52
implying that only the vertical velocity is modied within the damping zone, corresponding to S(x) 6 0 (see
Fig. 7). Owing to denition of one body force, this term (
~
F
d
) is added to the right-hand side of the NS equa-
tions (see Eq. (5) for creating the negative contribution against incident waves. According to the slightly dif-
ferent formulation adopted in [47], the fringe function S dened is given by
0
0.04
0.08
0.12
3 4 5 6 7 8 9 10 11 12

(
m
)
x (m)
t = 2.0 s 1.5 1.0
0.4
0.2
0.1
U
c
= -0.5 ms
-1
-0.45
-0.3
-0.1
0
0.1
0.2
0 5 8 10 14
z

(
m
)
x (m)
t = 33.9 (s) U
c
= -0.5 ms
-1
wave absorber
0.1C m/s
Fig. 6. Spatial variation of tidal elevation from starting initially and the instantaneous velocity elds over grids (291 56) for Case 3
(regular waves), by using the internal source function.
878 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
Sx; t S
max
t F
x x
s
LD
rise
_ _
F
x x
t
LD
fall
1
_ _ _ _
53
with the aid of the length (L) of the computational domain and two coecients (D
rise
, D
fall
), where S
max
(t) is
the maximum strength of the fringe function and the subscripts (s, t) stand for the position of two ends (start-
ing and terminal) in the damping zone, respectively. To avoid partial reection, S(x) is smoothly varied along
the damping zone, especially near the entering point. For this reason, F(r) is assigned by
F r
0 r 6 0;
1
1e
1
r1

1
r
0 < r < 1;
1 r P1;
_

_
54
as shown in Fig. 7, where the curve plotted in the damping zone (total length of about 1.3k) consists of three
parts, starting from x = 12.8 (m) to nal step constant value S
max
= 50.0 s
1
near the end wall (D
rise
= 0.6 and
D
fall
= 0.1 in this case).
In fact, one often couples the sponge layer technique with a piston-like Neumann condition as the hybrid
algorithm, which is very ecient in a broader spectrum of waves [6]. Only requirement is to extend the extra
length (normally several wavelengths) at the end of a tank so that energy of outgoing waves is eectively dis-
sipated. Consequently, the performance depends heavily on the ratio of the beach length to wavelength k [5].
Note that such restriction also holds for an internal wavemaker [35], by adding a length at both ends of the
tank in general.
2.7.2. A breaking-type wave absorber
We attempt to apply a simple wave absorber, i.e. an articial sloping ridge is equipped in the extremity of
the NWT. It can be considered as one type of physical techniques related with the physical model (see [26]
with a perforated metal sheet or [12] for a similar device). In our study, the smooth impermeable structure
with constant slope (1/6 or less) is connected to a region with constant water depth, while the top of the
structure lies just slightly above the still water level (see Fig. 9, for example). One possible explanation is
that the wave energy can be eectively dissipated by breakwaters without the need of adding an extra
length, before waves reach its extremity. Accordingly, a simple local boundary condition at the outow
boundary can yield a well-posed solution, especially the wave reection can be restricted to the minimal
level.
A typical example, Case 3 (see Table 1), is selected for testing the performance of a breaking-type wave
absorber in 2D. By visualizing the typical velocity elds (see Figs. 1416), the high speed ows are created,
while an expanding jet overows along the crest of the absorber, leading to localized high turbulence. Conse-
quently, dissipation of energy due to waves breaking is the dominant process in the proximity of this absorber,
and as a result of a loss of wave energy in breaking waves. For this situation, the reected waves are expected
to be small, which can be convinced in Fig. 8: the general trend of the wave proles (averaged over 10 periods)
and the mean velocity proles (by interpolation at x = 7.8 m) follows closely that of experiments. It reveals
that considerable turbulence is rapidly damped out, under the dynamic SGS model. For this reason, our
NWT can be made small, which helps to reduce the CPU time.
0
10
20
30
40
50
0 5 10 15 20
S
(
x
)
x (m)
damping zone
Fig. 7. The damping zone described by the fringe function S(x) for Case 3.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 879
Based on our analysis, the high computational eciency is achieved, by using a breaking-type wave absor-
ber. This speculation can be further conrmed, by comparison of two types of traditional approaches applied
at the outow boundary: a homogeneous Neumann boundary and the open condition (see Eq. (37). Denitely,
such comparison is still of considerable interest for nonlinear water-wave propagation problems, especially
when fully turbulent ows are caused by waves breaking.
2.7.3. Comparison of two classes of wave absorbers
As expected, a more interesting reference case is concerned with the homogeneous Neumann boundary (i.e.
ow
ox
0 and u = 0, implying a free-slip rigid wall), since there exists a perfect reecting wall articially placed at
the outlet, where the return ow (undertow) occurs. In contrast of the Neumann boundary, the weakly reect-
ing condition (see Eq. (37) allows radiating long waves passing through the outlet. Owing to possible gener-
ation of the fully (or partially) reected waves from the boundary, a long domain is useful for both
approaches, which should cover a damping zone. Fig. 9 displays reection and absorption eciency in terms
of the wave train, instantaneous velocity elds and velocity proles (throughout this paper, the symbol C in
the velocity elds drawn stands for the linear shallow wave celerity, unless otherwise stated), when compared
with measurements available.
By observation, the perturbations caused by the wall in the early stages are negligible almost everywhere.
Probably, there is no serious problem with the simple treatment at the outlet, when modeling of a solitary
wave or dam-break ow [23]. After that, the degree of wave reection varies with the Neumann boundary
to breaking-type wave absorber, being strong for the former and practically negligible or weak for the latter
(see Fig. 9). This is one example that has been turned out to be very ecient for wave absorption, by using the
breaking-type wave absorber (also see Fig. 8). The wave trains, for instance, exhibit typically nonlinear fea-
tures: high crests and smaller troughs. The fair agreement between predicted and measured results gives partial
conrmation: the perfect energy dissipation on a beach with a shorter scale, implying only mildly reecting.
With the Neumann boundary condition, however, the results are contaminated due to occurrence of appre-
ciable and successive wave reection, leading to badly-behaved solution. One reason for that is because the
incident waves propagate against the return ow to occur under troughs or interact with the reected waves
from the vertical wall or both. Additionally, introduction of the articial damping zone cannot eectively dis-
sipate the energy of outgoing waves, although the waves progressively lose some energy in this area (see Fig. 9,
where the damping zone starts from x = 16.0 m). Probably, the waves may be entirely absorbed but less e-
cient, provided that the length of damping zone is sucient (the several wavelength or more). The same con-
sequence also holds with the open condition, as in this case an impulsive velocity is generated at the outow
boundary, resulting in the articially reected waves there. Consequently, high absorption eciency of our
physical absorbing device and its cost are highlighted by comparison with the two reference simulations men-
-0.1
-0.05
0
0.07
0.14
0 0.25 0.5 0.75 1

(
m
)
t /T
a
U
c
= 0.0 ms
-1
Waves alone
Exp
-0.45
-0.3
-0.1
0
0.1
-0.1 -0.05 0 0.05 0.1
z

(
m
)
u (ms
-1
)
U
c
= 0.0 ms
-1
Waves alone
Exp
Fig. 8. Comparison between computations and measurements on grids (291 56) for Case 3 (regular waves alone): the mean surface
elevation and the vertical proles of mean horizontal velocity, by using the breaking-type wave absorber.
880 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
tioned. An additional benet is that it can apply to various wave conditions without the requirement in
knowledge of incident waves, except for extra cost adding one simple geometry (a 2D mountain ridge) to
our solver.
3. Test cases
To validate the accuracy of our design wavecurrent generators, ve test cases (see Table 1) under consid-
eration are applied, which involve three types of dierent applications of interest:
the fully coupled wavecurrent motions in the open channel over constant water depth (labelled Case 1 to
Case 3), for which measurements are available [26,27,12];
overtopping of violent waves over a vertical barrier, responsible for impact of waves onto the structure,
leading to strong reection at a seawall (referred to as Case 4);
and overtopping of waves over a seadike, characterized by various physical phenomena in water of varying
depth (named Case 5).
-0.1
-0.05
0
0.06
0.12
0 10 20 30 40

(
m
)
t (s)
U
c
= 0.222 ms
-1
WG: x = 4.0 (m)
-0.45
-0.3
-0.1
0
0.1
0.2
0 3 6 10 15 20
z

(
m
)
x (m)
t = 31.0 (s) damping + Neumann 0.1C m/s
-0.45
-0.3
-0.1
0
0.1
0.2
0 3 6 10 15 20
z

(
m
)
x (m)
t = 31.0 (s) damping + open 0.1C m/s
-0.45
-0.3
-0.1
0
0.1
0.2
0 3 6 10 14
z

(
m
)
x (m)
t = 31.0 (s) U
c
= 0.222 ms
-1
wave absorber
0.1C m/s
-0.45
-0.3
-0.1
0
0.1
0 0.1 0.2 0.3 0.4
z

(
m
)
u (ms
-1
)
U
c
= 0.222 ms
-1
Fig. 9. Comparison of three types of wave absorbers over grids (291 56) for Case 3 (the regular following waves): the wave trains at one
certain position, subsequent velocity elds and vertical proles of horizontal velocity (by interpolation at x = 2.9 m). the breaking-
type wave absorber; a Neumann condition plus the damping zone; the open condition plus the damping zone; } measurements.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 881
When simulating waves following or opposing a turbulent current (Case 1 to Case 3), various theory mod-
els are applied and there has been signicant progress in this area. Based on the linear wave approximation,
for example, Huang and Mei [22] proposed an analytical boundary-layer method for Case 1 plus a little in
Case 2. For the latter, more detailed study was achieved by Groeneweg and Battjes [21], according to the gen-
eralized Lagrangian mean formulation. Both also involve computations of waves on a negative current, while
Fredse et al. [12] borrowed the Reynolds-averaged NavierStokes (RANS) equations with a kx model for
solution of the near-bottom motions in Case 3. But the limited values are given in the turbulent regions near
the free surface (such as the instantaneous velocity elds), especially the lack of the wave train (importantly, as
such assessment helps to observe whether mass conservation holds for a large number of wave periods). Addi-
tionally, no one attempts to address grid resolution, which is essential on uncertainty analysis, including the
study of convergence issues.
In our present study, all ve test cases are demonstrated with our own solver. Denitely, this is very valu-
able in assessing our design wavecurrent generators. As illustrated below, we rst concern the open channel
ow and then study waves on beaches. It is shown that our solver captures most of intrinsic features that
would arise in the wavecurrent interaction, especially it provides more detailed quantitative analysis of the
velocity elds plus waves. For example, the surface proles averaged over 10 periods are given (see Fig. 11)
and the mean horizontal velocity proles at a single station are easily obtained by interpolation at
x = 7.0 m (see Fig. 12). The computations are performed over a cut-cell grid, which are identical to physical
conditions in experiments, except at the bed boundary, where the rough bed in experiments is against the
-0.04
-0.02
0
0.02
0.04
0 0.25 0.5 0.75 1

(
m
)
t /T
a
a U
c
= -0.11 ms
-1
411 X 80
291 X 56
Exp
-0.1
-0.05
0
0.07
0.14
0 0.25 0.5 0.75 1

(
m
)
t /T
a
b
U
c
= 0.222 ms
-1
291 X 80
291 X 56
Exp
Fig. 11. The mean surface proles over two meshes. (a) the opposing waves for Case 1 (regular waves), by the internal generator; (b) the
following waves for Case 3 (regular waves), by the external generator.
-0.1
-0.05
0
0.07
0.14
2 4 6 8 10

(
m
)
x (m)
U
c
= 0.222, 0.0, -0.222 ms
-1
t = 33.5 (s)
Following waves
Pure waves
Opposing waves
Fig. 10. The spatial surface proles under three dierent current speeds (positive, zero, negative) over grids (291 56) for Case 3 (regular
waves), where more details are given in Fig. 16.
882 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
smooth one in simulations. Two Cartesian grid levels with varying cell sizes are employed for the grid rene-
ment study. The results are represented in terms of the instantaneous velocity elds in 2D (see Fig. 1(b)) and
the time history of the surface elevation g (m) at the selected wave gauge (WG), as these variables are directly
given, by solving the NavierStokes equations with the free surface.
3.1. Wavecurrent coupling in water of constant depth
As incident waves propagate into the fully developed shear ow, the adverse current is strongly coupled
with the waves, forming the opposing waves (see Figs. 1416). Very encouragingly, our results give the com-
prehensive comparison for waves on three types of the current states varying with the signed number for Case
1 to Case 3. Denitely, they also provide the physical insight for the interaction of waves with a current. On
these gures, a major feature is essentially that waves on an adverse current are shorter than waves on a posi-
tive one, as the current compresses the wavelength somehow (see Fig. 10), which is consistent with the point
drawn in Eq. (49) for the given values of T
a
and d. An expected phenomenon is that the phase-averaged sur-
face proles display typically nonlinear features, especially in the opposing waves. For the latter, the dynamics
of waves can vary signicantly with current-induced wave breaking, in case the waves are blocked by the
strong opposing currents [7,8], where the change in the gross topology will undergo the processes of merging
and breakup. For some particular positive current problems (Case 1 to Case 2), the attenuation of wave
heights relative to pure waves is obvious, as the waves are travelling on the current.
Interestingly, the proles of the longitudinal velocity (u) may t the sequence near the free surface (see
Fig. 12, where the corresponding data are not available in measurements), which display a notable reduction
in following waves, as observed in the physical model available (see Fig. 9). According to our computations,
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
-0.3 -0.2 -0.1 0 0.1 0.2
z

(
m
)
u (ms
-1
)
U
c
= -0.16, 0.16 ms
-1
Opposing waves
Exp
Following waves
Exp
Current alone
Exp
Fig. 12. The vertical proles of mean velocity under two dierent currents (positive, negative) over grids (291 56) for Case 2 (regular
waves), including current-only situation, by using the internal generator.
-0.04
-0.02
0
0.02
0.04
0 10 20 30

(
m
)
t (s)
U
c
= -0.11 ms
-1
WG: x = 1.0 (m)
-0.2
-0.1
0
0.1
-0.2 -0.15 -0.1 -0.05 0 0.1 0.2
z

(
m
)
u (ms
-1
)
U
c
= -0.11 ms
-1
Finite-volume source
Point source
Exp
Fig. 13. Comparison of two dierent source sizes over grids (291 56) for Case 1 (the regular opposing waves): the wave train at a certain
location and the mean velocity proles (by interpolation at x = 1.25 m), by using the nite-volume source () and the point source
( ), respectively.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 883
the predicted velocity diminishes continuously to zero towards two ends of the curve, whether or not there
exists a current (also see Fig. 8). One possible explanation is the presence of a transition zone, due to the eects
of the surface boundary layer both at the seabed and free surface. For situations where the velocity proles are
predicted well, for example, in the case of following waves, simulation of the wave height tends to be high
accurate (see Fig. 11): the shape of the wave proles may be resolved adequately, especially an excellent agree-
ment between theory and experiment is achieved with increased mesh renements. In this sense, the discrep-
ancy attributed to the grid eects is relatively small for variable of the wave height (in regular waves). This
behavior also holds in the case of the opposing waves but for the velocity proles the discrepancy is apparent
in the proximity of the free surface, where the waves experience an increase in speed towards the free surface
(also see Fig. 13). The reason for that is probably attributed the turbulent structures in this area, where the
turbulent properties are highly anisotropic [49]. On the other hand, the wave boundary layer is likely to extend
up into the current-induced logarithmic layer. A local adaptive mesh adjacent to interfaces is attractive for
capture of the fully developed turbulent boundary layer. In the near-bottom region, the dierent local features
are partly due to dierent seabed conditions regardless of resolution that is employed.
-0.2
-0.1
0
0.1
0 2 4 6 8
z

(
m
)
x (m)
t = 15.5 (s) U
c
= 0.11 ms
-1
wave absorber
0.1C m/s
-0.2
-0.1
0
0.1
0 2 4 6 8
z

(
m
)
x (m)
t = 15.5 (s) U
c
= 0.0 ms
-1
wave absorber
0.1C m/s
-0.2
-0.1
0
0.1
0 2 4 6 8
z

(
m
)
x (m)
t = 15.5 (s) U
c
= -0.11 ms
-1
wave absorber
0.1C m/s
-0.04
-0.02
0
0.02
0.04
0 5 10 15 20

(
m
)
t (s)
U
c
= 0.11,0.0, -0.11 ms
-1
WG: x = 1.0 (m)
-0.02
0
0.02
0.04
0 0.25 0.5 0.75 1

(
m
)
t /T
a
U
c
= 0.11,0.0, -0.11 ms
-1
Following waves
Pure waves
Opposing waves
Fig. 14. Following waves, pure waves and opposing waves over grids (291 56) for Case 1 (regular waves): the instantaneous velocity
elds, the wave trains at one certain position and the mean surface proles. U
c
= 0.11; 0.0; 0.11 ms
1
.
884 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
This is our rst results in the study of waves opposing or following a turbulent current, based on our
internal or external generator, respectively. In the case of the former, the strength of the source distributed
spatially is massless, which allows the ow to pass through the source region, represented with a small box
in our test case. For example, the height (l
z
) of this box and its length (l
x
) for Case 1 are 0.093 and 0.09 m,
respectively, covered 11 and 3 cells (i.e. the nite-volume source) on the given grid. By using the point
source (only one cell in height and length), it is found that the results are not much sensitive to variation
of the size in the source region (see Fig. 13), whereas signicant errors are produced in some transient
time, which seems less ecient.
3.2. Waves on beaches
Our solver describes well waves following (or opposing) a turbulent current over constant water depth for
Case 1 to Case 4. The reminder is that we would address the aspect in 2D and 3D for waves on beaches (Case
5), which is a more challenging case, because of the impressive natural phenomenon. We rst investigate irreg-
ular waves in 2D and then discuss the 3D eects in regular waves. The results for waves following a positive
current are summarized below, based on our external generator.
-0.5
-0.3
-0.1
0
0.1
0.2
0 5 10 14
z

(
m
)
x (m)
t = 32.4 (s) U
c
= 0.16 ms
-1
wave absorber
0.1C m/s
-0.5
-0.3
-0.1
0
0.1
0.2
0 5 10 14
z

(
m
)
x (m)
t = 32.4 (s) U
c
= 0.0 ms
-1
wave absorber
0.1C m/s
-0.5
-0.3
-0.1
0
0.1
0.2
0 5 10 14
z

(
m
)
x (m)
t = 32.4 (s) U
c
= -0.16 ms
-1
wave absorber
0.1C m/s
-0.1
-0.05
0
0.05
0.1
0 10 20 30 40

(
m
)
t (s)
U
c
= 0.16, 0.0, -0.16 ms
-1
WG: x = 6.1 (m)
-0.1
-0.05
0
0.05
0.1
0 0.25 0.5 0.75 1

(
m
)
t /T
a
U
c
= 0.16, 0.0, -0.16 ms
-1
Following waves
Pure waves
Opposing waves
Fig. 15. Following waves, pure waves and opposing waves over grids (291 56) for Case 2 (regular waves): the instantaneous velocity
elds, the wave trains at a certain position and the mean surface proles. U
c
= 0.16 ; 0.0; 0.16 ms
1
.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 885
By visualizing the wave-induced motions in regions of interest with several stages (see Fig. 18), it is found
that ows become highly turbulent under waves breaking, often subjected to identiable jet-splash cycles,
when waves overtop over the seadike. After that, the breaking waves transform into a turbulent bore, char-
acterized by a steep turbulence front and an area of recirculating ow. In the wave front, the particle velocity
exceeds the speed limit underneath the region of ows, leading to the rapid increase in high turbulent kinetic
energy (t = 32.4 s, for example). Interestingly, there exists two possible sources of turbulence generation for
waves on beaches. One is the initial instability that leads to breaking, because of continuous increase of the
wave steepness (t = 4.5 s); the other is the proceeding of the downwash, which also induces waves breaking
(t = 33.4 s). Denitely, the presence of a current modies wave characteristics, especially in regions subjected
to strong turbulent ows, where a large hydraulic jump is created nearby at x = 4.0 m (t = 32.2 s), under cer-
tain circumstance of irregular waves alone.
In contrast to this, the positive current stretches the wavelength, leading to reduction in the steepness of
waves. Consequently, the wavecurrent interaction tends to be less severe, probably similar to computa-
tions in deep water waves, characterized by a thick sheet of water on the beach (dened as the so-called
swash), in which the nonlinearity is weak. Such eects tend to be considerable in general, as a current
increases.
-0.45
-0.3
-0.1
0
0.1
0.2
0 5 10 14
z

(
m
)
x (m)
t = 33.5 (s) U
c
= 0.222 ms
-1
wave absorber
0.1C m/s
-0.45
-0.3
-0.1
0
0.1
0.2
0 5 10 14
z

(
m
)
x (m)
t = 33.5 (s) U
c
= 0.0 ms
-1
wave absorber
0.1C m/s
-0.45
-0.3
-0.1
0
0.1
0.2
0 5 10 14
z

(
m
)
x (m)
t = 33.5 (s) U
c
= -0.222 ms
-1
wave absorber
0.1C m/s
-0.08
-0.04
0
0.06
0.12
0 10 20 30 40

(
m
)
t (s)
U
c
= 0.222, 0.0, -0.222 ms
-1
WG: x = 6.1 (m)
-0.06
-0.03
0
0.06
0.12
0 0.25 0.5 0.75 1

(
m
)
t /T
a
U
c
= 0.222, 0.0, -0.222 ms
-1
Following waves
Pure waves
Opposing waves
Fig. 16. Following waves, pure waves and opposing waves over grids (291 56) for Case 3 (regular waves): the instantaneous velocity
elds, the wave trains at a certain location and the mean surface proles. U
c
= 0.222; 0.0; 0.222 ms
1
.
886 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
On the other hand, the wave pattern develops rapidly after the transient phase and highly sharp-crest waves
(i.e. one typical feature of irregular waves) are well captured (see Fig. 17), which always occur in real sea. As
expected, the phase shift is the result of the Doppler eects due to the presence of a current. Interestingly, the
high turbulence activity occurs at phase mostly within the surf zone, subjected to rapid deformation of the free
surface, where the eddy viscosity m
t
dominates over the molecular viscosity m (see Fig. 17 for varying with m
t
).
3.3. 3D eects
To observe the 3D eects, the source of the wavemaker in 2D (ranging from the bottom to the free surface)
is uniformly distributed along the inlet section extended to the width of 1.0 m. Geometrically, this shape con-
structs a cubic tank of edge length 6.3 1.0 1.0 m
3
separately in the longitudinal, lateral and vertical direc-
tions (see Fig. 1(a)). For the present 3D seadike problem, any longitudinal prole is identical with the 2D case
(see Fig. 1(b), referred to as the symmetry plane of y = 0, for convenience). Based on the usual CFD validation
procedure, the convergence performance is rst discussed for regular waves alone, and then computations with
dierent grids are conducted for investigation of the asymptotic grid convergence. Accordingly, an interesting
comparison between 3D and 2D cases is described. Finally, the asymmetry of ows at the cross-sections (i.e.
the yz plane) is detected for various horizontal positions in the presence of a current or not. Summarily, our
results demonstrate that the gravity-wave breaking is essentially a three-dimensional process only in the surf
zone (where the ow shape loses its symmetry and becomes fully 3D), as surface waves move up the slope,
resulting in a 3D convective instability. In this sense, the wave propagations dier substantially from strong
internal waves that may cause intense shear in 3D [3].
3.3.1. Convergence history
By examining the L
2
norm of the residuals (Res), Res
u
f

N
m1
u
n
i;j;k
u
n1
i;j;k

2
=Ng
1
2
with the total number
of cells N and u = (U, V, W, P), Fig. 19 exhibits a set of the convergence history in iterations for the momen-
tum equations (U, V, W) in the (x, y, z)-directions plus the pressure (P). In this case, computations are con-
ducted over one single mesh (205 21 40) with varying cell sizes, except in the y-direction, where the grid
spacing is uniform.
As expected, the fast convergence rate about the pressure is achieved due to the computationally ecient
algorithm developed in our solver. This is divided into two regimes according to the trend of curves: rst to
decay nearly exponentially with a relatively short time and then turns to the at until a certain time level spec-
ied by a user (see Fig. 19). Based on the fact that solution of the pressure-Poisson equation is essentially the
most time-consuming part for any incompressible ow problems, therefore, our solve exhibits the well overall
performance in computations. Furthermore, mass conservation ensures that our simulations are stable even
during lengthy computations. This is because errors are eliminated exponentially with the iterative number,
in case the errors in mass are mainly caused by the residuals from the continuity equation (in fact, these
are principally composed of high-frequency error components [36]). According to a decoupled approach, it
0
0.1
0.2
0.3
0 5 10 15 20 25 30 35

(
m
)
t (s)
WG: x = 5.2 (m)
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)

t
/ = 92.3
20.1
free surface
U
c
= 0.14 m/s
Fig. 17. Current eects over grids (251 40) for Case 5 (irregular waves): the time history of the surface elevation at one wave gauge
selected and the eddy viscosity contours (m
t
/m) at t = 32.2 s (magnied view: 5 to 100 with interval 10). U
c
= 0.14 ms
1
; U
c
= 0.0.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 887
is seen that the residuals of the velocities (U, V, W) almost follow the similar trend (see Fig. 19), whereas some
local peaks over these three curves are clearly magnied, probably corresponding to the steep gradients of the
velocity in waves breaking.
3.3.2. Grid renement study
Three dierent meshes are applied for study of the eects of increasing grid density on waves, as illustrated
in Figs. 20 and 21 for the wave train and velocity elds.
First, all of the simulations in 3D start from rest as one kind of the most general situations in the NWT,
which obviously dier from the previous studies [32,13]: for the former, the steady-state solution on the 3D
coarse grid is assigned as an initial guess for the ne solution; on the other hand, the calculated results in
2D are interpolated onto the 3D domain for the latter. Hence, no white noise (or an initial disturbance) is
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 3.0 (s)
U
c
=0.0 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 3.0 (s)
Uc = 0.1 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 3.0 (s)
U
c
=0.14 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 4.5 (s)
Uc = 0.0 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 4.5 (s)
Uc = 0.1 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5

6
z

(
m
)
x (m)
t = 4.5 (s)
U
c
=0.14 m/s
0.1C m/s
0.3
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
6 5 4 3 2 1 0
z

(
m
)
x (m)
U
c
=0.0 m/s
0.1C m/s t = 32.2 (s)
0.3
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
6 5 4 3 2 1 0
z

(
m
)
x (m)
U
c
=0.1 m/s
0.1C m/s t = 32.2 (s)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 32.2 (s)
U
c
=0.14 m/s
0.1C m/s
0.3
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
6 5 4 3 2 1 0
z

(
m
)
x (m)
U
c
=0.0 m/s
0.1C m/s t = 32.4 (s)
0.3
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
6 5 4 3 2 1 0
z

(
m
)
x (m)
U
c
=0.1 m/s
0.1C m/s t = 32.4 (s)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 32.4 (s)
U
c
=0.14 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 33.2 (s)
Uc = 0.0 m/s
0.1C m/s
0.3
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
6 5 4 3 2 1 0
z

(
m
)
x (m)
U
c
=0.1 m/s
0.1C m/s t = 33.2 (s)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 33.2 (s)
U
c
=0.14 m/s
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 33.4 (s)
Uc = 0.0 m/s
0.1C m/s
0.3
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
6 5 4 3 2 1 0
z

(
m
)
x (m)
U
c
=0.1 m/s
0.1C m/s t = 33.4 (s)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 33.4 (s)
U
c
=0.14 m/s
0.1C m/s
Fig. 18. Waves on the beach over grids (251 40) for Case 5 (irregular waves): the instantaneous velocity elds, where the current speed
increases gradually from zero to U
c
= 0.14 ms
1
.
888 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
-7
-5
-3
-1
0
2
0 5 10 15 20 25 30 35
L
o
g
1
0

(
R
e
s
)
t (s)
U
V
W
P
Fig. 19. Convergence history in 3D on the number of grids 172,220 for Case 5 (regular waves alone): the L
2
norm of the residual
(U, V, W,P).
-0.12
-0.06
0
0.06
0.12
0 5 10 15 20 25 30

(
m
)
t (s)
WG: x=1.0 (m)
-0.16
-0.08
0
0.08
0.16
0 5 10 15 20 25 30

(
m
)
t (s)
WG: x=3.81 (m)
Fig. 20. The eects of three dierent grids on waves (3D) for Case 5 (regular waves alone): the wave trains at two certain points on the
symmetry plane of y = 0. 205 21 40; 205 21 28; 145 21 28; } Exp.
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.0 s 1.0 m/s
145X21X28 (x=2.02 m)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s 1.0 m/s
145X21X28 (x=2.02 m)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.0 s 1.0 m/s
205X21X40 (x=2.02 m)
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s 1.0 m/s
205X21X40 (x=2.02 m)
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
0.5 0.25 0 -0.25 -0.5
z

(
m
)
y (m)
205X31X40 (x=2.02 m)
1.0 m/s t = 11.0 s
0.2
0.1
0
-0.1
-0.3
-0.5
-0.7
0.5 0.25 0 -0.25 -0.5
z

(
m
)
y (m)
205X31X40 (x=2.02 m)
1.0 m/s t = 11.5 s
Fig. 21. The eects of three dierent grids on waves (3D) for Case 5 (regular waves alone): the instantaneous velocity elds at a certain
cross-section.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 889
imposed on our three-dimensional turbulent computation, until the rst breaking event occurs, when the ini-
tial development of the free surface is driven by our external generator, characterized by the capture of the
transient in waves (see Fig. 20). By observation, a common knowledge for the wave pattern is that more useful
information in the motions can be represented with increased mesh renement, especially close to the surf
zone, where the ows are essentially in 3D (see Fig. 20). Secondly, an excellent agreement with measurements
available is achieved (whereas there exists some discrepancy at one farther position, x = 3.81 m, from the
inlet), as a mesh increases. This is the expected behavior, as physically a higher resolution allows the numerical
schemes to dissipate on a broader wavenumber range. Finally, the grid resolution is sensitive for the velocity
elds (see Fig. 21), which probably describe the obliquely descending eddies (t = 11.5 s) entraining sediments.
Such coherent turbulent structures descend from the surface towards the bottom. Interestingly, this phenom-
enon causes asymmetry of the ow in the lateral direction, closely related with the longitudinal vortices
(dened by x
x

ow
oy

ov
oz
). For our 3D problems, therefore, the conclusion of the mesh renements is almost
identical with the case in 2D, as stated in our study [30] for the comments.
3.3.3. Comparison between 2D and 3D
On the assumption that comparison is made on the central plane (y = 0), Fig. 22 provides an illustrative
observation for this one in terms of the wave trains selected. One feeling is that the results in 3D remain very
similar to those in 2D within the non-breaking region (WG: x = 2.02 m) but far from satisfactory at the posi-
tion adjacent to the wave-breaking region (WG: x = 3.81 m). It is clear that these results directly give reliable
evidence: waves in breaking are a rapidly dissipative process in the localized region. Hence, the use of a 3D
solver coupled with SGS models helps to capture the wave-induced small-scale turbulent ows, in case the
conguration in 2D becomes unstable due to waves breaking. In particular, 3D gives more promising results,
as compared to 2D (see Fig. 22), where the development of waves reaches a fairly steady pattern after the tran-
sient from rest. Of course, the shape of waves exhibits some temporal uctuations, which may be indicative of
the generation of high-order harmonics in computations, resulting from a slight asymmetry for the surface
proles.
Before breaking, a particular situation exhibits that the wave dynamics is still in 2D, as the pre-breaking
simulations are almost identical but quite dierent after breaking. In 3D, the breaking wave-induced second-
ary circulation usually emerges as counter-rotating vortex (pairs) in conned channels (see Fig. 21). It moves
with a strong central ow, occupying the most of the cross-section, especially in increasing the mesh. But the
-0.12
-0.06
0
0.06
0.12
0 5 10 15 20 25 30

(
m
)
t (s)
WG: x=2.02 (m)
-0.16
-0.08
0
0.08
0.16
0 5 10 15 20 25 30

(
m
)
t (s)
WG: x=3.81 (m)
0
0.05
0.1
0.15
0.2
0 5 10 15 20 25 30

(
m
)
t (s)
WG: x = 5.2 (m)
Fig. 22. Comparison of 3D with 2D for Case 5 (regular waves alone): the wave trains at three positions on the symmetry plane (y = 0).
3D (205 21 40); 2D (205 1 40); } Exp.
890 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
predicted secondary current is rather weak, based on the results in 2D and 3D at the same position (see
Fig. 22: x = 2.02 m), indicating that small disturbances in 3D are not sucient to aect the two-dimensional
results (although the results in 3D look better). This is consistent with the study in [21], as the secondary cir-
culations tend to be one or two orders of magnitude smaller than the longitudinal velocity.
Owing to stretch and intensication in shear, the three-dimensional dynamic characteristic of the motions
will grow up rapidly in the surf zone during waves breaking, and as a result of an initial two-dimensional con-
vective instability that leads to the transition to three-dimensional vortex structure. Physically, the convective
instability causes the longitudinal rolls (i.e. spanwise instability), responsible for dissipative energy in the
broader range, implying that dissipation for the two-dimensional ow is substantially limited (due to only
one cell in the y-direction). In this sense, such mechanism may be considered as the principal dierence
between turbulence dynamics in 2D and 3D. Based on the above analysis, it can be thought that the three-
dimensional dynamics is dominated by the secondary cross-stream convective rolls, that account for a great
portion of the energy loss [13]. For this reason, the velocity in the shear ow is overestimated [53], resulting
from the larger surface proles in this area (see Fig. 22: x = 5.2 m), due to the lack of physical dissipation in
the 2D case.
3.3.4. Sidewall eects
Figs. 23 and 24 display the eects of sidewall on the ow in the (y, z)- and (x, y)-planes, respectively. Both
describe the growth in the 3D motions with the vortex patterns in the presence of the positive current or not,
when waves propagate in tanks with the nite width.
Under waves breaking, the rapid growth of the spanwise perturbations causes the counter-rotating vortices,
leading to the three-dimensional characteristics of the breaking dynamics. But in the non-breaking region
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s
U
c
= 0.0 (x=1.0 m)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s
U
c
= 0.17 ms
-1
(x=1.0 m)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s
U
c
= 0.0 (x=3.81 m)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s
U
c
= 0.17 ms
-1
(x=3.81 m)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s
U
c
= 0.0 (x=4.5 m)
0.1C m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
-0.5 -0.25 0 0.25 0.5
z

(
m
)
y (m)
t = 11.5 s
U
c
= 0.17 ms
-1
(x=4.5 m)
0.1C m/s
Fig. 23. Sidewall eects on the lateral motions in 3D over grids (205 21 40) for Case 5 (regular waves): the instantaneous velocity elds
varying with three cross-sections (a xed time level).
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 891
there is no obvious sign for the strong coupling between the transverse and longitudinal dynamics, due to the
weak solution of the lateral motions (see Fig. 23). Consequently, the eect of the longitudinal scale on the hor-
izontal velocity is dominant over that of the cross-sectional part. As a highly turbulent bore generated prop-
agates along the beach, the obliquely descending eddies, found by Nadaoka et al. [41] in experiments, can also
be described by the velocity vectors in the plane (see Fig. 24: t = 12.0 s). Such turbulent structures in 3D
behave a pair of counter-rotating vortices, forming the three-dimensional vorticity elds but decay with time
or destroy. Interestingly, the inuence of walls may not be only the three-dimensional forcing from which the
instability develops [19], as for waves on beaches high turbulence in 3D can be induced by the breaking waves.
Physically, the near-wall motions are substantially suppressed by the sidewall boundary-layer thickness,
while varying with the addition of even the small current (see Fig. 24). Suciently far from the walls (perfectly
in the center of plane, y = 0), the signicant inuence of the sidewall does not detected. Anyway, this eects
could be further reduced by use of a wider ume, especially the width of a tank becomes innite. Note that the
ows become three-dimensional (although weak) at toe of the dike (see Fig. 23: x = 1.0 m), due to the longi-
tudinal modulation of the bottom topography. In fact, only within the breaking region, the motions are sig-
nicantly three-dimensional (see Fig. 23: x = 4.5 m), where the surface prole shows pronounced local lateral
variability, and as a result of advection of a larger-scale three-dimensional feature from the swash zone.
4. Conclusions
By design of the dierent wavecurrent generators applied to the NWT, we study interactions of breaking
waves with currents over a structure, based on a dynamic Smagorinsky model. Computations are conducted
on Cartesian cut-cell grids with our developed NavierStokes solver. By superimposing a steady current onto
waves, our external generator technique provides well-dened wavy inow conditions, which have been
0.5
0.25
0
-0.25
-0.5
3 2 1 0
y

(
m
)
x (m)
Uc=0.0 m/s t = 11.5 (s) 1.0 m/s
0.5
0.25
0
-0.25
-0.5
3 2 1 0
y

(
m
)
x (m)
Uc=0.17 m/s t = 11.5 (s) 1.0 m/s
0.5
0.25
0
-0.25
-0.5
3 2 1 0
y

(
m
)
x (m)
Uc=0.0 m/s t = 12.0 (s) 1.0 m/s
0.5
0.25
0
-0.25
-0.5
3 2 1 0
y

(
m
)
x (m)
Uc=0.17 m/s t = 12.0 (s) 1.0 m/s
0.5
0.25
0
-0.25
-0.5
3 2 1 0
y

(
m
)
x (m)
Uc=0.0 m/s t = 12.5 (s) 1.0 m/s
0.5
0.25
0
-0.25
-0.5
3 2 1 0
y

(
m
)
x (m)
Uc=0.17 m/s t = 12.5 (s) 1.0 m/s
Fig. 24. Sidewall eects on the longitudinal motions in 3D over grids (205 21 40) for Case 5 (regular waves): the instantaneous velocity
elds varying with three time levels (a xed vertical position of z = 0.28 m).
892 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
addressed in the study of waves following a current. For waves propagating against strong currents, the inter-
nal source function introduced is used as one type of a net driving force for creation of the desirable opposing
current, forming our internal generator.
Our solver has be extensively illustrated with ve test cases studied. These are of particular values in dis-
tinctly dierent applications, especially in the area of nonlinear shallow water-wave propagation problems,
signicantly characterized by overtopping of waves. Given nite amplitude surface water waves, we concern
a direct comparison of the numerical simulation with the experimental ows varying with dierent currents
(positive, zero, negative). By validation, the results in 2D or 3D agree well with measurements available in
terms of the wave or horizontal velocity proles. Hopefully, our study can provide valuable base for the anal-
ysis of wavecurrent-structure coupling, induced by regular and irregular waves. The major conclusions are
drawn as follows:
The grid dependency is still a crucial aspect for our VOF-based solver. Near the free surface, for example,
the minimum grid space (in the vertical direction) varies with case by case. When validated the results, visu-
alization of the wave pattern plus the instantaneous velocity elds are both important, based on our
experience.
By using a breaking-type wave absorber, the reected waves can be eectively dissipated without loss of the
ow characteristics of interest, which diers substantially from various theoretical schemes proposed for
absorbing waves reected at the outow boundary. In this sense, our computations are without the insta-
bility problems even for a long duration. This is very useful for situations involving some actual cases, as
the lengthy computations are necessary, especially such information helps to conrm whether the global
mass and energy are conserved.
Under the test for waves against the adverse current, the dynamics of ows could more dramatically vary,
as in this case the waves shorten the wavelength, relative to waves superimposed on the positive current.
Owing to existence of a transition zone (related with two dierent roughness values caused by the actual
bed and free surface), our results reveal that the ow is increased in the outer layer after a drop in the log-
arithmic description of the velocity proles, vanishing at the free surface. This is true for waves following a
current but in the situation of opposing waves it appears that growth of the mean horizontal velocity is
underestimated, as compared to measurements available.
Owing to the lack of the spanwise convective rolls, the physical dissipation for the two-dimensional ow is
substantially limited in the case of waves breaking. On the other hand, wave dynamics remains two-dimen-
sional in regions of non-breaking waves but the dominant motions in the surf zone are obviously a three-
dimensional convective instability that leads to spanwise convective rolls. It is clear that the disturbance
caused by 3D enhances the strength of this instability at the expense of waves breaking.
Ongoing and future work will incorporate the particular adverse current case for the study of wave
blocking.
Acknowledgments
The present work is supported by the Funds for Scientic Research Flanders, Belgium (Project No.
G.0199.06). Additionally, we also thank the two referees for their helpful suggestions.
Appendix A.
A.1. Eects of a dynamic SGS model on MILES in breakwaters
More recently, irrotational ow models are still widely applied for wavecurrent interactions, in case the
eects of viscosity and vorticity are weak, as studied in [8,55] with the Boussinesq-type equations. Under
breaking waves, however, the ows become complex unsteady ows. To confront the scale-complexity prob-
lem inherent to turbulent ows, the technique of LES has emerged as an alternative to the RANS approach
[14]. In the context of this theory, MILES probably provides a promising practical tool due to without any
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 893
further modelling: it lets the numerical viscosity dampen the numerical simulation, by introducing high-reso-
lution locally monotone algorithms [15]. More specically, an elementary approach in MILES is to discretize
the convective terms, based on the schemes like Jameson, MUSCL and ENO [9,16]. But a hybrid approach
(called MILES + SGS in our study), constructed by MILES that is coupled with a dynamic SGS model, is
essential for accurate estimation of periodic breaking water-wave propagation problems, especially for the
lengthy computations, as the SGS model incorporates a model bringing a turbulent viscosity. The discussion
regarding the issues is given as follows.
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 3.5 s (MILES)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 3.5 s (MILES+SGS)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 15.4 s (MILES)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 15.4 s (MILES+SGS)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 32.0 s (MILES)
1.0 m/s
-0.7
-0.5
-0.3
-0.1
0
0.1
0.2
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 32.0 s (MILES+SGS)
1.0 m/s
-0.15
-0.08
0
0.08
0.15
0 5 10 15 20 25 30 35

(
m
)
t (s)
WG: x=2.02 (m)
-0.16
-0.08
0
0.08
0.16
0 5 10 15 20 25 30 35

(
m
)
t (s)
WG: x=3.81 (m)
MILES+SGS MILES Exp.
Fig. 25. Comparison between MILES and MILES + SGS over grids (251 40) for Case 5 (pure regular waves): the instantaneous velocity
elds and the wave trains at two certain positions.
894 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
Flows behave as localized high turbulence in waves breaking, as given with two typical case studies related
to coastal structures (see Figs. 25 and 26), leading to the rapid increase in high turbulent kinetic energy. But
the energy of waves accumulated at the high wavenumbers can be eectively dissipated by turbulence, when
using a dynamic SGS model. Otherwise, the spurious currents are probably generated only with MILES, and
as a result of the damage to true solution at the same grid. Hopefully, this provides one example: MILES can-
not resolve LES properly [46], in case the principal mechanism of turbulence diusion is through the action of
viscous on the breaking waves.
When the ows become self-sustaining turbulence, the wave-induced motions evolve into the following two
distinct parts: (i) In the early stage, the viscous eects may be regarded as negligible, and the results look very
-0.6
-0.4
-0.2
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 4.1 (s)
MILES
1.0 m/s
-0.6
-0.4
-0.2
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 4.1 (s)
MILES+SGS
1.0 m/s
-0.6
-0.4
-0.2
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 4.9 (s)
MILES
1.0 m/s
-0.6
-0.4
-0.2
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 4.9 (s)
MILES+SGS
1.0 m/s
-0.6
-0.4
-0.2
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 12.9 (s)
MILES
1.0 m/s
-0.6
-0.4
-0.2
0
0.1
0.2
0.3
0 1 2 3 4 5 6
z

(
m
)
x (m)
t = 12.9 (s)
MILES+SGS
1.0 m/s
-0.4
-0.2
0
0.15
0.3
0 5 10 15 20 25 30 35 40

(
m
)
t (s)
WG: x=3.81 (m)
-0.25
-0.15
0
0.1
0 5 10 15 20 25

(
m
)
t (s)
WG: x=5.2 (m)
MILES+SGS
MILES
Fig. 26. Comparison between MILES and MILES + SGS over grids (251 40) for Case 4 (pure irregular waves): the instantaneous
velocity elds and the wave trains at two certain locations.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 895
close except for a small phase shift (also see [30] with our shock-capturing Euler solver), as in this case the
leading-order truncation error (growing with increasing wavenumber) acts as the intrinsic SGS model that dis-
sipates the energy of waves. It means that our MILES approaches capture well growth rates of the initial per-
turbation, caused in waves breaking. Consequently, the total contribution of the turbulence increase can be
almost balanced under the numerical dissipation mechanism. (ii) Nevertheless, dierent features are apparent
in waves continuously breaking, and viscous diusion plays an important role in the ow regime eventually
dominated by turbulence. After the initiation of progressive waves breaking, the performance of MILES is
far from satisfactory in terms of the surface proles (as a function of time), especially for the oweld. A
major reason is that the numerical viscosity in MILES cannot restrain the increase of energy uctuated on
a broad range of length and time scales, partly due to insucient grid resolution or the frequent occurrence
of small-scale structures or both. On the contrary, the hybird scheme maintains a transfer of energy in the
nearly equilibrium state of turbulence eld, by modelling of small-scale turbulent uctuations, responsible
for the kinetic energy dissipation. MILES + SGS is probably a more viable approach for breaking water-wave
problems, therefore, based on the fact that SGS models provide an additional physics mechanism, by which
dissipation of a rapid build-up of kinetic energy at high wavenumber can occur [16]. This behavior is partic-
ularly important in waves breaking.
References
[1] K. Anthony, P. Ugo, B. Elias, A priori and a posteriori tests of inow conditions for large-eddy simulation, Phys. Fluids 16 (3) (2004)
46964712.
[2] J.A. Battjes, Surf-zone dynamics, Ann. Rev. Fluid Mech. 20 (1988) 257293.
[3] S. Belcher, J. Hunt, Turbulent ows over hills and waves, Ann. Rev. Fluid Mech. 30 (1998) 507538.
[4] N. Booij, R.C. Ris, L.H. Holthuijsen, A third-generation wave model for coastal regions, Part I, model description and validation, J.
Geophys. Res. 104 (C4) (1999) 76497666.
[5] M.S. Celebi, M.H. Kim, R.F. Beck, Fully nonlinear 3-D numerical wave tank simulation, J. Ship Res. 42 (1998) 3345.
[6] A. Clement, Coupling of two absorbing boundary conditions for 2D time-domain simulations of free surface gravity waves, J.
Comput. Phys. 126 (1996) 139151.
[7] A. Chawla, J.T. Kirby, Monochromatic and random wave breaking at blocking points, J. Geophys. Res. 107 (2002) 119.
[8] Q. Chen, Per A. Madsen, H.A. Schaer, D.R. Basco, Wavecurrent interaction based on an enhanced Boussinesq approach, Coast.
Eng. 33 (1998) 1139.
[9] F.K. Chow, P. Moin, A further study of numerical error in large-eddy-simulation, J. Comput. Phys. 184 (2003) 366380.
[10] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulations of waves, Math. Comput. 31 (1977) 629651.
[11] J.H. Ferziger, Large eddy simulation, in: M.Y. Hussaini, T. Gatski (Eds.), Simulation and Modelling of Turbulent Flows, Cambridge
University Press, New York, 1995.
[12] J. Fredse, K.H. Andersen, B.M. Sumer, Wave plus current over a ripple-covered bed, Coast. Eng. 38 (1999) 177221.
[13] O.B. Fringer, R.L. Street, The dynamics of breaking progressive interfacial waves, J. Fluid Mech. 494 (2003) 319351.
[14] U. Frisch, Turbulence, Cambridge University Press, Cambridge, UK, 1995.
[15] C. Fureby, F.F. Grinstein, Large eddy simulation of high-Reynolds-number free and wall-bounded ows, J. Comput. Phys. 181
(2002) 6897.
[16] E. Garnier, M. Mossi, P. Sagaut, P. Comte, M. Deville, On the use of shock-capturing schemes for large-eddy simulation, J. Comput.
Phys. 153 (1999) 273311.
[17] M. Germano, U. Piomelli, P. Moin, W. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids 3 (1991) 180200.
[18] S. Ghosal, T.S. Lund, P. Moin, K. Akselvoll, A dynamic localization model for large-eddy simulation of turbulent ows, J. Fluid
Mech. 286 (1995) 229255.
[19] F. Gheusi, J. Stein, O.S. Ei, A numerical study of three-dimensional orographic gravity-wave breaking observed in hydraulic tank, J.
Fluid Mech. 410 (2000) 6799.
[20] W.D. Grant, O.S. Madsen, Combined wave and currents interaction with a rough bottom, J. Geophys. Res. 84 (1979) 17971808.
[21] J. Groeneweg, J.A. Battjes, Three-dimensional wave eects on a steady current, J. Fluid Mech. 478 (2003) 325343.
[22] Z. Huang, C.C. Mei, Eects of surface waves on a turbulent current over a smooth or rough seabed, J. Fluid Mech. 497 (2003) 253
287.
[23] K.M.T. Kleefsman, G. Fekken, A.E.P. Veldman, B. Iwanowski, B. Buchner, A volume-of-uid based simulation method for wave
impact problem, J. Comput. Phys. 206 (2005) 363393.
[24] M. Israeli, S.A. Orszag, Approximation of radiation boundary condition, J. Comput. Phys. 41 (1981) 115135.
[25] F.J. Kelecy, R.H. Pletcher, The development of a free surface capturing approach for multidimensional free surface ows in closed
containers, J. Comput. Phys. 138 (1997) 939980.
[26] P.H. Kemp, R.R. Simons, The interaction between waves and a turbulent current: waves propagating with the current, J. Fluid Mech.
116 (1982) 227250.
896 T. Li et al. / Journal of Computational Physics 223 (2007) 865897
[27] G. Klopman, Secondary circulation of the ow due to waves and currents: laser-Doppler ow measurements for waves following or
opposing a current, Technical Report Z2249, Delft Hydraulics, 1997.
[28] M. Lesieur, O. Metails, New trends in large eddy simulations of turbulence, Annu. Rev. Fluid Mech. 28 (1996) 145.
[29] D.K. Lilly, A proposed modication of the Germano subgrid-scale closure method, Phys. Fluids A 4 (1992) 633635.
[30] T. Li, P. Troch, J. De Rouck, Wave overtopping over a sea dike, J. Comput. Phys. 198 (2004) 686726.
[31] T. Li, P. Troch, J. De Rouck, D. Goossens, Numerical simulation of water wave impacts using a NavierStokes solver, in:
Proceedings of the 29th International Conference on Coastal Engineering, vol. 4, 2004, pp. 41004112.
[32] T. Li, Computation of turbulent free-surface ows around modern ships, Int. J. Numer. Meth. Fluids 43 (2003) 407430.
[33] T. Li, A LVOF user guide, Technical Report, Department of Civil Engineering, Ghent University, Belgium, 2006.
[34] T. Li, P. Troch, J. De Rouck, On the study of wavecurrent-structure interactions in a numerical wave tank, in: Proceedings of the
16th International Oshore and Polar Engineering Conference, vol. 3. 2006. pp. 356362.
[35] P. Lin, P.L.F. Liu, Internal wave-maker for NavierStokes equations models, J. Waterway Port Coast Ocean Eng. 4 (1999) 207215.
[36] D. Lo rstad, L. Fuchs, High-order surface tension VOF-model for 3D bubble ows with high density ratio, J. Comput. Phys. 200
(2004) 153176.
[37] M.S. Longuet-Higgins, R.W. Stewart, The changes of amplitude of short gravity waves on steady non-uniform currents, J. Fluid.
Mech. 10 (1961) 529549.
[38] P. Lubin, S. Vincent, S. Abadie, J.P. Caltagirone, Three-dimensional large eddy simulation of air entrainment under plunging
breaking waves, Coast. Eng. 52 (2006) 631655.
[39] C. Meneveeau, T.S. Lound, W.H. Cabot, 1A Lagrangian dynamic subgrid-scale model of turbulence, J. Fluid Mech. 319 (1996) 353
385.
[40] F.M. Najjar, D.K. Tafti, Study of discrete test lters and nite dierence approximations for the dynamic subgrid-scale model, Phys.
Fluids 8 (1996) 10761088.
[41] K. Nadaoka, M. Hino, Y. Koyano, Structure of the turbulent ow eld under breaking waves in the surf zone, J. Fluid Mech. 204
(1989) 359387.
[42] I. Orlanski, A simple boundary condition for unbounded hyperbolic ows, J. Comput. Phys. 21 (1976) 251269.
[43] J.C. Park, M.H. Kim, H. Miyata, Three-dimensional numerical wave tank simulations on fully nonlinear wavecurrentbody
interaction, J. Mar. Sci. Technol 6 (2001) 7082.
[44] D.H. Peregrine, Interaction of water waves and currents, Adv. Appl. Mech. 16 (1976) 9117.
[45] U. Piomelli, J. Liu, Large-eddy simulation of rotating channel ows using a localized dynamic model, Phys. Fluids 5 (7) (1995) 839
848.
[46] P. Sagaut, Large Eddy Simulation for Incompressible Flows, third ed., Springer, 2005.
[47] P. Schlatter, N.A. Adams, L. Kleiser, A windowing method for periodic inow/outow boundary treatment of non-periodic ows, J.
Comput. Phys. 206 (2005) 505535.
[48] J.U. Schlu ter, H. Pitsch, P. Moin, Large eddy simulation inow conditions for coupling with Reynolds-averaged ow solvers, AIAA
42 (3) (2004) 478484.
[49] L. Shen, D.K.P. Yue, Large-eddy simulation of free-surface turbulence, J. Fluid Mech. 440 (2001) 75116.
[50] J. Smagorinsky, General circulation experiments with the primitive equations, Part I: the basic experiment, Mon. Weather Rev. 91
(1962) 99152.
[51] O.V. Vasilyev, T.S. Lund, P. Moin, A general class commutative lter for LES in complex geometries, J. Comput. Phys. 146 (1998)
82104.
[52] B. Vreman, B. Geurts, H. Kuerten, On the formulation of the dynamic mixed subgrid-scale model, Phys. Fluids 6 (12) (1994) 4057
4059.
[53] Y. Watanabe, H. Saeki, Three-dimensional large eddy simulation of breaking waves, Coast. Eng. 41 (1999) 281301.
[55] S.B. Yoon, P.L.F. Liu, Interaction of currents and weakly nonlinear water waves in shallow water, J. Fluid Mech. 205 (1989) 397419.
T. Li et al. / Journal of Computational Physics 223 (2007) 865897 897