Modeling nonlinear waves in a numerical wave tank with localized meshless RBF method

Computers & Fluids 44 (2011) 221–228

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Modeling nonlinear waves in a numerical wave tank with localized meshless RBF method Utku Senturk ⇑ Mechanical Engineering Department, Faculty of Engineering, Ege University, 35100 Bornova, Izmir, Turkey

a r t i c l e

i n f o

Article history: Received 27 August 2009 Received in revised form 10 December 2010 Accepted 5 January 2011 Available online 9 January 2011 Keywords: Meshless Radial basis functions Local Nonlinear Wave Simulation

a b s t r a c t This paper presents the numerical simulation of free surface waves in a 2D domain which represents a wave tank, using a localized approach of the meshless radial basis function (RBF) method. Instead of global collocation, the local approach breaks down the problem domain into subdomains, leading to a sparse global system matrix which is particularly advantageous in tackling the time consuming simulation process. Mixed Eulerian–Lagrangian approach is adopted for free surface tracking and fourth order Adams– Bashforth–Moulton scheme (ABM4) is used for time stepping. Both linear and nonlinear Stokes waves are simulated and compared to analytical solutions. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Many problems in the field of engineering and science are introduced by partial differential equations. Extensive amount of research was handled in order to develop approximate solution methods in the past decades, such as FDM [1], FEM [2], BEM [3] and FVM [4], some reached a mature stage and commercialized. Lately, meshless methods which are based on the use of scattered nodal points rather than polygonisation, have undergone intense study as an alternative method. Amongst a variety of meshless methods, using multiquadrics (MQ) for radial basis function (RBF) networks, as originally presented by Kansa [5,6], offers a continuously differentiable and integrable grid free scheme for representing surfaces and partial derivative estimation where the function that is to be approximated is decomposed into a set of RBFs and the derivatives can be calculated directly by differentiating. Alternatively, the indirect RBF method was introduced by MaiDuy et al. [7] which is based on the integration of decomposed derivatives to approximate the original function in a least squares sense. In this manner, several test problems, including Navier– Stokes benchmarks, yielding a superior accuracy to the direct approach were studied with IRBF [8–10]. Despite the accuracy, global meshless methods result in a fully populated system matrix that is computationally expensive to deal and this is the major drawback. To tackle this problem, a localized scheme was suggested [11,12], ⇑ Tel./fax: +90 232 388 85 62. E-mail address: [email protected] 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.01.004

where the problem domain is broken into subdomains leading to a sparse global system matrix. Several studies were presented concerning a range of problems from diffusion to thermo-fluid problems, including implementation to Navier–Stokes equations [13– 16]. Computation of free surface flows is a well studied area where much research has been being carried out using various numerical methods. A major group of studies deals with the wave propagation simulation problem as well as wave–structure interactions [17–19] using boundary element methods. Grilli et al. [20,21], investigated the numerical wave tank problem using high order boundary element method. Zhang et al. [22], studied the numerical wave tank problem using desingularized boundary integral equation method. Additionally, several researches [23–26] gather around the background of VOF method, introduced by Hirt and Nichols [27]. Wu et al. [28–30], and Xiao et al. [31] presented a meshless numerical model for nonlinear free surface water waves by collocation of boundary points and using the fundamental solution of Laplace equation as the RBF. Notably, most of these studies have the Mixed Eulerian–Lagrangian approach for free surface tracking that is introduced by the pioneer work of Longuet-Higgins and Cokelet [32]. A more thorough review on theory and technology of numerical wave tanks is presented by Tanizawa [33]. The numerical simulation of waves in a 2D numerical wave tank subjected to nonlinear boundary conditions on the free surface, using localized meshless RBF method is the scope of this study. The paper is organized as follows. In Section 2, the numerical wave tank as a boundary value problem is formulated and the method-

222

U. Senturk / Computers & Fluids 44 (2011) 221–228

ology is introduced. Section 3 gives detailed information about the theory and the strategy of implementing localized meshless RBF method. Numerical results and analytical solutions are presented in Section 4. 2. Numerical wave tank simulations

