Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle-hydrodynamics formulations in simulating free-surface flows

The theoretical formulation of the smoothed particle hydrodynamics (SPH) method deserves great care because of some inconsistencies occurring when considering free-surface inviscid flows. Actually, in SPH formulations one usually assumes that (i) surface integral terms on the boundary of the interpolation kernel support are neglected, (ii) free-surface conditions are implicitly verified. These assumptions are studied in detail in the present work for free-surface Newtonian viscous flow. The consistency of classical viscous weakly compressible SPH formulations is investigated. In particular, the principle of virtual work is used to study the verification of the free-surface boundary conditions in a weak sense. The latter can be related to the global energy dissipation induced by the viscous term formulations and their consistency. Numerical verification of this theoretical analysis is provided on three free-surface test cases including a standing wave, with the three viscous term formulations investigated.

Figure 1

Configurations of the kernel support $\Omega \left(\mathit{r}\right)$ with respect to the fluid domain boundary.

Figure 2

Configuration (left) and SPH setup and initial pressure for the test case in Sec. 6a.

Figure 3

Evolution of a circular patch of fluid with a nonuniform initial vorticity distribution. Left: initial angular velocity $\omega$. Right: angular velocity $\omega$ at time $t=100\text{s}$ predicted by the SPH simulation using the MGF. The contour levels are representative of the intensity of the vorticity field.

Figure 4

Evolution of a circular patch of fluid with a nonuniform initial vorticity distribution. Comparison between the velocity component $u_{\theta }$ predicted by the SPH simulation using the MGF (solid line), and the analytical solution (dashed line) at different times.

Figure 5

Evolution of a circular patch of fluid with a nonuniform initial vorticity distribution. Left: time histories of the kinetic energy dissipation $d{E}_k/dt$ evaluated using ${\langle \nabla ·\mathbb{V}\rangle }^{PVW}$, ${\langle \nabla ·\mathbb{V}\rangle }^{\mathrm{MG}}$, ${\langle \nabla ·\mathbb{V}\rangle }^{\mathrm{MEA}}$. The dashed-dotted line represents the analytical solution. Right: relative error on the kinetic energy evolution between the analytical solution and SPH solutions for ${\langle \nabla ·\mathbb{V}\rangle }^{\mathrm{PVW}}$ and ${\langle \nabla ·\mathbb{V}\rangle }^{\mathrm{MG}}$.

Figure 6

Stretching of a square patch of viscous fluid. Left: free surface configuration at time $t=0$ and $t{A}_0=1$. Right: enlarged view of the $y$-component of the vector ${\langle \nabla ·\mathbb{V}\rangle }^{MG}$ evaluated along the symmetry axis AA’ at initial time $t=0$ for decreasing smoothing lengths $h$.

Figure 7

Stretching of a square patch of viscous fluid. Enlarged views of the $y$ component of the vectors ${\langle \nabla ·\mathbb{V}\rangle }^{\mathrm{PVW}}$ (left) and ${\langle \nabla ·\mathbb{V}\rangle }^{\mathrm{MEA}}$ (right) evaluated along the symmetry axis $A{A}^{\prime }$ at initial time $t=0$ for decreasing smoothing lengths $h$.

Figure 8

Stretching of a square patch of viscous fluid. Left: time histories of the kinetic energy dissipation $d{E}_k/dt$ for the three formulations. The dashed-dotted line represents the solution obtained through a finite-difference scheme solver. Right: relative error between the kinetic energy evaluated with the SPH MGF formulation and the one obtained by the FDM solver.

Figure 9

Standing wave: sketch of the problem (top left) and the decay of the kinetic energy as predicted by using the PVWF (top right), the MGF (bottom left), and the MEAF (bottom right).