A Numerical Model Based on Navier-Stokes Equations to Simulate Water Wave Propagation with Wave-Structure Interaction

1. Introduction

The interaction between waves and structures is a very important subject of the coastal engineering. Many numerical models have been developed over recent decades to analyze these types of cases, which involve phenomena that combine reflection, shoaling, refraction, diffraction, breaking and wave-wave interaction. These non-linear effects provide harmonic generation, including energy transferences with high complexity.

Models based on the Laplace equation assume the potential flow, in which the movement is irrotational and the flow is incompressible. Models based on the boundary element technique [1,2,3] and spectral methods [4,5] are some examples. This theory applies neither to viscous flows nor to situations in which there are flow separations, vortex generations and turbulences.

Other models, called depth-integrated models [6], based on a Boussinesq-type equation for variable depth, consider polynomial approximations for the vertical velocity distribution and vertical integration in the resulting equations at a certain depth. The simplified hypotheses, that include slight non linearity and dispersion, limit the applicability of this type of models to shallow and intermediate waters. Several researches have been developed to extend the applicability of these equations, including high order terms, to deeper water and strong non linearity cases in the last decades. Wave propagation phenomena, such as breaking, bottom friction and run-up, have also been included in these extended Boussinesq equation [7,8,9,10,11,12,13,14]. The accuracy of these type of models has recently been improved by the implementation of the multi-layer concept, in which the water column is divided into layers and a velocity profile is adopted to each one [15,16]. Although the accuracy has been improved significantly, the simplified hypotheses, related to the vertical integration along each layer, limit the use of these models to depths without strong variations.

Many efforts have been carried out to develop non hydrostatic models to prevent the difficulties of the Boussinesq models [17]. The non hydrostatic models capture the free surface movement using a function on the horizontal plane, requiring lower vertical discretization by comparison with those that use classic methods to describe the free surface. In some of these models, the pressure and velocity fields are decomposed by hydrostatic and non-hydrostatic pressure to improve their efficiency.

The numerical solution of the fully Navier-Stokes equations to determine the tridimensional velocity and pressure fields and the free surface position demands high computational cost, due to the large horizontal scale of many coastal engineering problems. However, in cases in which there are flow separation, vortex shedding and turbulence, these models provide more real results. There are several methods to capture the free surface movements, such as the arbitrary lagrangian eulerian (ALE) [22,46], the marked and cell [25], the volume of fluid [26,27] and the level-set methods [28].

This text describes a code (in Fortran 90 language) that integrates the Navier-Stokes equations using a fractional method to simulate 3D incompressible flow problems with free surface, named FLUINCO [29]. The model employs a semi implicit two-step Taylor Galerkin method to discretize the Navier-Stokes equations in time and space; uses the ALE method and a mesh velocity distribution technique to deal with free surface movement.

To show the applicability of the code, two study cases are analyzed: the wave propagation over a submerged horizontal cylinder and submerged trapezoidal breakwaters. In Section 2, the numerical model is described. Section 3 presents study cases, their results and discussion. Finally, Section 4 concludes the analysis.

2. Numerical model

2.1. Governing equations

The algorithm is based on the continuity equation, given by:

$$\frac{\partial \rho }{\partial t}=\frac{1}{{{\displaystyle c}}^2}\frac{\partial p}{\partial t}=-\frac{{{\displaystyle \partial U}}_i}{{{\displaystyle \partial x}}_i}\left(i=\text{1},\text{2},\text{3}\right),$$E1

and momentum equations, that are represented by the following equations according to the ALE formulation:

$$\frac{{{\displaystyle \partial U}}_i}{\partial t}=-\frac{{{\displaystyle \partial f}}_{ij}}{{{\displaystyle \partial x}}_j}+\frac{{{\displaystyle \partial \tau }}_{ij}}{{{\displaystyle \partial x}}_j}-\frac{\partial p}{{{\displaystyle \partial x}}_i}+{{\displaystyle \text{w}}}_j\frac{{{\displaystyle \partial U}}_i}{{{\displaystyle \partial x}}_j}\left(i,j=\text{1},\text{2},\text{3}\right),$$E2