On the upstream boundary (C1), the fluid motion is generated either by a prescribed wave maker motion, or by specifying the wave properties according to a chosen wave theory at the boundary. Because of its simplicity, latter is implemented. Therefore, C1 is a Dirichlet type boundary, which specifies the value of the unknown potential function (/) that will be given analytically in the subsequent section for each problem that is under consideration.

2.1. Boundary value problem 2.2. Analytical solutions The fluid considered in the problem is assumed inviscid, incompressible and initially at rest. Using Navier–Stokes equations without viscous terms and Kelvin’s theorem of the conservation of circulation will yield,

In the small amplitude wave theory, the linearized solution for the velocity potential and the free surface elevation can be derived as

r~ V ¼0

/ðx; z; tÞ ¼

ð1Þ

which states that the motion of such an ideal fluid is irrotational [34]. For an irrotational flow, the entire flow field can be described by ~ V ¼ r/, evaluating into continuity gives the Laplace equation,

r2 / ¼ 0

ð2Þ

where /(x, z, t) represents the velocity potential in 2D space. Eq. (2) is the governing equation for the water wave boundary value problem which is elliptic and boundary conditions must be applied over the entire boundary [35]. Rectangular problem domain (x, z) e [0, 30]  [1, 0] and boundary conditions are given in Fig. 1. On the free surface boundary (C4), fully nonlinear, kinematic and dynamic free surface boundary conditions are,

@ g @/ @ g @/  on z ¼ gðx; tÞ ¼ @z @x @x @t "   2 # 2 @/ 1 @/ @/ on z ¼ gðx; tÞ ¼ g g  þ @t 2 @x @z

gðx; tÞ ¼

H g cosh kðh þ zÞ sinðkx  xtÞ 2 x coshðkhÞ

H cosðkx  xtÞ 2

ð7Þ ð8Þ

where the wave properties H, x, k are wave height, circular frequency and wave number, respectively. On the other hand, the perturbation approach of Stokes and the governing equation with nonlinear boundary conditions lead to Stokes waves [37]. The potential function and the free surface elevation function for Stokes’ 2nd order wave theory are defined by,

/ðx; z; tÞ ¼

H g cosh kðh þ zÞ sinðkx  xtÞ 2 x coshðkhÞ þ

ð3aÞ

gðx; tÞ ¼ ð3bÞ

3 2 cosh 2kðh þ zÞ sin 2ðkx  xtÞ H x 4 32 sinh kh

H H2 k cosh kh ð2 þ cosh 2khÞ cosðkx  xtÞ þ 2 16 sinh3 kh  cos 2ðkx  xtÞ

ð9Þ

ð10Þ

where g (x,t) is the displacement of the free surface about the horizontal axis, z = 0. In the small amplitude wave theory, these boundary conditions simplify into their linearized form,

In addition, Stokes’ 3rd order solutions for a given wave parameter (a) are defined by Ma and Li [38],

@ g @/ ¼ @z @t

/ðx; z; tÞ ¼

on z ¼ 0

ð4aÞ

@/ ¼ g g on z ¼ 0 @t

ð4bÞ

gðx; z; tÞ ¼

ð11Þ

3 X

ei cos iðkx  xtÞ

ð12Þ

where the coefficients are, 2

E1 ¼

ð5Þ

where C denotes phase velocity of the waves. The bottom boundary (C3) is assumed to be horizontal, so a no-flow condition applies,

@/ ¼0 @z

Ei cosh ikðh þ zÞ sin iðkx  xtÞ

i¼1

i¼1

On the radiation boundary (C2), outgoing waves are allowed to leave freely without reflection [36]. Thus, the Sommerfeld–Orlanski radiation boundary condition reads

@/ @/ þC ¼0 @t @x

3 X

2

2

aC ka C cosh khð1 þ 5cosh khÞ  5 sinh kh 8sinh kh

ð13Þ

2

E2 ¼

3ka C

ð14Þ

4

8sinh kh 2

ð6Þ

E3 ¼

k a3 Cð11  2 cosh 2khÞ 7

64sinh kh

e1 ¼ a

ð15Þ ð16Þ

2

e2 ¼

ka cosh khð2 þ 2 cosh 2khÞ 3

4sinh kh 2

e3 ¼

ð17Þ

3

3k a3 ð1 þ 8cosh khÞ 6

64sinh kh

ð18Þ

