The smoothed particles hydrodynamics (SPH) technique has been used to treat some problems occurred in the free surface flow particularly dam-breaking flow using the open-source SPH code through the NAG FORTRAN compiler and the result presentation was done by using the Para view visualization software. Several methods are employed in order to solve the experienced problems. The most important features have been achieved from this work are: the density filtering method is used to treat pressure fluctuation and the operation time has decreased. Moreover, changing the boundary condition has provided two important things; preventing the particles from penetration the walls and also decreasing the operation time even more. Furthermore, changing the smoothing length from the optimum case (h=1.3 Δx) to (h=1.5 Δx), it shows a few pressure fluctuations and an increase in operation time with no wave propagation. The second change is from (h=1.3 Δx) to (h=1.0 Δx) the pressure noise rose a little bit further but the operation time decreased. The third change is from (h=1.3 Δx) to (h=0.5 Δx) which led to increasing the pressure oscillations and the solution accuracy has been decreased. In spite of decreasing the influence radius, the operation time has increased. Finally, there is no difference has been achieved from changing the pressure formulation in the shape of frames except the operation time has reduced.

Keywords: SPH, the dam-breaking flow, treatments

The usage of SPH in modeling the fluid motion has been increased gradually. This is owing to several advantages over classical methods of grids includes; the mass conservation is guaranteed without paying further computations since the mass is presented by particles themselves. In spite of originally developing the SPH for astrophysical applications, it had been employed in a successful path in a wide range of engineering applications.^1,2 Moreover, the SPH has an ability to compute pressure from weighted contributions of neighboring particles instead of solving linear equation systems. Furthermore, SPH is able to produce a free surface for two-phase fluids where the denser fluid (usually water) can be presented by the particles and the lighter fluid (usually air) can be presented by the empty space. These reasons afforded the SPH an opportunity to simulate fluid motion in real time, in the next sections, some SPH features are discussed in details like so.

Starting from a mesh free method: the stable and accurate numerical solutions can be achieved for the partial differential equations or integral equations with respect to all types of boundary conditions with a set of particles which are distributed arbitrarily without employing any mesh which joins these particles. Basically, the mesh free methods have been divided into three main types; methods based on weak form formulation, methods based on strong form formulation and particles methods. However, the major disadvantage of weak form is that it is said not to be truly mesh free and the main disadvantages of strong form are unstable and inaccurate. The SPH (particles method) has provided a stable solution with a large deformation for many different problems and its accuracy is dependent so much on the selection of smoothing function.³ Lagrangian grid-based method: in this method, the grid has been attached on or fixed to the materials in the computation process and thus it flows with the movement of the materials. It has several advantage includes; no convection term in relevant partial differential equations and therefore, its code is simpler and faster. Moreover, complicated geometries can easily be modeled. Furthermore, the data of time history can be easily obtained at any point.³ In the present work, the structure of this experiment is divided into several parts. Firstly, description a summary of the SPH method. Secondly, presenting the fluid mechanics for the basic simulation. Thirdly, using the density filter and define the main change. After that, change the boundary condition in order to obtain a better accuracy in solutions. Furthermore, the smoothing length is changed and presenting how this change influence on the numerical solution and the computation time. Eventually, the pressure gradient formulation is changed in order to investigate how the different formulas can agree or not with the fluid limits and the SPH cods. The objectives of the work are to gain familiarity with the behavior of the SPH by changing several variables in the SPH code. Furthermore, understand the flow physics for the dam-breaking under some different effects.

Smoothed particle hydrodynamics (SPH) can be introduced as a computational technique based on the mesh free, lagrangian methods used for simulating fluid flows. It was firstly used for solving the astrophysical problems in 3D and then developed and employed in a variety of engineering applications such as elastic flows, multiphase flow, flow through porous media etc.³ The fluid in SPH method is divided into a set of separate elements (called particles) which have a spatial distance called smoothing length (h). The heart of SPH is so-called smoothing kernel function (W) where its summation interpolation is used to estimate the density and then calculate the rest of SPH equations.

The smoothing particle hydrodynamics (SPH) approach has been basically in order to solve the problems related with the hydrodynamic problems, where the conservation equations of variables (density, velocity, and momentum) are given in forms of the partial differential equations. However, two very important steps should be achieved in order to obtain the SPH formulation which are; Kernel approximation and particle approximation. The first approximation uses the derivatives of a smoothing function so as to show the derivatives in continued representation as an integral illustration. While the computational domain in a particle approximation can be discretized with some particles. The field function A(x) can be given by the formulas below:

$\langle A(x)\rangle \approx {\displaystyle \sum _{j=1}^{N}A(x_j})W(x_i,h)\frac{m_j}{{\rho }_j}$ …………. (1a)

