- Email: [email protected]
Numerical simulations of water wave dynamics based on SPH methods
843
2009,21(6):843-850 DOI: 10.1016/S1001-6058(08)60221-8
NUMERICAL SIMULATIONS OF WATER WAVE DYNAMICS BASED ON SPH METHODS* ZHENG Kun, SUN Zhao-chen, SUN Jia-wen State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China, E-mail: [email protected] ZHANG Zhi-ming, YANG Guo-ping, ZHOU Feng CCCC Water Transportation Consultants Co. Ltd, Beijing 100007, China
(Received July 29, 2009, Revised October 29, 2009)
Abstract: A numerical model was established for simulating water wave dynamic problems by adopting the Smoothed Particle Hydrodynamics (SPH) methods of iterative solution of Poisson’s equation for pressure field, and meanwhile the sub-grid turbulence model was applied in the simulation so as to more accurately describe the turbulence characteristics at the time of wave breaking. In this article, simulation of the problem of the dam collapsing verifies the compoting accuracy of this method, and its results can be identical with the results of VOF method and the experimental results by comparison. Numerical simulations of the course of solitary wave and cnoidal wave run-up breaking on beaches were conducted, and the results are basically consistent with experimental results. This indicates that the SPH method is effective for the numerical simulation of the complex problems of water wave dynamics. Key words: Smoothed Particle Hydrodynamics (SPH), Poisson’s equation, water wave dynamics
1. Introduction Wave breaking in near-shore has been a hot issue in coastal engineering. The strong nonlinearity of wave breaking arouses the complex movement of water body. As a result, numerical simulation of wave breaking process on beaches has become a difficult problem. Due to the restriction of the grid scale, in the past, grid-based approach has limited accuracy in dealing with multiple free-surfaces in the process of breaking, while the particle method is not subject to restrictions of the grid-scale and thus can simulate more accurately the problem of multiple free-surfaces of fluid motion on beaches. In this article numerical simulation of the run-up breaking process of solitary wave on beaches is * Project supported by the National High Techology Research and Development Program of China (863 Program, Grant No. 2007AA11Z130) Biography: ZHENG Kun (1977-), Male, Ph. D. Candidate
conducted by using a grid-free method, the Smoothed Particle Hydrodynamics (SPH)[1-4]. The standard SPH method is to achieve the incompressibility of water body by adopting artificial compressing technique. Pressure is taken as a Lagrange factor, and the pressure field is expressed and solved by using the state equation. This method has the advantage of directly giving the relationship between density and pressure. But its computation of the pressure field is still accurate enough. In order to compute the pressure field more accurately, pressure is solved by employing the SPH methods of iterative solution of Poisson’s equation for pressure field, and the sub-particle equation turbulence model is adopted[5-7]. The numerical results indicate that it is more accurate and stable compared with the solution by using state equation.
844
2. Basic equations for SPH 2.1 Governing equations The N-S equation described can be written in the form
S = (2 Sij Sij )1/ 2 , and ' stands for the initial spacing of particles. The final form of Eqs.(3) and (4) can be written as
1 dU +
(1)
dU i = ¦ m j U ij iWij dt j
dU 1 = P + F +Q 0 2U dt U
(2)
§ P P · dU i = ¦ m j ¨ j2 + i 2 ¸ iWij + g + ¨U dt U j ¸¹ j © i
where P stands for isotropic pressure, U the fluid density, U is the fluid velocity, and Q 0 the viscosity coefficient. For turbulent flow the RANS equation and a sub-particle scale turbulence model are used. Equations (1) and (2) in the sub-grid scale can be written as[8,9]:
1 dU +
dU 1 1 = P + g +Q 0 2 U +
2 3
2 3
j
j
§ W j Wi · ¨¨ 2 + 2 ¸¸ iWij © Ui U j ¹
(7)
where W is the kernel function, expressed usually with the B-spline kernel function
(3)
1 · §2 W R, h = D d ¨ R 2 + R 3 ¸ , 0 d R 1 (8a) 2 ¹ ©3
(4)
3º ª1 W R, h = D d « 2 R » , 1 d R 2 ¬6 ¼
(8b)
W R, h = 0 , R t 2
(8c)
where W stands for the tensor of sub-particles scale, and its various variables are expressed as
W ij = U (2Q t Sij kG ij ) U CI ' 2G ij
¦m
(6)
(5)
in which Ci = 0.0066 , Q t stands for the viscosity coefficient of turbulence, k the turbulent kinetic energy of the sub-grid-scale, Sij the stress tensor of subgrid-scale :
1 § wu wu · Sij = ¨ i + j ¸ 2 ¨© wx j wxi ¸¹ The turbulence viscosity coefficient is calculated by using the standard Smagorinsky model:
Q t = (Cs ') 2 S in which Cs = 0.12 is the Smagorinsky constant, the local strain-rate is expressed in the form of
For two-dimensional problems, D d = 15 / 7ʌh2 , h stands for the smooth length, and R = r / h standing for the relative length. 2.2 Poisson’s equation for giving pressure In the solution process of pressure P , we will divide each time-step into two steps, for the first step, using the source item at the time of t = n to work out the velocity and displacement variables u , r generated by the non-pressure terms, and then making amendment, according to the continuity equation and momentum equation. Thus we can obtain the following:
u = u n + 'u
(9)
'u = f 't
(10)
r = r n + u 't
(11)
For the second step, assume that u n +1 = u + 'u , where 'u is the change of speed generated by the pressure term:
845
'u =
1
U
P n +1't
(12)
To introduce this speed term into the continuity equation, Poisson’s equation for pressure can be derived as follows:
§ 1 · U U <¨ P n +1 ¸ = 2 ©U ¹ U't
(13)
Then, using Poisson’s equation to work out the pressure and inserting it into the above formula, we can obtain the speed generated by pressure term and thus derive u n +1 and r n +1 at the time of t = n 1 .
r n +1 = r n + u n +1't
(14)
Here we use the improved form of the pressure gradient term proposed by Shao et al. [10]:
where the exponents n1 and n2 are assumed to be 4 and 2, D is a constant, taking different values for different problems, r0 refers to the influencing scope of the boundary forces (also known as cut-off length). We assume that the initial spacing of particles is around 2/3. In addition, N layers of virtual particles are pre-arranged outside the computational domain, and the outside layers of virtual particles have the same initial density of particles as the particles within the domain, and their density and speed always remain unchanged, zero speed is used for eliminating the truncation error of the internal particles near the boundary. In order to get accurate solution of the pressure at a certain point on the boundary, an estimation method proposed by Oger[13] is applied. As is shown in Fig.1, the pressure at the Point M can be estimated on the basis of the pressure on the water particle point in the domain around the Point M. S sensor can be seen as the location of pressure sensor in the experiment, and then the pressure at the Point M can be written as
Pij rij <Wij § 1 · 8 <¨ P n +1 ¸ = ¦ m j 2 2 j ©U ¹i Ui + U j rij +K 2 (15) where
Pij = Pi Pj
and
rij = ri rj
are a small
quantity that guarantees the denominator not being zero in the middle steps of calculation, generally taking 0.1 times of the smooth length. The coefficient matrix of this equation is symmetric positive definite and the main diagonal is dominant. Therefore, it is convenient for using a high-efficiency numerical method for solution. The bi-conjugate gradient method[11,12] is employed herein. 2.3 Boundary conditions Virtual particles on the boundary are arranged. When the particles in the computational domain draw near to the boundary, the virtual particles will apply a repulsive force upon these particles so as to ensure that these particles can not penetrate the solid boundary. The boundary force model in the form of the Lennard-Jones force describing molecular interaction is used.
ª§ r PBij = D «¨ 0 «¨© rij ¬ §r PBij = 0 , ¨ 0 ¨r © ij
n1 n2 · § r0 · º xij § r0 · ¸¸ ¨¨ ¸¸ » 2 , ¨¨ ¸¸ d 1 (16a) ¹ © rij ¹ »¼ rij © rij ¹
· ¸¸ ! 1 ¹
(16b)
Fig.1 Sampling area
PM =
³ ³
S sensor
P dS
S sensor
dS
(17)
If S sensor is very small, the pressure PM at the Point M can be approximately given. As is shown in Fig.2, first select all particles (except for the boundary particles) with the distance to the boundary smaller than d (the value of d is related to the smooth length) in the projection domain of S sensor . Each particle can be projected onto the boundary surface S sensor and sequenced in accordance with its coordinate q along the boundary line. Then the surface element dSic can be given by following equation:
dS i c = where
qi '+1 qi '1 2
(18)
qi c is the coordinate of i c . Hence, pressure at
846
the Point M can be approximately obtained through the pressure on a series of element dSic :
PM =
¦ P dS ¦ dS i'
i'
i'
i'
experimental data and VOF results.
(19)
i'
As the particle is projected onto the boundary, there is a small amount of difference between pi and pi c , resulting in a certain amount of error. When the particles are projected along the direction of gravity, pi c should be written as
Pi ' = Pi Ui g ( yi ' yi )
(20)
where y stands for the vertical coordinate.
Fig.2 Particle sampling
3. Numerical results 3.1 Problems of dam collapsing Dam collapse is a classic problem for free-surface flow. In this example it is assumed that the water body height is H = 0.2m , width B = 0.1m , and water flows along the horizontal solid boundary under the action of gravity. Numerical simulation takes 800 fluid particles and 0.005 m as the initial spacing of particles. Full-slip condition is applied at solid wall boundary. Viscosity coefficient is 10-6 kg/ms-1. The time-step is taken as a constant 10-4 s. The kernel function is cubic spline function. Figure 3 shows the snapshots of water particles distribution at the instants t = 0.05s, 0.1s, 0.15s and 0.18s after dam collapsing. The results obtained by the SPH are compared with the results by the VOF method (see Fig.3). These results are basically consistent. In the flow front and wet dry transformation zone on solid boundary, the front does not “raise its head” by detaching from the boundary, which indicates that the pressure on the fluid particles in this region is accurately given. The relation between non-dimensional time T = t ( g / H )1/ 2 and non-dimensional front position X = x / H is shown in Fig.4. The SPH results are consistent with
Fig.3 Snapshots of water particle at the instants t = 0.05s, 0.10s, 0.15s and 0.18s
Fig.4 Relation between the standardization time and front position after the dam collapsing
847
Fig.5 Snapshots of water particle at the time of t = 0.22s, 0.30s
existence of the vertical wall and thus form the slamming role upon the vertical wall, resulting in a peak value of impacting pressure. 3.2 Numerical simulation of breaking process of solitary wave on beaches Wave breaking in near shore has been a hot issue in coastal engineering. The strong nonlinearity of wave breaking arouses the complex movement of water body. As a result, numerical simulation of wave breaking process on beaches has become a difficult problem. Due to the restriction of the grid scale, grid-based approaches have limited accuracy in dealing with multiple free-surfaces in the process of breaking, while the particle method is not subject to the restrictions of the grid-scale, and thus can be used to simulate the problem of multiple free-surfaces of fluid with higher accuracy. In this study, a numerical wave tank model was developed, and used to simulate the wave run-up and breaking process of solitary wave on beaches. 3.2.1 Layout of numerical wave tank The layout of wave tank is shown in Fig.7. The length of the tank is 8.5 m, and the length of horizontal section is 3.5 m. A piston type wave-maker is located at 0.5 m from the left end. At the other end there is a slope with the gradient 1:15. The water depth is 0.2 m. The solitary wave height is 0.08 m. Particle spacing is taken as 0.005 m. The number of the particles in the entire wave tank is 2.5175×104. Three points are selected at the positions 1.0 m, 1.5 m and 2.0 m from the wave-maker to detect the transmission and wave height modulation of solitary wave in the wave tank.
Fig.6 Evolution of impact pressure
If there is a vertical wall downstream of the breaking dam, the dam-collapsing flow will impose impact on the wall, as shown in Fig.5. The vertical wall is assumed to be located at the position x = 0.4m . The water front will reach the wall at the time of t = 0.22s , after the dam breaking. As t = 0.30s , the water front initially runs up along the wall (see Fig.5). The water pressure at the position y = 0.025m on the vertical wall are computed. The pressure evolution is shown in Fig.6. It can be seen from the figure that dam-collapsing flow imposes impact pressure on the wall, which is much larger than static water pressure. The impact pressure is zero before the water front reaches the wall. Horizontal movement of water particles are blocked due to the
Fig.7 Layout of numerical wave tank
3.2.2 Boundary conditions of wave-maker To generate solitary wave, the following displacement[14] is given to the wave maker peddle.
[ (t ) =
T 4H 3H [C (t 0 ) d{1+ tanh 3 2 3d 4d S ([ )]} 2
where d is water depth, H is wave height, C
(21) is
848
the wave velocity:
C = g ( H + d ) , T0 =
2 C
4d 3 § H· ¨ 3.8 + ¸ 3H © d ¹
is acting time of the solitary wave, and S =
into the water body below, and a large number of bubbles being wrapped, and then wave splashes once again, presenting a clear characteristics of wave breaking in the form of a roller.
16 H d 3d
refers to the stroke of wavemaker plate. Figure 8 shows the comparison between the numerical and theoretical results for the solitary wave. It can be seen that they are quite consistent. As the method of iterative Poisson’s equation for the solution of pressure field is applied in this article, the decrease in wave peak does not appear comparing contrast with the standard SPH method by solving the state equation. As id shown in Fig.9, there is almost no wave height attenuation in the process of solitary wave passing through three monitoring points, which further indicates that the iterative Poisson equation for the solution of pressure field can effectively simulate the position of wave peak.
Fig.8 Comparison between the simulated and analytical results
Fig.9 Propagations of solitary wave in numerical flume
Figure 10 shows the comparison of flow fields given by the numerical method and physical model test at typical instants in the course of solitary wave breaking on beaches. It can be seen from the figure that the wave peak is at rather high speed, leading to the vertical shape of wave surface and the wave peak rushing towards in the form of a tongue until falling
Fig.10 Snapshots of water particle and flow field
3.3 Simulation of running-up of cnoidal wave on beaches Ting and Kirby[15] carried out a physical model test of the movement process of cnoidal wave on beaches in a wave tank. The layout of the tank in the experiment is shown in Fig.11. The length of the tank is 24 m with a slope gradient of 1:35 at the right end. The starting point of the slope is at 0.7 m. The water
849
Fig.11 Initial arrangement plan of the wave flume
depth d is 0.4 m. Incident wave height is 0.12 m, and wave period T is 2.0 s. From the experimental data, we can see that the wave breaking point is at 6.4 m. The water depth of the breaking point is 0.199 m. In the numerical simulation of this paper, all physical conditions are the same as the model tests. In the numerical simulation, the particle spacing is 'x = 'y = 0.01m , and the total number of particles is 5.7196×104. To validate the wave breaking effect, six points on the slope are selected to monitor the variation of the water level. The comparison between the numerical results and the test data is given in Fig.12. The change of wave-surface indicates that the wave crest is steep while the wave trough is relatively flat, presenting clear properties of shallow-water wave. Before waves reach the breaking point, the water level rises gradually, and wave peak becomes steeper until the wave breaks.
Fig.13 Distributions of wave crests, wave troughs and mean water levels for the spilling breaker on the slope
Figure 13 shows the changes of wave crest, trough and mean water level along the wave tank. The numerical results agree with the data of the physical model test. The numerical model can describe the whole process of waves running-up on beaches. 4. Conclusions A numerical model for simulating water wave problems has been developed based on the SPH methods. A sub-grid turbulence model has been employed to describe the turbulence characteristics more accurately for wave breaking. Iterative solution of Poisson’s equation for pressure field is used for pressure simulation. Good agreement between the numerical results and physical model test data shows that the model can be used to simulate accurately and efficiently the dam collapse and wave running-up and breaking on beaches. References [1] Fig.12 Free surface elevations at different positions
MONAGHAN J. J., KOS A. Solitary waves on a Cretan beach[J]. J. Waterw. Port Coast. Ocean Eng., 1999,
850
[2] [3] [4]
[5] [6] [7]
[8]
[9]
125(3): 145-154. MONAGHAN J. J., KOS A. and ISSA N. Fluid motion generated by impact[J]. J. Waterw. Port Coast. Ocean Eng., 2004, 129(6): 250-259. LIU G. R., LIU M. B. Smoothed Particle Hydrodynamics: A meshfree particle method[M]. Singapore: World Scientific, 2003. ROGERS B. D., LEDUC J. and MARONGIU J. C. et al. Comparison and evaluation of multi-phase and surface tension models[C]. Proc. 4th SPHERIC Int. Workshop. Nantes, France, 2009, 30-37. SHAO S. D., JI C. and GRAHAM D. I. et al. Simulation of wave overtopping by an incompressible SPH model[J]. Coastal Engineering, 2006, 53(9): 723-735. DALRYMPLE R. A., ROGERS B. D. Numerical modeling of water waves with the SPH method[J]. Coastal Engineering, 2006, 53(9): 141-147. 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. SHAO S. D., JI C. M. SPH computation of pluning waves using a 2-D sub-particle scale turbulence model[J]. International Journal for Numerical Methods in Fluids, 2006, 51(8): 913-936. CUI Yan,WU Wei and GONG Kai et al. Numerical simulation of sloshing in two dimensional rectangular tanks with SPH[J]. Chinese Journal of
Hydrodynamics, 2008, 23(6): 618-624(in Chinese). [10] 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. [11] ZHANG Yong-jie, SUN Qin and LI Jiang-hai. An improved ICCG method for large scale sparse linear equations[J]. Chinese Journal of Computational Physics, 2007, 13(5): 581-584(in Chinese). [12] YUAN Wei-ran, CHEN Pu and LIU Kan-xin. High performance sparse equations with solver for unetrical linear out of core strategies and its application on meshless methods[J]. Applied Mathematics and Mechanics, 2006, 27(10): 1339-1348. [13] OGER G., DORING M. and ALESSANDRINI B. et al. Two-dimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics, 2006, 213(2): 803-822. [14] ZOU Zhi-li, QIU Da-hong and WANG Yong-xue. Numerical simulation of nonliner wave generated in wave flume by VOF technique[J]. Journal of Hydrodynamics, Ser. A, 1996,11(1): 93-103(in Chinese). [15] TING F. C. K., KIRBY J. T. Dynamics of surf-zone in a spilling breaker[J]. Coastal Engineering, 1996, 27(3): 131-160.
2009,21(6):843-850 DOI: 10.1016/S1001-6058(08)60221-8
NUMERICAL SIMULATIONS OF WATER WAVE DYNAMICS BASED ON SPH METHODS* ZHENG Kun, SUN Zhao-chen, SUN Jia-wen State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China, E-mail: [email protected] ZHANG Zhi-ming, YANG Guo-ping, ZHOU Feng CCCC Water Transportation Consultants Co. Ltd, Beijing 100007, China
(Received July 29, 2009, Revised October 29, 2009)
Abstract: A numerical model was established for simulating water wave dynamic problems by adopting the Smoothed Particle Hydrodynamics (SPH) methods of iterative solution of Poisson’s equation for pressure field, and meanwhile the sub-grid turbulence model was applied in the simulation so as to more accurately describe the turbulence characteristics at the time of wave breaking. In this article, simulation of the problem of the dam collapsing verifies the compoting accuracy of this method, and its results can be identical with the results of VOF method and the experimental results by comparison. Numerical simulations of the course of solitary wave and cnoidal wave run-up breaking on beaches were conducted, and the results are basically consistent with experimental results. This indicates that the SPH method is effective for the numerical simulation of the complex problems of water wave dynamics. Key words: Smoothed Particle Hydrodynamics (SPH), Poisson’s equation, water wave dynamics
1. Introduction Wave breaking in near-shore has been a hot issue in coastal engineering. The strong nonlinearity of wave breaking arouses the complex movement of water body. As a result, numerical simulation of wave breaking process on beaches has become a difficult problem. Due to the restriction of the grid scale, in the past, grid-based approach has limited accuracy in dealing with multiple free-surfaces in the process of breaking, while the particle method is not subject to restrictions of the grid-scale and thus can simulate more accurately the problem of multiple free-surfaces of fluid motion on beaches. In this article numerical simulation of the run-up breaking process of solitary wave on beaches is * Project supported by the National High Techology Research and Development Program of China (863 Program, Grant No. 2007AA11Z130) Biography: ZHENG Kun (1977-), Male, Ph. D. Candidate
conducted by using a grid-free method, the Smoothed Particle Hydrodynamics (SPH)[1-4]. The standard SPH method is to achieve the incompressibility of water body by adopting artificial compressing technique. Pressure is taken as a Lagrange factor, and the pressure field is expressed and solved by using the state equation. This method has the advantage of directly giving the relationship between density and pressure. But its computation of the pressure field is still accurate enough. In order to compute the pressure field more accurately, pressure is solved by employing the SPH methods of iterative solution of Poisson’s equation for pressure field, and the sub-particle equation turbulence model is adopted[5-7]. The numerical results indicate that it is more accurate and stable compared with the solution by using state equation.
844
2. Basic equations for SPH 2.1 Governing equations The N-S equation described can be written in the form
S = (2 Sij Sij )1/ 2 , and ' stands for the initial spacing of particles. The final form of Eqs.(3) and (4) can be written as
1 dU +
(1)
dU i = ¦ m j U ij iWij dt j
dU 1 = P + F +Q 0 2U dt U
(2)
§ P P · dU i = ¦ m j ¨ j2 + i 2 ¸ iWij + g + ¨U dt U j ¸¹ j © i
where P stands for isotropic pressure, U the fluid density, U is the fluid velocity, and Q 0 the viscosity coefficient. For turbulent flow the RANS equation and a sub-particle scale turbulence model are used. Equations (1) and (2) in the sub-grid scale can be written as[8,9]:
1 dU +
dU 1 1 = P + g +Q 0 2 U +
2 3
2 3
j
j
§ W j Wi · ¨¨ 2 + 2 ¸¸ iWij © Ui U j ¹
(7)
where W is the kernel function, expressed usually with the B-spline kernel function
(3)
1 · §2 W R, h = D d ¨ R 2 + R 3 ¸ , 0 d R 1 (8a) 2 ¹ ©3
(4)
3º ª1 W R, h = D d « 2 R » , 1 d R 2 ¬6 ¼
(8b)
W R, h = 0 , R t 2
(8c)
where W stands for the tensor of sub-particles scale, and its various variables are expressed as
W ij = U (2Q t Sij kG ij ) U CI ' 2G ij
¦m
(6)
(5)
in which Ci = 0.0066 , Q t stands for the viscosity coefficient of turbulence, k the turbulent kinetic energy of the sub-grid-scale, Sij the stress tensor of subgrid-scale :
1 § wu wu · Sij = ¨ i + j ¸ 2 ¨© wx j wxi ¸¹ The turbulence viscosity coefficient is calculated by using the standard Smagorinsky model:
Q t = (Cs ') 2 S in which Cs = 0.12 is the Smagorinsky constant, the local strain-rate is expressed in the form of
For two-dimensional problems, D d = 15 / 7ʌh2 , h stands for the smooth length, and R = r / h standing for the relative length. 2.2 Poisson’s equation for giving pressure In the solution process of pressure P , we will divide each time-step into two steps, for the first step, using the source item at the time of t = n to work out the velocity and displacement variables u , r generated by the non-pressure terms, and then making amendment, according to the continuity equation and momentum equation. Thus we can obtain the following:
u = u n + 'u
(9)
'u = f 't
(10)
r = r n + u 't
(11)
For the second step, assume that u n +1 = u + 'u , where 'u is the change of speed generated by the pressure term:
845
'u =
1
U
P n +1't
(12)
To introduce this speed term into the continuity equation, Poisson’s equation for pressure can be derived as follows:
§ 1 · U U <¨ P n +1 ¸ = 2 ©U ¹ U't
(13)
Then, using Poisson’s equation to work out the pressure and inserting it into the above formula, we can obtain the speed generated by pressure term and thus derive u n +1 and r n +1 at the time of t = n 1 .
r n +1 = r n + u n +1't
(14)
Here we use the improved form of the pressure gradient term proposed by Shao et al. [10]:
where the exponents n1 and n2 are assumed to be 4 and 2, D is a constant, taking different values for different problems, r0 refers to the influencing scope of the boundary forces (also known as cut-off length). We assume that the initial spacing of particles is around 2/3. In addition, N layers of virtual particles are pre-arranged outside the computational domain, and the outside layers of virtual particles have the same initial density of particles as the particles within the domain, and their density and speed always remain unchanged, zero speed is used for eliminating the truncation error of the internal particles near the boundary. In order to get accurate solution of the pressure at a certain point on the boundary, an estimation method proposed by Oger[13] is applied. As is shown in Fig.1, the pressure at the Point M can be estimated on the basis of the pressure on the water particle point in the domain around the Point M. S sensor can be seen as the location of pressure sensor in the experiment, and then the pressure at the Point M can be written as
Pij rij <Wij § 1 · 8 <¨ P n +1 ¸ = ¦ m j 2 2 j ©U ¹i Ui + U j rij +K 2 (15) where
Pij = Pi Pj
and
rij = ri rj
are a small
quantity that guarantees the denominator not being zero in the middle steps of calculation, generally taking 0.1 times of the smooth length. The coefficient matrix of this equation is symmetric positive definite and the main diagonal is dominant. Therefore, it is convenient for using a high-efficiency numerical method for solution. The bi-conjugate gradient method[11,12] is employed herein. 2.3 Boundary conditions Virtual particles on the boundary are arranged. When the particles in the computational domain draw near to the boundary, the virtual particles will apply a repulsive force upon these particles so as to ensure that these particles can not penetrate the solid boundary. The boundary force model in the form of the Lennard-Jones force describing molecular interaction is used.
ª§ r PBij = D «¨ 0 «¨© rij ¬ §r PBij = 0 , ¨ 0 ¨r © ij
n1 n2 · § r0 · º xij § r0 · ¸¸ ¨¨ ¸¸ » 2 , ¨¨ ¸¸ d 1 (16a) ¹ © rij ¹ »¼ rij © rij ¹
· ¸¸ ! 1 ¹
(16b)
Fig.1 Sampling area
PM =
³ ³
S sensor
P dS
S sensor
dS
(17)
If S sensor is very small, the pressure PM at the Point M can be approximately given. As is shown in Fig.2, first select all particles (except for the boundary particles) with the distance to the boundary smaller than d (the value of d is related to the smooth length) in the projection domain of S sensor . Each particle can be projected onto the boundary surface S sensor and sequenced in accordance with its coordinate q along the boundary line. Then the surface element dSic can be given by following equation:
dS i c = where
qi '+1 qi '1 2
(18)
qi c is the coordinate of i c . Hence, pressure at
846
the Point M can be approximately obtained through the pressure on a series of element dSic :
PM =
¦ P dS ¦ dS i'
i'
i'
i'
experimental data and VOF results.
(19)
i'
As the particle is projected onto the boundary, there is a small amount of difference between pi and pi c , resulting in a certain amount of error. When the particles are projected along the direction of gravity, pi c should be written as
Pi ' = Pi Ui g ( yi ' yi )
(20)
where y stands for the vertical coordinate.
Fig.2 Particle sampling
3. Numerical results 3.1 Problems of dam collapsing Dam collapse is a classic problem for free-surface flow. In this example it is assumed that the water body height is H = 0.2m , width B = 0.1m , and water flows along the horizontal solid boundary under the action of gravity. Numerical simulation takes 800 fluid particles and 0.005 m as the initial spacing of particles. Full-slip condition is applied at solid wall boundary. Viscosity coefficient is 10-6 kg/ms-1. The time-step is taken as a constant 10-4 s. The kernel function is cubic spline function. Figure 3 shows the snapshots of water particles distribution at the instants t = 0.05s, 0.1s, 0.15s and 0.18s after dam collapsing. The results obtained by the SPH are compared with the results by the VOF method (see Fig.3). These results are basically consistent. In the flow front and wet dry transformation zone on solid boundary, the front does not “raise its head” by detaching from the boundary, which indicates that the pressure on the fluid particles in this region is accurately given. The relation between non-dimensional time T = t ( g / H )1/ 2 and non-dimensional front position X = x / H is shown in Fig.4. The SPH results are consistent with
Fig.3 Snapshots of water particle at the instants t = 0.05s, 0.10s, 0.15s and 0.18s
Fig.4 Relation between the standardization time and front position after the dam collapsing
847
Fig.5 Snapshots of water particle at the time of t = 0.22s, 0.30s
existence of the vertical wall and thus form the slamming role upon the vertical wall, resulting in a peak value of impacting pressure. 3.2 Numerical simulation of breaking process of solitary wave on beaches Wave breaking in near shore has been a hot issue in coastal engineering. The strong nonlinearity of wave breaking arouses the complex movement of water body. As a result, numerical simulation of wave breaking process on beaches has become a difficult problem. Due to the restriction of the grid scale, grid-based approaches have limited accuracy in dealing with multiple free-surfaces in the process of breaking, while the particle method is not subject to the restrictions of the grid-scale, and thus can be used to simulate the problem of multiple free-surfaces of fluid with higher accuracy. In this study, a numerical wave tank model was developed, and used to simulate the wave run-up and breaking process of solitary wave on beaches. 3.2.1 Layout of numerical wave tank The layout of wave tank is shown in Fig.7. The length of the tank is 8.5 m, and the length of horizontal section is 3.5 m. A piston type wave-maker is located at 0.5 m from the left end. At the other end there is a slope with the gradient 1:15. The water depth is 0.2 m. The solitary wave height is 0.08 m. Particle spacing is taken as 0.005 m. The number of the particles in the entire wave tank is 2.5175×104. Three points are selected at the positions 1.0 m, 1.5 m and 2.0 m from the wave-maker to detect the transmission and wave height modulation of solitary wave in the wave tank.
Fig.6 Evolution of impact pressure
If there is a vertical wall downstream of the breaking dam, the dam-collapsing flow will impose impact on the wall, as shown in Fig.5. The vertical wall is assumed to be located at the position x = 0.4m . The water front will reach the wall at the time of t = 0.22s , after the dam breaking. As t = 0.30s , the water front initially runs up along the wall (see Fig.5). The water pressure at the position y = 0.025m on the vertical wall are computed. The pressure evolution is shown in Fig.6. It can be seen from the figure that dam-collapsing flow imposes impact pressure on the wall, which is much larger than static water pressure. The impact pressure is zero before the water front reaches the wall. Horizontal movement of water particles are blocked due to the
Fig.7 Layout of numerical wave tank
3.2.2 Boundary conditions of wave-maker To generate solitary wave, the following displacement[14] is given to the wave maker peddle.
[ (t ) =
T 4H 3H [C (t 0 ) d{1+ tanh 3 2 3d 4d S ([ )]} 2
where d is water depth, H is wave height, C
(21) is
848
the wave velocity:
C = g ( H + d ) , T0 =
2 C
4d 3 § H· ¨ 3.8 + ¸ 3H © d ¹
is acting time of the solitary wave, and S =
into the water body below, and a large number of bubbles being wrapped, and then wave splashes once again, presenting a clear characteristics of wave breaking in the form of a roller.
16 H d 3d
refers to the stroke of wavemaker plate. Figure 8 shows the comparison between the numerical and theoretical results for the solitary wave. It can be seen that they are quite consistent. As the method of iterative Poisson’s equation for the solution of pressure field is applied in this article, the decrease in wave peak does not appear comparing contrast with the standard SPH method by solving the state equation. As id shown in Fig.9, there is almost no wave height attenuation in the process of solitary wave passing through three monitoring points, which further indicates that the iterative Poisson equation for the solution of pressure field can effectively simulate the position of wave peak.
Fig.8 Comparison between the simulated and analytical results
Fig.9 Propagations of solitary wave in numerical flume
Figure 10 shows the comparison of flow fields given by the numerical method and physical model test at typical instants in the course of solitary wave breaking on beaches. It can be seen from the figure that the wave peak is at rather high speed, leading to the vertical shape of wave surface and the wave peak rushing towards in the form of a tongue until falling
Fig.10 Snapshots of water particle and flow field
3.3 Simulation of running-up of cnoidal wave on beaches Ting and Kirby[15] carried out a physical model test of the movement process of cnoidal wave on beaches in a wave tank. The layout of the tank in the experiment is shown in Fig.11. The length of the tank is 24 m with a slope gradient of 1:35 at the right end. The starting point of the slope is at 0.7 m. The water
849
Fig.11 Initial arrangement plan of the wave flume
depth d is 0.4 m. Incident wave height is 0.12 m, and wave period T is 2.0 s. From the experimental data, we can see that the wave breaking point is at 6.4 m. The water depth of the breaking point is 0.199 m. In the numerical simulation of this paper, all physical conditions are the same as the model tests. In the numerical simulation, the particle spacing is 'x = 'y = 0.01m , and the total number of particles is 5.7196×104. To validate the wave breaking effect, six points on the slope are selected to monitor the variation of the water level. The comparison between the numerical results and the test data is given in Fig.12. The change of wave-surface indicates that the wave crest is steep while the wave trough is relatively flat, presenting clear properties of shallow-water wave. Before waves reach the breaking point, the water level rises gradually, and wave peak becomes steeper until the wave breaks.
Fig.13 Distributions of wave crests, wave troughs and mean water levels for the spilling breaker on the slope
Figure 13 shows the changes of wave crest, trough and mean water level along the wave tank. The numerical results agree with the data of the physical model test. The numerical model can describe the whole process of waves running-up on beaches. 4. Conclusions A numerical model for simulating water wave problems has been developed based on the SPH methods. A sub-grid turbulence model has been employed to describe the turbulence characteristics more accurately for wave breaking. Iterative solution of Poisson’s equation for pressure field is used for pressure simulation. Good agreement between the numerical results and physical model test data shows that the model can be used to simulate accurately and efficiently the dam collapse and wave running-up and breaking on beaches. References [1] Fig.12 Free surface elevations at different positions
MONAGHAN J. J., KOS A. Solitary waves on a Cretan beach[J]. J. Waterw. Port Coast. Ocean Eng., 1999,
850
[2] [3] [4]
[5] [6] [7]
[8]
[9]
125(3): 145-154. MONAGHAN J. J., KOS A. and ISSA N. Fluid motion generated by impact[J]. J. Waterw. Port Coast. Ocean Eng., 2004, 129(6): 250-259. LIU G. R., LIU M. B. Smoothed Particle Hydrodynamics: A meshfree particle method[M]. Singapore: World Scientific, 2003. ROGERS B. D., LEDUC J. and MARONGIU J. C. et al. Comparison and evaluation of multi-phase and surface tension models[C]. Proc. 4th SPHERIC Int. Workshop. Nantes, France, 2009, 30-37. SHAO S. D., JI C. and GRAHAM D. I. et al. Simulation of wave overtopping by an incompressible SPH model[J]. Coastal Engineering, 2006, 53(9): 723-735. DALRYMPLE R. A., ROGERS B. D. Numerical modeling of water waves with the SPH method[J]. Coastal Engineering, 2006, 53(9): 141-147. 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. SHAO S. D., JI C. M. SPH computation of pluning waves using a 2-D sub-particle scale turbulence model[J]. International Journal for Numerical Methods in Fluids, 2006, 51(8): 913-936. CUI Yan,WU Wei and GONG Kai et al. Numerical simulation of sloshing in two dimensional rectangular tanks with SPH[J]. Chinese Journal of
Hydrodynamics, 2008, 23(6): 618-624(in Chinese). [10] 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. [11] ZHANG Yong-jie, SUN Qin and LI Jiang-hai. An improved ICCG method for large scale sparse linear equations[J]. Chinese Journal of Computational Physics, 2007, 13(5): 581-584(in Chinese). [12] YUAN Wei-ran, CHEN Pu and LIU Kan-xin. High performance sparse equations with solver for unetrical linear out of core strategies and its application on meshless methods[J]. Applied Mathematics and Mechanics, 2006, 27(10): 1339-1348. [13] OGER G., DORING M. and ALESSANDRINI B. et al. Two-dimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics, 2006, 213(2): 803-822. [14] ZOU Zhi-li, QIU Da-hong and WANG Yong-xue. Numerical simulation of nonliner wave generated in wave flume by VOF technique[J]. Journal of Hydrodynamics, Ser. A, 1996,11(1): 93-103(in Chinese). [15] TING F. C. K., KIRBY J. T. Dynamics of surf-zone in a spilling breaker[J]. Coastal Engineering, 1996, 27(3): 131-160.