For the linear wave theory and the 2nd order wave theory, the wave length (k) is found using,

k¼ Fig. 1. Sketch of the problem domain and boundaries on the xz plane.

gT 2 tanh kh 2p

ð19Þ

223

U. Senturk / Computers & Fluids 44 (2011) 221–228

On the other hand, the wave length and the wave height for 3rd order wave theory are [38],

k¼

" # 2 2 gT 2 k a2 ð7 þ 2cosh 2khÞ tanh kh 1 þ 4 2p 8sinh kh "

H ¼ 2a 1 þ

# 2 3 3k a2 ð1 þ 8cosh khÞ 6

64sinh kh

ð20Þ

ð21Þ

The wave number is,

/ðt þ DtÞ ¼ /ðtÞ þ

tþDt

f2 dt

ð28Þ

t

These integrations are calculated numerically using the fourth order Adams–Bashforth–Moulton scheme (ABM4). This is a two step method that requires solving the algebraic equation system that is constructed using the meshless method, twice at each time step. At the predictor step or the Adams–Bashforth scheme, following equations are used:

gðt þ DtÞ ¼ gðtÞ þ

2p k¼ k

Z

ð22Þ /ðt þ DtÞ ¼ /ðtÞ þ



 Dt ð55f1t  59f1tDt þ 37f1t2Dt  9f1t3Dt Þ 24



for all calculations.

 Dt ð55f2t  59f2tDt þ 37f2t2Dt  9f2t3Dt Þ 24

ð29Þ

ð30Þ

The corrector step or the Adams–Moulton scheme involves, 2.3. Mixed Eulerian–Lagrangian (MEL) approach

gðt þ DtÞ ¼ gðtÞ þ Initiated by Longuet-Higgins and Cokelet [32], the MEL approach enabled the development of several studies on nonlinear wave propagation modeling in numerical wave tanks as well as wave body interaction simulations. In this method, the boundary value problem is solved using the Eulerian form of equations. Obtained solutions are used to track the free surface nodes using the Lagrangian form of the free surface boundary conditions with an appropriate time stepping integration scheme. Velocity potential values on the free surface are also found using same integration procedure so that these will appear as a Dirichlet type boundary condition in the next time step. The material derivative is utilized to define the time variations for the Lagrangian frame of reference, given by,

D @ dx @ dz @ ¼ þ þ ¼~ Vp  r Dt @t dt @x dt @z

ð23Þ

where ~ V p is the velocity of the moving nodes on the free surface. MEL formulation can be implemented in two ways, namely, the material node approach and the Semi-Lagrangian approach. In this study, the Semi-Lagrangian approach is adopted where the nodes on the free surface boundary are allowed to move vertically only. Thus, nodal velocity on the free surface becomes [18,22],

/ðt þ DtÞ ¼ /ðtÞ þ



 Dt ð9f1tþDt þ 19f1t  5f1tDt þ f1t2Dt Þ 24



 Dt ð9f2tþDt þ 19f2t  5f2tDt þ 9f2t2Dt Þ 24

2.5. Sponge layer and Sommerfeld–Orlanski condition Time domain simulation of the numerical wave tank needs an appropriate modification at the downstream boundary to ensure that outgoing waves are not reflected and do not influence the domain. In this study, a sponge layer scheme, alternatively known as the artificial beach or the damping zone is used along with the Sommerfeld–Orlanski radiation condition. In the prescribed section of the domain, the waves gradually vanish by adding extra terms to the free surface boundary conditions such as [22,33],

f1 ¼

Dg @/ @/ @ g  ¼  v ðxÞg @z @x @x Dt

Using this definition, the Lagrangian form of the dynamic and the kinematic free surface boundary conditions are as follows:

f2 ¼

D/ Dt

Dg @/ @/ @ g  ¼ @z @x @x Dt "   2 #   2 D/ 1 @/ @/ @/ @/ @/ @ g þ ¼ g g   þ f2 ¼ Dt 2 @x @x @z @z @x @x

f1 ¼

ð25Þ

ð26Þ

2.4. Time stepping integration Variations in free surface node positions and the velocity potential values of fully nonlinear wave propagation simulations can be calculated treating Eqs. (25) and (26) as first order ordinary differential equations and integrating,

Z t

tþDt

