A NEW NUMERICAL WAVE FLUME COMBINING THE 0-1 TYPE BEM AND THE VOF METHOD*

GUO Li-dong, SUN Da-peng, WU Hao

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China, E-mail: gldfirst@mail.dlut.edu.cn

(Received October 18, 2011, Revised May 10, 2012)

A NEW NUMERICAL WAVE FLUME COMBINING THE 0-1 TYPE BEM AND THE VOF METHOD*

GUO Li-dong, SUN Da-peng, WU Hao

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China, E-mail: gldfirst@mail.dlut.edu.cn

(Received October 18, 2011, Revised May 10, 2012)

A new coupling numerical wave model, based on both the Boundary Element Method (BEM) and the Volume Of Fluid (VOF) method, is established by taking advantages of the both methods to solve the wave-structure interaction problems. In this model, the wave transformation in front of structures is calculated by the 0-1 type BEM, and the intense wave motions near the structures are calculated by the VOF method. In this paper, the characteristics of the BEM and the VOF method are discussed first, and then the coupling treatments are described in detail. In the end, the accuracy and the validity of the coupling model are examined by comparing the numerical results with experiment results and other numerical results available for the interactions between regular waves with a monolayer horizontal plate.

Boundary Element Method (BEM), Volume Of Fuid (VOF) method, coupling, numerical wave model, horizontal plate

Introduction

Interactions between waves and maritime structures are an important issue in ocean and coastal engineering. The development of efficient and accurate numerical models is an important way to study the wave-structure interactions and the Boundary Element Method (BEM) and the Volume Of Fluid (VOF) method are widely used in the numerical simulation of the interactions between waves and structures.

Many non-linear wave models, established by using the BEM, can deal with both arbitrary incident waves and complex bottom topographies[1-3]. Sun et al.[4]established a numerical wave model by the 0-1 type BEM. In the model, the constant element method is used for the normal derivatives of the wave potential function, and the value of the function at the geometric center of any element is adopted for the valueof the element. Meanwhile, the linear element method is used for the wave potential function.

The VOF method is an efficient way to simulate the intense motion of waves, such as turbulence, sloshing, breaking and wave-structure interactions[5,6]. Zhao[7]carried out numerical simulations of the extreme wave generation by using the VOF method. Li et al.[8]established a numerical model based on the VOF method to investigate the regular waves interacting with a submerged horizontal twin-plate breakwater.

Though the BEM has prominent features of high accuracy and fast calculation speed, difficulties will arise when it is used to simulate intense wave motions as the BEM is based on the potential flow theory. Meanwhile, the VOF method is very suitable to capture the real-time free surface of fluid and thus can be used to describe the wave breaking and the water particle turbulence well, but this method requires a huge numerical resource in solving the N-S equations. Furthermore, the VOF models suffer from the numerical dissipation over a long distance of propagation. Hence, to simulate the wave and structure interactions efficiently and accurately, a new model should be established.

In many cases, the wave motion does not change greatly far away from structures, and the intense wave motion occurs only in regions near the structures.Therefore, it is a good choice to construct a model by combining the BEM and the VOF method together.

Based on this idea and considering the merits of the both methods mentioned above, a new numerical coupling model is established in this paper. In this model, the wave transformation in front of structures is calculated by the 0-1 type BEM, and the intense wave motions near the structures are calculated by the VOF method. This paper will be arranged as follows. The basic equations and boundary conditions are discussed in Section 1. The comparison between the BEM and the VOF method is presented in Section 2. The coupling treatments of the new model are described in Section 3. The verification of the coupling model is given in Section 4. Finally, the conclusions are presented in Section 5.

1. Basic equations and boundary conditions

1.1 The BEM

In the coordinate system for a two-dimensional potential non-linear wave problem, the computational domain boundaries can be divided into the incident boundary Sc, the free surface boundary Sf, the downstream boundary Sqand the bottom boundary Sv. The potential function of velocity φ(x,z,t) is defined, which satisfies the Laplace equation

On the incident boundary, the outside-normal velocity of the potential function can be written as

