Modelling Dam Break Evolution over a Wet Bed with Smoothed Particle Hydrodynamics : A Parameter Study

When investigating water flow in spillways and energy dissipation, it is important to know the behavior of the free surfaces. To capture the real dynamic behavior of the free surfaces is therefore crucial when performing simulations. Today, there is a lack in the possibility to model such phenomenon with traditional methods. Hence, this work focuses on a parameter study for one alternative simulation tool available, namely the meshfree, Lagrangian particle method Smoothed Particle Hydrodynamics (SPH). The parameter study includes the choice of equation-of-state (EOS), the artificial viscosity constants, using a dynamic versus a static smoothing length, SPH particle spatial resolution and the finite element method (FEM) mesh scaling of the boundaries. The two dimensional SPHERIC Benchmark test case of dam break evolution over a wet bed was used for comparison and validation. The numerical results generally showed a tendency of the wave front to be ahead of the experimental results, i.e. to have a greater wave front velocity. The choice of EOS, FEM mesh scaling as well as using a dynamic or a static smoothing length showed little or no significant effect on the outcome, though the SPH particle resolution and the choice of artificial viscosity constants had a major impact. A high particle resolution increased the number of flow features resolved for both choices of artificial viscosity constants, but at the expense of increasing the mean error. Furthermore, setting the artificial viscosity constants equal to unity for the coarser cases resulted in a highly viscous and unphysical solution, and thus the relation between the artificial viscosity constants and the particle resolution and its impact on the behavior of the fluid needed to be further investigated.


Introduction
In Sweden, hydropower plays a significant role in the supply of energy and it generates roughly 45% (66.0 TWh in 2011) of the total electricity power production [1].The first larger hydropower plant in Sweden commenced its operation for almost 100 years ago and the latest was put into service in the 1970s.Hence, most of the hydropower plants in Sweden are old, which has initiated a number of projects on refurbishments.Due to the deregulated energy market, hydropower plants are more and more run at off design conditions.Consequently, to have an optimal hydropower that, for instance, can run at the best efficiency point when the winnings are high, the plants must be modernized and adapted to the energy market.Other key aspects are new requirements on renewable and "green" energy sources.This makes the prospect of redesigning old hydropower stations even higher.Finally, the demand on hydropower as a short term regulating source is increasing to take care of more frequent and larger fluctuations in the energy system, mostly due to the on-going expansion of wind power.Regarding the Swedish hydropower sector as a whole focus is set on the service and administration as well as refurbishments of existing plants, to meet the new demands.To assess hydraulic properties in these projects, physical scale models are often used.Performance, accuracy and reliability of mathematical models are in these cases considered to be too poor.Hydropower hydraulics is characterized by large scales and flow rates.The geometry is usually partly formed by nature, i.e. to some extent chaotic at a scale of roughness.Furthermore, highly aerated and disturbed free surface flows are frequently encountered, for instance, in spillway channel flows and in energy dissipaters.This constitutes significant mathematical simulation challenges.Thus, there is a need of computationally robust and reliable numerical methods to handle these fundamentally complex problems.Modelling of highly disturbed aerated free surface flows is complex when grid based method is used [2].Severe problems with mesh entanglement and determination of the free surface have been encountered [3].The meshfree, Lagrangian particle based method Smooth Particle Hydrodynamic (SPH), has shown to be a good alternative to grid based methods to overcome such problems [4] [5].The technique was first proposed independently by [6] [7] in the late seventies to solve astrophysical problems in a three-dimensional open space.Today, the SPH method has been applied to a number of fields and problems [5] and the maturity of the method has increased significantly.Thus, SPH was chosen as the computational method in the present work.In the SPH method, the fluid domain is represented by a set of non-connected particles which possess individual material properties such as density, velocity and pressure [8].Besides representing the problem domain and acting as information carriers, the particles also act as the computational frame for the field function approximations.As the particles move with the fluid, their material properties change as a function of time and spatial co-ordinate due to interactions with neighboring particles.A major advantage of SPH is that the method is mesh-free, and thus considerable time is saved as compared with methods that need a predefined mesh.However, as the SPH is relatively unexplored as compared with traditional finite volume method (FVM) and finite element method (FEM) methods, there are some areas that still need considerable attention.For instance, wall boundary conditions are generally difficult to set in SPH as they do not appear in a natural way within the SPH formalism.Furthermore, the SPH method is typically slower computational wise since the time step depends on the speed of sound and the explicit integration techniques used.The idealized dam break problem, where a given volume of initially stationary water is released onto a dry channel bed by a sudden removal of a gate, is a well-studied problem (e.g.[9]- [11]).After the gate has been removed, a negative surge propagates upstream and a dam break wave moves rapidly downstream.If, however, the dam break wave propagates over a wet bed (an initial fluid layer at rest), the flow field will become considerably different, see for instance [12] and [13].Reference [13] studied the influence of adding a polymer in order to reduce the turbulence drag using a dam break setup.However, clean water test was conducted as well and the experimental results had been used by many authors as a validation test case such as [14]- [16].The SPH European Research Interest Community (SPHERIC) has adopted the work by [13], and digitized wave profiles, wave velocity data and snapshots are available at wiki.manchester.ac.uk/spheric.The aim of this paper is to study the effects of several test parameters when applied to the SPHERIC test case based on the work by [13] of dam break evolution over a wet bed.The test parameters are the equation-of-state, the artificial viscosity constants, using a dynamic or a static smoothing length, the SPH particle resolution and the FEM mesh scaling as thin finite shell elements are used as solid wall boundaries.