f1 dt

1 ¼ g g  2

"   2 #   2 @/ @/ @/ @/ @/ @ g þ  þ  v ðxÞ/ @x @z @z @z @x @x

ð33Þ

ð34Þ

where m is the damping coefficient that is calculated using,

Using meshless methods for such moving boundary problems is particularly advantageous since only nodal positions need to be updated rather than remeshing the entire domain at each time step during the solution process.

gðt þ DtÞ ¼ gðtÞ þ

ð32Þ

where the values obtained from the first solution are directly used in the corrector step. After updating the required derivatives and defining the boundary values of the velocity potential gathered from the corrector step, the algebraic system of equations are solved again to arrive at the final solution for the current time step. Calculated nodal values of the free surface elevation function will be used to define the free surface profile of the next time step.

ð24Þ

Dg ~ ~ Vp ¼ k Dt

ð31Þ

ð27Þ

VðxÞ ¼

8 <

as x

:

0

h

xðLbkÞ bk

i2

x  ðL  bkÞ

ð35Þ

x < ðL  bkÞ

where b and as are the sponge layer length factor and the tuning factor, respectively. These parameters can be used to adjust the damping effect and their values are given in the results section. The implementation of Sommerfeld–Orlanski type radiation boundary condition which is introduced by Eq. (5) is handled by replacing time derivative term with a forward finite differencing scheme as,

/tþDt ¼ /t  DtC



t @/ @x

ð36Þ

Though the phase velocity is defined by k/T for steady waves, it’s quite different here because of the large transient effects since water is the initially at rest. To reduce possible numerical errors,

224

U. Senturk / Computers & Fluids 44 (2011) 221–228

the dynamic free surface boundary condition given in Eq. (3b) is used to calculate phase velocity at each time step by [33],

h  2 @/ 1

C¼

gn þ 2

@x

@/

þ

@/2 i @z

ð37Þ

This procedure is applied to several nodes on the free surface close to the radiation boundary, and their arithmetic mean is assigned to C in Eq. (36) unless it’s negative. C is set to zero for negative arithmetic mean. 2.6. Ramp function Whenever the analytical solutions given in Eqs. (7)–(21) are used for wave generation, applying a ramp function to these solutions becomes crucial. This ramp function gradually increases to unity as the simulation proceeds, avoiding the numerical error that may be experienced at early steps since the simulation starts with a still water surface. For this purpose, the following function is used, 1 2

RampðtÞ ¼

1  cos

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx  xi Þ2 þ ðz  zi Þ2 þ c2

Ri ðxÞ ¼

@x

( h

Multiquadric RBF is given by,

 i pt

Tm

1

t  Tm

ð38Þ

t > Tm

where Tm controls the time interval for the ramping. Tm is set to 2  T throughout the study.

ð44Þ

for 2D geometry where c is the free shape parameter. Polynomials can be generated for desired m values using Pascal triangle of monomials such as,

PT ðxÞ ¼ f1 x z x2 xz z2 . . .g |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð45Þ

m number of monomials

The idea behind building this augmented structure is to avoid the singularity problem that is experienced in the shape function construction [12]. In Eq. (42), RBF moment matrix (R0) is constructed using multiquadric RBFs,

2

3 Rn ðx1 Þ .. 7 7 . 5 R1 ðxn Þ    Rn ðxn Þ

R1 ðx1 Þ 6 . R0 ¼ 6 4 ..

 .. .

ð46Þ

and polynomial moment matrix (Pm) is,

2 6 6 6 6 T Pm ¼ 6 6 6 6 4

2.7. Smoothing scheme

1

1



x1

x2



z1 ...

z2 ...

 .. .

1

3

7 xn 7 7 7 zn 7 7 .. 7 . 7 5

ð47Þ

pm ðx1 Þ pm ðx2 Þ    pm ðxn Þ

As reported in several articles, a numerical instability in the calculated free surface profile is experienced with a saw-toothed appearance. In order to remove this instability, a smoothing procedure which is given by Longuet-Higgins and Cokelet is applied to the wave profile and the velocity potential values on the free surface such as,

In addition, since there are n + m unknowns in Eq. (42), the following m equations can be added as constraint conditions [12],

1 gj ¼ ðgj2 þ 4gj1 þ 10gj þ 4gjþ1  gjþ2 Þ 16

as,

/j ¼

1 ð/j2 þ 4/j1 þ 10/j þ 4/jþ1  /jþ2 Þ 16

ð39Þ ð40Þ

n X

pj ðxi Þai ¼ PTm a ¼ 0 j ¼ 1; 2; . . . ; m

ð48Þ

i¼1

Combining these equations, the system of equations are formed

 



a ~ s ¼ Us ¼ R0 Pm U T Pm 0 b 0 |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}