where ρ is the specific mass, p is the pressure, $${{\displaystyle U}}_i=\rho {{\displaystyle \text{v}}}_i$$, $${{\displaystyle f}}_{ij}={{\displaystyle \text{v}}}_j\left(\rho {{\displaystyle \text{v}}}_i\right)={{\displaystyle \text{v}}}_j{{\displaystyle U}}_i$$, v_i are the velocity components, w_i the reference system velocity components and τ_ij is the viscous stress tensor (i,j.=1,2,3).

2.2. Semi-implicit two-step Taylor-Galerkin method

Basically, the algorithm consists of the following steps [30]:

1. Calculate non-corrected velocity at Δt/2, obtained by time discretization of Eq. (2), where the pressure term is at t instant, according to Eq. (3).

$${{\displaystyle \tilde{U}}}_i^{n+1/2}={{\displaystyle U}}_i^n-\frac{\Delta t}{2}\left(\frac{{{\displaystyle \partial f}}_{ij}^n}{{{\displaystyle \partial x}}_j}-\frac{{{\displaystyle \partial \tau }}_{ij}^n}{{{\displaystyle \partial x}}_j}+\frac{{{\displaystyle \partial p}}^n}{{{\displaystyle \partial x}}_i}-{{\displaystyle \rho g}}_i-{{\displaystyle \text{w}}}_j^n\frac{{{\displaystyle \partial U}}_i^n}{{{\displaystyle \partial x}}_i}\right)\left(i,j=\text{1},\text{2},\text{3}\right),$$E3

where g_i are the gravity acceleration components.

1. Update the pressure p at t+Δt, obtained by time discretization of Eq. (1), given by the Poisson equation:

$$\frac{1}{{{\displaystyle c}}^2}\Delta p=$$

$$-\Delta t\left[\frac{{{\displaystyle \partial \tilde{U}}}_i^{n+1/2}}{{{\displaystyle \partial x}}_i}-\frac{\Delta t}{4}\frac{\partial }{{{\displaystyle \partial x}}_i}\frac{\partial \Delta p}{{{\displaystyle \partial x}}_i}\right]$$E4

where $$\Delta p={{\displaystyle p}}^{n+1}-{{\displaystyle p}}^n$$and i = 1,2,3.

1. Correct the velocity at t+Δt/2, adding the pressure variation term from t to t+Δt/2, according to the equation:

$${{\displaystyle U}}_i^{n+1/2}={{\displaystyle \tilde{U}}}_i^{n+1/2}-\frac{\Delta t}{4}\frac{\partial \Delta p}{{{\displaystyle \partial x}}_i}$$ (i =1,2,3),

1. Calculate the velocity at t+Δt using variables updated in the previous steps as follows:

$${{\displaystyle U}}_i^{n+1}={{\displaystyle U}}_i^n$$ $$-\Delta t\left(\frac{{{\displaystyle \partial f}}_{ij}^{n+1/2}}{{{\displaystyle \partial x}}_j}-\frac{{{\displaystyle \partial \tau }}_{ij}^{n+1/2}}{{{\displaystyle \partial x}}_j}+\frac{{{\displaystyle \partial p}}^{n+1/2}}{{{\displaystyle \partial x}}_i}-{{\displaystyle \text{w}}}_j^{n+1/2}\frac{{{\displaystyle \partial U}}_i^{n+1/2}}{{{\displaystyle \partial x}}_i}-{{\displaystyle \rho g}}_i\right)$$ (i,j =1,2,3),

2.3. Space discretization

The classical Galerkin weighted residual method is applied to the space discretization by using a tetrahedron element. In the variables at t+Δt/2 instant, a constant shape function PE is used, and in the variables at t and t+Δt, a linear shape function N is employed. By applying this procedure to Eq. (3), (4), (5) and (6), the following expressions in the matrix form are obtained [29]:

$${{\displaystyle \Omega }}_E^{n+1/2}{{\displaystyle \tilde{\overline{U}}}}_i^{n+1/2}={{\displaystyle C}}_{}^{}{{\displaystyle \overline{U}}}_i^n$$

