NUMERICAL WAVE FLUME WITH IMPROVED SMOOTHED PARTICLE HYDRODYNAMICS*

ZHENG Jin-hai

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China, E-mail: jhzheng@hhu.edu.cn

SOE Mee Mee, ZHANG Chi

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China

HSU Tai-Wen

Department of Hydraulic and Ocean Engineering, National Cheng Kung University, Tainan 701, China

NUMERICAL WAVE FLUME WITH IMPROVED SMOOTHED PARTICLE HYDRODYNAMICS*

ZHENG Jin-hai

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China, E-mail: jhzheng@hhu.edu.cn

SOE Mee Mee, ZHANG Chi

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China

HSU Tai-Wen

Department of Hydraulic and Ocean Engineering, National Cheng Kung University, Tainan 701, China

(Received April 25, 2010, Revised August 8, 2010)

This article presents an improved Nearest Neighboring Particle Searching (NNPS) technique for numerical modeling of water waves with the Smoothed Particle Hydrodynamics (SPH) method. The proposed technique differs from others by introducing the concept of Inner and Outer Particle Searching (IOPS) and shifting most of advanced CPU operations into simple addition operations. The IOPS method is shown to significantly improve the computational efficiency and reduce the CPU time especially for large number of particles, based on comparisons with other two NNPS methods. This method is implemented in a 2DV numerical wave flume conducted by the SPH method. Three test cases are examined, including generations and propagations of dam-breaking induced waves, solitary wave and irregular wave. Calculated results are in good agreements with experimental data and theoretical solutions with fairly satisfactory CPU time-consuming. The wave motions observed in physical facilities are successfully reproduced by the SPH numerical wave flume, revealing its robust capability of modeling realistic wave propagation and substantial potential for a wide variety of hydrodynamic problems.

particle searching technique, Smoothed Particle Hydrodynamics (SPH), numerical wave flume, Navier-Stokes equation

1. Introduction

Smoothed Particle Hydrodynamics (SPH) method was originally developed for astrophysical applications, and has later been extended to a number of other fields. The basic idea of SPH is to use the collective motions of large number of particles to represent a flow in a Lagrangian way rather than Eulerian way. In a particle approach, the governing equations are discretized and solved with respect to the individual particles filled within the computational domain. This method is conceptually simple without adding new physics and high accuracy can be achieved by increasing number of particles.

To date, the SPH becomes increasingly popular and finds wide applications in computational fluid mechanics, including free surface flow[1], porous flow[2], landslide-induced flow[3]and fluid-structure interactions[4]. In particular, the SPH method has been adapted in a variety of numerical wave flumes for studying green water overtopping[5,6], wave breaking in the surf zone[7,8], and wave-structure interaction[9-11]. The most attracting feature of SPH for wave modeling is that it naturally needs no special approach for dealing with the free surface. It is meshfree and every particle at the free surface can be easily tracked in the Lagrangian frame. When conducting a numerical wave flume involved with flow separation and large deformation, this appears a significant advantage over other traditional Eulerian methods, including those

2. Numerical method

The basic concept of SPH method is to treat a flow as a sum of moving particles. Each particle has its own physical quantities such as mass, velocity, density and pressure gradient. During flow motions, all particles change their positions and corresponding properties with time. The SPH model then solves the trajectory of each particle. Through the use of integral interpolations, the field variables are expressed by integrals approximated by summation interpolations over neighboring particles. Firstly, the functions of variables at an arbitrary location are transformed into integral forms, referred to the kernel approximation. Secondly, those integral expressions are approximated by summations of variables of relevant scattered particles.

2.1 Kernel approximation

In a continuous variable field, any function can be expressed in terms of its values at a set of particles by use of a weighting function

where A is the function with respect to the spatial tensor, Ω denotes the computational domain, and δ(r?r′) is Dirac delta function. Replacing the Dirac delta function by a kernel W, leads to

where h is the smoothing length which controls the size of the area around a given particle where the contribution from other particles can not be neglected, and W is a cubic spline kernel. Representing the approximation process with a sign of〈〉and introducing the gradient of A(r), we get

Equation (4) can be written as