Governing Equations
The superscripts α and β are used to denote coordinate directions and the subscripts i and j denotes particle indices.The continuity and momentum equations can be written in the SPH formalism according to, where i ρ is the density, m j is the particle mass, g is the acceleration due to gravity and x i and ij i j = − v v v is the position and velocity vector respectively. .
In Equation (3), dim is the number of space dimensions and is the derivative of the function used to define the kernel function W, i.e.
( ) ( ) and h is the smoothing length.The commonly used cubic B-spline kernel is obtained by choosing θ as, ( ) ( ) where C is a constant of normalization defined for one-, two-and three-dimensional spaces as: The NULL material model will be applied that defines the deviatoric viscous stress as [17], , , where , α β ε is deviatoric strain rate.By setting the dynamics viscosity µ to zero the total stress tensor reduc- es to the isotropic pressure which is defined by an equation-of-state (EOS).In literature, several EOSs can be found but in the present work the Gruneisen, Linear Polynomial and the Morris EOS are considered only.The Gruneisen EOS employs the cubic shock velocity-particle velocity and defines the pressure for compressed materials as [17] ( ) ( ) ( ) and for expanded materials as, ( ) where S 1 , S 2 and S 3 are coefficients of the slope of the u s -u p curve, u s and u p are the chock and particle velocity respectively.The variable c 0 is the intercept of the curve corresponding to the adiabatic speed of sound, γ 0 is the Gruneisen gamma, α is the first volume correction to γ 0 , E 0 is the initial internal energy and ρ is the initial density.Properties and constants for water given by [18] are applied, see Table 1.The Linear Polynomial EOS is defined as [17], ( ) .
References [19] and [20] used the following expressions for the constants in Equation ( 10), where c 0 in Equation ( 11) is the adibatic speed of sound and , where 0 ρ is the initial density.
Reference [20] used k = 2 and 5 . The Morris EOS [21] defines the pressure as, ( ) which can be obtained by setting the constants 1 in the Linear Polynomial EOS.The artificial viscosity proposed by [22] is implemented as, q 1 and q 2 are constants and ij i j = − r r r and ij i j = − v v v are the distance and velocity vector, respectively.Over bars denote average values.If the smoothing length h is dynamic an average value is used, i.e.