$$-\frac{\Delta t}{2}\left[{{\displaystyle L}}_j^{}\left({{\displaystyle \overline{f}}}_{ij}^n-{{\displaystyle \overline{\tau }}}_{ij}^n\right)+{{\displaystyle L}}_i^{}{{\displaystyle \overline{p}}}_{}^n-{{\displaystyle T}}_{}^{}{{\displaystyle \overline{U}}}_i^n-{{\displaystyle \Omega }}_E^{n+1/2}{{\displaystyle \overline{\rho }\overline{g}}}_i\right]$$E5

$$\left(\tilde{M}+\frac{{{\displaystyle \Delta t}}^2}{4}H\right)\Delta \overline{p}=\Delta t\left({{\displaystyle L}}_i^T{{\displaystyle \tilde{\overline{U}}}}_i^{n+1/2}+{{\displaystyle f}}_a\right)$$E6

$${{\displaystyle \overline{U}}}_i^{n+1/2}={{\displaystyle \tilde{\overline{U}}}}_i^{n+1/2}-\frac{\Delta t}{4{{\displaystyle \Omega }}_E}{{\displaystyle L}}_i^{}\Delta \overline{p}$$E7

$${{\displaystyle M}}^{n+1}{{\displaystyle \overline{U}}}_i^{n+1}={{\displaystyle M}}^n{{\displaystyle \overline{U}}}_i^n+\Delta t{{\displaystyle L}}_j^T\left({{\displaystyle \overline{f}}}_{ij}^{n+1/2}-{{\displaystyle \overline{w}}}_j^{n+1/2}{{\displaystyle \overline{U}}}_i^{n+1/2}\right)$$E8

$$-\Delta t{{\displaystyle Q}}_j^{}{{\displaystyle \overline{\tau }}}_{ij}^n+\Delta t{{\displaystyle Q}}_i^{}\left({{\displaystyle \overline{p}}}_{}^n+\Delta \overline{p}/2\right)+\Delta t{{\displaystyle S}}_{bi}-\Delta t{{\displaystyle C}}^T{{\displaystyle \overline{g}}}_i$$E9

where variables with upper bars at n and n+1 instants indicate nodal values, while those at n+1/2 instant represent constant values in the element. The matrices and vectors from Eq. (7) to (10) are volume and surface integrals that can be seen in detail in [30].

Equation (8) is solved using the conjugated gradient method with diagonal pre-conditioning [31]. In Eq. (10), the consistent mass matrix is substituted by the lumped mass matrix, and then this equation is solved iteratively.

The scheme is conditionally stable and the local stability condition for the element E is given by

$${{\displaystyle \Delta t}}_E\le \beta {{\displaystyle h}}_E/\left|u\right|$$E10

where h_E is the characteristic element size, β is the safety factor and u is the fluid velocity.

2.4. Mesh movement

The free surface is the interface between two fluids, water and air, where atmospheric pressure is considered constant (generally the reference value is null). In this interface, the kinematic free surface boundary condition (KFSBC) is imposed. By using the ALE formulation, it is expressed as:

$$\frac{\partial \eta }{\partial t}+{{\displaystyle {}^{(s)}\text{v}}}_i^{}\frac{\partial \eta }{{{\displaystyle \partial x}}_i}={{\displaystyle {}^{(s)}\text{v}}}_3^{}$$E11

where η is the free surface elevation, $${{\displaystyle {}^{(s)}\text{v}}}_3^{}$$ is the vertical fluid velocity component and $${{\displaystyle {}^{(s)}\text{v}}}_i^{}$$ (i=1,2) are the horizontal fluid velocity components in the free surface. The eulerian formulation is used in the x₁ and x₂ directions (horizontal plane) while the ALE formulation is employed in the x₃ or vertical direction.

The time discretization of KFSBC is carried out in the same way as the one for the momentum equations as presented before. After applying expansion in Taylor series, the expressions for η at n+1/2 (first step) and n+1 (second step) instants are obtained:

$${\eta }^{n+1/2}={\eta }^n+\frac{\Delta t}{2}(\begin{array}{c}\\ \end{array}{}^{(s)}\text{v}{}_3-{}^{(s)}\text{v}{}_1\frac{\partial \eta }{\partial x_1}-{}^{(s)}\text{v}{}_2\frac{\partial \eta }{\partial x_2}{\begin{array}{c}\\ \end{array})}^n$$E12

