Abstract

Simulation of irregular waves over submerged obstacle on a NURBS potential numerical wave tank

Arash Abbasnia^I; Mahmoud Ghiasi^II

^IDepartement of Maritime Engineering, Amirkabir University of Technology, Tehran, Iran, 15875-4413, a_abbasnia@aut.ac.ir

^IIDepartement of Maritime Engineering, Amirkabir University of Technology, Tehran, Iran, 15875-4413, mghiasi@aut.ac.ir

ABSTRACT

In this paper, a fully non-linear three-dimensional Numerical Wave Tank (NWT) is developed for studying propagation and scattering of non-linear random sea wave over bottom submerged bars. The simulation of fully non-linear free-surface is based on Non-Uniform Rational B-Spline formulation (NURBS) as a novel approach and Mixed Eulerian-Lagrangian method (MEL). High-order boundary integral equation is used to solve the Laplace equation in the Eulerian frame. To update the free-surface, time marching approach including material node method and fourth order Runge-Kutta time integration scheme is used. To obtain appropriate numerical solutions for wave propagation problem, damping zone is set at the downstream. Also, the NURBS approximation is employed to evaluate the velocity of the free-surface particles. Propagation of regular and irregular waves in a NWT is investigated and compared with the available experimental and numerical data. Transmission of the random sea wave over submerged bars is also compared with the experimental and prior numerical studies.

Keywords: Irregular wave, Numerical wave tank, Potential flow, Mixed Eulerian-Lagtrangian, Non-Uniform rational B-spline

1 INTRODUCTION

Description of free-surface in the numerical simulation of free-surface flow is important to obtain accurate solutions for the problem. Numerous remedies have been developed to find accurate approximation of free-surface. They sometimes are complicated to implement and some of them are time consuming. Polynomial interpolation functions have been widely used as shape function to define marine structures and boundary geometry. Eight-node quadrilateral elements were used by Ning and Teng (2006) and Ning et al. (2009) to simulate free surface with biquadratic shape function in time marching scheme. To achieve more accurate results, higher order shape function can be used by more nodal points. Numerical procedure becomes time consuming with increment of nodal points and by decrement of the mesh size, numerical instability may occur. Also, influence coefficient matrix of boundary integral may become near singular by use of the fine mesh size.

A two dimensional potential wave tank was developed by Tang and Huang (2008) for propagation of the second-order wave and simulation of Bragg bottom effect on the free surface evolution. They employed the linear boundary element method to solve the non-linear problem. Third order shape function with twelve-node quadrilateral curvilinear elements was employed by Shao (2010) in weakly non-linear wave-body interaction. Flow field around an offshore monopile was computed using FEM method by Li et al. (2011) where eight-node quadrilateral elements with second-order shape function were used to describe the hull boundary geometry.

In the two past decades, Non-Uniform Rational B-Spline surface (NURBS) has been widely applied to define the complex shape of marine structures. Diffraction problem of floating body was solved by Datta & Sen (2006) in which hull shape was approximated with NURBS. Review of the past studies shows that using of B-spline surface for free-surface modeling is a new approach in time domain simulation. Desingularized boundary integral equation, in both direct and indirect formulations was given by Cao et al. (1991). Three-dimensional NURBS indirect Boundary Integral Equation (BIE) was proposed by Gao and Zou (2008). Two dimensional B-spline boundary element formulations for different degree of continuity of the geometric boundary and variables are developed by Cabral et al. (1990, 1991). A two dimensional super-parametric boundary element method was formulated by Damanpack et al. (2013) to solve Poisson equation for bending analysis of thin functionally graded plates based on Green second identity. Two dimensional potential numerical wave tank based on NURBS boundary integral equation is developed by Abbasnia and Ghiasi (2014) to simulate interaction between non-linear wave and truncated breakwater.

