Sådhanå (2022)47:45  Indian Academy of Sciences

[Link]
Sadhana(0123456789().,-volV)FT3](012345
6789().,-volV)

Study of water wave in the intermediate depth of water using second-

order Stokes wave equation: a numerical simulation approach
DEEPAK KUMAR SINGH and PRADIP DEB ROY*

Department of Mechanical Engineering, NIT, Silchar, Assam 788 010, India

e-mail: pdroy02@[Link]

MS received 14 June 2020; revised 20 August 2021; accepted 8 December 2021

Abstract. The paper presented for simulating the complete physics of the second-order Stokes wave equation
using CFD code ANSYS-FLUENT software and volume of fluid (VOF) method in a 2-D numerical wave tank
(NWT). The main objective of this work is to study the effect of wave steepness on ocean waves in the
intermediate depth of water at low Ursell numbers (Ur \ 18). The inlet velocity method generates the second-
order water wave in the numerical model for achieving this objective. Generally, three methods are used to
generate the wave in a numerical wave tank: (1) inlet velocity method; (2) flap-type method; and (3) piston type
method. I have used the inlet velocity method to generate the second-order water wave in the simulation. In this
study, the simulation results of the model are compared with the analytical results. It shows that the accuracy of
the numerical results is in good agreement with the analytical results. The MATLAB code solves the second-
order Stokes wave equation. Six tests under different wave heights have been conducted at two regions in an
NWT to analyze the wave theory. It also shows that high energy flow is concentrated near the free surface. It
also shows that the nonlinearity effect increases with the increase of wave steepness and found horizontal and
vertical velocity increases with wave steepness. This water wave study improves the coastal engineering
application for further investigation and extracts wave energy from the ocean.

Keywords. CFD; second-order stokes wave theory; surface elevation; horizontal and vertical velocities;
velocity vector.

1. Introduction performed using numerical models. However, some diffi-

culties exist in a numerical model for generating the wave.
Ocean wave is an abundant source of renewable energy. These difficulties are moving boundaries on the free sur-
According to the World Energy Council [1], the global face, wavemaker boundary conditions, and the selection of
ocean wave power has been estimated at 2-TW, and the appropriate non-reflecting far-field boundary conditions.
energy density (energy per unit length) of the waves is very Many authors studied various types of wave theories in
high. Despite the energetic advantages, wave energy con- different water depths under different wave steepness and
version into mechanical or electrical power is developed by relative depth conditions using ANSYS-FLUENT software
different types of wave energy converters (WEC). The run and the volume of fluid (VOF) method. Such wave theories
of WEC in the ocean is not accessible due to the complex are linear wave theory, stokes wave theory, etc. Machado
behavior of the ocean. To date, about one thousand proto- et al [3] extended regular wave theory to the second-order
types WEC are known, with only two hundred having Stokes theory by CFD ANSYS CFX software. The waves
reached the test phase. None of them is still economically have been tested in two ways: the inlet velocity method and
viable for commercialization and robust design due to the the piston-type wavemaker theory. Most of the other
complex behavior of the ocean (Mustapa et al [2]). researchers also investigated the second-order wave theory
Studying water wave properties and their impacts on (Flick and Guza [4]; Moubayed and Williams [5] etc.).
coastal structures and related coastal phenomena like WEC Maâtoug and Ayadi [6] numerically simulate the second-
is often performed in physical wave tanks and flumes, order Stokes wave theory by the finite difference method.
where the inlet velocity method produces the desired water The crank–Nicholson method and an implicit approach
waves. A numerical wave tank is an alternative to physical were used for the solution. The second-order wave theory
modeling because studying different wave conditions and leads to an unusual bump in the wave trough for large
implementing the modifications are more conveniently waves (Dean and Dalrymple [7]). Therefore, higher-order
solutions were proposed, such as third-order by Borgman
*For correspondence and Chappelear [8] and fifth-order by Fenton [9]. But
45 Page 2 of 15 Sådhanå (2022)47:45