$${\eta }^{n+1}={\eta }^n+\Delta t(\begin{array}{c}\\ \end{array}{}^{(s)}\text{v}{}_3-{}^{(s)}\text{v}{}_1\frac{\partial \eta }{\partial x_1}-{}^{(s)}\text{v}{}_2\frac{\partial \eta }{\partial x_2}{\begin{array}{c}\\ \end{array})}^{n+1/2}$$E13

Linear triangular elements coincident with the face of the tetrahedral elements on the free surface are used to the space discretization by applying the Galerkin method.

The mesh velocity vertical component w₃ is computed to diminish element distortions, keeping prescribed velocities on moving (free surface) and stationary (bottom) boundary surfaces. The mesh movement algorithm adopted in this paper uses a smoothing procedure for the velocities based on these boundary surfaces. The updating of the mesh velocity at point i of the finite element domain is based on the mesh velocity of the points j that belong to the boundary surfaces, and is expressed in the following way [32]:

$${{\displaystyle \text{w}}}_3^i=\frac{{\displaystyle \sum _{j=1}^{ns}{{\displaystyle a}}_{ij}{{\displaystyle \text{w}}}_3^j}}{{\displaystyle \sum _{j=1}^{ns}{{\displaystyle a}}_{ij}}}$$E14

where ns is the total number of points belonging to the boundary surfaces and a_ij are the influence coefficients between the point i inside the domain and the point j on the boundary surface given by the following expression:

$${{\displaystyle a}}_{ij}=\frac{1}{{{\displaystyle d}}_{ij}^4}$$E15

with d_ij being the distance between points i and j. In other words, a_ij represents the weight that every point j on the boundary surface has on the value of the mesh velocity at points i inside the domain. When d_ij is low, a_ij has a high value, favouring the influence of points i, located closer to the boundary surface containing point j.

The free surface elevation, the mesh velocity and the vertical coordinate are updated according to the following steps:

1. Calculate η ^n+1/2 and $${{\displaystyle \tilde{U}}}_i^{n+1/2}$$, Eq. (13) and Eq. (3), respectively.

2. Calculate $$\Delta p$$, Eq. (4).

3. Calculate $${{\displaystyle U}}_i^{n+1/2}$$, Eq. (3).

4. Calculate $${{\displaystyle U}}_i^{n+1}$$, Eq. (6).

5. Calculate η ⁿ⁺¹, Eq. (14).

6. Update the mesh velocity w₃ and the vertical coordinate x₃:

7. Calculate the mesh velocity in the free surface at t + Δt: $${{\displaystyle {}^{(S)}w}}_3^{n+1}$$ = $$\left({{\displaystyle \eta }}^{n+1}-{{\displaystyle \eta }}^n\right)/\Delta t$$.

8. Calculate the mesh velocity in the interior of the domain at n+1 e n+1/2 by using Eq. (13) and $${{\displaystyle \text{w}}}_3^{n+1/2}=\frac{{{\displaystyle (\text{w}}}_3^{n+1}+{{\displaystyle \text{w}}}_3^n)}{2}$$, respectively.

9. Update the vertical coordinates in the interior of the domain: $${{\displaystyle x}}_3^{n+1/2}={{\displaystyle x}}_3^n+{{\displaystyle \text{w}}}_3^n\frac{\Delta t}{2}$$, $${{\displaystyle x}}_3^{n+1}={{\displaystyle x}}_3^n+{{\displaystyle \text{w}}}_3^{n+1/2}\Delta \text{t}$$.

2.5. Wave generation and radiation conditions

The wave generation is considered imposing the free surface elevation and the fluid velocity components to each time step directly, considering the linear wave equations [33].

The Flather’s radiation condition [34] is used to deal with open boundaries. In this method, the Sommerfeld condition to free surface elevation is combined with one-dimensional version of the continuity equation. Then, the normal velocity of the boundary can be expressed by:

$$u=\eta \sqrt{\frac{g}{h}}$$E16

where g is the gravity acceleration and h is the depth.

3. Study cases

3.1. Submerged cylinder