ð49Þ

G

These smoothing formulas are used every 10 time steps and the saw-toothed appearance is effectively removed.

where G is a symmetrical matrix which is constructed using symmetrical RBF moment matrix (R0). Thus, the vector of unknown coefficients is,

3. Numerical procedure

 a

In this study, Localized Meshless Method with radial basis functions (RBF) is implemented for simulating nonlinear water waves. Using this method, the approximation of a function u at any point x can be written as a linear combination of n multiquadric RBFs (Ri(x)) and m polynomials (pj(x)),

uðxÞ ¼

n X

Ri ðxÞai þ

i¼1

m X

pj ðxÞbj ¼ RT ðxÞa þ pT ðxÞb

ð41Þ

j¼1

where xT = [x, z] and a, b are unknown coefficients [12]. At this point, primary goal is to find ai and bj. In order to achieve this, a subdomain is formed around the point of interest at x including n number of neighbor points in this subdomain. Then ai and bj coefficients can be determined by enforcing Eq. (41) to be satisfied at these n nodes around the point of interest. This will lead to linear equation systems, one system for each node and each having n number of linear equations. So, in matrix form, we have,

Us ¼ R0 a þ Pm b

ð42Þ

where Us collects the function values for each subdomain such as,

Us ¼ fu1

u2

. . . un g

T

ð43Þ

b

~s ¼ G1 U

ð50Þ

Using this information, Eq. (41) can be re-written as,

es uðxÞ ¼ RT ðxÞa þ PT ðxÞb ¼ fRT ðxÞ PT ðxÞgG1 U

ð51Þ

Augmented shape functions are defined by,

f T ðxÞ ¼ fRT ðxÞ pT ðxÞgG1 SF ¼ fsf1

sf2

. . . sfn

sfnþ1 . . . sfnþm g

ð52Þ

So the problem has evolved into building shape functions rather than directly attempting to solve for ai and bj. Using nodal values, the value of the function u and its derivatives can be approximated rapidly e.g.,

~ ¼ SFT Us u

ð53Þ

@ @ ~ ¼ SFT Us u @x @x

ð54Þ

@2 @2 ~¼ u SFT Us @x2 @x2

ð55Þ

225

U. Senturk / Computers & Fluids 44 (2011) 221–228

z subdomains

x

Δz

z=-h Δx Fig. 2. Subdomains are allowed to overlap. From left to right, three subdomains are shown with a1 < a2 < a3.

where SF is the vector of shape functions whose entries are the first n entries of the augmented shape functions given in Eq. (52), i.e.,

SFT ¼ fsf1

sf2 . . . sfn g

ð56Þ

in which n is the number of nodes used in the subdomain [12]. The problem domain which is represented by uniformly (Dx, Dz) distributed, N number of nodes is broken into overlapping subdomains, leading to small systems of equations rather than a

Set RBF parameters (α,c) and problem constants (h, g, k, t, ω etc.) Build nodes, generate subdomains Generate shape functions and their derivatives

dense global system matrix. This is achieved by assigning circular zones centered on each node of the problem domain, each of those with a radius of,

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r s ¼ a Dx2 þ Dz 2

where a is a predefined parameter that controls the size of the subdomain. Thus, appropriately modifying a changes the number of points covered by the subdomain boundary, in this case, a circle (Fig. 2). Accounting into appropriate boundary conditions, given boundary value problem is reduced to a set of algebraic equations. For an internal node xT = [xI, zI], the governing equation (Eq. (2)) can be discretized by,

! @2 @2 T T U ¼ |{z} SF þ SF 0 @x2 @z2 TI |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð58Þ

KI

Predict

Similarly, for a Dirichlet boundary node, say xT = [xb, zb] with boundary value ub,