higher-order wave theory is highly complicated to explain development of three-dimensional viscous wave flume. The
analytically. Therefore, higher-order wave theories need to other researchers also used this method (Huang and Dong
study numerically. [22]; Karim et al [29]) in various coastal applications.
The nonlinear wave theory of higher-order was first The second category of the VOF models is the internal
introduced numerically by Chappelear [10] and then further wave generation method; the mass source function is
developed by (Dean [11]; Chaplin [12]). Other numerical introduced in a particular region inside the computational
methods based on internal wave generation models devel- domain (Lin and Liu [30]). In this method, the fluid is
oped for generating waves. These models have some alternatively injected or sucked into the computational
advantages for avoiding interference with the boundary domain such that it produces the same physical effect as the
conditions. The source line method with the Boussinesq desired wave. Various types of waves were generated using
equations was used by Larsen and Dancy [13] to make short this method through the proper definition of the source
waves. Several researchers developed this approach to function, such as linear and nonlinear; irregular, Stokes
generate linear and nonlinear waves for regular waves, second and higher orders, solitary and cnoidal waves (Lin
irregular waves, and multi-directional waves (Brorsen and and Liu [30]). This model has been successfully used by
Larsen [14]; Li et al [15]). Kawasaki [31] to study breaking wave over the submerged
Inviscid fluid and irrotational flow are applicable in non- breakwater by Hieu and Tanimoto [32] to simulate wave-
viscous flow, and it is not acceptable in the studies structure interactions, and by Hur and Mizutani [33] to
mentioned above. Therefore, viscous fluid assumptions are determine the transverse wave forces that act on 3D
inevitable in numerical models. Marker and Cell (MAC) asymmetric structures on a submerged permeable break-
method is a modified version used by Harlow and Welch water. Zhi and Zhan [34] investigated numerical modeling
[16] and Park et al [17] on the free-surface flows to study of wave evolution, breaking, and run-up in shallow water
the solitary wave in shallow water. Anbarsooz et al [18] by the VOF method. After modification, Hafsia et al [35]
and Wu and Hsiao [19] also developed nonlinear viscous introduced this method in the one-dimensional region by
waves numerically to study the solitary waves in a two- reducing the source domain.
dimensional numerical wave tank in the intermediate depth The two categories mentioned above are producing
of water under low and high steepness conditions. After artificial free surface wave profiles, and they are different
that, a modified version of the accurate technique, called from the natural physical wave generation phenomenon. In
the Stanford University Modified MAC (SUMMAC), was the first category, the velocity components at the inflow
used to determine the velocity components at the free boundary are set according to an analytical solution of the
surface cells. Tang et al [20] used a mathematical expres- desired wave. Therefore, the velocity profile at this
sion that becomes more accurate than previous for the boundary is a function of the vertical direction. For
dynamic boundary conditions on a free surface, including example, a piston-type wavemaker moves horizontally with
the viscosity and surface tension effects. The SUMMAC a velocity that has no vertical variation. The added mass
method and exact mathematical expression simultaneously source term is computed according to a prescribed free
were also used by Huang et al [21] to investigate the surface profile in the other category. It shows that the flow
nonlinear viscous wavefields generated by a piston-type pattern close to the source region is different from the flow
wavemaker. Huang and Dong [22] and Dong and Huang pattern around an actual physical wave generator paddle.
[23] used this dual numerical scheme to generate small and It shows from the literature that the resultant wavelength
finite-amplitude waves and solitary waves in a two-di- in a natural wave generation mechanism is a function of
mensional wave flume. Wang et al [24] also applied the wavemaker period, stroke, and still water depth. Predicting
same method for simulating a three-dimensional numerical the wavelength is performed using wavemaker theories
viscous wave tank equipped with a piston-type wavemaker. (Dean and Dalrymple [7], Madsen [36]). The theories are
The Volume of Fluid (VOF) method was first introduced linear wavemaker theory and the second-order wavemaker
by Hirt and Nichols [25] in 1981 for treating the flows with theory. However, the wavelength is obtained, either by the
a free surface and then developed many numerical wave wave speed or by the mean velocity at a point in the fluid. If
tanks based on this method. The literature study shows the none of these are known, then Fenton [9] stated that the
VOF method is mainly classified into two categories for theory is irrational and likely to be in error at first order.
wave generation. The first category is the VOF-inflow That is why, for the steeper waves, the discrepancy between
method; here, the boundary conditions of the desired wave the experimental results and those of the wavemaker theory
are set based on the free surface elevation and the velocity becomes more pronounced (Ursell et al [37]).
components obtained from the analytical solution of the Therefore, developing a numerical wave tank that sim-
desired wave. This technique was first used by Lin and Liu ulates the real physical process of wave generation is of
[26] for the wave generation in a 2D wave flume and then great interest. This goal is achieved by suitable modeling
used by Troch and De Rouck [27] to develop an active the motion of the fluid in the tank. In this case, the wave
wave generating-absorbing boundary condition. Li and period, the wavelength, and the wave height calculated
Fleming [28] used the VOF-inflow technique for the through the complete solution of the Navier–Stokes
Sådhanå (2022)47:45 Page 3 of 15 45

