In the literature one of the most used turbulence
models for the simulations of the breaking waves [1],[2],[3],
[4],[5] and turbulent phenomena is the Smagorinsky
model. This model is characterized by the presence of
the Smagorinsky coefficient, which is defined (in the
previous paper) without distinction in the shoaling zone, in
the region around the breaking point and in the surf zone.
In the models present in the literature the turbulent
phenomena are underestimated by the Smagorinsky model
and the dissipation of the averaged motion is due to the
numerical scheme.
Another commune turbulence model used in the
literature in the context of wave breaking [6],[7] is the
model ( is the turbulent kinetic energy and is the mixing
length), in which there is no distinction of the representation
of the production and dissipation of the turbulent kinetic
energy before and after the breaking point.
In the simulations presented in the literature, both
the abovementioned models do not take into account the
turbulent phenomena at the bottom boundary layer (buffer
layer and turbulent core), because the first calculation grid
cell is outside the boundary layer, where these phenomena
are dominant.
In this paper a new is proposed that is able to
dissipate the average motion energy and overcomes the
limitations of the abovementioned turbulence model present
in the literature. To take into account the different behavior
of the turbulence in the wave propagation, the mixing length
is a function of the first and second spatial derivatives of the
maximum water surface elevation.
In this numerical model the Navier-Stokes
equations, written in integral contravariant form,
expressed in a generalized time-dependent curvilinear
coordinate system in
which the vertical coordinate moves following the
free surface, are numerically solved. The numerical scheme
uses 5th order WTENO (Wave Targeted Essentially
Non-Oscillatory) reconstruction technique (from the cell
average value of the variables to the point ones) and an
Exact Riemann
Solver. The WTENO technique is collocated in the context
of the TENO and ENO reconstruction technique [8],[9], but
applicated to the waves. The Poisson equation and the motion
equations are expressed in terms of the conserved variables
( e ).
The three-dimensional Navier-Stokes equations without
the Christoffel symbols and the equation relative to the
movement of the free surface are expressed in integral
contravariant form in a time-dependent curvilinear coordinate
system, that adopt conserved variables e ,
and are given by
(1)
(2)
The transformation from the Cartesian coordinate system
to the curvilinear one is given by
(3)
Breaking wave simulations by a new 𝑘−𝑙 turbulence model
FRANCESCO GALLERANO, BENEDETTA IELE, FEDERICA PALLESCHI, GIOVANNI CANNATA
Department of Civil, Constructional and Environmental Engineering Sapienza University of Rome Rome, ITALY
Abstract: The three-dimensional motion equations are used to simulate the wave and velocity fields. These equations are
written in integral contravariant form on a time-dependent curvilinear coordinate system. In this paper a new 𝒌−𝒍 turbulence
model in contravariant form is proposed for three-dimensional simulation of breaking waves. In this model the mixing
length is defined as a function of the first and second spatial derivatives of the maximum water surface elevation.
Keywords: Wave breaking; turbulence model; turbulent kinetic energy; boundary layer; boundary conditions.
Received: May 29, 2021. Revised: June 18, 2022. Accepted: July 17, 2022. Published: September 13, 2022.
1. Introduction
2. Motion Equations
WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS
DOI: 10.37394/232011.2022.17.15
Francesco Gallerano, Benedetta Iele,
Federica Palleschi, Giovanni Cannata
E-ISSN: 2224-3429
113
Volume 17, 2022
In (3) in the undisturbed water depth,
is the total water depth and is the
free surface elevation with respect to the undisturbed water
depth. In (1) and (2)
and
are the cell average values of
the conserved variables. and are
respectively the contravariant components of the fluid
velocity and the velocity of the moving coordinates.
and
are the contravariant and
covariant base vectors. , in which
,
indicates the vertical unit vector and
indicates the vector product is the Jacobian of the
transformation. is the covariant base vector defined at the
center of the control volume. The volume
and the areas
are not time dependent. The
boundary surfaces of the control volume , on which the
coordinate is constant, that surfaces are in correspondence
of the larger and smaller value of (, and are cyclic),
are
and
. the contour lines of the surface element
, on which is constant and are located
respectively in correspondence of the larger and smaller value
of , are e (with . are the
contravariant components of the stress tensor that contains the
dynamic pressure term, is the fluid density and is the
gravity acceleration.
For obtaining a divergence-free velocity field the
advancing in time of the numerical solution a predictor-
corrector procedure is used. To an approximated not
divergence-free velocity field is added a corrector
divergence-free velocity field, , that takes into account
the dynamic pressure.
(4)
is the potential scalar, obtained by the solution of the
Poisson equation given by
(5)
(6)
is the dimension of the grid size and is
the abovementioned Smagorinsky coefficient.
In the lower part, close to the bottom, and in the higher part
of the boundary layer there is the dominance respectively of
viscous stresses and turbulent stresses. In the middle part of
the boundary layer viscous and turbulent stresses coexist.
The motion equations, in which turbulent stress tensor is
modelled by Smagorinsky, are solved in the turbulent core. In
this model in the turbulent boundary layer the turbulent eddy
viscosity is given by
(7)
in which is the von Kàrmàn constant.
In the model, in which is the turbulent kinetic energy
and is the mixing length, the turbulent kinetic energy
equation is solved. This equation is expressed in integral
contravariant form in a time-dependent curvilinear
coordinate system.
(8)
in which
is the dissipation of turbulent kinetic
energy, is the production of turbulent kinetic energy. The
turbulent eddy viscosity is given by
with
(9)
In the new model the mixing length is a function of
the first and second spatial derivatives of the maximum water
surface elevation,
and
.
Let be the maximum
wave height point by point in the time;
be the maximum total water depth point by
point in the time;
and
be respectively the maximum and minimum
water surface elevation.
(10)
in which the coefficients
;
;
;
(11)
3. Turbulence Model
In the following paragraph the Smagorinsky model is
presented and the new proposed model is explained.
3.1 Smagorinsky Turbulence Model
Many authors [1],[2],[3],[4],[5] use Smagorinsky model
to represent the turbulence in the breaking waves. In
turbulent stress tensor, , where is the
strain rate tensor, the turbulent eddy viscosity is given by
3.2 𝑘−𝑙 Turbulence Model
WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS
DOI: 10.37394/232011.2022.17.15
Francesco Gallerano, Benedetta Iele,
Federica Palleschi, Giovanni Cannata
E-ISSN: 2224-3429
114
Volume 17, 2022
are the multiplier of the undisturbed water depth .
The first derivative, , allows to
differentiate the behavior of the turbulence model before, after
and around the breaking point. Before the breaking point this
derivative is positive and after is negative. The effects of the
diffusive terms in the motion direction of the momentum can
be reduced in the zone immediately after the breaking point
by this new mixing length.
Let define two different configurations for the boundary
conditions of this model.
The main difference between the two configurations is that
the second configuration solves the motion equations and the
turbulent kinetic energy equation also in the buffer layer
unlike the first configuration. In the first configuration these
equations are solved in the turbulent core. Between the buffer
layer and the turbulent core, it is known that the balance
between production and dissipation of turbulent
kinetic energy is valid [10],[11]. The boundary conditions
for the first configuration are given by the (7) for the
turbulent eddy viscosity and for the turbulent kinetic
energy outside the buffer layer is
(12)
The turbulent phenomena and the distribution of the
turbulent kinetic energy is influenced by the production of
turbulent kinetic energy, in the buffer layer and the turbulent
core of the bottom boundary layer this production is high.
The second configuration is characterized by the fact that
the motion equations and the turbulent kinetic energy equation
is solved also in the buffer layer (near the viscous sublayer),
as already said, to correctly represent the strong variability of
the production of turbulent kinetic energy closer to the bottom.
The boundary condition for the turbulent kinetic energy is zero
to the bottom, the turbulent eddy viscosity and the mixing
length in the boundary layer are given respectively by (7) and
(13).
(13)
Fig. 1. Computational domain for simulation of Ting and Kirby’s test case
[12].
It has been considered three configurations.
In the first configuration, named C1, Smagorinsky
turbulence model is adopted. The first calculation grid cell for
the motion equations is in the turbulent core. The velocity
boundary condition and the eddy viscosity in (7) are placed at
the border between the buffer layer and the turbulent core,
( is the adimensionless wall distance, , and
are the vertical distance away from the wall, the bottom
friction velocity calculated by a logarithmic law [11] and the
kinetic viscosity coefficient).
In the second configuration, named C2, the new
turbulence model is adopted. In this configuration, the first
calculation grid cell in which the motion equations are solved,
is in the turbulent core, and the lower face of this cell is at the
border between the buffer layer and the turbulent core,
. The velocity boundary condition is assigned identically to
the ones in C1. In addition, the boundary condition for the
turbulent kinetic energy (12) is placed at the upper part of the
turbulent core, .
The third configuration, named C3, uses the new
turbulence model, but the lower face of the first calculation
grid cell is placed at the border between the viscous sublayer
and the buffer layer, . The first calculation grid cell
for the motion equations and for the turbulent kinetic energy
equation is in the buffer layer.
The numerical results obtained with the configuration C1
are shown in Fig. 2
Fig. 2. Configuration C1: maximum water surface elevation.
Numerical results (dotted line , solid line , dashed
line ) and experimental measurements (square) [12].
The breaking point has various positions in the three
simulations obtained with different Smagorinsky coefficients;
it is postponed in the blue line and it is anticipated in the red
one. The maximum water surface elevation is overestimated
in the simulation obtained with (blue line) and it is
underestimated in the one obtained with (red line).
The dissipation of the energy of the averaged motion is greater
in the simulations with a high value of the Smagorinsky
coefficient (, red solid line) and is lower in the ones
with a small value of the same coefficient (, blue
solid line). The increase of the Smagorinsky coefficient
produced an increase of the turbulent eddy viscosity that
reduces the wave height, because the diffusion in the motion
direction of the momentum is greater. It is evident that the
abovementioned coefficient influences a lot the numerical
simulations and the prediction of the turbulent phenomena in
the surf zone.
To overcome the limitation of the Smagorinsky model a
new turbulence model is presented, the numerical
simulations of which produced the following results. In Fig. 3
. 5HVXOWV
In this paragraph the numerical results obtained with two
different turbulence models are compared with the
experimental ones obtained by Ting and Kirby [12] to
validate the new model. The case study consists in a spilling
breaker wave on a slope channel. The channel is represented
in Fig. 2. The parameters of the cnoidale wave reproduced
are: undisturbed water dept of , initial wave
height and period . The simulations are
made by grid points with and layers in
vertical direction (Fig. 1).
WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS
DOI: 10.37394/232011.2022.17.15
Francesco Gallerano, Benedetta Iele,
Federica Palleschi, Giovanni Cannata
E-ISSN: 2224-3429
115
Volume 17, 2022
the numerical results obtained in the configuration C2 are
compared with the experimental ones.
Fig. 3. Configuration C2: maximum water surface elevation.
Numerical results solid line and experimental measurements
(square) [12].
The new model overcomes the limitations of the
Smagorinsky model. By comparing the results in Fig. 2 and
Fig. 3, it can be noticed that the new 𝑘 − 𝑙 model well predicts
the breaking point and slightly underestimates the maximum
water surface elevation around the breaking point. The
steepness of the black line immediately after the breaking
point is in good accordance with the experimental results,
because the mixing length in this zone is reduced and so
reduces the diffusion in the motion direction of the
momentum.
The time mean turbulent kinetic energy profiles at 𝑥 =
7.27𝑚 e 𝑥 = 7.88𝑚, obtained in the configuration C2, are
shown in Fig. 4.
(a)
(b)
Fig. 4. Configuration C2: vertical distribution of the time mean
turbulent kinetic energy. Numerical results (triangle) and
experimental results (square) [12]. (a) 𝑥 = 7.27𝑚, (b) 𝑥 = 7.88𝑚 .
The vertical coordinate is 𝑍 = (𝑧 − 𝜂)⁄𝐻 and the horizontal one is
𝐾 = √𝑘/(𝑔𝐻).
This numerical simulation does not produce results in
good agreement with the numerical results in terms of
time mean turbulent kinetic energy profiles, as shown in
Fig. 4.
The production of turbulent kinetic energy, that is
generated in the buffer layer and in the turbulent core, is
not
well predicted by this numerical model (configuration C2),
because the first calculation grid cell for the turbulent kinetic
energy is in the upper part of the turbulent core. For this reason
the variability of the production of turbulent kinetic energy
along the vertical direction cannot be taken into account.
The configuration C3 solves the turbulent kinetic energy
equation and the motion equations also in the buffer layer to
take into account the variability of the vertical distribution of
turbulent kinetic energy and so the turbulent phenomena.
Fig. 5. Configuration C3: maximum water surface elevation.
Numerical results solid line and experimental measurements
(square) [12].
(a)
(b)
Fig. 6. Configuration C3: vertical distribution of the time mean
turbulent kinetic energy. Numerical results (triangle) and
experimental results (square) [12]. (a) 𝑥 = 7.27𝑚, (b) 𝑥 = 7.88𝑚 .
The vertical coordinate is 𝑍 = (𝑧 − 𝜂)⁄𝐻 and the horizontal one is
𝐾 = √𝑘/(𝑔𝐻).
The maximum water surface elevation and the time mean
turbulent kinetic energy distribution along the vertical
direction are in good agreement with the experimental
results as it can see in Fig. 5 and Fig. 6. By placing the
lower face at the first calculation grid cell at 𝑦+ = 10 and by
solving the al the equations in the buffer layer, the turbulent
phenomena are well represented.
In Fig. 7 an instantaneous wave and velocity field is shown
(𝑇 = 100𝑠). One vector of each three is shown. In the zone
after the breaking point the turbulent eddy viscosity, shown by
the contour, is reduced by the use of the new formula for the
mixing length.
WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS
DOI: 10.37394/232011.2022.17.15
Francesco Gallerano, Benedetta Iele,
Federica Palleschi, Giovanni Cannata
E-ISSN: 2224-3429
116
Volume 17, 2022
Fig. 7. Configuration C3: Instantaneous wave fields. Contour of turbulent eddy viscosity. Vector field in which one vectors of every three are shown.
.
In this paper a new turbulence model in
contravariant form is proposed for three-dimensional
simulation of breaking waves. In this model the mixing length
is defined as a function of the first and second spatial
derivatives of the maximum water surface elevation. It has
been demonstrated that the new turbulence model, in which
the motion equations and the turbulent kinetic energy equation
are solve also in the buffer layer, produces results in good
agreement with the experimental ones. The new model is able
to represent the maximum water surface elevation and the time
mean turbulent kinetic energy distribution along the vertical
direction
REFERENCES
[1] G. Ma, F. Shi, J.Y. Kirby, “Shock-capturing non-hydrostatic
model for fully dispersive surface wave processes”, Ocean Model, vol.
43-44, pp. 22-35, 2012.
[2] G. Cannata, C. Petrelli, L. Barsi, F. Gallerano, “Numerical integration
of the contravariant integral form of the Navier–Stokes equations in
time-dependent curvilinear coordinate systems for three-dimensional
free surface flows”, Contin. Mech. Thermodyn. vol. 31, pp. 491–519,
2019.
[3] B. Iele, F. Palleschi, “Boundary conditions for the simulation of wave
breaking”, WSEAS Transac. Fluid Mech., vol. 15, pp. 41-53, 2020.
[4] K.Z. Fang, Z.B Liu, “Modeling breaking waves and wave-induced
currents with fully nonlinear Boussinesq equations”, WSEAS Transac.
Fluid Mech., vol. 9, pp. 131-143, 2014.
[5] F. Palleschi, B. Iele, F. Gallerano, “Integral contravariant form of the
Navier-Stokes equations”, WSEAS Transac. Fluid Mech., vol. 14, pp.
101-113, 2019.
[6] S. F. Bradford, “Numerical simulation of surf zone dynamics”, J.
Waterw. Port, Coast. Ocean Eng., vol. 126, pp. 1-13, 2000.
[7] F. Gallerano, G. Cannata, L. Barsi, F. Palleschi, B. Iele, “Simulation of
wave motion and wave breaking induced energy dissipation”, WSEAS
Transac. Fluid Mech., vol. 14, pp. 62-69, 20191.
[8] G.Jiang, C. Shu, “Efficient Implementation of Weighted ENO
Schemes”, J. Comput. Phys., vol. 126, pp. 202-228, 1996.
[9] J. Peng, S. Liu, S. Li, K. Zhang, Y. Shen, “An efficient targeted ENO
scheme with local adaptive dissipation for compressible flow
simulation”, J. Comput. Phys., vol. 425, pp. 1-25, 2021.
[10] F. Gallerano, G. Cannata, “Noll’s axioms and formulation of the
closure relations for the subgrid turbulent tensor in Large Eddy
Simulation”, WSEAS Transac. Fluid Mech., vol. 15, pp. 85-90, 2020.
[11] P.L.F. Liu, P. Lin, “A numerical model for breaking waves: the volume
of fluid method”, Res Report No. CACR-97-02, pp. 1-54, 1997.
[12] F.C.K. Ting, J.T. Kirby, “Observation undertow and turbulence in a
wave period”, Coast. Eng., vol. 24, pp. 177–204, 1995.
. Conclusion
WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS
DOI: 10.37394/232011.2022.17.15
Francesco Gallerano, Benedetta Iele,
Federica Palleschi, Giovanni Cannata
E-ISSN: 2224-3429
117
Volume 17, 2022
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The author(s) contributed in the present research,
at all stages from the formulation of the problem
to the final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
No funding was received for conducting this study.
Conflict of Interest
The author(s) declare no potential conflicts of
interest concerning the research, authorship, or
publication of this article.
Creative Commons Attribution
License 4.0 (Attribution 4.0
International, CC BY 4.0)
This article is published under the terms of
the Creative Commons Attribution License
4.0
https://creativecommons.org/licenses/by/4.0/deed.en_US
