Pressure evaluation during dam break using weakly com- pressible SPH

This paper presents a solution of a dam break problem in two dimensions obtained with smoothed particle hydrodynamics (SPH) method. The main focus is on pressure evaluation during the impact on the wall. The used numerical method and the way of pressure evaluation are described in detail. The numerical results of the kinematics and dynamics of the flow are compared with experimental data from the literature. The abilities and limitations of the used methods are discussed.


Introduction
A Dam-break problem is a relatively widely studied freesurface problem. In its simplest configuration, a liquid column collapses in a gravity field, forming a surge on a flat dry bed of a tank. The surge then impacts a vertical wall of the tank, shaping a vertical jet along the wall and a rolling wave afterwards. The relative simplicity of the problem makes it suitable for both experimental and numerical investigation.
The experimental works dealt mainly with kinematics and dynamics of the flow. Martin et al. investigated the position of the surge front and the height of the column as a function of time [1]. Lobovský et al. focussed primarily on pressure loads on the vertical wall during the surge impact; pressure was measured in various heights above the bottom of the tank, and the results were statistically evaluated [2].
Works dealing with a numerical solution of the dam break problem were often combined with an experiment, which proved the numerical method to be capable of simulating free-surface flows. Such works are for example Koshizuka and Oka (moving particle semi-implicit) [3], Cruchaga et al. (finite element method) [4], and Hu and Sueyoshi (constrained interpolation profile and moving particle semi-implicit) [5]. These works deal only with the kinematics of the flow.
Some works presenting numerical solution of freesurface flows using SPH have also been published. These which dealt with dynamics of the flow are for example Colagrossi and Landrini [6], Ferrari et al. [7], and Laigle and Labbe [8]. All these authors compared their results with experimental data published by Zhou et al. [9].
The aim of this paper is to present a method of pressure evaluation in SPH and compare the obtained results * e-mail: petr.jancik@fs.cvut.cz * * e-mail: tomas.hyhlik@fs.cvut.cz with experimental data published in [2]. Two main issues to address in order to reach this goal are a correct formulation of wall boundary condition, and pressure oscillations occurring in weakly compressible SPH solution. The basic kinematics of the flow is compared with the experiment as well to validate the method.