equations. Wood et al [38] used the Fluent software for the based on Neumann type and Dirichlet type boundary con-
piston-type wavemaker, and Finnegan and Goggins [39] ditions. Figure 1 represents the geometry of the numerical
used the Ansys CFX commercial software for the flap-type wave tank with its dimensions and boundaries.
wavemaker. Dao et al [40] used Open FOAM software to The rectangular grid has been constructed for computa-
investigate the ocean’s wave characteristics by modeling a tional domain discretization. Grid with approximately
physical wave tank with a flap paddle and porous beach. 200901 nodes has been tested for investigating the
The other researchers were also using ANSYS, Fluent numerical results (figure 2). To analyze the free surface
software (Zhao et al [41]; Prasad et al [42]). more accurately, mesh refinement has needed near the
The past literature shows the different numerical meth- vicinity of the mean water level, as shown in figure 3. The
ods used to solve the higher-order problem in an NWT, and mesh refinement in a numerical simulation performed on
the simulation results validated with the mathematical the Mesh 3 system supports the grid convergence test. The
model. It also shows from the past literature that most finer meshes are observed near the free surface at the
researchers have not mentioned the importance of Ursell working zone, while coarser meshes have been found
numbers in higher-order problems and the effect of wave towards the right wall of the domain. The mesh size
steepness on the ocean wave. Velocity contour and the increased gradually towards the end of the NWT to provide
velocity vector of water particles are the essential factors in a numerical damping effect. Nonlinear waves of wave
the coastal engineering sciences not mentioned in their heights H = 0.1 m, 0.15 m, and 0.2 m generated in an
study. However, a numerical method is applied in the NWT for T = 3.23 s.
present study, which simulates the wave generation phe- The specification of the system at which simulation has
nomenon in an NWT using CFD code ANSYS-FLUENT run is an Intel(R) Xeon(R) [email protected] GHz
software. The VOF method is implemented accurately in 3.31 GHz processor and a 64.0GB RAM. A simulation was
the model for the wave generation performed by the inlet performed based on RANSE with the viscous model. A first-
velocity method. This study mainly discussed the effect of order implicit scheme and a second-order upwind
wave steepness on ocean waves at low Ursell numbers scheme are used here to discretize. The time step size has
(Ur \ 18) and plot velocity contour and the velocity vector been set at 0.01 s. The number of time steps has been set at
of water particles, which are the essential factors in the 4000. The maximum iteration has been considered 20 to
coastal engineering sciences. The presented model is cap- complete the solution. The Courant number has been set at
able of producing a second-order Stokes wave equation. 0.25 for stability. Based on these parameters and the com-
The accuracy of the proposed numerical model was verified putational capabilities, the simulation took 12 hrs to solve.
by comparing the results of the analytical data. MATLAB
code is used to obtain the analytical data.
2.1 Geometry of damping
2. Numerical model Dissipate the wave energy and avoid secondary wave
reflection; a numerical damping zone is conducted in lx2, as
The Numerical Simulation in the computational domain is a shown in figure 1. Different authors propose different
rectangle of size (lx and lz). Where lx is the length of the tank, methods in their literature for the elimination of wave
and lz is the height of the tank. One end of the tank is a wave reflections, namely:
generating zone (lx1), and the other is the damping zone (lx2).
• controlling the total time of simulation to ensure that
Calm water depth (d) is 1 m. Cartesian coordinates are
the first wave generated does not reach the end of the
employed in the x-z plane, vertical to the fluid surface. The
tank (Silva et al [43]; Finnegan and Goggins [39]).
z-axis is directed vertically upward from the still water level
• mesh size increases at the end of the tank to dissipate
(SWL), which is positive. The x-axis is along the direction
wave energy (Kim et al [44]; Silva et al [43]).
of the propagation of waves. The inlet velocity method
• implementing a beach (Elangovan [45]; Saincher and
generates the wave at a constant velocity at the inlet (x = 0)
Banerjee [46]).
of the NWT and flows in the x-direction. Consider the free
• size of the tank length increase for reducing the wave
surface is oscillating in the z-direction. The effect of flow in
reflection from the affecting area. (Morgan et al [47];
the y-direction, is not considered. The top portion of the
Prasad et al [42]).
NWT is open to the atmosphere and is assumed to be zero
• creating a damping zone (Koo and Kim [48];
pressure. The bottom and right walls of the numerical wave
Westphalen et al [49]).
tank (NWT) are solid walls where the no-slip boundary
conditions are applied. The no-slip boundary condition In this study, some of the above methods are imple-
means that the velocity component normal to the solid mented to dissipate the wave energy: (1) mesh size grad-
boundary will have zero velocity relative to the wall. The ually increases at the damping zone, (2) length of the
tangential component is unaffected by the boundary con- numerical wave tank increase, and (3) the simulation time
dition. Therefore the present numerical model is solved is controlled to prevent the waves from reaching the end of
45 Page 4 of 15 Sådhanå (2022)47:45

Figure 1. Definition sketch of a numerical wave tank.

the tank. Furthermore, the height of the numerical wave

tank above the SWL is increased to 2.5 m to prevent air
circulation at the top boundary from interfering with the
wave elevation.

2.2 Governing equation

The governing equations are the mass continuity equation
and the Navier-Stokes (N-S) equations for incompressible
fluid flow, given by:
ou ow
þ ¼0 ð1Þ
ox oz

Figure 2. Time histories of surface elevation (g) at H/L = 0.0210

for the convergence test on various mesh sizes.

Figure 3. Mesh refinement near the free surface.

Sådhanå (2022)47:45 Page 5 of 15 45
 2  condition. The left wall boundary is a wavemaker, while
ou ou ou 1 op o u o2 u
þu þw ¼ þm þ ð2Þ the right part of the domain is a damping zone.
ot ox oz q ox ox2 oz2
The boundary conditions are:
 2 