On the bottom boundary, the full reflective boundary condition is assumed

On the downstream open boundary, the Sommerfeld radiation condition is used

On the free surface boundary, the waves satisfy the kinematics and dynamics boundary conditions as follows:

The boundary conditions of the free surface Eqs.(5a), (5b) can be written in Lagrangian form and with time-stepping:

where DXF(t)=[xF(t),zF(t )]is the position vector of any fluid particle in the instantaneous free surface.

The semi-mixed Eulerian-Lagrangian method[9]is adopted in the model, so the expression of the free surface boundary conditions expressed by Eqs.(6a), (6b) can be written as:

where δ/δt=?/?t +Vp?, Vp=(0,δn/δt)is the speed of fluid particles on the free surface in the semimixed Eulerian-Lagrangian method.

P(x,z) and Q(x,z) are two arbitrary points in the calculation domain and on the boundary of the calculation domain, respectively. According to the Green integration theorem, Eq.(1) can be expressed as

By introducing the boundary conditions Eqs.(2) to (5a), (5b) and (7a), (7b), the integral equation under the conditions of the steady boundaries can be established as follows

Next, the BEM integral equation Eq.(9) is discretized by the 0-1 combined element method[4]. In addition, the potential function and the normal derivative of the potential function are expanded in the time domain, and then the linear equation system with φΔ and ηΔ as unknown quantities can be obtained as

The evaluations of the equation at the time step n are based on the iterative loop after the wave surface subdivision at the time step 1n-. This process of iteration continues until the error of computation of the increment of the wave surface at the time step n is within the permitted limit, that is,

where ε is the given allowable error, m is used for the thm cycle at the time step n.

1.2 The VOF method

The governing equations for the VOF method are the continuity equation, the Reynolds time-averaged equation and the two-equation kε- model[10]for an incompressible fluid:

where x and z are the horizontal and vertical coordinates, respectively, u and v are the mean velocity components in the x and z directions, respectively, g is the acceleration of gravity, p is the pressure, ρ is the fluid density, ν is the coefficient of the kinematic viscosity, and νt=Cuk2/ε is the coefficient of the turbulent viscosity,1Cε,2Cε,kσ and εσ are the parameters in the kε- model, and in the present study the following standard values are used: Cε1=1.43, Cε2=1.92, σk=1.0, σε=1.30.

The free surface is described by means of a fluid fraction function F(x,z,t), which has a value between one to zero, representing the volume fraction of a cell being occupied by the fluid. Thus, F=1 denotes a cell full of fluid, while an empty cell will have F=0. A cell with F value between one and zero is either intersected by a free surface or contains voids that will be partially filled with fluid. Further, a free surface cell can be identified as a cell containing a non-zero value of F and having at least one neigh-boring cell where F=0. The time variation of the function F(x,z,t) is governed by

The finite difference method is adopted to carry out the numerical calculations. The staggered cell method is used to generate the computational grids and then the values of all variables in each cell can be obtained.

A damping layer is employed to absorb the wave energy gradually on the downstream boundary, and the non-slip and impermeable conditions are adopted on both the obstacle surface and the bottom boundary and the bottom boundary.

In the model, the time intervaltΔshould satisfy the Courant condition so that the stable numerical calculations are completed as soon as

2. Comparison between the BEM and the VOF method

2.1 Verification of the BEM and the VOF method

To validate the both numerical models for simulating the propagation of waves, a comparison is made between the numerical results obtained by the models and the existing experimental data for the transformation of regular waves over a submerged breakwater. The physical experiment[11]was conducted in a flume of 23 m long and 0.8 m wide and with a water depth of 0.4 m (Fig.1). A submerged breakwater with a slope β=1/10 for both incident and leeside slopes was installed in the wave flume. The regular wave period is T=1.6s , and the wave height is H= 0.2m.

Fig.1 Experiment layout of wave propagation over a submerged breakwater

The comparison between the experimental data and the numerical results is shown in Fig.2, which shows that the discrepancy is very small and can be ignored, indicating that both of the BEM and the VOF methods can simulate the propagating waves well.