(
) is used to avoid the denominator to go to zero.According to [22], the linear term associated with q 1 in the artificial viscosity expression produces both bulk and shear viscosity.The second quadratic term is incorporated to handle high Mach number shocks and is intended to suppress particle penetrations [22] and [8].Furthermore, the values of the constants q 1 and q 2 are according to [14] problem dependent.In literature, many authors suggest setting q 2 to zero and q 1 between 0.01 and 0.1 [23] and [24].However, [8] states that the values of q 1 and q 2 are usually set around unity.The impact of the choice of the constants will be another parameter to be studied in the present work.With the reduction of the total stress tensor, the inclusion of the artificial viscosity term together with the symmetric properties of the ( ) term, the momentum equation (Equation ( 2)) reduces to the familiar expression, ( ) A first order time integration scheme is applied and the time step is determined according to [25] as, Table 1.Gruneisen EOS properties and constants for water [18].
where C CFL is the Courant number usually set to unity.The smoothing length h is either static or dynamic.For the dynamic case h is obtained from, ( ) For numerical reasons, a minimum and maximum value is imposed on the smoothing length, i.e.H min h 0 < h < H max h 0 .In the present study, H min and H max were set equal to 0.5 and 1.5 respectively while h 0 is the initial smoothing length given by 0 , where d 0 is the initial distance between particles.The impact on the results from using a static or dynamic smoothing length will also be studied in the current study.

Boundary Conditions
Solid wall boundaries are modelled as rigid shell finite elements and the coupling between the boundaries and the SPH-particles are governed by a penalty based "node to surface" contact-algorithm [26].This is a so-called one-way contact where the SPH particles are defined as the "slave side" and the FEM elements as the "master side".In applying the penalty method, the slave nodes are checked for penetration through the master surface.If a slave node has penetrated, an interface spring is placed between the master surface and the node.The spring stiffness is chosen approximately in the order of magnitude of the stiffness of the interface element normal to the interface.The resultant force applied to the SPH particle in the normal direction of the FEM element is proportional to the amount of penetration.The wall boundaries used here can be compared with the Lennard-Jonestype boundary condition [24] where a central force is applied to a fluid particle.However, the normal component is considered only.Hence, the boundary condition used here is truly frictionless.For more information regarding the contact algorithm see [17].

Geometrical Setup
To validate the SPH model experiment, the SPHERIC Benchmark test case of dam break evolution over a wet bed is used.Numerous authors have adopted this test case for validation proposes, e.g.[14]- [16].The test case is based on the experimental work by [13] where water initially stored in a tank is released through a gate into a prefilled rectangular channel.The schematic arrangement and geometry is illustrated in Figure 1.
In the experiment by [13], several tank and channel depths were investigated.However, in the present work only one combination of tank and channel depth were considered, namely the tank depth L TD = 0.15 m and channel depth L CD = 0.018 m.The tank length L TL = 0.38 m was set as in the experiment and the channel length L CL = 1.62 m was smaller than in the experiment (L CL = 9.55 m) in order to reduce the overall number of particle and hence computational time.According to [13], a propagating bore develops after the removal of the gate.The stationary water in the channel resists the oncoming water and, as a consequence, a unstable "mushroom-like" wave forms and breaks in both forward and backward directions.In the experiment, the kinetic energy was usually enough to build up several breaking events until a smooth propagating bore was formed.Two CCDcameras recorded the events, thus digitized snapshots, wave profile data and average wave front velocities is available for comparison.Here, wave profile is compared only.According to [2], a general difference of approximately 0.043 s is observed between experiments and SPH simulations, thus this time shift is accounted for in the result section.The gate velocity was constant and measured to 1.5 m/s.According to the Benchmark test case instruction, it is vital to include the gate movement in simulations.

