Introduction

Many analytical one-dimensional hydraulic models that find practical application involve flow lengths that are much greater than its depth or width, such as that in channels, rivers and piping systems. For cases where vertical fluid movement within the flow or variations in depth are small but not negligible, analytical solutions for the shallow water equations (SWE) have been developed [1] to provide a good approximation if the boundaries are such that they have an analytical solution. For more complex boundaries, SWE can be used to generate a set of differential equations to be solved numerically at a much lower computational cost than applying a general purpose two-dimensional numerical model of fluid flow.

However, precision is just one of many properties that a numerical model may exhibit, and the particular case of flow with a high aspect ratio can give greater insight as to the impact of highly non-symmetric resolution when using a general-purpose numerical fluid flow model such as WSPH. This is especially attractive for cases with an open boundary since no dynamic remeshing algorithms or surface detection procedures are needed, as is the case of the Finite Volume Method (FVM), among many others. The free movement of particles does mean solid or elastic boundaries need coherent meshless formulation through equations that have numerical parameters that must be chosen carefully to avoid deficient modeling of flow close to discharge areas, like those at a sluice gate or pipe outlet.

The transient case of a sluice gate opening at the end of a long channel was chosen to set up the models for this study, because it serves as an excellent showcase for potentially insightful phenomena. The results were then compared directly with previously published data, including an analytical solution [1], and the vena contracta locus and depth from physical experiments [2], which were for a dammed channel with a 50:1 length to depth ratio and constant rectangular cross section.