ow ow ow 1 op o w o2 w (a) Dynamic free surface boundary condition (DFSBC)
þu þw ¼ þm þ ð3Þ
ot ox oz q oz o x2 o z2 based on Bernoulli’s equation
"    #
Where u and w are the velocity components of the x and
ou 1 ou 2 ou 2 patm
z-axis, respectively [m s-1]; t is the time [s]; p is the þ þ þ þ gg ¼ 0 z ¼ g
ot 2 ox oz q
pressure [Pa]; q is the density [kg m-3], which varies over
the air and water interface, and m is the kinematic viscosity ð8Þ
[m2 s-1].In the present study, the interface of the two fluids
has been modeled by the VOF method (Hirt and Nichols where / is the velocity potential. g = g(x,z,t) is the
[25]). This method determines the volume fraction of each water surface elevation, and g is the acceleration due to
fluid (air and water) in each control volume. The volume gravity, and patm is the atmospheric pressure on the free
fraction function (aq) is the ratio of the volume occupied by surface and usually taken as gage pressure, patm = 0
the qth phase in a cell to the cell’s total volume. Accord- (Dean and Dalrymple [7]). Hence Eq. (8) reduced to:
ingly, there can be three possible conditions. Such that "    #
ou 1 ou 2 ou 2
aq = 1 means that the cell filled with the fluid of the qth þ gg þ þ ¼ 0 z ¼ g ð9Þ
phase, aq = 0 means that the cell does not contain any fluid ot 2 ox oz
of the qth phase, and when 0 \ aq \ 1, the cell is called an
(b) Kinematic free surface boundary condition (KFSBC):
interface cell. For water wave problems, q = 1 or 2 since
there are only two phases: air and water. The following ou og ou og
equations determined aq: ¼ þ z¼0 ð10Þ
oz ot ox ox
oaq   (c) Bottom boundary condition (BBC):
þ r  aq V ¼ 0 ð4Þ
ot
ou
X ¼ 0 z ¼ d ð11Þ
2
oz
aq ¼ 1 ð5Þ
q¼1
The velocity potential of second-order Stokes wave
Where V is the velocity vector. The density of the mixture equation (u) and the free surface profile (g) is obtained by
is computed based on the volume fraction, as follows: substituting the boundary conditions Eqs. (8)–(11), into
  Laplace equation Eq. (7) and then expanding through the
q ¼ aq qw þ 1  aq qa ð6Þ Taylor series expansion (Dean and Dalrymple [7]).

Where qw is the water density, and qa is the air density. In ag cosh kðz þ dÞ sinðk x  x tÞ
u ðx; z; tÞ ¼
the present work, the density of air and water is defined as x coshðk dÞ
qa = 1.225 kg/m3, qw = 998.2 kg/m3, respectively. 3 x a2
þ coshð2kðz þ dÞÞ sinð2ðk x  x tÞÞ
8 sinh4 ðk dÞ
ð12Þ
3. Theoretical formulation
k a2 coshðk dÞ
The governing equation is the Laplace equation for irrota- g ¼ a cosðk x  x tÞ þ ð2
4 sinh3 ðk dÞ
tional wave motion: þ cosh 2k dÞ cos 2ðk x  x tÞ ð13Þ
r u¼0
2
0  x  lx ; 0  z   d ð7Þ Where ’a’ is the wave amplitude and xð¼ 2p=TÞ is the
angular frequency of the wave, and the wavenumber
To solve the governing equations, it is necessary to
kð¼ 2p=LÞ satisfies the dispersion relation given in
specify appropriate boundary conditions at all domain
Eq. (14). The dispersion relation is obtained by expanding
boundaries. The boundary conditions which should be
the KFSBC through the Taylor series expansion (Dean and
satisfied are as follows: (1) kinematic boundary conditions
Dalrymple [7]) by considering the linear terms at z = 0.
and dynamic boundary conditions at the free surface. The
free surface opened to the atmosphere is the top boundary
x2 ¼ kg tanhðkdÞ ð14Þ
of the NWT (2) No-slip boundary condition at the tank base
of the computational domain. The no-slip means that the Eq. (15) and Eq. (16) give the local velocity in the lon-
velocity component normal to the wall equals 0, and the gitudinal and vertical directions, respectively, and these
tangential component is unaffected by the boundary velocities are obtained by partially differentiating the
45 Page 6 of 15 Sådhanå (2022)47:45

potential function / (Eq. 12) with respect to x and z, theories by an Ursell number (Ur = 40). The theory is
respectively. applicable at intermediate depths (0.05Bd/LB0.5) and
deep water depths (d/L C 0.5). The wave must be steep
ag 1
u¼ k coshðkðz þ dÞÞ cosðk x  x tÞ
x coshðk dÞ
3 x a2 Table 2. Mesh size parameters.
þ 2k coshð2kðz þ dÞÞ cosð2ðk x  x tÞÞ
8 sinh4 ðk dÞ
Mesh Dx (m) Dz (m) Nodes
ð15Þ
1 0.20000 0.00666 90601
ag 1 2 0.15789 0.00526 145161
w¼ k sinh ðkðz þ dÞÞ sinðk x  x tÞ
x coshðk dÞ 3 0.12000 0.00500 200901
3 x a2 4 0.10909 0.00444 248501
þ 2k sinh ð2kðz þ dÞÞ sin ð2ðk x  x tÞÞ
8 sinh4 ðk dÞ
ð16Þ

4. Region of possible waves and the validity

of the theory

Fenton proposed a range over which periodic solutions

for different waves theory can occur, shown in figure 4,
which shows the limit between cnoidal and Stokes wave

Figure 4. Region of possible Stokes wave theory (Fenton [50]). Figure 5. Limits of validity for various wave theories (B. Le
Méhauté [53]).

Table 1. Wave parameters were used in this study for the six tests.

Tank Time Water Wave Wave Ursell No.