To keep the stability of solutions, free-surface particle velocity has to be evaluated accurately. To obtain tangential velocity for isoparametric linear elements, double node approach was developed by Grilli and Svendsen (1990). A polynomial formulation was presented by Grilli et al. (2001) for biquadratic and high-order curvilinear elements. For time marching scheme, first-order and second-order finite difference formulation in time was employed in Numerical Wave Tank (NWT) by Wu et al. (2005) and Xiao et al. (2009) as low-order time integration methods. Fourth-order Runge-Kutta method was used by Koo (2003), and fifth-order Runge-Kutta-Gil and fourth-order Adams-Bashforth-Moulton methods were used by Zhang et al. (2005) as high-order integration scheme to update the free-surface.

Different types of wave-makers on the inflow boundary were reviewed by Tanizawa (2000) and Newman (2010). In the opposite side of the wave-maker, the wave absorber is adopted to prevent wave reflection from the end wall. Artificial damping zone has also been applied by Cointe (1991) on the free surface boundary to minimize the wall effect on computational domain. During the free-surface simulation of nonlinear waves, the non-physical saw-tooth instability may occur. Instabilities may also occur due to variable mesh size or natural singular treatment at the intersection of the wave-maker and the free-surface. To treat the so-called saw-tooth instability, different smoothing schemes can be used such as Chebyshev five-point smoothing scheme presented by Koo and Kim (2004) and the B-spline smoothing scheme applied by Tanizawa (2000).

This paper is mainly focused on the development of three-dimensional potential NWT based on NURBS approximation. Desingularized direct boundary integral coupled with NURBS formulation is developed to solve the boundary value problem in the Eulerian frame. In the material node approach based on Mixed Eulerian-Lagrangian (MEL) method presented by Longuet-Higgins and Cokelet (1976), nodal points on the free surface are allowed to move freely with the free-surface particles and are traced in the Lagrangian frame. Derivatives of B-spline basis function are used to obtain tangential derivations of potential. The fourth-order Runge-Kutta time integration scheme is applied into the fully non-linear free-surface boundary condition to obtain instant position of the moving boundary for the next time step. To avoid non-physical saw-tooth instability, five-point Chebyshev smoothing scheme is applied on the obtained instantaneous variables for every few time steps.

Series of tests are performed to verify the present numerical procedure. Superposition of several first-order and second-order regular waves at a certain spot is measured as focused wave. Three random waves are propagated in the present NWT and their propagations are compared with the experimental measurements and available numerical results. Effect of bottom profile on the free surface elevation is also investigated. For a submerged bar on the bottom the wave transformations are investigated and compared with the experimental data and the numerical results.

2 NUMERICAL MODEL

Consider a three-dimensional numerical wave tank with depth d , width B and length L . A layout of the computational domain is illustrated in Figure 1. A wave generator is placed on the upstream wall of the tank. A right-hand Cartesian coordinate system (Oxyz) is defined on the mean water surface where the origin is on the corner of the tank, the x -axis is laid on the length of the tank and the positive z -axis directed vertically upwards. An artificial damping zone is defined at adjacent of the end tank wall.

2.1 Governing equation and boundary conditions

It is assumed that the fluid is homogeneous, incompressible and inviscid and the flow is irrotational and surface tension on the free-surface is neglected. Therefore, velocity potential ϕ(x, y, z, t) can be introduced. Velocity potential satisfies the Laplace equation in the fluid domain Ω:

There are fully non-linear boundary conditions, Kinematic Free Surface Boundary Condition (KFSBC) and Dynamic Free Surface Boundary Condition (DFSBC) on the free-surface (S_f) which can be written as the following (Tanizawa, 2000):

where, η(x, y, t) is the wave elevation measured from still water level, t shows the time of simulation and g is the gravitational acceleration. Both boundary conditions are satisfied on the exact free-surface.

Impermeable condition is applied on the rigid bottom boundary (S_b) and sides and downstream end walls boundaries (S_sw, S_d) so that,

where, is the normal vector directed outward of the fluid. Inflow boundary condition can be written as (Tang and Huang, 2008):