Fig.2 Comparison between computed wave surface elevation and experimental data

2.2 Computational stability, accuracy and numerical dissipation

The numerical simulation of the regular wave propagation along a horizontal bottom is carried out (the water depth D=0.48m , the flume length L= 18m) by the BEM and the VOF method, respectively. The regular wave period is T=1.4s , and the wave height is H=0.06m .

Fig.3 The time history of wave profiles at =3mx

Figures 3 and 4 are the time history of the wave profile at x=3m and x=13.8m calculated by the BEM and the VOF method, respectively.

Fig.4 The time history of wave profiles at x=13.8m

At=3mx, the amplitude of the waves obtained by the BEM (Fig.3(a)) increases gradually from 0 m to 0.03 m, and then the wave profile is almost stable, while, in the case of the VOF method, before the wave amplitude reaches a stable value, in the process of the wave amplitude increasing gradually (Fig.3(b)), there is an unstable wave with the wave height greater than H. That shows that the stability and the accuracy of the BEM method are better than those of the VOF method in the simulation of wave propagation.

Atx=13.8m (about 5 wavelength distance away from the wave board), the wave amplitude of the time history of the wave profile obtained by the BEM (Fig.4(a)) increases from 0 m to about 0.03 m. However, the maximum wave height of the time history of the wave profile obtained by VOF is less than 0.06 m (Fig.4(b)). In view of the stable wave height of about 0.05 m, it is indicated that the numerical dissipation in the VOF method is very serious. Thus, the VOF method is not suitable for large-scale or long-range numerical calculations.

2.3 Free surface tracking

With the BEM method, based on the potential flow theory, it is difficult to track accurately the surface of the waves near the stage of breaking, such as in the case of studying the overturning phenomenon of the water surface. In contrast, the VOF method, based on the viscous theory, is well known for its capacity to simulate the free surface of fluid. The VOF models were used to model the wave breaking and the postbreaking[12,13]. Figure 5 shows that the process of the wave breaking is simulated by the VOF method.

Fig.5 The process of wave breaking obtained by VOF method

2.4 Computational time and storage use

With the BEM method, only the boundaries of the computational domain are subdivided and solved, which means that not too much computational time and storage are required.

On the other hand, the VOF method is based on the finite difference method or the finite element method, where the whole computational domain is to be subdivided and solved. The refinement and the diversity of the computational grids require enormous computation time and storage. Therefore, the VOF method is far more expensive than the BEM in terms of computational time and storage requirements to obtain computational results of the same resolution.

Fig.6 Schematic diagram of the coupling BEM+VOF meathod

3. Coupling of the BEM and the VOF method

3.1 Construction of the coupling model

With the merits of the 0-1 type BEM method and the VOF method in mind, a new hybrid model is developed in this paper by combining the 0-1 type BEM with the VOF method. In the model, the intense wave motions at the proximity of the structures are calculated by the VOF method and the rest of the computational domain is modeled by the BEM (Fig.6).

The new model has the following features:

(1) The model enjoys the merits of stability and precision of the BEM in the wave propagation and deformation, and reduces the numerical dissipation of the VOF method. Thus the accuracy of the model is superior to the VOF method, and it can be used forlarge-scale or long-distance numerical calculations.

(2) This model can be used for strong non-linear problems, such as wave breaking problems, which can not be adequately tackled by the BEM, and it makes full use of capability of the VOF method in tracking the free surface.

(3) The model can considerably reduce the timeconsuming VOF domain, so it is superior to the VOF method in terms of computational time and storage requirements.

3.2 Boundary matching conditions

As shown in Fig.7, the regions marked by 1 and 2 represent the front-domain (the BEM computational domain) and the back-domain (the VOF computational domain), respectively.qS is the downstream boundary of the front-domain,oS is the incident boundary of the back-domain. The domain between Soand Sqis the Transmission Domain (TD). The length of TD is chosen as one wavelength in this paper. The boundariesoS andqS are also known as the matching boundary of the transmission domain.