Numerical Setup
As stated above, the three EOSs including different values of the artificial viscosity constants q 1 and q 2 are to be tested and compared with experimental results.Thus, six cases is set up, see Table 2. Furthermore, for each case additional parameters were considered.Firstly, the spatial resolution of the SPH particles which was defined by dividing the channel depth L CD , here defined as the characteristic length scale of the problem, by 4, 5, 6 and 10.All SPH particles were initially placed on a structured grid with zero initial velocity.Secondly, the scaling of the FEM boundaries was tested where the side of a FEM element was set equal to 2x, 1x and 0.5x the initial interparticle spacing.Finally, the dynamic (D) and static (S) smoothing length was compared.All parameters are summarized in Table 3.
The initial density 0 ρ of the SPH particles was set equal to 1000 kg/m 3 and all simulations were stopped at 0.488 s due to the time shift proposed by [2].As the tank and channel depths was fixed at 0.15 m and 0.018 m respectively, the total number of SPH particles depend on the resolution only and ranges from about 4550 to 28600.A total of 144 simulations were conducted using the nonlinear finite element code LSTC LS-DYNA v. 971 R7.0.0 on multicore Linux workstations.

Post Processing
The post processing method used in this work incorporates the detection of the free surface as well as interpolation of data.The method is divided into two steps.Firstly, all nodes representing the boundary must be detected and secondly a Delaunay triangulation and Barycentric interpolation is performed, for a detailed description of Delaunay triangulation and Barycentric interpolation see [27] [28] and [29].The boundary node identification technique is based on the "ARC" method presented in [30].Thus, a boundary node is defined as a node with a non-overlapped part of its circumference.For these nodes the non-overlapped circumference is seeded with "dummy" which can be considered as the free-surface.The boundary node identification and especially the seeding of "dummy" nodes is crucial when arbitrary free-surfaces are to be detected.This is due to the Delaunay triangulation method which impose that the outer boundary is always convex, better known as the convex hull.The dummy nodes and the triangulation can be seen in Figure 2. Finally, the Barycentric interpolation using the triangulation onto a structured grid is performed.Where Barycentric interpolation f for an arbitrary interpolation point is defined as where b x is the barycentric coordinate and f x are the function values such as velocity at the triangle vertices.As

EOS Artificial viscosity constants
Gruneisen (q 1 = 0.1; q 2 = 0) and (q 1 = 1; q 2 = 1) Linear polynomial (q 1 = 0.1; q 2 = 0) and (q 1 = 1; q 2 = 1) Morris (q 1 = 0.1; q 2 = 0) and (q 1 = 1; q 2 = 1)  data is mapped to a structured grid, traditional techniques for data visualization such as streamline plotting can then be used.Apart from data visualization, the post processing method can be used to statistically quantify the difference between numerical and experimental results.By imposing both the numerical and experimental results on the same grid the non-overlapped grid points can be identified.The error is here defined as number of non-overlapped grid points divided by the total number of points in the experiment for each time step.The mean error is then the average of all time step errors.Non-overlapped grid points can be seen as an area and hence when the non-overlapped area goes to zero the mean error goes to zero as well.