provided that the filtering range kh is inside the computational domain . Ω

Note that, if the filtering range extends beyond the computational domain, Eq.(5) will break down for the lack of particles and a special treatment should be applied at the boundary[19].

2.2 Particle approximation

If the number of particles constructing the whole domain is N, Eq.(5) can be approximated by a summation, given as

where r, m and ρ represent the spatial vector, mass and density of an individual particle, respectively. Subscript b represents the index number.

It can be seen from Eq.(6) that the variable A at the point r is affected by all N particles in the computation domain. In practice, the influence of kernel is restricted to a radial distance of an order of kh. That means only the particles within this distance (so-called neighboring particles) are considered to contribute to the summation in Eq.(6). While particles keep moving with flow motions, the relative positions of all particles should be calculated to find the effective particles for each point at every time step. This work is done by the NNPS method.

2.3 IOPS method

This section describes in detail the concept and implementation of the IOPS method proposed in this study. The IOPS method saves the computational time by shifting most of advanced CPU operations (e.g., multiplication, division) into simple addition operations. Firstly, each individual particle is marked by the so-called inner and outer grids. Secondly, for a given particle inside the domain, the particles in its filtering control unit in the inner grid is marked. Finally, the distances between this particle and other particles in the filtering control unit are computed to find the nearest neighboring particles.

If the number of particles is N and the number of inner grids is m for the whole domain, the number of particles in each inner grid would be N/m . Figure 1 shows the grid structures. Overlaps can be seen between boundaries of neighboring outer grids. The number of particles in each outer grid is larger than that in each inner grid. When total number of particles is relatively large, however, that difference can be neglected in the estimation of computational effort. Therefore, the computation order of marking all particles by inner and outer grids is O(2×N×m).

For a given particle in an inner grid, a square filtering control unit is defined with its side length of kh. Then, particles belonging to both the filtering unit and the outer grid are marked, with a computation order of O(N2/m2). For the whole domain, that is O(N2/m). Till now, only particles inside the filtering control unit need to be considered to search the nearest neighboring ones. The computation order for this process appears a constant and can no longer be reduced. Eventually, it is easy to find that the smallest total computation order can only be obtained when 2×N×m+N2/m is smallest. That is whenwe get the minimum CPU time consuming.

2.4 Numerical wave flume

In Lagrangian form, the conservation of mass and the conservation of momentum for a Newtonian incompressible fluid can be expressed as

where u is the particle velocity tensor, p the fluid pressure,ν the kinematic viscosity, F the body force per unit mass.

The pressure of each particle is obtained from the equation of state. By artificially enhancing the fluid compressibility, the following form is proposed

then the particle is regarded as a free-surface particle (constant β=0.8-0.98) for which a density equal to ρ0is imposed. This treatment is based on the fact that the calculated particle density on the free surface drops abruptly for the lack of particles in the outer region of the free surface.

3. Model results

where H is the wave height, d the water depth, and C the wave celerity.

The time-evolution of surface elevation at threemeasuring points x=1.5m , x=2.5m and x=3.5mare collected to verify the numerical results, as shown in Fig.8. The three surface profiles show the good consistency. For each point, the surface elevation remains stationary before wave passing through it and the peak value of surface displays no significant deviation from other two points. That means a stable pressure field is obtained in computation. The observed time lag between two adjacent points is 0.917 s, and the corresponding calculated wave celerity is 1.09 m/s, which agrees well with the analytical value of 1.08 m/s.

Figure 9 presents the simulated free surface at t=4s, as well as the first-order Boussinesq solution. The overall agreement is satisfactory. The calculated profile is steeper than the theoretical result in front of wave crest while milder behind it, which is due to the kinematic viscosity included in the numerical model but not presented in the analytical solution.

Figure 10 illustrates the velocity field at t=3s. Realistic result is shown. The horizontal velocity increases with decreasing depth with its peak value locating at the wave crest, where the vertical velocity is negligible. In the water column where wave has passed through, some minor disturbances are found. Nevertheless, generally realistic distributions of both value and direction of velocity are maintained for the whole domain, justifying the ability of the numerical wave flume developed in the present study.

3.3 Irregular waves

