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

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). Highorder 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.

Simulation of irregular waves over submerged obstacle on a NURBS potential numerical wave tank 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 func-Latin American Journal of Solids and Structures 11 (2014) 2308-2332 tion 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. (1990Cabral et al. ( , 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).
Latin American Journal of Solids and Structures 11 (2014) 2308-2332 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.

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 adja- cent of the end tank wall.

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 f can be introduced.Velocity potential satisfies the Laplace equation in the fluid domain W : There are fully non-linear boundary conditions, Kinematic Free Surface Boundary Condition (KFSBC) and Dynamic Free Surface Boundary Condition (DFSBC) on the free-surface ( f S ) which can be written as the following (Tanizawa, 2000): where, ( ) where, n  is the normal vector directed outward of the fluid.Inflow boundary condition can be written as (Tang and Huang, 2008): where, I f 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): Latin American Journal of Solids and Structures 11 (2014) 2308-2332 where, ( ) c q equals zero when the point is outside of the fluid domain and for a point on and inside the boundary domain ( ), it is solid angle i.e., 2p for a point on the smooth boundary and 4p for a point inside the boundary.For three-dimensional problems, ( ) , G q p is given as (Brebbia and Dominguez, 1992): where, q c  and p c  are the source and field points location, respectively.

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, ,  (Piegl and Tiller, 1996): where, i u is the knot given by Piegl and Tiller (1996).
 and tangential vectors s  and t  for an arbitrary point on the surface can be obtained as follows: Latin American Journal of Solids and Structures 11 (2014) 2308-2332   ( where, u are partial derivatives of NURBS surface in u and v directions shown in Figure 2. Location of a nodal (field) point p on the surface is expressed as (Piegl and Tiller, 1996) ( ) ( ) While source point q c  is hold outside of W and p c  on the boundary integration surface.Equation 6can be written as: Solids and Structures 11 (2014) 2308-2332   ( Equation 15 can be rearranged on the boundaries as (Cao et al. 1991): where, N G represents the boundaries in which the normal flux of potential is known and D G indi- cates boundary in which the known potential is applied, whereas Since q c  and p c  are not coincident, integrand singularity is removed.To find the location of q c  , the desingularization distance is used which is a function of local grid sizes on the boundary surface as shown in Figure 3. Distance of a source point from the boundary has been proposed as (Cao et al. 1991): ( ) ( ) Solids and Structures 11 (2014) 2308-2332 where, d L and u are constant and equal to 1.0 and 0.5, respectively, and m D is the square root of the panel area.
The boundary integral is computed by Gaussian quadrature scheme.Gaussian points are distributed over ( ) ) and field points p c  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 ( ) . First term of Equation 16, N b S G Î , can be written as: Similarly, other Dirichlet boundaries can be evaluated.Suppose that the bottom surface is divided into ( ) b M ´N = M segments in u and v directions, respectively.Then, Equation 19 can be written as:  ,  , , , Latin American Journal of Solids and Structures 11 (2014) 2308-2332 where, M and N are the number of segments in u and v directions and n d is the normal deriva- tive of potential on the free surface.Other integrals for the remaining boundaries can be discretized similarly.

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, c  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.

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 Q is potential on the control points and evaluated by potential of nodal points which do not participate in the basis function.Directional derivatives of f with respect to u and v can be written as (Piegl and Tiller, 1996): where, ( ) The particle velocity represented on the boundary can be represented by in which n f is the normal derivative of potential.

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, ( ) h correspond to first-order potential and wave profile, respectively, and h 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, m T is the modulation time which equals to the incident wave period in the present study.

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: Latin American Journal of Solids and Structures 11 (2014) 2308-2332 where, the subscript e corresponds to the reference configuration for the fluid.The function ( ) in which, w is a characteristic wave frequency and k is a characteristic wave number.The param- eters a and b control the strength and extent of the damping zone, respectively.

, 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 ( ) , which is continuous and continuously differentiable.

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.45When 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 40 T .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 ´  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 a and b 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 secondorder stokes wave for 2 b = and different a 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 a = , increment of the wave height is happen due to wave reflection.For 1 a = and 1.5 a = , the wave reflection is reduced from the wave absorber.For different length of the damping zone, the wave elevation for 1 a = is given along the numerical wave tank in Figure 7.For 1 b = , the wave absorber do not dissipate the incident wave and for long simulation the wave reflection is happened severely.For 2 b = and 3 b = , the incident wave is fully dissipated within the wave absorber and open water condition is kept.

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 repre- sents wave components with different amplitudes and periods which interact with each other and consequently make extreme wave.Amplitude of each wave component i a is obtained as: where, A is the general linear focused wave amplitude given by Ning et al. (2009) ´ 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 is used in high-order time integration scheme.The maximum linear and nonlinear wave elevations max.h 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 CV-FE 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 m and the water depth of 0.4 m .Spectrum of the fully propagated incident wave with significant wave height of 1 3 0.03 H m = and the significant period of where, the peak period p T is related to the significant period, and the coefficients are  merical 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.

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 freesurface.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

Figure 1 :
Figure 1: A definition sketch and coordinates system.
of Solids and Structures 11 (2014) 2308-2332 the surface, and , m n are the numbers of control points in u and v directions, respectively.k and l are the orders of the B-spline basis functions

Figure 2 :
Figure 2: Partial derivatives and unit normal vector.

Figure 3 :
Figure 3: Distribution of source points and collocation points.Ox are the Gaussian points located on the bottom sur- face.For the free surface boundary ( D f S G Î ), the second integral on the LHS of Equation16is discretized as: derivatives of B-spline basis functions.Derivatives of free-surface bound- ary value in tangential directions are computed as follows: Latin American Journal of Solids and Structures 11 (2014) 2308-2332 0 x and 1 x indi- cate the edges of the damping zones in horizontal plane on the free-surface.The terms e f and e c  are reference values.When reference values are set for calm water condition (

Figure 4 :
Figure 4: Comparison of analytical and numerical modeling of linear wave, A h vs. dimensionless time.
increasing of RMS is infinitesimal.Effect of time step on the numerical solution is given in Figure5.In this Figure the solutions of input second-order stokes wave with 0.104 kA = 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 Table2.For time step90 T, fully development of the input wave is postponed due to the coarse time steps.For 10 t T  , the wave is fully developed for three time steps and the numerical solutions are independent from the time step.
of Solids and Structures 11 (2014) 2308-2332 , ( ) i S f is the spectral density of each wave component i and f D is the frequency step depending upon the number of wave components n .For a wave flume with water depth of presented in Figure8.Characters of each wave spectrum are depicted in Table3based on Westphalen et al. (2012) in order to simulate the wave in which p T is the peak period of wave spectrum and p l is the wave length corresponding to the peak period.In this paper, an artificial wave generator is located at 0 x = 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 provided and the focused point and time are set to 0 case 3 and case 4, respectively.The phase angle i e is taken equal to zero. 4 36

Figure 5 :
Figure 5: Comparison of numerical modeling of linear wave for different time steps, A h vs. dimensionless time.

Figure 6 :
Figure 6: Performance of the end wall damping zone for different strength of the damping zone.

Figure 7 :
Figure 7: Performance of the end wall damping zone for different length of the damping zone.

Figure 8 :
Figure 8: Wave spectrum of focused wave and properties of simulated cases.

Figure 9 :
Figure 9: Comparison of time history of wave elevation of the first-order focused wave cases 2-4.

Figure 10 :
Figure 10: Comparison of time history of wave elevation of the second-order focused wave cases 2-4.

Figure 11 :
Figure 11: Comparison of JONSWAP spectrum recorded in physical and numerical wave tank.It seems that there are trivial mismatch at the peak spectrum inHsu et al. (2007) computations which is improved by the present model.Furthermore, the deviation between theoretical and nu-

Figure 10 :
Figure 10: Comparison of irregular wave transformation when propagating over the submerged bar.

Table 1 :
RMS error and CPU time of linear wave simulation.

Table 2 :
RMS error and CPU time of linear wave simulation.

Table 3 :
Properties of wave cases for simulation.

Table 4 :
Maximum wave elevation of three focused wave.