and

SFT Us ¼ ub |{z} |{z}

Assemble K, T matrices and solve for

Kb

Calculate derivatives of

and

Correct

and Calculate phase velocity C

ð59Þ

Tb

These small matrices are stacked together row by row to form the global system matrix K and global vector T so that discretized system equations are obtained as,

KðNxNÞ UðNx1Þ ¼ TðNx1Þ Assemble K, T matrices and solve for Update derivatives of

ð57Þ

ð60Þ

The global system matrix K, is now a sparse matrix, constructed using appropriate shape functions and their derivatives for each subdomain center and T is constructed using the partial differential equation and related boundary conditions where both K and T are problem dependent. The algorithm of the code development process that introduces the steps of the numerical procedure is given in Fig. 3.

Correct

4. Results Smooth

and

Rebuild nodes, generate subdomains Generate shape functions and their derivatives

Fig. 3. Flowchart of the developed code.

Simulations are carried out in a rectangular domain that represents a wave tank having a water depth h = 1 m and length L = 30 m. Linear, Stokes’ 2nd order and Stokes’ 3rd order waves are generated by assigning the analytical solutions given by Eq. (7)–(21). The results are presented for 5  150 nodes with parameters: c = 0.6, a = 1, Tm = 2  T, as = 1 and b = 1 unless otherwise stated. In the first attempt, linear waves are generated. Because the linearized free surface boundary conditions are defined on z = 0, the nodes in the domain stay fixed throughout the entire simulation which eventually simplifies the procedure. Next, fully nonlinear

226

U. Senturk / Computers & Fluids 44 (2011) 221–228

0

z (m)

-0.2 20

22

24

26

28

30

26

28

30

-0.4 -0.6 -0.8 -1

z (m)

x (m) 0.2 0 -0.2 20 -0.4 -0.6 -0.8 -1

22

24

x (m) Fig. 4. Nodal positions: at t = 0 and t = 25 s.

0.08

η (m)

0.04 0 -0.04 Linear Theory

-0.08

0

5

10

Numerical

15

20

25

20

25

30

x (m) 0.15

η (m)

0.075 0 -0.075 nd

2 Order Theory

-0.15

0

5

10

Numerical

15

30

x (m) 0.15

η (m)

0.075 0 -0.075 rd

3 Order Theory

Numerical

-0.15 0

5

10

15 x (m)

20

25

30

Fig. 5. Free surface elevations versus position at t = 20 s. (a) Linear, (b) Stokes’ 2nd, (c) Stokes’ 3rd.

wave propagation is simulated using Stokes’ 2nd and 3rd order theories, so that both problems are handled following the steps given in Fig. 3. The nodal distributions for linear and nonlinear cases are shown in Fig. 4 for 1/3 portion of the tank. Calculated free surface wave elevations throughout the tank for all investigated problems overlayed with analytical solutions are presented in Fig. 5. The wave parameters are: k = 2, H = 0.1 m for linear waves, H = 0.2 m for 2nd order waves and a = 0.1 for 3rd order waves where spanned time interval is 30 s. The sponge layer parameters are set to b = 1 and as = 1. Here, numerical simulation data are observed to be in phase with the linear theory, however 2nd and 3rd order waves tend to travel faster compared to theory which may be interpreted as numerical error in the celerity. Wave elevation history recorded by a wave probe positioned right in middle of the tank for all three

problems are shown in Fig. 6. Here the wave number is set to k = 1 and the time step is 0.05 s. Propagation of first group of waves are shown in Fig. 7 for 3rd order theory. These wave profiles are recorded at 5, 10, 15 and 20 s respectively. Leading wave has a noticeably smaller amplitude since a ramp function is used where this effect can be modified by tweaking Tm. The effect of the sponge layer length factor b is presented in Fig. 8. All three profiles have same parameters except having different length factors as b = 1, 2 and 3.The tuning factor is fixed, as = 1 as given in [33]. Finally, effects of the RBF parameters are investigated. Particularly, the free shape parameter c has a peculiar influence on the solutions. For c values smaller than 0.1 and larger than 1 lead to divergence, conversely, 0.1 < c < 1 gives close results. Three simula-

227

U. Senturk / Computers & Fluids 44 (2011) 221–228