Test length lx (m) period T (s) depth d (m) length L (m) height H (m) (d/L) (H/L) d/gT2 H/gT2 Ur < 40
On the surface z = 0
1 60 3.23 1 9.4826 0.10 0.105 0.0105 0.01 0.0010 8.9919
2 60 3.23 1 9.4826 0.15 0.105 0.0158 0.01 0.0015 13.4879
3 60 3.23 1 9.4826 0.20 0.105 0.0210 0.01 0.0020 17.9839
On the bottom z = -1
4 60 3.23 1 9.4826 0.10 0.105 0.0105 0.01 0.0010 8.9919
5 60 3.23 1 9.4826 0.15 0.105 0.0158 0.01 0.0015 13.4879
6 60 3.23 1 9.4826 0.20 0.105 0.0210 0.01 0.0020 17.9839
Sådhanå (2022)47:45 Page 7 of 15 45

Figure 6. Comparison of horizontal velocity (u) between analytical and numerical results under different wave steepness (H/L) at z = 0
and z = -1 m.

crested, and flatter trough, with no wave breaking HL2

occurring. The range of applicability of different wave Ur ¼ \40 ð17Þ
d3
theories depends on the water depth (d), the wavelength
(L), and the wave height (H). For these criteria, second- where Ur, known as the Ursell number, was proposed by
order Stokes waves theory must satisfy the Ursell number Ursell [52]). It has shown in figure 4 that the three nonlinear
given in Eq. (17) (Fenton [50]; Hedges and Ursell [51] waves satisfied the conditions of second-order Stokes
and table 1). waves and their corresponding coordinates are (0.10,
45 Page 8 of 15 Sådhanå (2022)47:45

Figure 7. Comparison of vertical velocity (w) between analytical and numerical results under different wave steepness (H/L) at z = 0.

9.4826), (0.15, 9.4826), and (0.20, 9.4826) respectively for gradually. The node numbers of four different mesh
the water depth d = 1 m. systems are 90601, 145161, 200901, and 248501. For the
convergence test, surface elevation (g) at the wave steep-
ness (H/L = 0.0210) from the four different meshes were
5. Convergence test compared, shown in figure 2. It is observed from the con-
vergent results that there is no difference between Mesh 3
Four different meshes shown in table 2 are used for the and Mesh 4 as the size of the meshes decreases from Mesh
convergence test. ‘Mesh 1’ is the coarse mesh system, and 3 to Mesh 4. Therefore, all numerical simulations were
the size of the mesh decreases from ‘Mesh 1’ to ‘Mesh 4’ performed on Mesh 3 system.
Sådhanå (2022)47:45 Page 9 of 15 45

Figure 8. Velocity contour in the (a) x-direction and (b) the z-direction.

Figure 9. Velocity vector of water particles at x = 27 m.

6. Results and discussions vector of the water particle. Different parameters were used to
solve the problem: depth of water d = 1 m, wave heights
This section mainly discusses the flow behavior of the second- H = 0.1 m, 0.15 m, and 0.2 m, wavelength L = 9.4826 m, and
order Stokes wave theory in the intermediate depth of water the time period T = 3.23 s are shown in table 1. These
and the effect of wave steepness on ocean waves at low Ursell parameters are more appropriate for the test zone according to
numbers numerically. The study parameters are horizontal the limit of Ursell number Ur \ 18. According to the wave-
velocity, vertical velocity, surface elevation, and the velocity length of the ocean wave, tank length is considered 60 m.
45 Page 10 of 15 Sådhanå (2022)47:45

Figure 10. Comparison of surface elevation (g) between analytical and numerical at different wave steepness (H/L).

The diagram of B. Le Méhauté [53], represented in figure 5, steepness, which is a function of the relative depth. Corre-
shows approximately the ranges of validity and the test of the sponding coordinates of the red closed area where the test has
zone of the domain of the second-order Stokes wave theory by been conducted are d/gT2 = 0.01 and H/gT2 = 0.0010, 0.0015,
a red closed area in the intermediate water. The graph is and 0.0020. These coordinates shown in table 1 are considered
established for two-dimensional periodic waves, but it indi- within the domain for analyzing the present wave theory. It
cates any water waves. The graph is limited by breaking observed, as the values of H/gT2 increase, the nonlinearity
waves, which implies a maximum value for the wave effect of the wave increase.
Sådhanå (2022)47:45 Page 11 of 15 45

Figure 11. Evolution of the numerical free surface profile in an NWT at different time-step t for H = 0.1 m.

Figure 6 shows that the horizontal velocity (u) plotted Figure 6(c) shows that the velocity curve is more flatter in
against the time (t) at the free surface z = 0 and at the bottom trough shown by the black circle for the wave steepness
z = -1 m for T = 3.23 s and d/L = 0.105. The figures show H/L = 0.0210, and d/gT2 = 0.01, H/gT2 = 0.002 (refer to
the continuous line and the dotted line for analytical results figure 5) compared to figure 6(a), indicating that the nonlin-
and numerical results, respectively. The figures also show earity effect increases with the increase of wave steepness.
that the horizontal velocity depends on the wave steepness, The results obtained from the numerical simulation are in
and the velocity increases with H/L = 0.0105 to 0.0210. excellent agreement with the analytical solution.
45 Page 12 of 15 Sådhanå (2022)47:45

from CFD. The comparison illustrates under different wave