Previous studies using SPH to model sluice gate opening [3, ]4], successfully have concentrated on the formation of the hydraulic jump and the effects of the shape and roughness of the discharge area. Other studies [5, 6] deal with the addition of an obstacle close to the discharge area and how the shape of the hydraulic jump changes in time. Studies using the analytical SWE equations to calibrate SPH models [1, [4], do not mention the effect of low resolution on shallow depths or rely on relatively short lengths to model high aspect ratios.

This study focuses on the effect of using the entire length of the channel according to boundaries described in previously published studies, and to show the behavior of the shallower portion of flow during the first few instants after the opening of the sluice gate over a dry bed, to illustrate the impact of vertical resolution on the SPH model of this type of flow more specifically.

A Brief outline of the wsph method

The smoothed particle hydrodynamics (SPH) approach is a meshless discrete method where the volume of fluid is divided into a finite number of particles, usually uniform in size and nature. As with any other Lagrangian approach to discretization, since the particles are not fixed in space, their movement is modelled with the fluid dynamics equations for conservation of mass and conservation of momentum [7]. Any other applicable equations required to model the physical state of the fluid (changes in density, viscosity, temperature, for example) must, of course, be considered as well. The main difference of the SPH method in comparison with Newtonian particle methods¹, such as the Multiple Particle Method (MPM) and Particle In Cell (PIC) method [8], is that any single particle interacts both with immediately adjacent neighbors and those within a spatial radius equal to khas shown in Figure 1.

To avoid an unnatural distorted behavior of the particle, the farther away the neighboring particle is, the less it affects its neighbor. In the SPH method, this is done through a weighting function called the kernel functionW(x,h)[9], which can be any continuous, spatially symmetric function. However, to guarantee conservation of physical meaning, the value of the scalar or vector field that the function is used with must not introduce artificial alterations [7]. Also, it should be as computationally cheap as possible. This function smooths the distribution of the value of any variable associated to particle dynamics, such as speed, acceleration, density, enthalpy, entropy, or viscosity, among many others that may be necessary for the particular case. This smoothing helps to stabilize the computational process and is extensive to the values of any integrals or derivatives that need to be considered in the conservation of mass and the conservation of momentum equations applicable to the flow [10] With a larger spatial radius, more smoothing effect occurs, however, too large of a value will cancel out any local dynamics and the flow will exhibit numerical clumping. Therefore, a proper choice for smoothing length h (usually a value between 1.1 and 2.0 times the nominal particle size) is a matter that requires careful consideration [10,11].

Figure 1
A visual representation of the smoothing function W(x,h).

Since the central challenge for solving the fluid mechanics equations is obtaining the value f(x_i) of the dynamic and physical variables pertaining to the state of the flow, the SPH method approximates them through a numerical scheme known as the particle approximation[7]. The general statement is that for a particle i, located at position x_i an approximation of f(x_i) can be calculated as a weighted average of the value of f(x)for all N neighboring particles using (1), which can be readily transformed to computational source code:

According to this definition, the position of any particle must be known, making this scheme an explicit forward time-step method. The origin and mathematical development of (1) and the choice for the smoothing function and many other associated parameters can be studied in detail in monographs on the topic of SPH applied to fluid flow [10,11].

A complete description of all applicable equations of fluid flow is beyond of the scope of this article, so only the main numerical driver for particle dynamics will be mentioned in this section. This is the application of the particle approximation to the equation of conservation of momentum [7], shown as (2) for two-dimensional flow in orthogonal coordinates:

The usual terms for fluid flow are: the speed vector u, pressure P, viscosity m, plane strain rate tensor \epsilon^{xy}, and density \rho. However, particle approximation produces additional terms which include particle mass m, a second appearance of density \rho and partial derivatives of the smoothing function. The reason to mention this equation here, aside from its importance, is that the multiplied density terms \rho_i\rho_jnot only affect the value of acceleration or any other variable through (1) but also drive the variations. In the case of the weakly compressible SPH method (WSPH), this is due to the form of the equation of state, that gives the fluid an artificial compressibility to correlate the variation of pressure as a function of local density \rho_i. In this study, (3) was used and is known as the Monaghan equation [7]for local effective pressure affecting particle i:

The value of \rho_0 is a constant reference density of the fluid being modelled; the true speed of sound in the fluid c₀ can be used, but for many subsonic flow phenomena, a much lower value can be used without considerable loss of accuracy [10], so larger time steps can be used and therefore lower computational cost is required. However, if this is done, to avoid numerical instability and particle interpenetration, it is important that c₀ is at least an order of magnitude larger than the highest expected fluid speeds in the case being studied. Empirical constant g is conventionally set at 7, but can be adjusted by the user. Most importantly, the value of \rho_1 evolves according to (1) which, like all other evolving physical properties in SPH, is a weighted local average of the implied densities of surrounding particles. Values of \rho_i far enough from \rho_0 can produce artificially high values of P that can turn into small local numerical instabilities or even artificial numerical shock waves all across the modelled domain.

Therefore, these terms are key to the stable evolution of all variables calculated this way, if an incompressible fluid is to be modeled, one way to prevent numerical instabilities is to add artificial compressibility to the numerical scheme. The addition of such an equation of state for nearly incompressible fluids is generally known as the weakly compressible SPH method (WSPH) used for this study [7, 10,11].

This is the general framework of SPH as a numerical method applicable to solve the equations of fluid dynamics, shown as a starting point. However, it is a complex, well-evolved method as shown by the considerable number of publications that use SPH to model many different fluid flow phenomena, like those mentioned in comprehensive reviews of the state of the art such as [12] and up-to-date revisions of future challenges for the application of SPH methods in general [13].

Methodology

As mentioned above, available validation data drove to the selection of the particular dimensions and boundary conditions for the numerical case setup. References [1,2], published data for the case of a long (50:1), horizontal open flow channel with a completely absent gate, a conventional dam break. These same references also study this setup but with an open gate with a height equal to 20 % of the initial dam height. The available analytical solutions are for any instant, but vena contracta size is given for a time point t = 5.0 s after the gate opening for the following three cases: a dry bed, shallow wet bed and deep wet bed. In this study, the dry bed configuration is of interest and was used for the case setup. Figure 2 shows the rectangular arrangement of the two-dimensional fluid that replicates the proportion and size of the physical experiment [2].

Figure 2
Initial particle arrangement for 2D WSPH model of an open dam gate.

In all cases, WSPH two-dimensional Dx uniformly sized particles were used and arranged in an evenly spaced rectangular array. The Morris equation for the conservation of momentum [10] shown in (1) was used to evolve momentum, and the Monaghan equation of state [7] with an artificial speed of sound of c₀ = 80 m/s was used since fluid speeds were expected to be well under 10 m/s. Properties of clean water at T = 20 ºC were used (viscosity m = 0.000896 Pa·s and density r = 998 kg/m³). Cubic form [7], for smoothing function W(x,h) was used in all cases with a smoothing length of h = 1.1·Dx among the simplest typical value, important in terms of computational cost due to the amount of times it needs to be evaluated along with its spatial derivatives, at least once for each particle and must be recalculated at each time step.

All solid boundaries, including the static sluice gate when applicable, were modelled using a single line of virtual particles spaced evenly at Dx/2; their behavior was modeled with the Monaghan repulsive force equation [14], which affects nearing fluid particles only when within a certain range r₀ of the boundary and becomes higher in inverse proportion to the distance between the boundary particle and the fluid particle (this proportionality can be adjusted empirically by the user). No turbulence models or correction factors were applied. Time step was constant and chosen using the Verlet criteria [11], which is based on expected numerical propagation of data. Iterations were set to continue up to a physical time of t = 5.0 s. All other significant parameters varied according to particle size and boundary dimensions, and are specified in Table 1 for each of the three variants used for this study.

Table 1

Subcase stup for conventional dam break and 20% sluice gate

The first two setups were intentionally prepared with shorter boundaries to reduce computational time and to determine how much of a difference appeared between a 70 % length model (Open HR) and a 100 % length model (Gate HR) both in terms of computational time and model precision. The low-resolution model (Open LR) was used to verify evolution to determine whether an even higher resolution model (>20 000 particles) was necessary.

These models were processed using the authors own source code that implements the WSPH method, written in C programming language. The first version of the functional source code, named Gap-SPH is freely available at a public repository (https://sourceforge.net/projects/gap-sph/) under GNU GPL v.3 license terms.

Results and discussion

The first run, which is with few layers of particles in the vertical direction (only ten rows of particles, in contrast to 350 columns horizontally), displays a reasonable exact solution except for the leading front segment, which has noticeable lag; more than 15 %, in comparison with the analytical model [1]. This may be due to the fact that this area is made up of a single layer of particles, as shown in Figure 3, where a good number of them are not affected by enough neighbors to actually move at a better pace. Also, the repulsive forces from the boundary virtual particles are large enough to become dominant over any inertial or viscous terms, forcing an almost perfect, but artificial, alignment of the first layer of particles. By increasing the vertical resolution of the problem to twenty rows (a total of 14 000 particles), the quality of the resulting distribution is greatly improved, as shown in Figure 4 for the same general area at the chosen reference instant t = 5.0 s. In this higher resolution model, similar problems to those for the lower resolution model can be observed at the boundary.

Figure 3
Detail of the water front for dam break with 3500 particles at t = 5.0 s (model OpenLR, 1 m x 35 m).

In these images, it is very important to keep in mind the nature of the model, which is of a high (1:50) aspect ratio; as such, there is correspondingly high graphical distortion in the vertical scale in relation to the horizontal scale. Some particles may appear overlapping, but this is not a particle penetration problem, it is only a consequence of the high horizontal density of the plot, necessary to show such a long model in a single graph.

Figure 4
Detail of the water front for dam break with 14 000 particles at t = 5.0 s (model OpenHR, 1 m x 35 m).

The complete profile for that instant is directly compared in Figure 5F with the analytical solution for the SWE developed by [1] and applied to the dam failure with the specified dimensions according to [2] Since the analytical model implies a 1 m x 50 m fluid, but the model labeled OpenHR is for a 1 m x 35 m fluid, the data has been purposefully displaced -15 m in the horizontal direction.

Figure 5
Calibration of hydraulic profile for dam break with 14 000 particles at t = 5.0 s.

The apparent rounding at the point of discontinuity at the surface close to x = 20 m cannot be considered analytical model imprecision or a defect of the numerical model, but rather an irrelevant property of the analytical solution; a real linear fluid would not display a sharp edge, so it is not of interest that this is reproduced in the numerical model. The comparatively minor inaccuracies related to full channel depth (1 m) that need to be addressed are a small lead (a maximum of +4,8 % for position x = 27 m) in the upper half of the fluid and a comparable lag in the lower half (a maximum of -5,1 % for position x = 44 m).

However, the differences are not enough to reject the validity of the long channel dam failure model and additional resolution did not improve precision, so the dam with an opened gate was modeled with this same resolution but for the full 1 m x 50 m fluid. The gate opening measures 0.20 m starting from the bottom of the dam, which amounts to 20 % of the nominal dam height. For this situation, an analytical solution is also available [1], and the vena contracta depth and position is checked with experimental data [2]. Figures 6 7 8 show position fields for the area of interest of the model for three significant instants.

Figure 6
Detail of hydraulic profile for 20 % sluice gate with 20 000 particles at t = 0.5 s

Figure 7
Detail of hydraulic profile for 20% sluice gate with 20 000 particles at t= 2.5 s

Figure 8
Detail of hydraulic profile for 20% sluice gate with 20 000 particles at t = 5.0 s

During the first few instants (t < 1.2 s), the particles are very regularly aligned in a square array, almost as if they were static, with the exception of the open discharge section (the open gate) and the particles immediately adjacent to boundaries, as in Figure 6. Later, the distribution of particles is more irregular due to the complex dynamics near the gate, and fluid accumulates in the discharge zone, which drops in level at that location as in Figures 78.

The effect of vena contracta is definitely present, and its depth and position coincide almost exactly with the analytical solution as seen for x = 50 m in Figure 8, as the aperture is 0.2 m, while the discharge jet is initially 0.12 m deep (analytical solution gives a depth of 0.13 m), then expands to about 0.16 m deep for a short distance (it is not yet a hydraulic jump as there is not enough fluid downstream to provide such a boundary) and then at an average depth of about 0,11 m before reducing its height in a manner similar to a conventional water head flowing over a dry bed.

Figure 9
Calibration for the hydraulic profile for 20% sluice gate with 20 000 particles at t = 5.0 s

A comparison with the analytical solution for this dataset is shown in Figure 9 where there is a good degree of correspondence (2.8 % RMS relative error to full channel depth for the model between 30 m < x < 60 m, with a maximum error of -7.1 % at x = 42 m) except for the area with only one layer of particles available to model the water front (occurring between 62 m < x < 80 m at t = 5.0 s). As in the low-resolution dam break case with 3 500 particles, there are not enough numerical energy transferal processes between particles to properly model the continuum. The error in that area is large, both in qualitative and quantitative terms, and part of the reason is because of the virtual particle equations that model the boundary through a strong repulsive force. Since this ends up being a series of particles in a line, there is not much of a chance of energy transfer thorough viscous shear stress or the increment of local density which would in turn contribute to local pressure through the equation of state.

As in previous cases, this can be partially improved by increasing the resolution, which obviously has a high computational cost not only because of the number of particles, but the time step must be decreased as well. As a reference, the model with 3 500 particles took only two hours of computational time; the model with 14 000 particles, with the same size boundaries and equal physical time took more than thirty hours of computer time to complete. On the other hand, additional tuning of the boundary repulsive force could also improve precision somewhat without allowing particle leaks, because the first layer of particles has practically no vertical movement but makes up almost a quarter of the 0.20 m gate aperture, so the boundary trigger distance r₀ may introduce unwanted distortions to gate flow condition since, to be effective, it must be in the order of magnitude of the nominal particle size \Deltax.

These results show a good correlation and have readily illustrated many characteristics of applying a general-purpose numerical method to cases with extreme proportions. This paves the way for the proposal of future developments for this type of application and gives pointers to tune applicable numerical model parameters of the WSPH method when modeling long, shallow cases of open channel flow.

Conclusiones

The application of the WSPH numerical scheme to a classical two-dimensional hydraulics problem proved to be effective to model flow over a high aspect ratio according to the special cases explored in this study, where mean RMS error was well under 5 % with a maximum of 7 % in critical areas in relation to upstream channel depth. Local error was in the order of 10 % (1 % in relation to upstream channel depth) at the gate discharge where the vena contracta forms. Vertical resolution and the repulsive force that acts at the boundaries both affect precision, and this means it can be compensated for in similar models. Other particularly interesting aspects that came from this study were:

• The numerical model exhibited the phenomenon of vena contracta properly, and further along, a small hydraulic jump in at the gate outlet, both of which are expected according to the analytical model and the experimental observations. Although precision was acceptable considering how low vertical resolution was in that area (one to four layers of particles, but only 10 % local error in predicted depth in relation to analytical and experimental information for that locus), the repulsive forces from the virtual particles used to model the solid boundary, also contributed to diminish precision.

• In the discharge area close to the gate, artifacts similar to physical waves and vortices appeared, such as those that appear in open channel flow phenomena. At the entrance to the gate, since vertical resolution is that of at least twenty particle layers, wavelengths can be measured to determine whether they are numerical or physical in origin. At the discharge, since resolution is so low, any periodic irregularities will have a numerical origin and cannot be considered as a replica of vortexes that may occur in the real fluid. This strategy can be very useful to quickly validate the behavior of a numerical model, since wavelengths in channel flow cases can be readily calculated with Froude number corrected wavelength equations or measured in a relatively straightforward manner with physical experiments.

• There was excellent correlation between the results of the numerical model, even in the absence of stabilization algorithms or turbulence models, so it is likely that, for cases with a high length-to-depth ratio, future studies will not need to assess its impact or consider its incorporation. A more precise correlation should be obtained through good enough resolution, since the dynamic phenomena related to SWE are so gradual that correctors for vorticity or impact have little or no effect.

• Low resolution in the water front is the main problem in the models studied, and this can be a serious issue when modeling sediment carry-over. In these cases, interaction with the boundary layer is key to its accuracy, so a different type of boundary particle will have to be adopted or developed to properly interact with the fluid particles.

• Even when sediment transport is not relevant, to improve the accuracy of SPH models of shallow water problems without requiring computationally costly increments of the number of particles, future studies may consider creating a new type of boundary particle customized specially for SWE-type problems.

The SPH scheme is well-adapted to many different open-flow hydraulics problems with geometric discontinuities such as the general case of flow channels or partially filled piping systems, and the convenience of not having to necessarily elaborate surface detection algorithms or dynamic remeshing, develops great advantages in comparison with mesh-dependent methods such as FVM or VOF, and even hybrid approaches can be developed to attempt to meld the positive aspects of each methodology. However, the proposal and modification of different solid or elastic boundary equations have proven in many instances that it is one of the most critical aspects of this method that is continually evolving and is far from being a standard procedure applicable to any case of fluid flow.

Acknowlegements

The authors would like to thank the anonymous peer reviewers and editors for their insightful suggestions and corrections, which greatly improved readability, style and completeness. We appreciate the time they invested in reviewing this paper and the consideration of the personal mark and purpose with which it has been written.

References

L. Cozzolino, L. Cimorelli, C. Covelli, R. Della Morte, and D. Pianese, “The analytic solution of the shallow-water equations with partially open sluice-gates: the dam-break problem”, Advances in Water Resources, vol. 80, pp. 90–102, Jun. 2015, doi:10.1016/j.advwatres.2015.03.010.

A. Defina and F. M. Susin, “Hysteretic behavior of the flow under a vertical sluice gate”, Physics of Fluids, vol. 15, no. 9, pp. 2541–2548, Jul. 2003, doi:10.1063/1.1596193.

M. Khanpour, A. R. Zarrati, and M. Kolahdoozan, “Numerical simulation of the flow under sluice gates by SPH model”, Scientia Iranica, vol. 21, pp. 1503–1514, Oct. 2014.

C.P. Sun, D.L. Young, L.H. Shen, T.F. Chen, and C.C. Hsian, “Application of localized meshless methods to 2D shallow water equation problems”, Engineering Analysis with Boundary Elements, vol. 37, pp. 1339–1350, Nov. 2013, doi:10.1016/j.enganabound.2013.06.010.

S. Gu, X. Zheng, L. Ren, H. Xie, Y. Huang, J. Wei, and S. Shao, “SWE-SPHysics Simulation of Dam Break Flows at South-Gate Gorges Reservoir”, Water. vol. 9, art. 387, May 2017, doi:10.3390/w9060387.

S. Gu, F. Bo, M. Luo, E. Kazemi, Y. Zhang, and J. Wei, “SPH Simulation of Hydraulic Jump on Corrugated Riverbeds”, Applied Sciences. vol. 9, art. 436, 2019.

J. J. Monaghan, “Smoothed particle hydrodynamics”, Annual Review of Astronomy and Astrophysics, vol. 30, pp. 543–574, 1992, doi: 10.1146/annurev.aa.30.090192.002551.

G.H. Cottet, “Multi-physics and particle methods”, presented at the Second MIT Conference on Computational Fluid and Solid Mechanics, Cambridge, United States of America, June 17-20, 2003, pp. 1296-1298. doi:10.1016/B978-008044046-0.50319-5.

R. A. Gingold, and J. J. Monaghan, “Smoothed particle hydrodynamics: theory and application to non-spherical stars”. Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375–389, 1977, doi:10.1093/mnras/181.3.375.

G. R. Liu and M. B. Liu, Smoothed Particle Hydrodynamics: A Meshfree Particle Method. Upper Saddle River, N.J.: World Scientific, 2003.

D. Violeau, Fluid Mechanics and the SPH Method: Theory and Applications. Oxford: Oxford University Press, Oxford, 2012.

J. J. Monaghan, “Smoothed particle hydrodynamics and its diverse applications”, Annual Review of Fluid Mechanics, vol. 44, pp. 323–346, January 2012. doi: 10.1146/annurev-fluid-120710-101220.

R. Vacondio, et al., “Grand challenges for smoothed particle hydrodynamics numerical schemes”. Computational Particle Mechanics, 2020, doi:10.1007/s40571-020-00354-1.

J. J. Monaghan, “SPH without a tensile instability”, Journal of Computational Physics, vol. 159, pp. 290–311, April 2000, doi:10.1006/jcph.2000.6439.

Notes