Computational method
SPH is a mesh-free particle method. Its Lagrangian nature makes it suitable for free-surface and multi-phase flow problems because interfaces are tracked naturally. The governing equations of fluid motion, continuity and momentum equation, in the Lagrangian framework are where t, ̺, v, p, and f denote time, density, velocity, pressure, and external body force. Fluid is modelled as inviscid.
Since the fluid is considered to be compressible, an equation of state is needed to close the system of equations. It is common in weakly compressible SPH that an artificial equation of state is employed. In this case a simple equation is used. Variables c and ̺ 0 are numerical speed of sound and reference density. Numerical speed of sound should be chosen approximately ten times higher than the maximal flow velocity. This choice keeps the maximal Mach number in a reasonable range of values. If it was too high, compressibility would become an issue; if it was too low, time steps would have to be too small, and the computation would become time expensive [10]. Note that the equation of state is not a function of temperature or internal energy thus energy equation does not have to be solved. Using SPH techniques, partial differential equations (1) and (2) are spatially discretised and transformed into the following ordinary differential equations: Indices i and j denote particles, m is mass, and W i j is so called smoothing function. There are more smoothing functions. In this work, the truncated Gaussian smoothing function was employed, which can be written as where h is so-called smoothing length d is number of spatial dimensions, and non-dimensional particle distance x denotes a vector of spatial coordinates. The truncation of the function does not affect the solution, but it has a significant effect on computational performance because the number of interacting particles is reduced. There is an additional diffusive term in discretised continuity equation (4). In more detail, it can be written as where δ is a coefficient of diffusion. This term reduces pressure oscillations in the solution [11]. Another additional term Π i j in the discretised momentum equation (5) is an artificial viscosity term. It serves for stabilisation of the numerical solution, and it can be written as where α is a tuning parameter and term with ε prevents singularity. All variables with both indices i and j denote a mean value taken from particles i and j.
Fluid is considered to be inviscid, so wall boundary condition is modelled as a free-slip. Walls are made of three layers of so-called dummy particles, which take part in summations in equations (4) and (5). Using this approach, particles near walls are uniformly surrounded by neighbouring particles and field variables can be evaluated correctly. The pressure of a dummy particle is determined by the formula where indices w and f denote wall and fluid particle [12]. Density is then obtained from the equation of state (3). Dummy particle velocities are set to zero.
To further refine the obtained pressure field, density reinitialisation is conducted. New density values are obtained using the formula This reinitialization is performed every 20 time steps. It is similar to what used Colagrossi and Landrini [6] but with simpler zero-order Shepard function [13]. Pressure values at certain locations are evaluated using virtual sensors. These sensors obtain values of pressure from fluid particles in the vicinity using the formula where index S denotes sensor. The size of the sensor influence domain is dependent on the smoothing length h of a smoothing function W. Too small value of the smoothing length leads to fewer interacting influencing particles and noisy output. On the other hand, if the smoothing length is too great, distant particles are taken into account and the output is not correct. A good choice of smoothing length appears to be the same value as the one used in the simulation. For this choice, the sensor influence domain contains between 10 and 20 fluid particles. This method is not suitable for coarse spatial resolution simulations because, again, to distant particles influence the output.
To further suppress high-frequency noise on sensor output, filtering in the time domain is conducted. Since this filtering is done after the simulation, Gaussian kernel function (6) was used as a low-pass filter. The only changes are that number of dimensions d = 1 and a smoothing length h is replaced by a smoothing time period τ. The time filtering can be written where indices m and n denote time steps. A correct choice of smoothing period τ is problem dependent. It should be tuned to suppress high-frequency numerical noise, but it should not affect sudden changes caused by physical phenomena.
The ordinary differential equations (4) and (5) have to be integrated in time. In this case, the modified leap-frog algorithm is used, which can be written as v ρ x n+1 i = x n i + ∆tv where indices in superscript denote time step. Since the leap-frog algorithm is an explicit integration scheme, its numerical stability depends on a size of a time step. To determine a maximal time step, the Courant-Friedrichs-Lewy condition for SPH in the form was used.

Results
The numerical simulation was set to match the experiment carried out by Lobovský et al. [2]. The solution is performed in two dimensions, and the surrounding gaseous phase is omitted. Both these simplifications significantly improve computational performance. The initial configuration of the problem is in Fig. 1. Initially, the surge propagates towards the right wall of the tank (Fig. 2). After the liquid impacts the wall, a vertical jet along this wall forms. The top of the jet loses its momentum, while the bottom is still heading upwards. The outcome of this is a rolling wave (Fig. 3). This wave then impacts the liquid surface. Consequently, a jet heading towards the left wall of the tank forms (Fig. 4).  influence of resolution shows up more in the pressure evaluation because of the larger sensor influence domain.

Kinematics
To compare kinematics of the flow with experiments, nondimensional variables are used. Non-dimensional time, surge front position, and surface height are defined where g is the vector of gravitational acceleration, w is surge front position, and h is liquid surface height. Variables w 0 and h 0 denote initial width and height of the liquid column. First, the surge front position obtained from the simulation was compared with experimental data from [1-3, 5] (Fig. 5). The computed position of surge front is in a relatively good agreement with all the experimental data.
To further validate the kinematics of the flow, liquid surface height was measured at positions H1 -H4 (see Fig. 1) and compared with experimental data from [2]. The comparison is conducted only before the rolling wave and the following jet reach the positions of measurement. These phenomena also occur in the experiment, but they are quantitatively different. The differences are probably caused by viscosity and turbulence, which is not modelled  in the simulation. The absence of the surrounding gaseous phase also plays its role [6].
The simulation is in a good agreement at positions H1 (Fig. 6), H2 (Fig. 7), and H3 (Fig. 8). However, a difference in a shape of the propagating wave is noticeable at positions H2 and H3 near T = 1.7 and T = 2 respectively. This difference further expands, and at position H4 (Fig. 9), the difference between the experiment and the simulation is quite significant. Again, this is possibly caused by the absence of viscosity and turbulence modelling.