The interaction among regular waves and submerged circular cylinders, with their axes parallel to the crests of the incident waves, has been studied analytically, experimentally and numerically by many authors. The presence of an obstacle near the free surface may cause reflected and modified transmitted waves. These phenomena depend on the characteristics of the incident wave, the obstacle geometry and the depth. Many studies of this interaction are available to provide a good example to validate numerical codes.

The first study was developed by [35] and, after that, by [36]. Considering a linear behavior, these authors showed that (a) the cylinder does not reflect any energy, regardless of its ray, depth or wave frequency; (b) the transmitted waves are out of phase, but their amplitudes are the same. Chaplin [37] studied the nonlinear forces and characteristics of the reflected and transmitted waves experimentally. He showed that the reflection is negligible up to the third order. This author and Schonberg and Chaplin [38] presented many experimental and numerical studies concerning the nonlinear interaction among waves and submerged cylinders. A detailed review of analyses for this case can be found in Paixão Conde et al. [39].

This case considers a 5.2 m long and 0.425 m deep channel with a submerged cylinder of r = 0.025m positioned 1.60 m from the wave generator (Figure 1). The cylinder center is 0.075 m (3r) from the free surface. The frequency wave is f = 1.4Hz; its amplitude is a = 0.0119 m and its wavelength is L=0.796 m, characterizing a deep water case.

Figure 1.

Geometry of the horizontal cylinder case.

Table 1 shows the period, the frequency and the wavelength for the fundamental frequency and its 2nd, 3rd and 4th harmonics, according to the linear wave theory.

	Fundamental	2^nd harmonic	3^rd harmonic	4^th harmonic
T (s)	0.7143	0.3571	0.2381	0.1786
f(Hz)	1.4	2.8	4.2	5.6
L(m)	0.796	0.199	0.0885	0.0498

Table 1.

Period, frequency and wavelength of the fundamental frequency and its 2nd, 3rd and 4th harmonics in the horizontal cylinder case.

The mesh, with 173900 nodes and 515623 elements, has one layer of elements in the transversal direction. The average element size on the cylinder boundary is 0.0015.m (105 divisions in the circumference). The element size diminishes from the ends to the region near the cylinder and from the bottom to the free surface. The element sizes on the end where the wave generator is located and on the opposite end are 0.015 m (53 points per fundamental wavelength) and 0.02 m (40 points per fundamental wavelength), respectively. On the bottom, 0.0015m is also used.

The initial conditions are: null velocity components in all domain and hydrostatic pressure (null on the free surface). The wave is generated by imposing the surface elevation and the velocity components. The non-slip condition is imposed to the bottom and to the cylinder wall. The time step is 0.0002s, which satisfies the Courant condition.

Figure 2 shows the free surface elevation obtained by the code and experimental tests, where xc is the horizontal coordinate of the cylinder center. In general, there is agreement between numerical and experimental results [39]. We can notice the free surface disturbance downstream the cylinder. When (x-xc)/L is above 1.7, the numerical results are smoother than the experimental ones, showing the necessity of a refinement in this region.

Figure 3 shows a comparison among numerical and experimental results in terms of free surface elevation on four gauges located at (x-xc)/L equal to -0.503, 0.0692, 0.509 and 1.264 (there is only a numerical result on the first gauge). We can observe the similarity among numerical and experimental results.

Figure 4 shows the streamlines and the velocity modulus distribution at the same instant used in Figure 2. Recirculation and separation cannot be observed at downstream. Due to the oscillatory flow behavior, there is no time for recirculation productions. We can notice the flow acceleration near the cylinder due to the boundary layer effect. The viscous effects have only local influence, without modifying the velocity field far from it.

Figure 2.

Free surface elevation in the submerged cylinder case (Numerical  ; Experimental ■).

Figure 3.

Free surface elevations on the gauges located at (x-xc)/L equal to -0.503, 0.0692, 0.509 and 1.264 in the horizontal cylinder case (Numerical  ; Experimental ●).

In Figures 5 and 6, velocity component profiles, u and v, on the same gauge positions are presented. These profiles were constructed at the same instant as that used in Figure 2. According to the linear theory, the maximum value for both horizontal and vertical components is equal to 0.105 m/s. For horizontal components, these values occur on the crest and the trough, while for vertical ones, these values occur on upward and downward zero-crossings. When one component is the maximum, another is null, because the phase difference is 90 degrees.