where, ϕ_I is the theoretical input wave potential. For instantaneous free-surface, the initial conditions are defined as:

The direct boundary integral equation based on the Green's second identity is used to solve the boundary value problem in the Eulerian frame (Brebbia and Dominguez, 1992):

where, c(q) equals zero when the point is outside of the fluid domain and for a point on and inside the boundary domain (Φ ∈ S_f ∪ S_b∪ S_d∪ S_sw), it is solid angle i.e., 2π for a point on the smooth boundary and 4π for a point inside the boundary. For three-dimensional problems, G(q, p) is given as (Brebbia and Dominguez, 1992):

where,

2.2 Desingularized NURBS boundary integral equation

An arbitrary nodal point on the curved free surface can be described by the NURBS as (Piegl and Tiller, 1996)

where,

where, u_i is the knot given by Piegl and Tiller (1996). (v) is evaluated analogously. The unit normal and tangential vectors and for an arbitrary point on the surface can be obtained as follows:

where,

Location of a nodal (field) point p on the surface is expressed as (Piegl and Tiller, 1996)

where,

While source point

Equation 15 can be rearranged on the boundaries as (Cao et al. 1991):

where, Φ_N represents the boundaries in which the normal flux of potential is known and Φ_D indicates boundary in which the known potential is applied, whereas Φ ∈ Φ_N ∪ Φ_D in which Φ_NS_b ∪ S_d ∪ S_sw ∪ S_w with ϕ = φ and Φ_D ∈ S_f with ∂ϕ / ∂n = ψ.

Since

where, L_d and u are constant and equal to 1.0 and 0.5, respectively, and D_m is the square root of the panel area.

The boundary integral is computed by Gaussian quadrature scheme. Gaussian points are distributed over (u, v) plane ([0,1] ×[0,1]) and field points _p are located on each boundary surface.

Thus, the location of the source points cloud can be obtained as:

where, is Gaussian integration weight factor corresponding to Gaussian points and Λ(u, v) = |_u ×_v|. First term of Equation 16, Φ_N ∈ S_b, can be written as:

Similarly, other Dirichlet boundaries can be evaluated. Suppose that the bottom surface is divided into M × N (= M_b) segments in u and v directions, respectively. Then, Equation 19 can be written as:

where, u_i and v_j(i= 1, ..., M, j = 1, ..., N) are the Gaussian points located on the bottom surface. For the free surface boundary (Φ_D ∈ S_f), the second integral on the LHS of Equation 16 is discretized as:

where, and are the number of segments in u and v directions and d_n is the normal derivative of potential on the free surface. Other integrals for the remaining boundaries can be discretized similarly.

2.3 Time marching scheme

Material node approach is used to transform Equation 2 from the Eulerian form to the Lagrangian form. Nodal points are allowed to freely move with the particles of free surface. Then, Equation 2 can be written as (Tanizawa, 2000):

where, is the location of particles on the free-surface. The position of the particles is updated by the evaluation of particle velocity on the free-surface and Fourth-order Runge-Kutta time integration approach.

2.4 Tangential derivatives

Since, geometry of the free-surface is described by NURBS, the normal and tangential unit vector components are provided by Equations 10, 11 and 12. The potential variation on the free-surface can also be described by NURBS as:

where, Τ_i,j is potential on the control points and evaluated by potential of nodal points which do not participate in the basis function. Directional derivatives of ϕ with respect to u and v can be written as (Piegl and Tiller, 1996):

where, and are derivatives of B-spline basis functions. Derivatives of free-surface boundary value in tangential directions are computed as follows:

The particle velocity represented on the boundary can be represented by in which ϕ_n is the normal derivative of potential.

2.5 Artificial wave generator

Extreme waves described by linear and second-order Stokes wave theories are imposed on the inflow boundary. These are written as follows:

and also,

where, and correspond to first-order potential and wave profile, respectively, and and correspond to the second-order Stokes wave properties (Ning et al., 2009). A ramping function is engaged to avoid the impulse-like behavior and keep the stability of solutions and to reach to the steady state properly. In the present modeling, the ramping function given by Tang and Huang (2008) is used:

where, T_m is the modulation time which equals to the incident wave period in the present study.

2.6 Artificial damping zones

To obtain proper numerical solution in a numerical wave tanks, the artificial damping zone (sponge layer) is adopted at the end of wave tank. The energy dissipation scheme is used including adding an artificial damping term to the fully non-linear free-surface boundary condition over the region of the free-surface adjacent to the rigid walls boundaries and the inflow boundary. Modified non-linear free-surface boundary condition with damping coefficient presented by Cointe (1991) is as follows:

where, the subscript e corresponds to the reference configuration for the fluid. The function ν (x) is damping coefficient defined by

in which, ω is a characteristic wave frequency and k is a characteristic wave number. The parameters α and β control the strength and extent of the damping zone, respectively. x₀ and x₁ indicate the edges of the damping zones in horizontal plane on the free-surface. The terms ϕ_e and e are reference values. When reference values are set for calm water condition (ϕ_e = 0, η= 0), damping zone acts as simple absorber. If a propagating wave is used as reference value, then the damping zone allows only this wave to pass through. In practice, the damping coefficient equals zero except in the damping zone (x₀ < x < x₁), which is continuous and continuously differentiable.

3 NUMERICAL RESULTS

3.1 Convergence test and model verification

To perform the numerical simulation, a numerical wave tank of Ning and Teng (2006) is provided. The dimension of the wave tank is taken as 15.45 m length, 0.3 m width and 0.8 m depth. Linear wave with the wavelength of λ = 5.15 m and wave slope of kA =0.033 is specified for inflow boundary value. An artificial damping zone is extended along a wave length foreside of end wall tank. A numerical wave probe is adopted at (0.28L,0.5B) in horizontal plane. The solution of proposed NURBS NWT for 4×24 nodal points in u and v directions respectively, time step T/ 40 and B-spline basis function order κ × l = 3×3 are compared with analytical and the available numerical calculations in Figure 4.

When the wave is fully developed, analytical solution and numerical computations with different number of nodal points are tabulated on Table 1. Root mean square of wave elevation for a complete wave period at the raised numerical probe is defined as:

where, K is the number of the time steps to complete a wave period with time step T /40. NURBS NWT is run with a desktop PC (Intel Core 2 Quad CPU, 2.66 GHz, and 2 GB of RAM). It is shown that to achieve about half of RMS corresponds to 4×20 nodal points, more than twice CPU time is taken for 4×48 nodal points at each time step.

Accuracy of the present model for different wave slopes and wave frequencies of second-order stokes wave is examined and given in Table 2. Hence, three wave slopes are chosen and RMS of wave elevations is determined for 4×30 nodal points on free surface and time step T /40 . The order of B-spline basis function is chosen κ × l = 3×3. It is shown that increasing wave slope decreases the accuracy of solution. For kA =0.104 , the wave slope is greater than kA =0.033 about three time but, increasing of RMS is infinitesimal.

Effect of time step on the numerical solution is given in Figure 5. In this Figure the solutions of input second-order stokes wave with kA =0.104 and different time step with water depth 0.5 is compared. Mesh size and order of B-spline basis function is the same as for Table 2. For time step T/ 90 , fully development of the input wave is postponed due to the coarse time steps. For t/T > 10 , the wave is fully developed for three time steps and the numerical solutions are independent from the time step.

To obtain the proper solution of NWTs, the damping zones should be adopted. A part of wave energy will return back to the computational domain from downstream boundary if the damping of wave energy is too weak. On the other hand, if the absorbing strength is too powerful, the damping zone will act as a solid boundary and the waves reflect from outflow boundary. The length and strength of damping zone control the performance of damping zone. Therefore, for different parameters α and β the wave elevation of free surface is compared in Figure 6 and Figure 7, respectively.