Results and Discussion
The post processing method described in the section above was used to visualize the results.To the left in Figure 3, the absolute velocity field from one SPH case together with the digitized experimental data points (red) of the wave profile for the successive time steps is shown.Note that the time step given includes the time shift proposed by [2].To the right in Figure 3, the dummy nodes (red) representing the free surface together with 35 randomly distributed streamlines (blue) are shown.The data is taken from the D12 case using the Gruneisen EOS with artificial viscosity constants q 1 and q 2 set equal to unity.One of the most significant discrepancies when comparing the experimental and numerical results is that the wave front obtained from the SPH simulation is ahead of the experimental wave front, see Figure 3.This is evident for all studied cases and especially evident in the early stages when the "mushroom" jet mentioned by [13] starts to form and break.The numerical over prediction of the speed of the wave front is not isolated to the early stages as the predicted secondary splash in the later stages is ahead of the experimental results and the splash height is over predicted.The secondary splash is resolved in L CD /10 and to some extent in the L CD /6 cases only.One possible explanation that the SPH wave front is ahead of the experiments is due to the frictionless contact between the SPH particles and the FEM mesh, i.e. the free-slip boundary condition.There is a possibility to include a friction force in the contact formulation.However, the proposed friction formulation available is not applicable for these types of problems [17].However it may be possible to implement the dynamic particle boundary condition without major alteration to the code (LS-DYNA).Crespo et al. (2008) used the dynamic particle boundary condition and obtained good results with no tendency of the SPH wave front to be ahead of the experimental results.No tests using the dynamic particle in LS-DYNA has been conducted in present work.All boundary conditions suffer from one or several drawbacks and further research is needed.In order to better refer to the different flow features in the later stages of the dam break evolution, the domain is divided into four zones, see Figure 4.
Figure 8.Comparison of experimental (grey) and numerical results.Case D2 (L CD /10), Gruneisen, Linear Polynomial and Morris EOS and artificial viscosity constants; q 1 = 1, q 2 = 1 (blue) and q 1 = 0.1, q 2 = 0 (red).shaded grey area is the experimental results.Furthermore, the choice of setting q 1 = 1 and q 2 = 1 is shown in blue and q 1 = 0.1 and q 2 = 0 in red respectively.Few flow features are visible for the with lowest resolution (L CD /4), see Figure 5.The mushroom jet is almost not detected and zone four is not resolved at all.The simulations seem to behave too stiff to capture the significant flow features, regarding both the combinations of the artificial constants q 1 and q 2 .When increasing the resolution to L CD /5 (Figure 6), few flow features are seen for cases with artificial viscosity constants sett equal to unity.For the other choice, more features are visible such as a minor breaking wave.However, neither the secondary splash nor zone four is detectible with this resolution.Regarding the cases with resolution L CD /6 the mushroom jet is more pronounced than in the cases with lower resolution, see Figure 7. Breaking wave and secondary splash as well as the initial stages of zone four are visible for all EOSs using artificial constants set equal to q 1 = 0.1 and q 2 = 0.For the cases with finest resolution (L CD /10) most of the flow features shown in the experiment are now visible in the numerical results for both choices of the artificial viscosity constants, see Figure 8.Thus, the resolution of the SPH particles has a major impact on the outcome.Generally, a good agreement between simulations and experimental data for zone one being closest to the dam is shown, see Figures 5-8.The bump in the second zone stems from the mushroom jet that broke down in the upstream direction.Although the bump often appears in the simulations the horizontal position of the bump is generally poorly predicted.Cases with low resolution of the SPH particles and artificial constants set equal to unity seems to be unable to resolve the bump as exemplified in Figure 5.The third zone or the secondary splash zone show large differences between the different cases, compare Figure 5 and Figure 8.The final zone, zone four, is resolved for L CD /10 and to some extent in the L CD /6 cases, see Figure 7 and Figure 8.The dynamic versus a static smoothing length as well as the scaling of the FEM mesh showed no significant effects on the results.
As shown above there is a qualitative agreement between simulations and experiments if the set-up is fine enough.Using the method proposed in the post-processing section a quantitative comparison can also be carried out in the form of the mean error that is summarized for all cases in Figure 9.The error for all cases is within 5% of each other and showing a total of roughly 15% to 20% discrepancy between experimental and numerical results.Three significant trends are evident in Figure 9. Firstly, the major difference between the six cases is the choice of artificial viscosity coefficients q 1 and q 2 not the EOS.Selecting q 1 and q 2 equal to unity seems to produce the better result.However, as can be seen in Figure 5 and Figure 6 for the two lower resolutions, this choice of constants results in a highly viscous and unphysical solution, i.e. the artificial viscosity term is greatly over predicted.For the two finest resolutions, the contrary seems to be true, i.e. a more dynamic and fluid like behavior is obtained.In the authors opinion, the lower error obtained by setting q 1 and q 2 equal to unity is due to the highly viscous solution which reduces the wave front velocity and hence lowering the mean error.This makes it cumbersome to find a good relation of spatial resolution and artificial viscosity constants to obtain a trustworthy solution.Secondly, the results connected to the L CD /10 cases or the case with the highest resolution seems to produce the worst results.This is in sharp contrast to the expected outcome, as generally more nodes in any numerical method should produce a better result and also the qualitative investigation indicated this.One plausible explanation to this behavior is that more features are resolved using more particles such as the breaking wave and zone four.Also, the height of the secondary splash produced mainly in the L CD /10 cases are generally over predicted thus increasing the number of none overlapping grid points, i.e. increasing the mean error, as can be seen in Figure 8. Thirdly, all EOSs produce similar results which are expected as the density differences seen in the simulations are very small.However, the density differences are not small enough and unphysical pressure overestimation and large pressure oscillations are obtained ( 6 10 ∓ Pa), especially close to wall boundaries as well as near the free-surface.It is likely that the over prediction of pressure is due to the speed of sound c 0 which was set equal to the physical speed of sound of 1484 m/s.Usually, c 0 is set equal to at least 10 times the expected maximum velocity in the flow [24] [31].Initial tests with reduced speed of sound ( ) were performed, where H ref is the maximum water height in the tank (H ref = 0.15 m) as according to [14].However, severe particle collapse and contact instabilities were obtained and, hence, the physical speed of sound was used.Generally, the pressure oscillations near wall boundaries and free surfaces can be avoided by use of either density filters or correction of the kernel and/or kernel gradient.Considerable efforts have been devoted to overcome these problems in the SPH research community [4].Several corrections and filtering techniques have been proposed and are available in literature, see for instance [32].
In the zeroth order Shepard density filter, the simplest and quickest type of density filter, oscillations are smoothed out by filtering density and then re-assigning a new density to each particle [33]- [35].Kernel and kernel gradient correction techniques are generally included to handle the truncated kernel function close to freesurfaces and boundaries where conditions of consistency and normalization fail [32] [34] [36] [37].Furthermore, the development of a truly incompressible SPH (ISPH) solver have shown great promise in estimation of pressure and show less oscillations, see for instance [16].As with boundary conditions no filtering or correction techniques are currently implemented in LS-DYNA and the effects of such methods are not possible to test at this point.