Gauge 1 ((x-xc)/L = -0.503) is located upstream, near the wave crest; no significant disturbance in u and v profiles is observed. The horizontal velocity component is positive and its maximum value is similar to the theoretical value in the crest. The wave trough passes by gauge 2; the vertical velocity component presents low values and the horizontal velocity component has negative values, reaching the maximum absolute value close to the theoretical ones (0.105 m/s). Gauge 3 is located near the first crest upstream, resulting in high horizontal component values. Finally, gauge 4 is on a region between the trough and upward zero-crossing. Both component profiles are negative and the vertical component magnitude shows how close the gauge is to upward zero-crossing.

The non-slip boundary condition on the bottom does not change the general behavior of the wave propagation, because this case is considered a deep water one.

Figure 4.

Streamlines and velocity modulus at the instant in which the free-surface elevation was captured (Figure 2) in the horizontal cylinder case.

Figure 5.

Horizontal velocity components at the same instant used in Figure 2 along the depth on gauges located at (x-xc)/L equal to -0.503, 0.0692, 0.509 and 1.264 in the horizontal cylinder case.

Figure 7 shows the frequency spectra for these four gauges distributed along the channel. In all cases, the energy is concentrated on the fundamental frequency and its harmonic waves. On gauge 1, the fundamental frequency presents most energy and the second harmonic shows a little value. On gauges 2 and 4, located upstream, significant energy up to the third harmonic appears, similar to the experimental results.

Figure 6.

Vertical velocity components at the same instant used in Figure 2 along the depth on gauges located at (x-xc)/L equal to -0.503, 0.0692, 0.509 and 1.264 in the horizontal cylinder case.

Figure 7.

Frequency spectra on gauges located at (x-xc)/L equal to -0.503, 0.0692, 0.509 e 1.264 in the horizontal cylinder case.

3.2. Submerged trapezoidal breakwaters

Whatever the numerical model characteristics, the simulation of wave propagation over submerged breakwaters are important tests to validate wave propagation models. In these cases, the harmonic generation [40,41] and the vortex formation, depending on the geometry [42], also occur. When waves propagate in deep waters over a submerged obstacle, part of the wave energy is transferred from the primary wave component to their harmonics, contributing to increase non linearity. Harmonic generation phenomena that occur when waves propagate over obstacles, such as natural reefs, were studied theoretically [6], experimentally [43,44,45] and numerically [46,23,17,47,48,45,49,50]. In some situations, the correct simulation of the flow can only be figured out considering the viscosity effects [51]. Huang and Dong [42] studied the interaction between solitary waves and rectangular submerged breakwaters using a model based on 2D Navier-Stokes equations and concluded that the flow around the breakwater is laminar, without turbulence. The experimental studies carried out by Ting and Kim [51] and Zhuang and Lee [52] show that velocity fluctuations do not exist around the breakwater.

Two different configurations of the trapezoidal breakwaters, with different level of non-linearity, are used to test the behaviour of the numerical model. In the first case, the downstream and upstream slopes are 1:20 and 1:10, respectively [44]. In the second one, both slopes are 1:2 [45], where the non-linear effects are more significant.

3.2.1. Breakwater with slopes 1:20 and 1:10

Figure 8 shows the channel and the submerged breakwater geometries, and the position of the gauges. The channel is 23m in length, 0.4m and 0.1m are the maximum and the minimum depths, respectively. In the channel entrance, a monochromatic wave is generated with a period of 2.02s and an amplitude of 0.01m.

Table 2 presents some parameters for this case study, in which H is the wave length, h is the depth, k=2π/L is the wave number and Ur.=.gHT²/h² is the Ursell number, where T is the wave period. H/h, even on the platform, has small values in comparison with breaking limit of approximately 0.8 [33]. The case involves intermediate water for the channel (0.314.<.kh.<.3.142) and shallow water for the platform (kh < 0.314). Ursell numbers show that the non-linear effects on the platform are more intensive.

Figure 8.

Channel geometry for the 1:20 and 1:10 breakwater

	H/h	kh	Ur
Channel (h = 0.4m)	0.050	0.674	5.0
Platform (h = 0.1m)	0.259	0.318	103.6