In Figure 6, the free surface oscillations along the numerical tank (x/d) due to the input second-order stokes wave for β = 2 and different α are shown. When the strength of the damping zone is increased the damping zone acts as the rigid walls and some portion of incident wave is reflected. For α = 3 , increment of the wave height is happen due to wave reflection. For α = 1 and α = 1.5 , the wave reflection is reduced from the wave absorber.

For different length of the damping zone, the wave elevation for α = 1 is given along the numerical wave tank in Figure 7. For β = 1 , the wave absorber do not dissipate the incident wave and for long simulation the wave reflection is happened severely. For β = 2 and β = 3 , the incident wave is fully dissipated within the wave absorber and open water condition is kept.

3.2 Focused wave generation and propagation

In the experimental work, data measurement of non-linear waves is transformed in the Fourier domain as presented by Ning et al. (2009) and wave spectrum S(f) is manifested. Indeed, it represents wave components with different amplitudes and periods which interact with each other and consequently make extreme wave. Amplitude of each wave component α_i is obtained as:

where, A is the general linear focused wave amplitude given by Ning et al. (2009), S_i(f) is the spectral density of each wave component i and Δf is the frequency step depending upon the number of wave components .

For a wave flume with water depth of d = 0.5 m, spectrums are presented in Figure 8. Characters of each wave spectrum are depicted in Table 3 based on Westphalen et al. (2012) in order to simulate the wave in which T_p is the peak period of wave spectrum and λ_p is the wave length corresponding to the peak period.

In this paper, an artificial wave generator is located at x = 0 and an artificial damping zone with the length twice as great as the length of the wave is placed at the end of the tank. To verify proposed numerical procedure, a tank with the dimensions of 12 m ×1 m × 0.5 m (L × B × d) is provided and the focused point and time are set to x₀ = 3 m and t₀ = 9.2 s for case 2, x₀ = 3.27 m and t₀ = 10 s for case 3 and case 4, respectively. The phase angle ε_i is taken equal to zero. 4×36 Gaussian points are adopted on the free-surface and the rational B-spline basis function of order 3×3 is used to simulate the free-surface elevation. Other boundaries are flat and described by isoperimetric linear quadrilateral elements. Time step Δt = T_p / 40 is used in high-order time integration scheme.

The maximum linear and nonlinear wave elevations η^max for the cases given in Table 3 are computed and compared with the physical measurements and numerical results of Ning et al. (2009) and Westphalen et al. (2012). The comparisons are presented in Table 4. Meanwhile, numerical results include a three-dimensional potential numerical wave tank computations based on the boundary integral equation with isoparametric quadratic element developed by Ning et al. (2009). In addition, the numerical wave tank calculations based on commercial CFD packages with both FV and CVFE solvers were conducted by Westphalen et al. (2012).

For case 2, the highest crest of the first-order wave input with NURBS NWT is higher than the evaluations with CV-FE and FV solver. For the second-order wave components, NURBS NWT predicts the maximum wave elevation closer to the experimental measurement than the FV solver.

For case 3, the trend is similar. For case 4, the numerical results show good agreements with the experimental data. It should be mentioned that in this case, the focused wave almost broke in the real wave tank and the nonlinearity predominates on the simulation.

Time series of wave elevations of the cases 2-4 are depicted at focused point in Figure 9 and Figure 10, respectively. The comparison of computed crest and trough focused waves using linear and nonlinear theory with experimental results is presented. For case2, computational wave crest on the center of wave group is reached to the experimental measurement by NURBS NWT. When the input wave includes the second-order wave components, prediction of the central wave group trough coincides with physical wave tank measurement as shown in Figure 10. The results of the surrounding crests and troughs are generally improved from the first-order to the second-order. For both input wave cases, the surrounding trough elevations are slightly higher than the physical experiments with respect to the surrounding crests and somewhat, NURBS NWT decreases the small differences. For case 3, it can be said that the numerical evaluations are closer to the physical experiments better than case 2 and the slight difference of the surrounding troughs and crests are decreased. It shows that moving from input linear wave components to the second-order wave components improves trough elevations.