steepness (H/L = 0.0105, 0.0158, 0.0210) at the free sur-
face z = 0 for the time period T = 3.23 s. The wave profile
has been analyzed between the time t = 10 s to 25 s at a
distance x = 27 m from the wavemaker. The almost con-
stant levels of wave crests and troughs also indicated that
no significant wave reflections have been found due to the
damping effect. The results obtained from the numerical
simulation and analytical methods are in good agreement
with each other. Figure 10(c) observed that the surface
profile is more flatter in the trough shown by the black
circle for the wave steepness H/L = 0.0210 compared to
figure 10(a) indicates that the nonlinearity effect improves
with the increase of wave steepness.
Figure 12. Comparison between the numerical and analytical Figure 11 shows the volume fraction of the time evolu-
free surface elevation at x = 27 m for H = 0.1 m. tion of the numerical free surface wave profile for the wave
height H = 0.1 m, wavelength L = 9.4826 m, and time
period T = 3.23 s. The red and blue regions are represented
Figure 7 shows that the vertical velocity (w) plotted by water and air, respectively, and the yellow line at the air-
against the time (t) under different wave steepness at the water interface represents the free surface (color interpre-
free surface z = 0 for T = 3.23 s and d/L = 0.105. The tation in this figure, the reader, referred to the web version
figures show that the vertical velocity increases with wave of this article). Plus sign is shown in the figures that the
steepness increase, and it is more at the wave steepness wave travels from left to right of the tank at different time-
H/L = 0.0210. Figure 7(c) shows that the curve is more steps t = 0 s, 1 s, 2 s, and 3 s. Corresponding coordinates
tilted than figure 7(a), indicating that the nonlinearity effect of the plus sign in the x-z plane concerning the time shown
increases with wave steepness. The numerical simulation in figures. The time steps are calculated from the simulation
results are in good agreement with the analytical solution. of the Ansys-Fluent software. It shows that nearly at
Figure 8(a) shows the velocity contour in the x-direction, t = 12 s (figure 11e), the wave profile reaches a steady-state
while figure 8(b) shows the velocity contour in the shape between the tank length 12 m to 36 m, and after that,
z-direction. The velocity contour is a volume-fraction weighted at t = 23 s and x = 40 m (figure 11f), the wave profile starts
velocity. For example, a cell containing 100% of the water to reduced its wave height gradually due to damping action.
has its velocity equal to the velocity of water, and the Figure 12 shows the time evolution of the numerical
velocity of an air-cell containing 0% of water is zero. The water-free surface elevation (measured from the still water
velocity contour of the water is plotted as the volume height) at x = 27 m away from the wavemaker compared to
fraction of the water in the air region is minimal and has zero the analytical results. It also shows that after nearly three
velocity. If air velocity is plotted, the air region would have periods, the wave approaches its steady profile, and the
a very high velocity, which does not make any physical wave becomes matches almost with the analytical results at
sense (Prasad et al [41]). Water particles in region 1 are the time t = 12 s at the exact location shown in figure 11(e).
traveling in the wave propagation direction, as shown in Similarly, figure 13 shows the time evolution of the
figure 8(a), while in region 2, particles are traveling in the numerical free surface wave profile for H = 0.2 m,
opposite direction. Due to this, it has been observed that the L = 9.4826 m, and T = 3.23 s. Plus sign shown in the
water particles pushed upwards, as shown in figure 8(b) in figures travels according to the time-steps t = 0 s, 1 s, 2 s,
region 3. As shown in figure 9 in region D, the velocity and 3 s. Corresponding coordinates of the plus sign in the
vector of the particles is easily understood this fact. It has x-z plane concerning the time are also shown in figures. It
been observed that a slight movement of the particles and shows that nearly at t = 13 s (figure 13e), the wave profile
velocity are recorded in region C. The critical point is that reaches a steady-state shape between the tank length 15 m
high-energy flow is concentrated near the free surface. It is to 39 m, and after that, at t = 25 s and x = 43 m
also essential that water particles are transported more in the (figure 13f), the wave profile starts to reduced its wave
wave propagation direction than in the vertical direction. height gradually due to damping action.
The reason for this velocity of particles recorded is high in Figure 14 shows the time evolution of the numerical
the x-direction than in the z-direction. This explanation water-free surface elevation at x = 27 m away from the
shows that the characteristic of intermediate water waves. wavemaker compared to the analytical results. It also shows
Figure 10 shows the time history of the surface elevation that after nearly three periods, the wave approaches its
(g) against the time (t) at d/L = 0.105. The result is pre- steady profile, and the wave becomes matches almost with
sented in figure 10 between the analytical solutions the analytical results at the time t = 13 s at the exact
obtained from Eq. (13) and the numerical solution obtained location shown in figure 13(e).
Sådhanå (2022)47:45 Page 13 of 15 45

Figure 13. Evolution of the numerical free surface profile in an NWT at different time-step t for H = 0.2 m.

7. Conclusions steepness on ocean waves at low Ursell numbers (Ur \ 18).

The inlet velocity method is used to generate the second-
In this study, second-order Stokes waves are simulating order water wave for achieving this objective. It has shown
numerically in an NWT in the intermediate depth of water. that the nonlinearity effect increases with the increase of
The simulation results of the model compared with the wave steepness and found horizontal and vertical velocity
analytical results. Six tests have been conducted at two increases with the increase of wave steepness. It has also
regions in an NWT to understand the wave theory. The shown in figure 6(c) that the horizontal velocity is more
main objective of this work, to study the effect of wave flatter in the trough shown by the black circle for the wave
45 Page 14 of 15 Sådhanå (2022)47:45