0.08

η (m)

0.04 0 -0.04 Numerical

Linear Theory

-0.08

0

5

10

15

20

25

20

25

20

25

t (s) 0.15

η (m)

0.075 0 -0.075 nd

2 Order Theory

Numerical

-0.15 0

5

10

15 t (s)

0.15

η (m)

0.075 0 -0.075 rd

Numerical

3 Order Theory

-0.15

0

5

10

15 t (s)

η

Fig. 6. Free surface elevations versus time at x = 15 m. (a) Linear, (b) Stokes’ 2nd, (c) Stokes’ 3rd.

0

5

10

15

20

25

30

x (m) Fig. 7. Propagation of first group of waves.

0.15

η (m)

0.075 0 -0.075 β=1

-0.15 0

5

10

β=2

β=3

15

20

25

30

x (m) Fig. 8. The effect of the sponge layer length factor: b = 1, 2 and 3.

0.15

0.15

c=0.8 c=0.6

0.075

0.075

0 -0.075 -0.15

α=1 α=2

c=0.4

η (m)

η (m)

rd

3 Order Theory

0 -0.075

10

11

12

t (s) Fig. 9. The effect of the free shape parameter: c = 0.4, 0.6 and 0.8.

13

-0.15 10

11

12 t (s)

Fig. 10. The effect of the subdomain size parameter: a = 1 and 2.

13

228

U. Senturk / Computers & Fluids 44 (2011) 221–228

tion results are plotted in Fig. 9 for c = 0.4, 0.6 and 0.8 between 10–13 s. Secondly, simulations are carried out for subdomain size parameter a = 1 and 2 (Fig. 10). Increasing subdomain size populates the K matrix, consequently boosts up the computational time. For this reason, it’s found that a = 1 is sufficient for accurate solutions and no further investigations are needed. 5. Concluding remarks Increasing hardware capacities in the computer technology enable modeling complex multiphysics problems, yet, dealing with dense global system matrices will still serve as the bottleneck of global meshless RBF methods for larger scale problems. To reduce the computational cost, a localized scheme of meshless RBF method is implemented to a 2D nonlinear wave tank problem in this study. As stated in Ref. [12], problems including Neumann type boundaries often lead to amplified numerical errors when a direct collocation scheme is used. For this reason, further studies may focus on implementing fictitious points or hermite type collocations which may yield better accuracy. In order to conduct a comparative study in the presence of analytical results, the numerical wave tank problem is used as a benchmark problem, by directly specifying the wave properties on the upstream boundary. Alternatively, a prescribed wave maker motion may be implemented on C1, either by attaching analytical solutions according to an appropriate wave maker theory or by treating it as a moving boundary where update of nodal distribution, shape functions and their derivatives will be necessary for each time step. This study favors the flexibility of meshless methods particularly for moving boundary problems where only nodal positions need to be updated rather than remeshing of the entire domain at each time step during the solution process. References [1] Smith GD. Numerical solution of partial differential equations: finite difference methods. USA: Oxford University Press; 1986. [2] Zienkiewicz OC, Taylor RL. The finite element method. Bristol: Butterworth; 2000. [3] Brebbia CA, Teles JCF, Wrobel LC. Boundary element techniques. Berlin: Springer; 1984. [4] Patankar S. Numerical heat transfer and fluid flow. USA: Taylor & Francis; 1980. [5] Kansa EJ. Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics-I. Comput Math Appl 1990;19: 127–45. [6] Kansa EJ. Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics-II. Comput Math Appl 1990;19:147–61. [7] Mai-Duy N, Tran-Cong T. Approximation of function and its derivatives using radial basis function networks. Appl Math Model 2003;27:197–220. [8] Mai-Duy N, Tran-Cong T. Numerical solution of differential equations using multiquadric radial basis function networks. Neural Networks 2001;14: 185–99.