Fig.7 The TD

The BEM requires the surface elevation and velocity conditions on the matching boundaryqS before the calculation at each time step. The matching boundaryqS is the downstream boundary of the front-domain, and contains also the computational nodes in the back-domain. Thus, the matching boundary conditions for the front-domain can be used by the flow field of the back-domain, that is:

The VOF method requires the velocity, the pressure and the volume of fluid function on the matching boundaryoS, before the calculation at each time step. The matching boundaryoS is the incident boundary of the back-domain, and contains also the computational nodes in the front-domain. Thus, the matching boundary conditions for the back-domain can be used by the flow field of the front-domain:

where z(k-1/2) and z(k+1/2) represent the vertical coordinates of points below and above the kth layer cell on the boundary So(Fig.8).

Fig.8 Sketch of the function ()Fk on the boundaryoS

3.3 Connection of the water surfaces

In the TD, the water surfaces calculated by the BEM and the VOF method do not necessarily coincide due to the difference in their calculation accuracy. That would be a cause for the instability of computation. To deal with this problem, the connection method proposed by Kim et al.[14]is adopted here. The water surfaces are connected smoothly in the TD at each time-step by using the weighted means with respect to the horizontal distance to the TD. The BEM surface line in the TD is modified to agree with thereference surface line, while the VOF surfaces in the TD are kept intact. As the BEM surface is modified, φ andnφ should be adjusted according to the new water surface. The evolution of the water surface and the velocity potential are derived from Eqs.(7a), (7b) with respect to the time increment tΔ as:

where

The value ofnφ at the previous time step (1)m- is corrected by Eqs.(24a), (24b).

3.4 Solution process of the coupling model

The incident boundary condition is given in the front-domain, and the full reflective boundary condition or the absorbing boundary condition is given in the back-domain. The two-way synchronization solution process of the coupling model is as follows:

The back-domain VOF model is started to be used when the waves propagate to the matching boundaryoS, and the both models receive the fluid information and affect each other on the matching boundaries. In the time step n, the results of the front-domain for the BEM model are obtained first, and then the surface elevation and the velocity on the boundary Sqare made to match the boundary conditions Eqs.(17) to (19), obtained by the VOF model at the time step 1n-. Secondly, the velocity, the pressure and the volume of fluid function on the boundaryoS, satisfying the boundary conditions Eqs.(20) to (23), are calculated by the BEM method. Lastly, according to Eqs.(17) to (19), the fluid information on the boundaryqS is obtained by the VOF method. Then, the boundary conditions of the front-domain at the time step +1n are obtained.

4. Verification of the new model

4.1 Numerical verification for regular waves

The numerical verification is made for regular waves in the area shown in Fig.9. The right boundary condition of the VOF domain is an open or wall condition. The water depth D=0.48m , the width of the TD is set as 0.3L, the regular wave period is T= 1.4s, the wave height is H=0.06m , and the time step is Δt=0.02s .

Fig.9 Domain for the numerical verification of the model

Fig.10 Comparison of water profile between model and the theoretical value

Figure 10 shows the comparison for the spatial profiles of the water surface ()xη between the numerically calculated values and the theoretical ones. It is seen that the new model gives values in very good agreement with the theoretical ones, indicating that the coupling 0-1BEM+VOF model is successful for two way propagations of incident and reflected waves.

Moreover, the computational time ratio between the BEM domain and the VOF domain is 1:3 in the same time step (CPU: AMD Athlon(tm) 64×2, memories:1GB). Therefore, the reasonable reduction of the VOF domain by using 0-1 BEM+VOF model allows an efficient computation to simulate the interaction between the waves and the structures.

4.2 Experiment verification of the coupling 0-1BEM+

VOF model

Ren[15]carried out an experimental study of wave impacts on a horizontal plate. The problem is obviously strongly nonlinear, due to the water particle turbulence and the wave breaking. Therefore, the 0-1 BEM+VOF model is used to solve the problem.

Fig.11 The sketch of the experimental setup of Ren[15]

Fig.12 Sketch of pressure transducers on the subface of the horizontal plate