The numerical flume has a length of 5 m and a height of 1 m, as well as an initial water depth of 0.65 m, as shown in Fig.11. A 0.9 m high wavemaker is placed at the left boundary of flume to generate irregular waves, and an artificially absorbing sponge layer is configured in the right flume section from 13 m to 16 m. The particle setting in the SPH method is the same as in the case of solitary wave. The time-series of surface elevation are recorded at seven points located at 7.5 m, 8 m, 8.5 m, 9 m, 9.5 m, 10.0 m and 10.5 m along the horizontal axis.

where1smf and2smf are smoothing functions expressed as

The smoothed time-series of wavemaker position and velocity are displayed in Fig.13. With this method, a continuous signal is obtained for the numerical experiment.

Figure 14 presents the comparisons of free surface variations at Point 3 and Point 5 between experimental data of Cox and Ortega[27]for the case without a deck, the SPH simulation results from present study and Gomez-Gesteira et al. model[5](hereafter GGM). The simulated profiles are in generally good agreements with the GGM results and the measurements, both in amplitude and phase, although there are some slight discrepancies between numerical and experimental results for secondary peaks. These discrepancies are probably due to the smoothing effects of wavemaker signal. Particularly for the phase of the first wave and the height of the highest wave, the present model provides fitter results to the observation than GGM’s.

Figure 15 shows the particle snapshots in the righter region of flume at t=11s, t=13s and t=15s. The primary aim of this figure is to examine the efficiency of the sponge layer in the model. Wavesenter the sponge layer (x =13m-16m ) att=13s and begin to dissipate energy in this region. At t=15s, the surface variation is undistinguishable, indicating that most of wave energy has been effectively absorbed due to the viscous effect in the sponge layer.

It can be reasonably concluded that the SPH numerical wave flume developed in this study is capable of generating reliable irregular wave trains and reproducing their propagations with good stability and boundary descriptions. Numerical results show the improvements over that of GGM to a reasonable extent. Moreover, only 25 min is needed to complete the computation for this case. This does benefit from the more efficient inner and outer particle searching technique proposed in this paper.

4. Conclusions

In this paper, an IOPS technique has been proposed for the SPH simulation. This technique effectively reduces the time requirement of computation by introducing inner and outer grids and shifting most of advanced CPU operations into simple addition operations. It is particularly useful to improve the SPH applications to larger study areas, at least partly. Employing the IOPS method, a SPH numerical wave flume is conducted based on the Navier-Stokes equations. Three test cases are used to investigate its performance, including dam breaking, generation and propagation of solitary wave and irregular wave. Numerical results are in good agreements with theoretical solutions, experimental data and results provided in previous studies. It is concluded that the SPH numerical wave flume with the IOPS method presented in this study can reproduce the generation and propagation of realistic waves in a physical tank, which could act as the foundation for further studies involving more complicated flows.

[1] ZHENG Kun, SUN Zhao-chen and SUN Jia-wen et al. Numerical simulations of water wave dynamics based on SPH methods[J]. Journal of Hydrodynamics, 2009, 21(6): 843-850.

[2] CHU Yang, LU Wen-qiang. Porescale SPH simulation of flow in porous media[J]. Journal of Engineering Thermophysics, 2009, 30(3): 437-440(in Chinese).

[3] DU Xiao-tao, WU Wei and GONG Kai et al. Two dimensional SPH simulation of water waves generated by underwater landslide[J]. Journal of Hydrodynamics, Ser. A, 2006, 21(5): 579-586(in Chinese).

[4] XU Qing-xin, SHEN Rong-ying. Fluid-structure interaction of hydrodynamic damper during the rush into the water channel[J]. Journal of Hydrodynamics, 2008, 20(5): 583-590.

[5] GOMEZ-GESTEIRA M., CERQUEIRO D. and CRESPO C. et al. Green water overtopping analyzed with a SPH model[J]. Ocean Engineering, 2005, 32(2): 223-238.

[6] SHAO S. D., JI C. M. and GRAHAM D. I. et al. Simulation of wave overtopping by an incompressible SPH model[J]. Coastal Engineering, 2006, 53(9): 723-735.