Table 2.

Wave parameters for the 1:20 and 1:10 breakwater

Table 3 presents periods, frequencies and wavelengths concerning the fundamental frequency and the harmonic components that occur along the wave propagation. The wavelength was estimated according to the dispersion equation of the linear theory. These values are references to determine discretizations in time and space to be used in the modeling.

	Fundamental	2^nd harmonic	3^rd harmonic	4^th harmonic
Period (s)	2.02	1.01	0.67	0.50
Frequency (Hz)	0.50	1.00	1.50	2.00
Wavelength (m)	3.73	1.46	0.70	0.39

Table 3.

Period, frequency and wavelength concerning the fundamental frequency, and 2^nd, 3^rd and 4^th harmonics for the 1:20 and 1:10 breakwater

FLUINCO used a mesh with 88700 elements and 37296 nodes. Twenty layers of elements were used in vertical direction, where small elements are located near the bottom and the free surface. Along the channel, the element sizes vary from Δx.=.0.08m in the boundary to Δx.=.0.025m around the platform. In the transversal direction, only one layer of elements is used, because the behavior of the flow is bi-dimensional. In the entrance of the domain, the wave generation condition is imposed while at the end the radiation condition is imposed. The velocity components are null on the bottom and the KFSBC is imposed in the free surface. The velocity component perpendicular to the surface is null for lateral walls (symmetry condition). As an initial condition, the velocity field is null and the pressure one is hydrostatic. The time step is 0.003s, a fact that satisfies the Courant stability condition.

Figure 9 shows the free surface elevations in gauge 3, located downstream the breakwater (x=5.7m); in gauge 6, on the platform (x=13.5m); in gauge 8, in the middle of the upstream slope (x=15.7m); and in gauge 11, on the upstream and far from the breakwater (x=23m). Results obtained by numerical model are compared with the experimental ones presented by Dingemans [44].

In general, there is good agreement between numerical results and experimental ones in gauges 3 and 6. In gauge 6, FLUINCO presents slightly smooth surface deformation. In gauges 8 and 11, corresponding to downstream, the nonlinear effects are more significant. The deformations in gauge 8 are well represented by FLUINCO; although the results get closer to the experimental ones in some regions, there are difficulties in representing the deformations related to higher harmonics, possibly due to the lack of an appropriate discretization to capture the nonlinear phenomena.

Figure 9.

Free surface elevation of the 1:20 and 1:10 breakwater

Figure 10 shows the frequency spectra obtained by the model in the gauges and a comparison with the experimental results. The differences found in the free surface elevation are confirmed in Figure 10, which shows differences in the intensity of harmonic components, mainly in gauges located at the end of the channel. The numerical model adequately simulate the position of the peaks of the fundamental frequency and the harmonic components throughout the domain. However, there are some differences in the amplitude of these peaks, especially in gauges 8 and 11.

Figure 11 presents the streamlines around the upstream slope of the breakwater in eleven instants completing one wave period obtained by FLUINCO. We can observe that the flow separation and the vortex do not exist at all instants, due to the mild inclination of the upstream slope.

Figure 10.

Numerical and experimental frequency spectrum in the gauges of the breakwater 1:20 and 1:10

Figure 11.

Streamlines of the 1:20 and 1:10 breakwater

3.2.2. Breakwater with slopes 1:2

In this case, the length of the channel is 35m and the maximum and the minimum depths are 0.5m and 0.15m, respectively (See Figure 13). In the entrance of the channel, a monochromatic wave is generated with a period of 2.68s, related to a wavelength of 5.66m in the channel, and an amplitude of 0.025m. This problem is case 6 studied by Ohyama et al. [45] who analyzed six different types of waves experimentally. Table 4 shows some parameters that characterize the problem, calculated according to the linear theory. The Ursell number on the platform is 210, indicating the strong non-linearity in this region. Parameter H/h shows that breaking does not even occur on the platform.

Figure 12.

Geometry of the channel for the 1:2 breakwater

	H/h	kh	Ur
Channel (h = 0.5m)	0.100	0.555	14.1
Platform (h = 0.15m)	0.355	0.294	210.0

Table 4.