[6] Maâtoug M A and Ayadi M 2016 Numerical simulation of

the second-order Stokes theory using finite difference
method. Alex. Eng. J. 55: 3005–3013
[7] Dean R G and Dalrymple R A 1984 Water Wave Mechanics
for Engineers and Scientists. Prentice-Hall, Englewood
Cliffs
[8] Borgman L E and Chappelear J E 1957 The use of the
Stokes-Struik approximation for waves of finite height.
Coast. Eng. Proc. 6: 252–280
[9] Fenton J D 1985 A fifth-order Stokes theory for steady
waves. J. Waterw. Port Coast. Ocean Eng. 111: 216–234
[10] Chappelear J E 1961 Direct numerical calculation of wave
properties. J. Geophys. Res. 66: 501–508
[11] Dean R G 1965 Stream function representation of nonlinear
ocean waves. J. Geophys. Res. 70: 4561–4572
[12] Chaplin J R 1979 Developments of stream-function wave
Figure 14. Comparison between the numerical and analytical
theory. Coast. Eng. 3: 179–205
free surface elevation at x = 27 m for H = 0.2 m.
[13] Larsen J and Dancy H 1983 Open boundaries in short wave
simulations, a new approach. Coast. Eng. 7: 285–297
steepness H/L = 0.0210 compare to figure 6(a), proving [14] Brorsen M and Larsen J 1987 Source generation of nonlinear
that the second-order Stokes wave equation is more pro- gravity waves with the boundary integral equation method.
ductive with the increase of wave steepness. It has also Coast. Eng. 11: 93–113
[15] Li Y S, Liu S X, Yu Y X and Lai G Z 1999 Numerical
been observed that the horizontal velocity on the free sur-
modeling of Boussinesq equations by finite element method.
face (z = 0) is high compared to the bottom of the water
Coast. Eng. 37: 97–122
(z = -1). It has also shown that the water particles on the [16] Harlow F H and Welch J E 1965 Numerical calculation of
free surface are transported more in the wave propagation time-dependent viscous incompressible flow of fluid with
direction than in the vertical direction; hence high energy free surface. Phys. Fluids 8: 2182–2189
flow is concentrated near the free surface. [17] Park J C, Kim M H and Miyata H 2001 Three-dimensional
Furthermore, it has also been studied that nonlinearity on numerical wave tank simulations on fully nonlinear wave-
the surface elevation has increased with the increase of current-body interactions. J. Mar. Sci. Technol. 6: 70–82
wave steepness. Numerical results and analytical results are [18] Anbarsooz M, Passandideh-Fard M and Moghiman M 2013
within 1.3752%, highlighting the validity of the CFD code. Fully nonlinear viscous wave generation in numerical wave
This study of water waves based on second-order Stokes tanks. Ocean Eng. 59: 73–85
wave theory on the free surface improves the coastal [19] Wu Y T and Hsiao S C 2018 Generation of stable and
accurate solitary waves in a viscous numerical wave tank.
engineering application for further investigation and
Ocean Eng. 167: 102–113
extracting wave energy from the ocean. Applications of [20] Tang C T, Patel V C and Landweber L 1990 Viscous effects
WEC are Oscillating Water Columns (OWC) (Celik A and on propagation and reflection of solitary waves in shallow
Altunkaynak [54]), Hinged contour devices (Senol K and channels. J. Comput. Phys. 88: 86–113
Raessi M [55]), Overtopping devices (Goulart M M et al [21] Huang C J, Zhang E C and Lee J F 1998 Numerical
[56]). simulation of nonlinear viscous wavefields generated by
piston-type wavemaker. J. Eng. Mech. 124: 1110–1120
[22] Huang C J and Dong C M 2001 On the interaction of a solitary
wave and a submerged dike. Coast. Eng. 43: 265–286
References [23] Dong C M and Huang C J 2004 Generation and propagation of
water waves in a two- dimensional numerical viscous wave
[1] World Energy Council 2017 Ocean energy. [Link] flume. J. Waterw. Port Coast. Ocean Eng. 130: 143–153
[Link] [24] Wang H W, Huang C J and Wu J 2007 Simulation of a 3D
[2] Mustapa M A, Yaakob O B, Ahmed Y M, Rheem C, Koh K numerical viscous wave tank. J. Eng. Mech. 133: 761–772
K and Adnan F A 2017 Wave energy device and breakwater [25] Hirt C W and Nichols B D 1981 Volume of Fluid (VOF)
integration: a review. Renew. Sust. Energy Rev. 77: 43–58 method for the dynamics of free boundaries. J. Comput.
[3] Machado F M M, Lopes A M G and Ferreira A D 2018 Phys. 39: 201–225
Numerical simulation of regular waves: optimization of a [26] Lin P and Liu P L F 1998 A numerical study of breaking
numerical wave tank. Ocean Eng. 170: 89–99 waves in the surf zone. J. Fluid Mech. 359: 239–264
[4] Flick R E and Guza R T 1980 Paddle generated waves in [27] Troch P and De Rouck J 1999 An active wave generating–
laboratory channels. J. Waterw. Port Coast. Ocean Div. 106: absorbing boundary condition for VOF type numerical
79–97 model. Coast. Eng. 38: 223–247
[5] Moubayed W and Williams A 1993 Second-order bichro- [28] Li B and Fleming C A 2001 Three-dimensional model of
matic waves produced by generic planar wavemaker in a Navier-Stokes equations for water waves. J. Waterw. Port
two-dimensional wave flume. J. Fluids Struct. 8: 73–92 Coast. Ocean Eng. 127: 16–25
Sådhanå (2022)47:45 Page 15 of 15 45