[7] DARLYMPLE R. A., ROGER B. D. Numerical modeling of water waves with the SPH method[J]. Coastal Engineering, 2006, 53(2-3): 141-147.

[8] KHAYYER A., GOTOH H. and SHAO S. D. Corrected incompressible SPH method for accurate water-surface tracking in breaking waves[J]. Coastal Engineering, 2008, 55(3): 236-250.

[9] GOMEZ-GESTEIRA M., DALRYMPLE R. A. Using a 3D SPH method for wave impact on a tall structure[J]. J. Wtrwy. Port, Coastal Ocean Eng., 2004, 130(2): 63-69.

[10] SHAO S. D. SPH simulation of solitary wave interaction with a curtain-type breakwater[J]. Journal of Hydraulic Research, 2005, 43(4): 366-375.

[11] GONG Kai, LIU Hua and WANG Ben-long. Water entry of a wedge based on SPH model with an improved boundary treatment[J]. Journal of Hydrodynamics, 2009, 21(6): 750-757.

[12] LV Biao, JIN Sheng and AI Cong-fang. A conservative unstructured staggered grid scheme for incompressible Navier-Stokes equations[J]. Journal of Hydrodynamics, 2010, 22(2): 173-184.

[13] CHRISTENSEN E. D., DEIGAARD R. Large eddy simulation of breaking waves[J]. Coastal Engineering, 2001, 42(1): 53-86.

[14] ZOU Lin, LIN Yu-feng and LAM Kit. Large-eddy simulation of flow around cylinder arrays at a subcritical Reynolds number[J]. Journal of Hydrodynamics, 2008, 20(4): 403-413.

[15] SUN D. P., LI Y. C. Wave damping absorber in numerical wave tank and calculation of wave transformation[J]. China Ocean Engineering, 2000, 18(2): 46-50.

[16] MONAGHAN J. J. Particle methods for hydrodynamics[J]. Computer Physics Report, 1985, 3(2): 71-124.

[17] HERNQUIST L., KATZ N. Tree SPH-A unification of SPH with the hierarchical tree method[J]. The Astrophysical Journal Supplement Series, 1989, 70: 419-446.

[18] MIHAI B., MARTY L. and NATHAN Q. Grid-assisted particle search in smoothed particle hydrodynamics[R]. Galway faculty of Engineering Research Day, 2004.

[19] LIU G. R., LIU M. B. Smoothed particle hydrodynamics: A meshfree particle method[M]. Singapore: World Scientific, 2003, 6-27.

[20] KOSHIZUKAR S., NOBE A. and OKA Y. Numerical analysis of breaking waves using the moving particle semi-implicit method[J]. International Journal for Numerical Methods in Fluid, 1998, 26(7): 751-769.

[21] VIOLEAU D., ISSA R. Numerical modelling of complex turbulent free-surface flows with the SPH method: An overview[J]. International Journal for Numerical Methods in Fluids, 2007, 53(2): 277-304.

[22] SHAO S. D., LO E. Y. M. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[J]. Advances in Water Resources, 2003, 26(7): 787-800.

[23] MARTIN J. C., MOYCE W. J. An experimental study of the collapse of liquid columns on a rigid horizontalplane[J]. Philosophical Transaction of the Royal Society of London, Ser. A, 1995, 244(882): 312-324.

[24] HIRT C. W., NICHOLS B. D. Volume of fluid method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201-225.

[25] HARLOW F. H., WELCH F. J. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[J]. Physics of Fluids, 1965, 8(12): 2182-2189.

[26] XIAO Bo. The theory of solitary wave generating in laboratory and its validity[J]. Chinese Port Engineering, 1991, 2: 19-24(in Chinese).

[27] COX D. T., ORTEGA J. A. Laboratory observations of green water overtopping a fixed deck[J]. Ocean Engineering, 2002, 29(14): 1827-1840.

10.1016/S1001-6058(09)60115-3

* Project supported by the Program for New Century Excellent Talents in University of China (Grant No. NCET-07-0255), the Special Fund of State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University (Grant No. 2009585812).

Biography: ZHENG Jin-hai (1972-), Male, Ph. D., Professor