The sketch of the experimental flume is shown in Fig.11. There is a wave maker on the left side of the flume to generate waves, and on the right side, there is a 1:30 slope beach to absorb the wave energy. The water depth is 0.6 m, and the distance between the wave-maker and the horizontal plate is 9 m. The distance between the wave-maker and the wave gauge is 8 m. The vertical space between the still water level and the horizontal plate is s, which will be called the headroom in this paper. In our numerical model, the distance from the wave-maker to the point =7mx is defined as the front domain, and the back-domain is from the point =6mx to the end of the flume. Eleven pressure measuring points are located along the horizontal plate, as shown in Fig.12. The wave height is chosen as 0.1m. Regular waves with periods of 1.2 s, 1.5 s and 2.0 s are generated. Four different s (/=sH0.1, 0.2, 0.3 and 0.4) are considered for each wave period.

4.2.1 Verification of water surface elevations

Figures 13 to 15 show the comparison of the numerical and measured water surface elevations at the wave gauge. It can be seen that the calculated values agree well with the experimental results, indicating that this model can be used to simulate the strongly non-linear wave propagation and deformation at different wave conditions and different headrooms. According to the comparison of the wave profiles at three wave periods, it can be found that the non-linear wave surface is strongly influenced by the horizontalplate when the wave period T is small, as shown in Fig.13(b). When the wave period T is large, however, the non-linearity is less severe, with better accuracy than that of the cases with small T, as shown in Fig.15(c).

Fig.13 The comparison between the numerical and measured water surface elevations at different headrooms (H= 0.1m, T=1.2s )

Fig.14 The comparison between the numerical and measured water surface elevations at different headrooms (H= 0.1m, T=1.5s )

Fig.15 The comparison between the numerical and measured water surface elevations at different headrooms (H= 0.1m, T=2.0s )

4.2.2 Verification of impact pressure

In the time domain, the process of the regular wave impacting is obviously cyclical and strongly stochastic. In order to compare the calculated results withthe experimental results, the data analysis method of Ren[15]is adopted, that is, taking a third of the peak of the impact pressure,1/3cP, as the characteristic value in the statistical analysis.

Fig.16 Pressure distribution on the subface of the horizontal plate (H=0.1m , T=1.2s)

Figures 16 to18 show the impact pressure comparison between the numerical results and the experimental results. At the same time, the numerical results obtained by using the Smooth Particle Hydrodynamics (SPH) method by Zheng[16]are also given in the figures.

Fig.17 Pressure distribution on the subface of the horizontal plate (H=0.1m , T=1.5s)

By comparing the numerical results with the experiment results and the previous numerical results, the following remarks can be made:

(1) The results of the present model and those obtained by Zheng[16]are in good agreement with experimental results.

(2) For most cases of overall impact pressure distributions, this model gives results in better agreement with experimental data than those obtained by Zheng[16](with better curve fitting). However, for the case of T=1.2s , s/H=0.1 (Fig.16), the present results are less consistent, which may be due to the stronger non-linearity when T is small.

Fig.18 Pressure distribution on the subface of the horizontal plate (H=0.1m , T=2.0s )

(3) In most cases, the numerical results show better performance than those obtained by Zheng[16]with respect to the experimental data in a single measuring point, indicating that the present model enjoys better numerical calculation accuracy than the SPH method. For the case of T=1.2s , s/H =0.3 (Fig.16), there are differences between the experimental results and the both numerical results at the measuring point 4, while the numerical results of this model is closer to the experimental results than that of Zheng[16].

(4) The agreement between the numerical results and the experimental results becomes better with the increase of the wave period.

5. Conclusions

A new numerical model combining the 0-1 type BEM and the VOF method is established by considering the merits of the both methods. The new method enjoys the merits of the BEM in stability and precision in simulating wave propagation and deformation, and reduces the numerical dissipation in the VOF method. The present model makes a full use of the VOF method in tracking the free surfaces near the structures, and can simulate wave breaking and post-breaking. In addition, the model can considerably reduce the computation time and the storage requirements of the VOF method.