[29] Karim M F, Tanimoto K and Hieu P D 2009 Modelling and Numerical simulation of monochromatic wave generated in
simulation of wave transformation in porous structures using laboratory: validation of a cfd code. In: 23 Congresso
VOF based two-phase flow model. Appl. Math. Model. 33: Nacional de Transport Aquaviario Construcao Naval Off-
343–360 shore, pp. 1–12
[30] Lin P and Liu P L F 1999 Internal wave-maker for Navier- [44] Kim S Y, Kim K M, Park J C, Jeon G M and Chun H H 2016
Stokes equations models. J. Waterw. Port Coast. Ocean Eng. Numerical simulation of wave and current interaction with a
125: 207–215 fixed offshore substructure. Int. J. Nav. Archit. Ocean Eng. 8:
[31] Kawasaki K 1999 Numerical simulation of breaking and 188–197
post-breaking wave deformation process around a sub- [45] Elangovan M 2011 Simulation of irregular waves by CFD.
merged breakwater. Coast. Eng. J. 41: 201–223 World Acad. Sci. Eng. Technol. 5: 1379–1383
[32] Hieu P D and Tanimoto K 2006 Verification of a VOF-based [46] Saincher S and Banerjee J 2015 Design of a numerical wave
two-phase flow model for wave breaking and wave–structure tank and wave flume for low steepness waves in deep and
interactions. Ocean Eng. 33: 1565–1588 intermediate water. Procedia Eng. 116: 221–228
[33] Hur D S and Mizutani N 2003 Numerical estimation of the [47] Morgan G C J, Zang J, Greaves D, Heath A, Whitlow C D
wave forces acting on a three-dimensional body on sub- and Young J R 2010 Using the rasinterFoam CFD model for
merged breakwater. Coast. Eng. 47: 329–345 wave transformation and coastal modeling. Coast. Eng.
[34] Zhi D O N G and Zhan J M 2009 Numerical modeling of Proc. 1: 1–9
wave evolution and run-up in shallow water. J. Hydrodyn. [48] Koo W C and Kim M H 2007 Fully nonlinear wave-body
21: 731–738 interactions with surface- piercing bodies. Ocean Eng. 34:
[35] Hafsia Z, Hadj M B, Lamloumi H and Maalel K 2009 1000–1012
Internal inlet for wave generation and absorption treatment. [49] Westphalen J, Greaves D M, Williams C J K, Hunt-Raby A
Coast. Eng. 9: 951–959 C and Zang J 2012 Focused waves and wave-structure
[36] Madsen O S 1971 On the generation of long waves. J. interaction in a numerical wave tank. Ocean Eng. 45: 9–21
Geophys. Res. 76: 8672–8683 [50] Fenton J D 1990 Nonlinear wave theories. The Sea Ocean
[37] Ursell F, Dean R G and Yu Y S 1960 Forced small- Eng. Sci. 9: 1–18
amplitude water waves: a comparison of theory and exper- [51] Hedges T S and Ursell, 1995 Regions of validity of
iment. J. Fluid Mech. 7: 33–52 analytical wave theories. Proc. Inst. Civ. Eng. Water Marit.
[38] Wood D J, Pedersen G K and Jensen A 2003 Modelling of Energy 112: 111–114
run up of steep non-breaking waves. Ocean Eng. 30: [52] Ursell F 1953 The long-wave paradox in the theory of
625–644 gravity waves. Math. Proc. Camb. Philos. Soc. 49: 685–694
[39] Finnegan W and Goggins J 2012 Numerical simulation of [53] Le Méhauté B 1976 An Introduction to Hydrodynamics and
linear water waves and wave-structure interaction. Ocean Water Waves, vol. 22, pp. 974–975. Springer-Verlag, New York
Eng. 43: 23–31 [54] Çelik A and Altunkaynak A 2019 Experimental investiga-
[40] Dao M H, Chew L W and Zhang Y 2018 Modeling physical tions on the performance of a fixed-oscillating water column
wave tank with flap paddle and porous beach in OpenFOAM. type wave energy converter. Energy 188: 116071
Ocean Eng. 154: 204–215 [55] Senol K and Raessi M 2019 Enhancing power extraction in
[41] Zhao X Z, Hu C H and Sun Z C 2010 Numerical simulation bottom-hinged flap-type wave energy converters through
of extreme wave generation using VOF method. J. Hydro- advanced power take-off techniques. Ocean Eng. 182: 248–258
dyn. 22: 466–477 [56] Goulart M M, Martins J C, Junior I C A, Gomes M D N,
[42] Prasad D D, Ahmed M R, Lee Y H and Sharma R N 2017 Souza J A, Rocha L A O, Isoldi L A and Santos E D D 2015
Validation of a piston-type wave-maker using Numerical Constructal design of an onshore overtopping device in real
Wave Tank. Ocean Eng. 131: 57–67 scale for two different depths. Mar. Syst. Ocean Technol. 10:
[43] Silva M C, Vitola M de A, Pinto W T and Levi C A 2010 120–129