Conclusion
The two dimensional (2D) SPHERIC Benchmark test case of dam break evolution over a wet bed has been used to investigate the impact of several parameters when performing free surface flow simulations.The choice of EOS, the artificial viscosity constants, a dynamic versus a static smoothing length, SPH particle resolution and the FEM mesh scaling of the boundaries were investigated.The numerical results showed generally a tendency to be ahead of the experimental results, i.e. to have a greater wave front velocity.This was argued to be an effect of using a frictionless contact algorithm between the SPH particles and the FEM boundaries.The choice of EOS, FEM mesh scaling as well as using a dynamic or a static smoothing length showed little or no significant improvement to the result.However, all EOSs showed larger than physical pressure levels in the order of 10 6 Pa and large pressure oscillations.This is believed to be an effect of not reducing the speed of sound in the calculations, which is usually done.The SPH particle resolution and the choice of artificial viscosity constants had a major impact.The method used for comparing numerical and experimental results showed increased mean error for both highly resolved cases and artificial viscosity constants set equal to 0.1 and 0. However, by visual observation, it was noted that increasing the number of particles representing the system increased the number of flow features resolved and that a highly viscous solution was obtained by setting the artificial viscosity constant to unity.Thus, the relation between the artificial viscosity constants and the particle resolution and its impact on the behavior of the fluid needed to be further investigated.

Figure 1 .
Figure 1.Schemetic arrangement and geometry for the dam break test case.

Figure 3 .
Figure 3. Absolute velocity field [m/s] together with experimental data (left) and dummy nodes and streamlines (right).

Figure 4 .
Figure 4.The later stages of the dam break evolution divided into four zones.

Figure 9 .
Figure 9.The calculated mean error for each EOS and artificial viscosity constant setup.

Table 2 .
Test case setup, EOS and artificial viscosity constants.

Table 3 .
Additional parameters tested for each case, SPH resolution, FEM scaling and smoothing length behavior, static (S) and dynamic (D).