[9] Mai-Duy N, Tran-Cong T. Numerical solution of Navier Stokes equations using multiquadric radial basis function networks. Int J Numer Methods Fluids 2001;37:65–86. [10] Mai-Duy N, Tanner RI. A collocation method based on one-dimensional RBF interpolation scheme for solving PDEs. Int J Numer Methods Heat Fluid Flow 2007;17:165–86. [11] Lee CK, Liu X, Fan SC. Local multiquadric approximation for solving boundary value problems. Comput Mech 2003;30:396–409. [12] Liu GR, Gu YT. An introduction to meshfree methods and their programming. New York: Springer; 2005. [13] Šarler B, Vertnik R. Meshfree explicit local radial basis function collocation method for diffusion problems. Comput Math Appl 2006;51:1269–82. [14] Kosec G, Šarler B. Solution of thermo-fluid problems by collocation with local pressure correction. Int J Numer Methods Heat Fluid Flow 2008;18:868–82. [15] Divo E, Kassab AJ. Localized meshless modeling of natural-convective viscous flows. Numer Heat Transfer Part B 2008;53:487–509. [16] Divo E, Kassab AJ. An efficient localized radial basis function meshless method for fluid flow and heat transfer. ASME J Heat Transfer 2007;129:487–509. [17] Koo W, Kim MH. Freely floating-body simulation by a 2D fully nonlinear numerical wave tank. Ocean Eng 2004;31:2011–46. [18] Koo W, Kim MH. Fully nonlinear wave–body interactions with surfacepiercing bodies. Ocean Eng 2007;34:1000–12. [19] Chistou M, Hague CH, Swan C. The reflection of nonlinear irregular surface water waves. Eng Anal Boundary Elem 2009;33:644–53. [20] Grilli ST, Horillo J. Numerical generation and absorption of fully nonlinear periodic waves. J Eng Mech 1997;123:1060–9. [21] Guyenne P, Grilli ST. Numerical study of three-dimensional overturning waves in shallow water. J Fluid Mech 2006;547:361–88. [22] Zhang XT, Khoo BC, Lou J. Wave propagation in a fully nonlinear numerical wave tank: a desingularized method. Ocean Eng 2006;33:2310–31. [23] Chen G, Kharif C, Zaleski S, Li J. Two dimensional Navier–Stokes simulation of breaking waves. Physics of Fluids 1999;11:121–33. [24] Park JC, Kim MH, Miyata H. Fully nonlinear free surface simulations by a 3D viscous numerical wave tank. Int J Numer Methods Fluids 1999;29:685–703. [25] Hieu PD, Tanimoto K. Verification of a VOF based two-phase flow model for wave breaking and wave–structure interactions. Ocean Eng 2006;33:1565–88. [26] Guignard S, Marcer R, Rey V, Kharif C, Fraunié P. Solitary wave breaking on sloping beaches: 2-D two phase flow numerical simulation by SL-VOF method. Eur J Mech – B/Fluids 2001;20:57–74. [27] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201–25. [28] Wu NJ, Tsay TK, Young DL. Meshless numerical simulation for fully nonlinear water waves. Int J Numer Methods Fluids 2006;50:219–34. [29] Wu NJ, Tsay TK, Young DL. Computation of nonlinear free surface flows by a meshless numerical method. J Waterway, Port, Coastal Ocean Eng 2008;134:97–103. [30] Wu NJ, Tsay TK. Applicability of the method of fundamental solutions to 3-d wave–body interaction with fully nonlinear surface. J Eng Math 2009;63:61–78. [31] Xiao LF, Yang JM, Peng T, Li J. A meshless numerical wave tank for simulation of nonlinear irregular waves in shallow water. Int J Numer Methods Fluids 2009;61:165–84. [32] Longuet-Higgins MS, Cokelet ED. The deformation of steep surface waves on water I: a numerical method of computation. Proc Roy Soc Lond A 1976;350:1–26. [33] Tanizawa K. The State of the Art on Numerical Wave Tank. In: Proceeding of the 4th Osaka colloquium on seakeeping performance of ships; 2000. p. 95– 114. [34] Newman JN. Marine hydrodynamics. MIT Press; 1977. [35] Anderson JD. Computational fluid dynamics: the basics with applications. New York: McGraw-Hill; 1995. [36] Orlanski J. A simple boundary condition for unbounded hyperbolic flows. J Comput Phys 1976;21:251–69. [37] Dean RG, Dalrymple RA. Water wave mechanics for engineers and scientists. World Scientific; 1991. [38] Ma R, Li G. Spectral analysis of Stokes waves. Ocean Eng 2006;29(6):593–604.