By comparisons of the water surface elevations in front of the horizontal plate and the impact pressure on the sub-face of the plate obtained by the present model with the experiment results[15]and other numerical results[16], the validity and the accuracy of the coupling model for wave impact programs are shown.

[1] SUNG H. G., SHOON C. H. Implicit formulation with the boundary element method for nonlinear radiation of water waves[J]. Engineering Analysis with Boundary Elements, 2010, 34(5): 511-529.

[2] KOO W. C., KIM M. H. Numerical simulation of nonlinear wave and force generated by a wedge-shape wave maker[J]. Ocean Engineering, 2006, 33(8-9): 983-1006.

[3] NING D. Z., TENG B. and EATOCK TAYLOR R. et al. Numerical simulation of non-linear regular and focused waves in an infinite water-depth[J]. Ocean Engineering, 2008, 35(8-9): 887-899.

[4] SUN Da-peng, LI Yu-cheng and GUO Hai-bin. A new pursuit mode of nonlinear wave surface[J]. China Ocean Engineering, 2004, 18(2): 257-266.

[5] HSIAO Shih-Chun, LIN Ting-Chieh. Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling[J]. Coastal Engineering, 2010, 57(1): 1-18.

[6] ZHANG Q., LIU P. L.-F. A numerical study of swash flows generated by bores[J]. Coastal Engineering, 2008, 55 (12): 1113-1134.

[7] ZHAO Xi-zeng, HU Chang-hong and SUN Zhao-chen. Numerical simulation of extreme wave generation using VOF method[J]. Journal of Hydrodynamics, 2010, 22(4): 466-477.

[8] LI Jing-bo, ZHANG Ning-chuan and GUO Chuansheng. Numerical simulation of waves interaction with a submerged horizontal twin-plate breakwater[J]. China Ocean Engineering, 2010, 24(4): 627-640.

[9] BECK R. F., REED A. M. Modern computational methods for ships in a seaway[J]. Society of Naval Architects and Marine Engineers, 2001, 109: 1-51.

[10] SHEN Y. M., NG C. O. and ZHENG Y. H. Simulation of wave propagation over a submerged bar using the VOF method with a two-equation kε- turbulence modeling[J]. Ocean Engineering, 2004, 31(1): 87-95.

[11] ZOU Zhi-li, WANG Tao and ZHANG Xiao-li et al. One-dimensional numerical models of higher-order Boussinesq equations with high dispersion accuracy[J]. Acta Oceanologica Sinica, 2003, 22(2): 287-299.

[12] LUBIN P., GLOCKNER S. and KIMMOUN O. et al. Numerical study of the hydrodynamics of regular waves breaking over a sloping beach[J]. European Journal of Mechanics - B/Fluids, 2011, 30(6): 552-564.

[13] JIANG C. B., CHEN J. and TANG H. S. et al. Hydrodynamic processes on beach: Wave breaking, up-rush, and backwash[J]. Communications in Nonlinear Science and Numerical Simulation, 2011, 16(8): 3126-3139.

[14] KIM Sang-Ho, YAMASHIRO Masaru and YOSHIDA Akinori. A simple two-way coupling method of BEM and VOF model for random wave calculations[J]. Coastal Engineering, 2010, 57(11-12): 1018-1028.

[15] REN B., WANG Y. Experimental study of irregular wave impact on structures in the splash zone[J]. Ocean Engineering, 2003, 30(18): 2363-2377.

[16] ZHENG Kun. Study on wave impact on horizontal plate based on SPH method[D]. Ph. D. Thesis, Dalian: Dalian University of Technology, 2009(in Chinese).

10.1016/S1001-6058(11)60272-2

* Project supported by the National Natural Science Foundation of China ( Grant No. 50921001), the Foundation of State Key Laboratory of Coastal and Offshore Engineering, Dalian University on Technology (Grant No. LP0804).

Biography: GUO Li-dong (1981-), Male, Ph. D. Candidate

SUN Da-peng,

E-mail: dpsun@dlut.edu.cn