A three-dimensional high-order numerical model for the simulation of
the interaction between waves and an emerged barrier
FRANCESCO GALLERANO, FEDERICA PALLESCHI, BENEDETTA IELE, GIOVANNI
CANNATA
Department of Civil, Constructional and Environmental Engineering
Sapienza University of Rome
Via Eudossiana 18, 00184, Rome
ITALY
Abstract: - We present a new three-dimensional numerical model for the simulation of breaking waves. In the
proposed model, the integral contravariant form of the Navier-Stokes equations is expressed in a curvilinear
moving coordinate system and are integrated by a predictor-corrector method. In the predictor step of the
method, the equations of motion are discretized by a shock-capturing scheme that is based on an original high-
order scheme for the reconstruction of the point values of the conserved variables on the faces of the
computational grid. On the cell faces, the updating of the point values of the conserved variables is carried out
by an exact Riemann solver. The final flow velocity field is obtained by a corrector step which is based
exclusively on conserved variables, without the need of calculating an intermediate field of primitive variables.
The new three-dimensional model significantly reduces the kinetic energy numerical dissipation introduced by
the scheme. The proposed model is validated against experimental tests of breaking waves and is applied to the
three-dimensional simulation of the local vortices produced by the interaction between the wave motion and an
emerged barrier.
Key-Words: - breaking waves, wave-structure interaction, vortices, high-order scheme, exact Riemann solver,
three-dimensional simulation
Received: June 9, 2021. Revised: May 21, 2022. Accepted: June 19, 2022. Published: July 18, 2022.
1 Introduction
The three-dimensional numerical simulation of
gravity wave motion can be carried out by
numerically solving the Navier-Stokes Equations. In
the numerical simulation of the Navier-Stokes
Equations for free-surface flows, one of the most
challenging problems is the calculation of the free
surface elevation. In this context, some of the first
numerical models proposed in the literature [1,2] are
based on a technique called volume of fluid (VOF)
method. The VOF method is used to simulate the
motion of two fluids (air and water) by solving a
momentum balance equation for the mixture of the
two fluids on a fixed Cartesian grid. The volume
fraction of the gaseous phase is calculated by its
continuity equation, while the continuity equation of
the liquid phase takes the form of the divergence of
the velocity field equal to zero. In the computational
cells occupied only by water, the volume fraction of
the gaseous phase is zero, while the volume fraction
of the liquid phase is equal to one. In the
computational cells occupied only by air, the
volume fraction of the gaseous phase is equal to
one, while the volume fraction of the liquid phase is
null. The free surface is located in the computational
cells occupied by both air and water. The numerical
solution of the continuity equation for the gaseous
phase is used to calculate the volume fraction of the
air and, thus, to calculate the position of the free
surface. In the above method, the free surface is not
sharply defined and its position does not coincide
with the boundary of a computational cell.
Consequently, by this method it is not
straightforward to impose boundary conditions for
the pressure and the kinematic boundary condition
at the free surface [3]. Furthermore, the adoption of
fixed Cartesian grids entails high computational
costs that make the application of this method to
real-scale numerical simulations of wave motion
very expensive. A different kind of numerical
models that are used to simulate free-surface flows
is based on the Smoothed Particle Hydrodynamics
(SPH) [4-6]. This approach can produce very
accurate representations of the free-surface
hydrodynamics but is usually limited to small-size
flow problems because its high computational cost.
In an alternative class of numerical models for the
three-dimensional simulation of free-surface flows,
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
128
Volume 17, 2022
the Navier-Stokes Equations are written in a system
of coordinates where the horizontal ones are
Cartesian, while the vertical one is a moving
coordinate (called 𝜎−𝑐𝑜𝑜𝑟𝑑𝑖𝑛𝑎𝑡𝑒) obtained by a
time-dependent transformation [7-9]. The above-
mentioned numerical models have good dispersive
properties that allow simulating the propagation of
non-breaking waves by computational grids with
less than ten nodes in the vertical direction. In
Rouvinskaia et al. 2018 [10], a numerical model
based on this method is used to simulate the internal
wave impact on the pillars of hydraulic engineering
constructions. In Zhang et al. 2021 [11] the same
approach is used in a three-dimensional model for
tsunami generation on irregular bathymetry.
The generalization of such an approach to time-
varying physical domains characterized by an
irregular shape has been proposed by [12], in which
a contravariant integral contravariant formulation of
the Navier-Stokes is written in time-dependent
curvilinear coordinates. In [13-15] the above-
mentioned motion equations are numerically
integrated by a second-order accurate Total
Variation Diminishing (TVD) shock-capturing
scheme commonly used in the literature [7,8], which
adopts an approximate Harten, Lax and van Leer
(HLL) Riemann solver. By using the above shock-
capturing scheme, a breaking wave is represented as
a discontinuity in the numerical solution. It is
known that, close to the discontinuities, the second-
order TVD scheme reverts to first-order locally and,
thus, introduces significant numerical dissipation of
kinetic energy in the numerical solution.
Furthermore, the approximate Riemann solver
assumes a simplified wave structure of the solution
of the Riemann problem that can produce too
dissipative solutions. As a consequence, the
numerical simulations of breaking waves carried out
by a second-order accurate TVD shock-capturing
scheme and an approximate Riemann solver are
affected by some errors: the increase in wave
amplitude during the shoaling process and the
maximum height of the waves are poorly predicted;
the position of the wave-breaking point and the
reduction of the wave height that take place after the
wave breaking are not correctly simulated. The
application of such low-order dissipative schemes
for the numerical simulation of breaking waves and
wave-induced currents requires very fine
computational grids, which make them suitable
mainly for simulations of laboratory-scale tests.
In this paper, we propose a new conservative non-
hydrostatic shock-capturing numerical scheme for
the numerical integration of the contravariant
integral form of the Navier-Stokes equations
proposed by [12]. The elements of novelty of the
proposed numerical scheme are three. The state of
the system is given by the field of the conserved
variables, 𝐻 and 𝐻𝑢, in which 𝐻 is the total water
depth and 𝑢 are the contravariant component of the
three-dimensional flow velocity. The first element
of novelty is given by the fact that, in the proposed
numerical scheme, the time-advancing of the
numerical solution is carried out by a predictor-
corrector procedure according to which an
approximate field of the conserved variable 𝐻𝑢∗
is corrected by the gradient of a potential scalar
function 𝜑 that take into account the non-hydrostatic
pressure component. Differently from [12], in the
present scheme, the scalar function 𝜑 is calculated
by numerically integrating an equation of Poisson-
type in which only the conserved variables 𝐻𝑢∗
are used, instead of the primitive variables 𝑢∗.
By this strategy, in the presence of discontinuities,
the proposed conservative scheme can converge to
the correct weak solution. The second element of
novelty concerns the reconstruction of the point
values on the cell faces of the computational grid. In
the proposed scheme, such point values are carried
out by a high-order reconstruction procedure, based
on an original Targeted Essentially Non-oscillatory
scheme which is designed for the simulation of
breaking waves. The third element of novelty
consist on the proposal of an exact Riemann solver
for updating the conserved variables on the cell
faces. The abovementioned modifications
significantly reduce the kinetic energy numerical
dissipation introduced by the scheme. Consequently,
the proposed numerical scheme differs from the
ones present in the literature, because allows us to
correctly represent the complex fully three-
dimensional flow patterns, that take place around
the coastal defence structures.
The proposed scheme is validated against
experimental tests of breaking waves and is applied
to the simulation of the vortices produced by the
interaction between the wave motion and an
emerged barrier. The hydrodynamic phenomena that
occur due to this interaction can induce significant
modifications in the coastal sediment transport and
local scouring around the barriers.
2 The proposed numerical scheme
In the proposed numerical scheme, the equations of
motion are obtained by the contravariant integral
formulation of the mass conservation and
momentum balance equations written in a time-
dependent curvilinear coordinate system. The main
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
129
Volume 17, 2022
element of the mathematical procedure consists on
the integration of the mass conservation and
momentum balance equations in a moving control
volume and the adoption of a time-dependent
coordinate transformation by which the moving
irregular control volume is converted in a fixed
prismatic one. By omitting the mathematical details
of the abovementioned procedure (that can be found
in [12]), the resulting motion equations are
𝑑
𝑑𝑡𝑔
()∙𝑔()𝐻𝑢𝑔𝑑𝜉𝑑𝜉𝑑𝜉
∆+
𝑔
()∙𝑔()𝐻𝑢(𝐻𝑢𝐻
⁄−𝑣)+𝑔
()
∆
∙𝑔()𝐺𝜂𝐻𝑔𝑑𝜉𝑑𝜉
− 𝑔
()∙𝑔()𝐻𝑢(𝐻𝑢𝐻
⁄−𝑣)+𝑔
()
∆
∙𝑔()𝐺𝜂𝐻𝑑𝜉𝑑𝜉=
− 𝑔
()∙𝑔()𝜕𝑝
𝜕𝜉𝐻𝑔𝑑𝜉𝑑𝜉𝑑𝜉
∆
+∑∫𝑔
()∙𝑔()
𝐻𝑔𝑑𝜉𝑑𝜉
∆
−
∫𝑔
()∙𝑔()
𝐻𝑔𝑑𝜉𝑑𝜉
∆
(1)
𝑑
𝑑𝑡𝐻𝑔𝑑𝜉𝑑𝜉𝑑𝜉
∆
+∑∫(𝐻𝑢−𝐻𝑣)𝑔𝑑𝜉𝑑𝜉
∆
−
∫(𝐻𝑢−𝐻𝑣)𝑔𝑑𝜉𝑑𝜉
∆
=0 (2)
where 𝐻 and 𝐻𝑢 are the conserved variables, given
by the water depth 𝐻 and its product by the
contravariant components of the flow velocity, 𝑢;
𝜂=𝐻−ℎ is the free-surface elevation; ℎ is the still
water depth; 𝐺 is the acceleration due to gravity; 𝜌
is the water density; 𝑝 is the dynamic pressure; 𝑅
are the contravariant components of the tensor of the
stresses devoid of the pressure term. 𝜉 is the 𝑙−𝑡ℎ
curvilinear coordinate, which is expressed as a time-
dependent function of the Cartesian coordinates, 𝑥,
𝜏=𝑡, 𝜉=𝜉(𝑥,𝑥,𝑥), 𝜉=
𝜉(𝑥,𝑥,𝑥),𝜉=𝑥+ℎ(𝑥,𝑥)𝐻(𝑥,𝑥,𝑡)
⁄.
In such coordinate transformation, 𝜉 is a moving
coordinate that changes over time according to the
variations of the water depth and assumes values
between 0 (at the bottom) and 1 (at the free surface).
In Eq. (1) and (2), 𝑔()=𝜕𝜉𝜕𝑥
⁄ and 𝑔()=
𝜕𝑥 𝜕𝜉
⁄ are the contravariant and covariant base
vectors, respectively; 𝐻𝑔=𝐻𝑘
∙𝑔()⋀𝑔() is
the specific expression assumed by the coordinate
transformation’s Jacobian, which is obtained
multiplying 𝐻 by the scalar product between the
vertical unit vector 𝑘
and the first two covariant
base vectors; ∆𝑉 is the volume of a computational
cell in the transformed space; ∆𝐴
and ∆𝐴
indicate, respectively, the area of the face of the
computational cell on which coordinate 𝜉 is
constant that is placed at larger and lower values of
𝜉; 𝑣 is the contravariant component of the
velocity vector that express the movements of the
generic coordinate 𝜉.
In this paper, Equations (1) and (2) are numerically
integrated by an original predictor-corrector finite-
volume scheme where the state of the system is
defined by the conserved variables 𝐻𝑢. In the
predictor step, a simplified form of Eq. (1), where
the dynamic pressure is omitted, is discretized by a
shock-capturing finite-volume scheme in which, at
the faces of each computational cell, a couple of
point values of each conserved variable is
reconstructed by an original high-order targeted
essentially non-oscillatory (TENO) procedure; these
couples of point values of the conserved variables
are assumed as initial values of the Riemann
problem. In the present numerical scheme, the
solution of each local Riemann problem is
calculated by extending to the three-dimensional
equations of motion the exact Riemann solver
proposed by [16]. In this way, the complete
structure of the wave solution is taken into account,
removing the simplifications adopted by the
approximate Riemann solvers proposed in [12-15].
The exact Riemann problem solution is used to
calculate a predictor field of the conserved
variables, 𝐻𝑢∗, that correspond to an
approximated three-dimensional velocity field with
non-zero divergence.
In the corrector step, by using exclusively the
predictor field of the conserved variables, we define
a Poisson equation which is written in the
curvilinear coordinate system,
()∙
()
=−∗
(3)
in which the unknown variable is a scalar function
φ. Eq. (3) is solved by a numerical iterative solver
in which the convergence is accelerated by a
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
130
Volume 17, 2022
multigrid method. The final field of the conserved
variables is obtained by adding the gradient of φ
(multiplied by 𝐻) to the predictor field
𝐻𝑢=𝐻𝑔()∙𝑔()
+𝐻𝑢∗ (4)
After the correction of the field of the conserved
variables, the position of the free surface is updated
by numerically solving the continuity equation
integrated over the vertical coordinate 𝜉, between
0 and 1.
The second element of novelty of the proposed
numerical concerns the reconstruction procedure in
the predictor step. Following the approach proposed
by [17], we define a convex combination of three
second-order interpolant polynomials, a function of
the smoothness of the numerical solution and a
dynamic threshold for choosing how many
polynomials participate in the given reconstruction.
If all the three polynomials are used, the value of
𝐻𝑢 on a given the cell face is calculated by a high
order approximation that significantly limits the
dissipation of kinetic energy due to the numerical
scheme. If only one or two candidate interpolant
polynomials participate to the reconstruction, 𝐻𝑢
is calculated by a low order approximation that
introduces numerical dissipation of kinetic energy
and avoids spurious unphysical oscillations in the
numerical solution. In the proposed procedure,
differently from the existing TENO schemes [17],
the value of this dynamic threshold depends on two
factors: the smoothness of the numerical solution
and the steepness of the front of the wave. During
the numerical simulations, in the instants and in the
grid nodes in which the front of the wave becomes
so steep to be considered a breaking wave front, the
dynamic threshold assumes its minimum values and,
consequently, is reduced the tendency of the
procedure to exclude one or two polynomials from
the reconstruction. Thus, at the wave breaking fronts
we obtain the maximum order of accuracy of the
reconstructions and minimize the dissipation of
kinetic energy introduced in the numerical solution
by the numerical scheme. By this strategy, we
entrust mainly to the turbulence model the task of
representing the wave breaking energy dissipation
and we leave mainly to the numerical scheme the
task of avoiding the unphysical oscillations at the
tail of the wave and at non-breaking fronts (where
the turbulence model is less effective). Below, we
describe the main mathematical aspects of the
proposed procedure.
The point value of a given conserved variable on a
face of the computational cell, 𝐻𝑢, is calculated
by a convex combination of interpolant polynomials
of second order, 𝑓 (with 𝑝=−1,0,1),
𝐻𝑢= ω𝑓 +ω𝑓+ω𝑓 (5)
where ω are the so-called nonlinear weights of the
combination, which are defined as
ω=
∑
(6)
In Eq. (6), δ are coefficients that can be 0 or 1 and
are used to exclude or not polynomials 𝑓 from the
combination; L are the linear weights and are
calculated to obtain fifth-order accuracy if no
polynomial is excluded from the combination. We
indicate by χ (with 𝑝=−1,0,1) a normalized
function that provides a measure of the smoothness
of the second-order polynomial 𝑓 and indicate by
𝐶 a dynamic threshold value for χ. At every
instant of the numerical simulation, χ are compared
with 𝐶, in order to assign 0 or 1 to coefficient δ
of Eq. (6),
δ=0,if χ<𝐶
δ=1,if χ>𝐶 (7)
The values of the dynamic threshold 𝐶 are given
by
𝐶=10 (8)
where exponent 𝑛 depends on two fixed parameters,
𝐴 and 𝐴, and two functions θ and θ
𝑛=𝐴+(θ+𝜃)(𝐴−𝐴) (9)
In Eq. (9), 𝐴 and 𝐴 are, respectively, the
assigned minimum and maximum for exponent 𝑛;
functions θ and 𝜃 can range between 0 and 1, and
their sum can be at most 1. Function θ depends on
the smoothness of the numerical solution [17]. 𝜃 is
an original function that we express as a function of
the wave front steepness
𝜃=
∗
−1 , if
>∗
𝜃=0, if
<∗
(10)
where 𝜕𝜂 𝜕𝑡
⁄ is the rate of change of the free-
surface elevation and 𝜕𝜂∗𝜕𝑡
⁄ is a threshold value
behind which the wave is considered as a breaking
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
131
Volume 17, 2022
wave. By Eq. (10), function 𝜃 is different from 0
only at the front of a breaking wave and increases as
the steepness of the wave front increases. The
proposed original reconstruction ensures high-order
of accuracy at the front of the breaking waves,
where kinetic energy dissipation is mainly
demanded to the turbulence model. At the wave tails
and non-breaking fronts (where 𝜃 is null), exponent
𝑛 depends only on θ, that is a function of the
smoothness of the numerical solution. In these
regions, where the turbulence model is less
effective, the ability of the proposed numerical
scheme to reduce the spurious (unphysical)
oscillations is maximum.
Unlike the second-order TVD reconstructions
adopted by [12-15], the proposed numerical scheme
ensures good non-oscillatory properties, without
introducing excessive numerical dissipation in the
numerical solution.
3 Results
Emerged barriers are commonly used in coastal
engineering. The interaction between the incoming
waves and an emerged barrier can produce velocity
and pressure fluctuations and very complex
instantaneous three-dimensional flow patterns
characterized by quasi-periodic vortices. These
vortices usually occur downstream of the edge of
the barrier and can produce the resuspension of solid
particles from the sea bottom and local scour
phenomena. In this section we apply the proposed
model to the simulation of the quasi-periodic vortex
formations that take place downstream of a vertical
barrier that interacts with breaking waves. A
rectangular 15 𝑚 × 3 𝑚 coastal area is discretized
by a computational grid shown in Fig.1, where
incoming cnoidal waves are produced at x = 0,
propagate from left to right on a 1:35 bottom slope,
and interacts with a 2 𝑚 long and 0.5 𝑚 wide
emerged barrier that is placed at 𝑥=7.1 𝑚 (red line
in Fig. 1). In both the x and y directions, the size of
the computational cells ranges between 0.025 𝑚
and 0.05 𝑚; in the z-direction the water depth is
discretized by 9 moving layers. During the
simulation, the first calculation node is placed at an
average distance from the bottom equal to 𝑧=30
(where 𝑧=𝑧𝑢∗𝜈
⁄ is a dimensionless distance
calculated by the friction velocity 𝑢∗ and the water
viscosity 𝜈). The remaining grid nodes in the
vertical direction are uniformly distributed along the
water column. The boundary on which 𝑦 =1 𝑚 is
treated as a reflecting boundary, on which the
normal velocity component is assumed equal to
zero, while a zero normal derivative is assumed for
the remaining quantities. The opposite boundary
(𝑦=4 𝑚) is treated as an open boundary, on which
the normal derivative of every quantity is assumed
equal to zero. As input boundary conditions we
impose 0.125 𝑚 high cnoidal waves, whose wave
period is equal to 2 𝑠. The eddy viscosity is
expressed by a Smagorinsky turbulence model in
which the Samgorinsky coefficient 𝐶 ranges
between 0.05 and 0.2. On the left boundary, at x = 0
m, a cnoidal wave with period 𝑇 = 2 𝑠 and wave
height 𝐻 = 0.125 𝑚 is imposed. Such a wave is
equal to the one experimentally reproduced by [18]
on a channel with the same initial water depth and
1:35 sloping beach (without the emerged barrier). In
Fig. 2 we show the comparison between the
experimental result of [18] and the numerical results
obtained by the proposed model. The wave height
obtained with the numerical model is in very good
agreement with the experimental results, as shown
in Fig. 2. The breaking point is located around 𝑥 =
6.5 𝑚 as the one obtained by experimental
measurements. In order to quantify the agreement
between the numerical and the experimental results,
we compare the free-surface elevation by using the
mean absolute percentage error (MAPE) [19]
𝑀𝐴𝑃𝐸=∑
×100 (11)
where 𝑅𝑒𝑥𝑝 and 𝑅𝑛𝑢𝑚 are respectively the
experimental and numerical measurements.
MAPE
4.85%
Table 1 Mean absolute percentage error for the free-
surface elevation.
The very low relative error between the proposed
model and the experimental results (less than 5 per
cent) demonstrates that the proposed model is able
to correctly simulate the wave shoaling, the
maximum wave height, the wave breaking point and
the wave height reduction in the surf zone.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
132
Volume 17, 2022
Fig.1. Computational grid (one grid line out of every three is drawn). The red line indicates the position of the
emerged barrier.
Fig. 2. Breaking wave. Maximum free-surface elevation, mean water level and minimum free-surface
elevation.
Circles:
experimental
measurements
by
[
18
]; lines:
proposed model numerical results
.
Fig. 3 shows three instants of the simulated free-
surface elevation given by the interaction between
the cnoidal waves and the emerged barrier. In such
figures, the emerged barrier is not drawn and its
basis is represented by the blue color. From Figs.
3(a) - 3(c) it is possible to see the wave height
reduction produced by the breaking of the waves
that approach the barrier and the reflection caused
by the wave collision with the barrier. The spatial
variations of the wave height caused by the
interaction with the barrier produce currents which
can be highlighted by averaging over time the
instantaneous flow velocity components.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
133
Volume 17, 2022
(a)
(b)
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
134
Volume 17, 2022
(c)
Fig. 3
(a)
-
(c).
A s
equence of
three
instantaneous views of the
simulated wave field
Fig. 4 shows the numerical results of the currents
induced by the wave-structure interaction at three
different depths: a) close to the seabed; b) at middle
depth; c) close to the water surface. Fig. 4 shows
that the interaction between the incoming waves and
the barrier produces, at the onshore side of the
barrier, a large eddy which turns in a clockwise
sense and increases its intensity as the distance from
the bottom increases.
(a)
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
135
Volume 17, 2022
(b)
(c)
Fig. 4. Plan view of the simulated circulation patterns: (a) close to the seabed; (b) at middle depth; (c) close to
the water surface.
In Fig. 5 the quasi periodic vortices that take place
near the barrier and propagate in the direction of the
waves are visualized by the so-called Q-method for
the vortices identification [20], according to which
positive values of the second invariant of the
velocity gradient tensor (denoted by Q) indicate the
presence of a vortex. Fig. 5 shows that vortices
characterized by a vertical axis take place close the
barrier edges, starting from the seabed; as the
distance from the barrier edges increases, the
orientation of these vortices changes and tends to
the direction of the propagation of the waves. These
vortex structures can put into suspension the solid
particles from the seabed and produce local erosion
phenomena near the barrier.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
136
Volume 17, 2022
(a)
(b)
(c)
Fig. 5.
Vortex structures identified by
the
Q
-
method (
𝑄
=
3
.
5
)
.
4 Discussion
The proposed numerical model differs from the ones
present in the literature; this model significantly
reduces the kinetic energy numerical dissipation
introduced by the scheme and allows us to correctly
represent the complex fully three-dimensional flow
patterns, that take place around the coastal defence
structures. By the simulation of the vortices
produced by the interaction between the wave
motion and an emerged barrier, it can be notice that
the proposed numerical model correctly represents
the hydrodynamic phenomena that can induce
significant modifications in the coastal sediment
transport and local scouring around the barriers.
5 Conclusion
A numerical model for the three-dimensional
numerical simulation of wave-structure interactions
has been presented. The equations of motion at the
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
137
Volume 17, 2022
basis of the numerical model are expressed in a
moving coordinate system and are integrated by an
original conservative shock-capturing numerical
scheme. The main element of originality of the
proposed numerical scheme is given by an original
TENO scheme specifically designed to simulate the
breaking of the waves; the local Riemann problem
produced by the TENO reconstruction procedure is
solved by an exact Riemann solver. By the proposed
numerical scheme, the three-dimensional flow
structures produced by the interaction between
trains of breaking waves and a coastal defence
structure parallel to shoreline has been simulated.
The results obtained by the proposed numerical
scheme show that by this approach it is possible to
represent both the large scale hydrodynamic
phenomena, like the wave-induced circulation
patterns downstream the barrier, and small scale
flow structures, like the quasi periodic vortices that
take place near the barrier’s edges, close to the
bottom. Both the hydrodynamic phenomena are
fully three-dimensional and can induce significant
modifications in the coastal sediment transport and
local scouring around the barriers.
References:
[1] Bradford S.F., Numerical simulation of surf
zone dynamics, Journal of Waterway, Port,
Coastal, and Ocean Engineering, 126, 1, 2000,
pp. 1-13.
[2] Higuera P, Lara J.L., Losada I.J., Simulating
coastal engineering processes with
OpenFOAM®, Coastal Engineering, 71, 2013,
pp. 119-134.
[3] Lin P., Li C.W., A σ‐coordinate
three‐dimensional numerical model for surface
wave propagation, International Journal for
Numerical Methods in Fluids, 38, 11, 2002, pp.
1045-1068.
[4] Price D.J., Smoothed particle hydrodynamics
and magnetohydrodynamics, Journal of
Computational Physics, 231(3), 2012, pp. 759–
794.
[5] Yang X., Peng S., Liu M., A new kernel
function for sph with applications to free
surface flows, Applied Mathematical
Modelling, 38(15-16), 2014, pp. 3822–3833.
[6] Perea J., Cordero J., A Stable Hybrid Potential–
SPH Technique to Enforce the Fluid
Incompressibility, WSEAS Transactions on
Fluid Mechanics, 13, 2018, pp. 50-59.
[7] Bradford S.F., Nonhydrostatic model for surf
zone simulation, Journal of Waterway, Port,
Coastal, and Ocean Engineering, 137, 4, 2011,
pp. 163-174.
[8] Ma G., Shi F., Kirby J., Shock-capturing non-
hydrostatic model for fully dispersive surface
wave processes, Ocean Modelling, 43, 2012,
pp. 22-35.
[9] Iele B., Palleschi F., Gallerano F., Boundary
conditions for the simulation of wave breaking,
WSEAS Transaction on Fluid Mechanics, 15,
2020, pp.41-53.
[10] Rouvinskaia E., Kurkina O., Kurkin A.,
Zaytsev A., Modeling of internal wave impact
on hypothetical pillars of hydraulic engineering
constructions in the conditions of the Sakhalin
island shelf, WSEAS Transactions on Fluid
Mechanics, 13, 2018, pp. 101-107.
[11] Zhang C., Kirby J.T., Shi F., Ma G., Grilli,
S.T., A two-layer non-hydrostatic landslide
model for tsunami generation on irregular
bathymetry. 2. Numerical discretization and
model validation, Ocean Modelling, 160, 2021,
pp. 101769.
[12] Cannata G., Petrelli C., Barsi L., Gallerano F.,
Numerical integration of the contravariant
integral form of the Navier–Stokes equations in
time-dependent curvilinear coordinate systems
for three-dimensional free surface flows,
Continuum Mechanics and Thermodynamics,
31, 2, 2019, pp. 491-519.
[13] Tamburrino M., Cannata G., Simulation of
wave run-up by means of the exact solution of
the wet/dry Riemann problem, WSEAS
Transactions on Fluid Mechanics, 14, 2019,
pp. 114-123.
[14] Gallerano F., Cannata G., Tamburrino M.,
Ferrari S., Badas M.G., Querzoli G., Water
waves overtopping over barriers, WSEAS
Transactions on Fluid Mechanics, 14, 2019,
pp. 84-91.
[15] Gallerano F., Cannata G., Barsi L., Palleschi,
F., Iele B. Simulation of wave motion and
wave breaking induced energy dissipation,
WSEAS Transactions on Fluid Mechanics, 14,
2019, pp. 62-69.
[16] Toro E.F., Shock-capturing methods for free-
surface shallow flows, Wiley-Blackwell, 2001.
[17] Peng J., Liu S., Li S., Zhang K., Shen Y., An
efficient targeted ENO scheme with local
adaptive dissipation for compressible flow
simulation, Journal of Computational Physics,
425, pp. 109902.
[18] Ting F.C.K., Kirby J.T., Observation of
undertow and turbulence in a wave period,
Coastal Engineering, 24, 1994, pp. 51-80.
[19] Hyndman R. J., Koehler A. B. Another look at
measures of forecast accuracy, International
Journal of forecasting, 22, 2006, pp. 679-688.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
138
Volume 17, 2022
[20] Hunt J.C., Wray A.A., Moin P., Eddies,
streams, and convergence zones in turbulent
flows, Studying turbulence using numerical
simulation databases, 2. Proceedings of the
1988 summer program, 1988, pp. 193-208.
Contribution of individual authors to
the creation of a scientific article
(ghostwriting policy)
Francesco Gallerano and Giovanni Cannata
conceptualized the model and ideated the scheme.
Federica Palleschi and Benedetta Iele carried out the
validation and application of the proposed model.
Sources of funding for research
presented in a scientific article or
scientific article itself
This research received no external funding.
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
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2022.17.13
Francesco Gallerano, Federica Palleschi,
Benedetta Iele, Giovanni Cannata
E-ISSN: 2224-347X
139
Volume 17, 2022