In the physical measurement, wave breaking almost occurs due to the steepness of the wave in case 4. Hence, severe nonlinearity is attended in the wave simulation. Nevertheless, computational results for central crest reasonably agree with the experimental measurements for both input wave cases. There are substantial differences in surrounding crests and troughs and the wave trends found by Ning et al. (2009). Capability of NURBS to predict sharp mutation of the free surface accommodates computations to agree with the physical measurements. As well, symmetry around the maximum crest of wave group in case 2 is retained by NURB NWT and for steeper wave is lost while the nonlinearity dominates in the numerical modeling.

Experimental investigation on propagation of incident irregular waves over a submerged bar is carried out by Beji and Battjes (1993) and its numerical simulations based on Boussinesq equations and mild-slope equations are conducted by Hsu et al. (2007) and (2002), respectively. JONSWAP spectrum of Hsu et al. (2007) given in Equation 35 is chosen to pass over a submerged bar. Dimensions of numerical wave tank is equal to the experimental wave flume of Beji and Battjes (1993) with the length of 37.7 m, the breadth of 0.8 mand the water depth of 0.4 m.

Spectrum of the fully propagated incident wave with significant wave height of H_1/3 =0.03 m and the significant period of T_s = 2.5 s is compared with the theoretical result and numerical computation in Figure 11.

where, the peak period T_p is related to the significant period, and the coefficients are ℘ = 3., σ₁ = 0.2189 , σ₀ = 0.07 for f < f_p and σ₀ = 0.09 for f > f_p.

It seems that there are trivial mismatch at the peak spectrum in Hsu et al. (2007) computations which is improved by the present model. Furthermore, the deviation between theoretical and numerical results is remained on the higher frequency domain, but it seems that NURBS NWT model is closer to the proposed spectrum. To measure the irregular wave evolution due to a submerged trapezoidal bar, eight wave probes are arranged along the wave flume as shown in Figure 12.

It is worth mentioning that the wave breaking does not occur during the simulation and the observation is conducted at the first wave probe. Comparison of the experiments and the numerical results is given in Figure 13. At the peak point of energy density at every wave probes, NURBS NWT is better than the other numerical computations, whereas the mild slop simulation for wave probes 5-8 and the Boussinesq simulation for wave probe 2 and 4 are overestimated. In addition, these simulations do not touch peak measurement at wave probes 1 and 3. For higher frequency, substantial deviation occurs at wave probes 5-8 for numerical simulation showing that computations are partly improved by NURBS NWT.

4 CONCLUSIONS

Development of a fully non-linear 3D NWT is considered in this paper for investigation of propagation and scattering of non-linear random sea wave due to bottom submerged bars. The simulation of fully non-linear waves using NURBS in a potential three-dimensional numerical wave tank is successfully completed.

MEL method and desingularized boundary element method is employed for numerical simulation. To keep numerical accuracy and avoid instability in MEL, five points Chebyshev smoothing scheme was adopted in time marching. To obtain appropriate numerical solutions for wave propagation problem in a numerical wave tank artificial damping zones (sponger layer) is adopted. Perturbation sources are placed on the fixed inflow boundary to make the free-surface oscillating during the simulation. Also, the NURBS is used to evaluate velocity of the free surface particles accurately. It is a novel procedure to calculate the free-surface kinematics.

The stability and accuracy of NURBS NWT were examined and verified to model the free-surface. It is shown that the present approach gives accurate solution, whereas the computation time was saved and the computational nodes were reduced. Simulations of propagation of irregular wave spectrum and wave transformation of a random sea waves due to submerged bar are presented and compared with the experimental and available numerical results. It is shown that the appropriate solution with sufficient accuracy can be obtained by the proposed numerical procedure.

Received 05.08.2013

In revised form 29.07.2014

Accepted 03.08.2014

Available online 16.08.2014

Publication Dates

History