Wave parameters for the 1:2 breakwater

Table 5 shows periods, frequencies and wavelengths concerning the fundamental frequency and the harmonic components that occur along the wave propagation.

A mesh with 120200 elements and 50526 nodes was used for FLUINCO in this simulation. The element sizes along the channel vary between dx=0.08m at the ends and dx=0.01m on the platform. The boundary and the initial conditions are similar to the ones in the previous case, and 0.002s was the time step.

	Fundamental	2nd harmonic	3rd harmonic	4th harmonic
Period (s)	2.68	1.34	0.89	0.67
Frequency (Hz)	0.373	0.746	1.124	1.493
Wavelength (m)	5.66	2.42	1.22	0.70

Table 5.

Period, frequency and wavelength related to the fundamental frequency, and 2nd, 3rd, and 4th harmonics for the 1:2 breakwaters

Figure 13 shows the free surface elevations in gauges 3 and 5 (gauge positions are indicated in Figure 5). Numerical results are compared with the experimental ones presented by Ohyama et al. [45]. The FLUINCO model represents the surface deformation recorded in gauge 3 well. The deformations of gauge 5 indicate that the nonlinearity increases. In this case, FLUINCO captures the variation of the surface elevation more accurately.

Figure 14 shows frequency spectra obtained in gauges 3 and 5. The fundamental and the harmonic waves are well represented by the models, but their amplitudes differ. The FLUINCO results are closer for the two gauges.

Figure 13.

Free surface elevation for the 1:2 breakwater in gauges 3 and 5.

Figure 14.

Breakwater case. Frequency spectra in gauges 3 and 5.

Streamlines during one wave period obtained by FLUINCO are presented in Figure 15. Unlike the previous case, a vortex, located between the upstream slope and the bottom, occurred during part of the wave period.

Figure 15.

Streamlines of the 1:2 breakwater

4. Conclusions

In this text, we showed a model, named FLUINCO, capable of simulating flows on free surface. It is based on the semi-implicit two-step Taylor- Galerkin method to integrate Navier-Stokes equations in time and space. An ALE formulation is employed to describe the free surface movement. The methodology was validated in two study cases: the wave propagation over a submerged horizontal cylinder and submerged trapezoidal breakwaters. Both study cases showed the application of a Navier-Stokes based code, which considers accurately vertical flow effects, in the wave-submerged structure interaction problems.

In the case of the submerged horizontal cylinder, the free surface elevations and the velocity profiles obtained by the model were similar to experimental ones [39]. The numerical results presented a slight free surface deformation downstream, possibly because of the lack of refinement that caused numerical diffusion. In this case, the viscous effects influenced the flow behavior locally whereas the viscosity was not important far from the cylinder. The non-slip condition on the bottom did not modify the wave propagation significantly because it is a deep water case.

In the case of trapezoidal breakwater, two analyses were carried out for different upstream and downstream slopes. The first analysis deals with upstream and downstream slopes of 1:20 and 1:10, respectively. The results obtained by the model were compared with Dingemans’ experimental data [44]. A comparison of the surface elevations and the energy spectrum for some gauges along the channel showed that the model provided good results. Although the FLUINCO results have been somewhat smoothed, they were closer to the experimental ones, including the ones in the gauges placed downstream, where nonlinear effects are more significant. Streamlines over a wave period showed that there was no flow separation in this case.

The second analysis consists of two 1:2 breakwater slopes; it showed a strong influence of nonlinear effects on the results of the surface elevation and the energy spectrum. The numerical results were compared with experimental ones presented by Ohyama et al [45]. The vertical velocity field obtained by FLUINCO showed that a vortex of non-turbulent origin was formed in the flow. The model obtained results closer to the experimental ones, including the ones downstream of the breakwater, where the nonlinearity effects are more significant. Both breakwater analyses showed that FLUINCO captures the nonlinear effects of the flow accurately, due to the fact that this model considers the influence of the vertical circulation in the flow.

Acknowledgement

The author acknowledges Conceição Juana Fortes, Eric Didier and José Manuel Paixão Conde and the partnership between FURG and Laboratório Nacional de Engenharia Civil (LNEC). The author also thanks the support of Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq - project 303308/2009-5).