The derivative of the field function can be written by:

$\langle \nabla A(x)\rangle \approx {\displaystyle \sum _{j=1}^{N}A(x_j})(x_i)\nabla W_{ij}\frac{m_j}{{\rho }_j}$ ………… (1b)

According to the two above equations, the governing equations of fluid flow and whole physics can be derived, therefore; the governing equations of mass and momentum can be respectively written by:

$\frac{d{\rho }_i}{dt}={\displaystyle \sum _{j=1}^N{m}_j}{v}_{ij}.{\nabla }_i{W}_{ij}$ ……………. (1c)

$\frac{dv}{dt}\sum _{j=1}^N\left(\frac{p_i}{{\rho }_i{}^2}+\frac{p_j}{{\rho }_j{}^2}{V}_i{W}_{ij}\right)$ ………………… (1d)

As can be shown in the Figure 1, the case is presented by considering the 2D flow produced after the dam-breaking and the wall has been positioned vertically at a distance from the broken dam and the flow is following until hitting by the vertical wall. After that the flow has overturned backwards onto the underlying fluid. The frame (at time=0) shows the water block has a large pressure at the left bottom corner owing to water column height and then (at time= 20 sec) the water flowed freely downstream due to the pressure drop. It is clear that there so many fluctuations in the predicting pressure in the rest frames which progressively destroy the pressure field and; therefore, it is physically unacceptable. This problem has been caused by small errors in the density calculation which might be amplified through the equation of state which links pressure and density. Effectively what happened in this circumstance is the number of particle neighbors seen by a boundary particle has decreased. Therefore, the density calculated would be less than the original value and; as a result, the equation of state has estimated an incorrect pressure distribution. The negative pressure prediction such fluctuations can possibly lead to tension-instability phenomena. This might be reason behind stopping the computation and generally producing the unphysical fragmentation of the water jets during the later overturning and also causes numerical instability.⁴ In order to alleviate this problem, an effective treatment is made below to correct the density field which strongly affects the pressure field and; consequently, the capacity of SPH simulation would be improved considerably.

When the dynamics are in general factual from SPH simulation, the pressure field corresponding to the particles presents has been very large pressure fluctuations. One of the most effective ways can be employed in order to overcome this problem is density re-initialization (density filtering) as shown in Figure 2. It is not expensive computationally and used widely which aim to process a filter over the density of the particles and the re-assign a density to each particle.⁴ Two types of correction approaches are used (zeroth order and first order) and both of them can be applied at every 30 time step. The zeroth order is more simple correction to the density filed and quick approach. However, the first order correction is able to reproduce the exact linear variation of density. In general, several benefits are achieved from applying the filtering to the density field (as shown) include; the pressure distribution gain good regularity, the results are more believable, with using of artificial viscosity, the energy has been conserved better, and uniformity between density, mass and occupied area has been restored and the flow representation of a dam-break is good.

It is clear that the density filtering afforded the dam-break flow an opportunity to gain a good representation and a pressure field without fluctuations and also the operation time has been decreased from (192.58 sec) in case of no density filtering to (185.9 sec) in the case of density filtering. On the other hand, it left a bad effect in which more complexity has been experienced in terms of wave propagation. Previous investigations pointed out that after comparing the numerical results for wave impact on a vertical breakwater with experimental measurements, it was shown that the density filter allows achieving a stable pressure field without oscillations. Nevertheless, a numerical diffusion on wave propagation and the consequent reduction in wave height have been obtained that is not seen in the experimental measurements.

In spite of the validation of SPH formulations for all inner particles, it is no guarantee to be accurate for the particles near the boundaries because the boundaries have truncated the distribution of Kernel function. Consequently, the SPH approach has experienced a problem in applying the boundary conditions and this would lead to the fact that the SPH approximation is not able to create the accuracy in the second order. Therefore, the appropriate and right selection of boundary treatments make a continuing concern in order to obtain the success and accuracy in the solution of the engineering problems with bounded fields using SPH technique. The incorrect choice of boundary condition would give two crucial results; firstly, the fluid particles will have an ability to penetrate the walls and then leave the bounded domain. Secondly, the solution will have some errors due to the kernel truncation at the boundary. This is basically because of the Lagrangian nature of SPH particles, the performing boundary conditions at fixed spatial boundaries has been less straightforward than in mesh-based methods.⁵ This phenomena is very obvious in Figure 2 where particles penetrate the wall; however, what happened in the Figure 3 is the boundary condition is changed from the dynamic boundary condition to the repulsive force boundary condition which acts to create a highly repulsive force to the particles close to the boundary and; therefore, prevents these particles from penetrating the wall which is physically unacceptable. The changing of this boundary provided two important advantages; first the penetration of particles through the wall has been prevented due to the increasing of the particles close to the wall. The second advantage of this changing is that the operation time has decreased even more from the previous case to (180.02 sec) which is less than that of the preceding case since in the previous case the computations did not take the wall particles into account and that means the decreasing in the number of particles led to decrease the operation time. All the preceding cases have been performed at optimizing smoothing length (h/Δx) = 1.3. Nevertheless, other three cases at different smoothing lengths (as shown in figures 4, 5, 6) have been examined in order to study the effect of changing the smoothing length on the performance of numerical resolution and in general on the performance of SPH approximation.