Dynamics
Pressure was evaluated at four points P1 -P4 (see Fig. 1) placed on the right wall of the tank. Pressure is taken as a non-dimensional quantity and it is compared with experimental data from [2]. The record from sensor P1 is shown in Fig. 10. There is a good agreement in the very initial phase of the impact.   There is a sudden rise in pressure and the maximal value predicted by the simulation is about 10% higher than the experimental value. Then the pressure decreases and the simulation values remain above the experimental ones. At about T = 6, there is a slight pressure rise in the experimental data. This rise occurs in the simulation as well, but later -at about T = 6.5. It is so because this pressure rise is connected to the rolling wave impact, which occurs later in the simulation than in the experiment.
Pressure from sensor P2 is in Fig. 11. There is a peak in pressure value as in the previous case, but its amplitude is lower. The maximal value is very similar in the simulation and the experiment. However, the following drop of pressure values is lower in the simulation, and it stays about 20% above the values from the experiment. Just like in the previous case, the pressure rise appears later in the simulation than in the experiment, and it reaches slightly higher values.  The record from sensor P3 is shown in Fig. 12. The peak occurs again, but it is lower and wider than in the previous cases. The simulation predicts the maximal value about 25% lower than the maximal measured value, but it drops slower, and it stays up to 30% higher than the experimental values. The pressure rise connected with the rolling wave impact occurs like in the previous cases.
Sensor P4 (Fig. 13) is located considerably higher above the tank bottom then the previous ones. That makes the output signal qualitatively different. The simulation predicts a gradual growth in pressure values, but a sudden rise followed by an almost constant level values occurs in the experiment. Consequently, pressure level obtained from the simulation is up to 50% lower than the values from the experiment. What stays the same for all sensors is the rolling wave impact rise, which is quite significant this time.

Conclusion
In this work, a dam break problem was solved using a weakly compressible SPH method. The obtained results were compared with experimental data from the literature. The kinematics of the flow is in relatively good agreement with available experimental data, especially at the beginning of the problem solution. Later, an absence of viscosity and turbulence modelling becomes evident, particularly during the formation of the rolling wave and during its impact. After the rolling wave impact, a neglection of the gaseous phase is not appropriate, as well as twodimensional simplification because air remains trapped under the rolling wave in experiments, and it escapes in the form of bubbles.
The pressure was evaluated in four points on the right wall of the tank. In general, the deviation of simulation and experiment increased with the height above the bottom. Three lower placed sensors gave good quantitative agreement with the experiment with the maximal pressure value deviation up to 30%. The experimental and numerical data from the highest sensor differ considerably. A common feature for all sensors in both experiment and simulation is a pressure rise after the rolling wave impact; however, it occurs later in the simulation than in the experiments.
To suppress pressure field noise in the space domain, artificial diffusion and density reinitialization were implemented. Possible problems with the used artificial diffusion model have been reported for long-time hydrostatic simulation, but they did not occur in the solved problem [14]. For correct evaluation of the kinematics of the flow, a relatively coarse resolution is sufficient. However, domain of influence of a pressure sensor is resolution dependent and a finer resolution is necessary.
The output signal from the virtual pressure sensors was further filtered in the time domain. A filtering function which uses 'future' values is employed since the time filtering is done after the computation. If the filtered values were needed during the computation, different filtering function would have to be employed. The filtering function and smoothing time period seem to be chosen correctly because high-frequency noise was mostly suppressed while sudden changes in pressure values, especially steep rises after the impact, were captured properly.
Despite the simplicity of the model, the results are in relatively good agreement with the experiments. To obtain better results, viscosity and turbulence should be modelled. The two-dimensional simplification and neglection of the surrounding gaseous phase are not entirely correct as well, especially in the later phases of the simulation. Surface tension can also play its role when droplets and bubbles emerge.