The smoothing length plays a key role in the SPH method where it has a significant effect on the results accuracy and the computational efficiency. Thus, if the smoothing length has been very small, that means the number of particles is not sufficient in the support domain of dimension kh and; therefore, the accuracy would be very low. On the other hand, if the smoothing length is very large, that means all particle details or local properties might be smoothed out and; consequently, the accuracy also would be influenced. However, the number of neighboring particles should be 5 (including the particle itself) in 1D and 21, 57 in 2D and 3D respectively with smoothing length about 1.2 times k=2 (for cubic spline kernel function) in order to make SPH works well.The smoothing length is strongly dependent on the density of a particle and from the last the smoothing length can easily be evolved according to this law:

$h=h_o{\left(\frac{{\rho }_o}{\rho }\right)}^{1/d}$ ……….(2)

Where, ho and ρo are the initial smoothing length and initial density respectively and d represents the dimension number.³

It is very clear in Figure 4 (case h/Δx =1.5) with large smoothing length, the frame shows some pressure fluctuations return to the field where the flow has overturned backwards onto the underlying fluid. This occurred due to some sudden change in the local density which is effectively influenced by the changing the smooth length which generally leading to decrease the solution accuracy with some increasing in time which is (198.04 se) owing to an increasing the radius of influence of the kernel function leading to more neighbor particles need to be computed. However, in this frame, an advantage can be noticed, which not exists in the other two frames, there is no wave propagation in this frame.

In the frame of case (h/Δx=1.0), see Figure 5, the smoothing length has been decreased. In spite of the presence of pressure fluctuations which has increased a little bit from what its value was at smoothing length 1.3, the reduction in the radius of influence seems to be more effective than the presence of pressure noise which is leading to decrease the operation time to (154.27 sec). However, comparing with the frame of (h/Δx =1.5), the amount of pressure fluctuations in smoothing length of 1.0 is less than those of smoothing length due to sharp change in the local density in the case of 1.5. On the other hand, the case of 1.0 smoothing length still suffers from the wave propagation.⁶

In the frame of case (h/Δx=0.5), shown in Figure 6, the smoothing length has been further decreased which led to changing the kernel function and also the number of particles is not sufficient in the support domain and; therefore, the accuracy is very low which is represented in giving the pressure fluctuations in some locations which are physically rejected. Basically, owing to the smoothed nature of SPH approach, that generates a problem which is sharp change in density has been smoothed or blurred, out over a small distance due to the smoothing length which provides pressure oscillations to the solution and in general has affected the results accuracy and the computational efficiency, the less particles the more the effects could be noticed. Several disadvantages have been noted in this case; first, the pressure fluctuations have been increased which would leads to the tension-instability. Secondly, a numerical diffusion on wave propagation and the consequent reduction in wave height have been obtained. Thirdly, the operation time has reached to the highest value in this study which is (233.17 sec).⁷

The changing of smoothing length means the numerical resolution has been changed and the optimizing selection of smoothing length in SPH technique provides effectively several merits to the solution including; high accuracy, less computation time and consistent contours with those of the exact results. To sum up, all smoothing lengths which have been examined (1.5, 1.0 and 0.5) provided many disadvantages to the solution and; therefore, the optimum smoothing length of 1.3 present’s better performance in terms of accuracy and computation efficiency.

The conservation momentum equation in Lagrangian form is:

$\frac{dv}{dt}=-\frac{1}{\rho }\nabla P$…… (3)

The standard SPH expression of pressure gradient can be given by:

$\langle \frac{1}{{\rho }^{(x)}}\frac{\partial \rho (x)}{\partial x}\rangle =\sum _{j=1}\frac{m_j}{{\rho }_j}\frac{(p_j-p_i)}{{\rho }_i}\frac{(x-x_j)}{\partial x_i}$ ………… (4)

In order to conserve linear and angular momentum, it is necessary to ensure that the discretized formulation of interaction forces between two particles is able to conserve the local linear and angular momentum. However, the formulation of momentum in equation (4) shows a very crucial feature. In the circumstance where the pressure field has been spatially constant, that would give us zero pressure of the particles in any arrangement of a system (because of the minus in the numerator between the pressures) and it is not able to provide an equal and opposite reaction between each two particles. Therefore, this formulation will never conserve neither linear nor angular momentum. When it has been applied, an error message we have obtained from the FORTRAN program which tells us that this formulation has a part which is physically unacceptable and; moreover, from the numerical procedure point of view, the solution would not be converged.

By using the Calculus of Variations afforded the equation (3) an opportunity to be derived in a very rigorous derivation, and the new formulation can be expressed by:

$\langle \frac{1}{{\rho }^{(x)}}\frac{\partial \rho (x)}{\partial x}\rangle =\sum _{j=1}(p_j+p_i)\frac{m_j}{{\rho }_j}\frac{\partial W(x-x_j)}{\partial x_i}$………..(5)

According to equation (5), it is possible to conserve linear and angular momentum since there is an equal and opposite reaction between each two particles and also if the pressure field has been spatially constant, that would not give us a zero pressure of the particles in any configuration of a system. As shown in Figure 7, the equation (5) has been used for the pressure gradient with changing the boundary condition to repulsive force boundary condition. If a comparison has been made with the results in Figure 3 which depends on the equation (6) for the pressure gradient, obviously there is no difference between the results in all frames even the small features have been noticed in both of them such as the wave propagation in the last frame. However, the operation time for applying equation (5) has increased to (183.17 sec) while it was (180.02 sec) for applying the equation (6) for pressure gradient where all previous cases have been performed using the basic formulation which has conserved linear and angular momentum.

$\langle \frac{1}{{\rho }^{(x)}}\frac{\partial \rho (x)}{\partial x}\rangle ={\sum _{j=1}}^m\left(\frac{p_j}{{\rho }^2}+\frac{p_i}{{\rho }_i{}^2}\right)\frac{\partial W(x-x_j)}{\partial x_i}$…………. (6)

All operation times are listed in Table 1.

Figure 1 The SPH simulation of dam-breaking without any treatment.

Figure 2 The SPH simulation of dam-breaking with density filtering in different times.

Figure 3 The SPH simulation of dam-breaking with changing the boundary condition.

Figure 4 The SPH simulation of dam-breaking by changing the smoothing length ((h/Δx)=1.5).

Figure 5 The SPH simulation of dam-breaking by changing the smoothing length ((h/Δx)=1.0).

Figure 6 The SPH simulation of dam-breaking with changing the smoothing length ((h/Δx)= 0.5).

Figure 7 The SPH simulation of dam-breaking with changing the pressure gradient formulation.

The SPH method is an attractive technique especially for the applications where the density is the field variable. This method has been employed treat some problems experienced in the 2D dam-breaking flow. The implementation of density filtering has removed the density oscillations which have influenced the pressure field as well as surface profile. Therefore, Filtering approach provides more accurate solution than unfiltered approach but with more complexes for wave propagation phenomena. After that, the repulsive force was employed as a boundary condition which provided two advantages; preventing the particles from penetration through the wall and reducing the operation time. Then, the optimal smoothing length has been changed from 1.3 to (0.5 then 1.0 then 1.5). However, all the other smoothing lengths (0.5, 1.0 and 1.5) provided different disadvantages (pressure fluctuations and increasing operation time etc.…) and in general the accuracy is influenced and the computation efficiency is badly affected. Furthermore, the pressure gradient has been produced by three different expressions, two of them are able to conserve linear and angular momentum and there is no difference has been noticed in their frames shapes. On the other hand, the third one basically has an unphysical term which led to present an error message from the FORTRAN programme and; therefore, it is not able to conserve momentum.

The author declares that there is no conflict of interest.

1. Lucy L. Numerical approach to the testing of the fission hypothesis. Astron J. 1977;82:12:1013–1024.
2. Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society. 1977;181:375–389.
3. Liu GR, Liu MB. Smoothed particle hydrodynamics a mesh free particle method. Singapore: World Scientific; 2003.
4. Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. Journal of Computational Physics. 2003;191(2):448–475.
5. Hosseini SM, Feng JJ. Pressure boundary conditions for computing incompressible flows with SPH. Journal of Computational Physics. 2011;230(19):7473–7487.
6. Liu M, Shao J, Chang J. On the treatment of solid boundary in smoothed particle hydrodynamics. Science China Technological Sciences. 2012;55(15):244–254.
7. Didier E, Neves DRCB, Martins R, et al. (2012). Modelling OF hydrodynamics around an impermeable breakwater: comparison between physical AND SPH numerical modelling. Engenharia Térmica. 2012;11(1–2):68–76.