A new 3D parallel SPH scheme for free surface flows

We propose a new robust and accurate SPH scheme, able to track correctly complex three-dimensional non-hydrostatic free surface flows and, even more important, also able to compute an accurate and little oscillatory pressure field. It uses the explicit third order TVD Runge-Kutta scheme in time, following Shu and Osher [Shu C-W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys 1988;89:439-71], together with the new key idea of introducing a monotone upwind flux for the density equation, thus removing any artificial viscosity term. For the discretization of the velocity equation, the non-diffusive central flux has been used. A new flexible approach to impose the boundary conditions at solid walls is also proposed. It can handle any moving rigid body with arbitrarily irregular geometry. It does neither produce oscillations in the fluid pressure in proximity of the interfaces, nor does it have a restrictive impact on the stability condition of the explicit time stepping method, unlike the repellent boundary forces of Monaghan [Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399-406]. To asses the accuracy of the new SPH scheme, a 3D mesh-convergence study is performed for the strongly deforming free surface in a 3D dam-break and impact-wave test problem providing very good results.Moreover, the parallelization of the new 3D SPH scheme has been carried out using the message passing interface (MPI) standard, together with a dynamic load balancing strategy to improve the computational efficiency of the scheme. Thus, simulations involving millions of particles can be run on modern massively parallel supercomputers, obtaining a very good performance, as confirmed by a speed-up analysis. The 3D applications consist of environmental flow problems, such as dam-break flows and impact flows against a wall. The numerical solutions obtained with our new 3D SPH code have been compared with either experimental results or with other numerical reference solutions, obtaining in all cases a very satisfactory agreement.     

Shock capturing by anisotropic diffusion oscillation reduction

This paper introduces an anisotropic diffusion oscillation reduction (ADOR) scheme for shock wave computations. The connection is made between digital image processing, in particular, image edge detection, and numerical shock capturing. Indeed, numerical shock capturing can be formulated on the lines of iterative digital edge detection. Various anisotropic diffusion and super diffusion operators originated from image edge detection are proposed for the treatment of hyperbolic conservation laws and near-hyperbolic hydrodynamic equations of change. The similarity between anisotropic diffusion and artificial viscosity is discussed. Physical origins and mathematical properties of the artificial viscosity are analyzed from the point of view of kinetic theory. A form of pressure tensor is derived from the first principles of the quantum mechanics. Quantum kinetic theory is utilized to arrive at macroscopic transport equations from the microscopic theory. Macroscopic symmetry is used to simplify pressure tensor expressions. The latter provides a basis for the design of artificial viscosity. The ADOR approach is validated by using (inviscid) Burgers' equation, the gas tube problems, the incompressible Navier-Stokes equation and the Euler equation. Both standard central difference schemes and a discrete singular convolution algorithm are utilized to illustrate the approach. Results are compared with those of third-order upwind scheme and essentially non-oscillatory (ENO) scheme.     

Free-surface flows solved by means of SPH schemes with numerical diffusive terms

A novel system of equations has been defined which contains diffusive terms in both the continuity and energy equations and, at the leading order, coincides with a standard weakly-compressible SPH scheme with artificial viscosity. A proper state equation is used to associate the internal energy variation to the pressure field and to increase the speed of sound when strong deformations/compressions of the fluid occur. The increase of the sound speed is associated to the shortening of the time integration step and, therefore, allows a larger accuracy during both breaking and impact events. Moreover, the diffusive terms allows reducing the high frequency numerical acoustic noise and smoothing the pressure field. Finally, an enhanced formulation for the second-order derivatives has been defined which is consistent and convergent all over the fluid domain and, therefore, permits to correctly model the diffusive terms up to the free surface. The model has been tested using different free surface flows clearly showing to be robust, efficient and accurate. An analysis of the CPU time cost and comparisons with the standard SPH scheme is provided.     

基于SPH方法的滴水涟漪动画模拟

光滑粒子流体动力学（Smoothed Particle Hydrodynamics,SPH）方法是一种新近发展的可用于流体模拟的无网格数值方法。文中基于SPH方法的基本原理,利用SPH方法求解描述水流现象的二维浅水波方程,根据具体模型使用Mon-aghan人工粘性的变形形式,有效地防止了相互靠近粒子的穿透,消除了SPH方法在模拟流体动力学问题时产生的数值振荡。通过使用可变光滑长度,使邻近粒子的数量保持相对稳定,提高了求解的计算效率和精度。同时,对光滑长度进行了修正以获取对称光滑长度,保持了粒子间相互作用对称性。全面考虑了各种定解条件的设置,对水滴的运动进行了模拟,SPH模拟结果与有限差分法、有限体积法结果非常吻合,验证了方法的准确性,为SPH方法的进一步发展和广泛运用奠定了基础。     

Efficient design sensitivity analysis of incompressible fluids using SPH projection method

Using the direct differentiation method, a design sensitivity analysis method for time-dependent incompressible fluids is developed. The fluid behavior is described as the motion of particles involved by the SPH method. In the SPH projection method, instead of changing the fluid density, incompressibility is enforced by the pressure Poisson equation derived from pressure projection, which enable to use larger time steps. In spite of the additional pressure Poisson equation, the computational cost for the design sensitivity is not expensive since the factorized system matrix of pressure Poisson equation can be utilized. Aforementioned computational efficiency is very beneficial for the design sensitivity computation required for every time step in explicit time integration and updated Lagrangian schemes, for which an update scheme of design velocity field is developed using the velocity sensitivity. Through demonstrative numerical examples, the developed DSA method turns out to be efficient and shows excellent agreement with finite differencing.     

An Essentially Non-oscillatory Crank–Nicolson Procedure for the Simulation of Convection-Dominated Flows

The Crank–Nicolson (CN) time-stepping procedure incorporating the second-order central spatial scheme is unconditionally stable and strictly non-dissipative for linear convection flows; however, its numerical solution in practice can be oscillatory for nonsmooth solutions. This article studies variants of the CN method for the simulation of linear convection-dominated diffusion flows, in which the explicit convection part is approximated by an upwind scheme, to effectively suppress nonphysical oscillations. The second-order essentially non-oscillatory scheme incorporated in the CN procedure (ENO-CN) has been found effective for a non-oscillatory numerical solution of minimum numerical dissipation. A stability analysis is provided for ENO-CN, which turns out to be unconditionally stable for problems of nonzero diffusion. However, for purely convective flows, it is stable only when the CFL condition is satisfied. Numerical results are presented to demonstrate its stability and accuracy.     

Smoothed particle hydrodynamics: Applications to heat conduction

In this paper, we modify the numerical steps involved in a smoothed particle hydrodynamics (SPH) simulation. Specifically, the second order partial differential equation (PDE) is decomposed into two first order PDEs. Using the ghost particle method, consistent estimation of near-boundary corrections for system variables is also accomplished. Here, we focus on SPH equations for heat conduction to verify our numerical scheme. Each particle carries a physical entity (here, this entity is temperature) and transfers it to neighboring particles, thus exhibiting the mesh-less nature of the SPH framework, which is potentially applicable to complex geometries and nanoscale heat transfer. We demonstrate here only 1D and 2D simulations because 3D codes are as simple to generate as 1D codes in the SPH framework. Our methodology can be extended to systems where the governing equations are described by PDEs.     

Improved Incompressible Smoothed Particle Hydrodynamics method for simulating flow around bluff bodies

In this article, we present numerical solutions for flow over an airfoil and a square obstacle using Incompressible Smoothed Particle Hydrodynamics (ISPH) method with an improved solid boundary treatment approach, referred to as the Multiple Boundary Tangents (MBT) method. It was shown that the MBT boundary treatment technique is very effective for tackling boundaries of complex shapes. Also, we have proposed the usage of the repulsive component of the Lennard-Jones Potential (LJP) in the advection equation to repair particle fractures occurring in the SPH method due to the tendency of SPH particles to follow the stream line trajectory. This approach is named as the artificial particle displacement method. Numerical results suggest that the improved ISPH method which is consisting of the MBT method, artificial particle displacement and the corrective SPH discretization scheme enables one to obtain very stable and robust SPH simulations. The square obstacle and NACA airfoil geometry with the angle of attacks between 0° and 15° were simulated in a laminar flow field with relatively high Reynolds numbers. We illustrated that the improved ISPH method is able to capture the complex physics of bluff-body flows naturally such as the flow separation, wake formation at the trailing edge, and the vortex shedding. The SPH results are validated with a mesh-dependent Finite Element Method (FEM) and excellent agreements among the results were observed.     

Investigation of wall bounded flows using SPH and the unified semi-analytical wall boundary conditions

Semi-analytical wall boundary conditions present a mathematically rigorous framework to prescribe the influence of solid walls in smoothed particle hydrodynamics (SPH) for fluid flows. In this paper they are investigated with respect to the skew-adjoint property which implies exact energy conservation. It will be shown that this property holds only in the limit of the continuous SPH approximation, whereas in the discrete SPH formulation it is only approximately true, leading to numerical noise. This noise, interpreted as a form of “turbulence”, is treated using an additional volume diffusion term in the continuity equation which we show is equivalent to an approximate Riemann solver. Subsequently two extensions to the boundary conditions are presented. The first dealing with a variable driving force when imposing a volume flux in a periodic flow and the second showing a generalization of the wall boundary condition to Robin type and arbitrary-order interpolation. Two modifications for free-surface flows are presented for the volume diffusion term as well as the wall boundary condition. In order to validate the theoretical constructs numerical experiments are performed showing that the present volume flux term yields results with an error 5 orders of magnitude smaller then previous methods while the Robin boundary conditions are imposed correctly with an error depending on the order of the approximation. Furthermore, the proposed modifications for free-surface flows improve the behavior at the intersection of free surface and wall as well as prevent free-surface detachment when using the volume diffusion term. Finally, this paper is concluded by a simulation of a dam break over a wedge demonstrating the improvements proposed in this paper.     

Numerical Methods with Adaptive Artificial Viscosity for Solving Navier−Stokes Equations

A numerical method for solving two-dimensional problems of a viscous compressible gas based on Navier–Stokes equations with the introduction of adaptive artificial viscosity is presented. The proposed method is implemented for areas of the general form on triangular grids. The method of the adaptive artificial viscosity is taken as the basis of the proposed numerical method and ensures the monotonicity of the solutions, even in the presence of shock waves. The artificial viscosity (introduced into the difference scheme) is constructed in such a way that it is absent in the boundary layer where the dynamic viscosity acts. The viscosity is determined from the conditions of the fulfillment of the maximum principle. The series of calculations of an external flow around a cylinder for various Reynolds and Mach numbers is described.     

A novel method for modeling Neumann and Robin boundary conditions in smoothed particle hydrodynamics

We present a novel smoothed particle hydrodynamics (SPH) method for diffusion equations subject to Neumann and Robin boundary conditions. The Neumann and Robin boundary conditions are common to many physical problems (such as heat/mass transfer), and can prove challenging to implement in numerical methods when the boundary geometry is complex. The new method presented here is based on the approximation of the sharp boundary with a diffuse interface and allows an efficient implementation of the Neumann and Robin boundary conditions in the SPH method. The paper discusses the details of the method and the criteria for the width of the diffuse interface. The method is used to simulate diffusion and reactions in a domain bounded by two concentric circles and reactive flow between two parallel plates and its accuracy is demonstrated through comparison with analytical and finite difference solutions. To further illustrate the capabilities of the model, a reactive flow in a porous medium was simulated and good convergence properties of the model are demonstrated.     

Particle packing algorithm for SPH schemes

Using some intrinsic features of the Smoothed Particle Hydrodynamics (SPH) schemes, an innovative algorithm for the initialization of the particle distribution has been defined. The proposed particle packing algorithm allows a drastic reduction of the numerical noise due to particle resettlement during the early stages of the flow evolution. Moreover, thanks to its structure, it can be easily derived starting from whatever SPH scheme and applies under the hypotheses that the fluid is weakly-compressible or incompressible as well. A broad range of numerical test cases proved this tool to be fast, robust and reliable also for complex geometrical configurations.     

A Lattice BGK Scheme with General Propagation

The standard Lattice BGK (LBGK) scheme often encounters numerical instability in simulation of fluid flow with small kinematic viscosity or as the nondimensional relaxation time is close to 0.5. In this paper, based on a time-splitting scheme for the Boltzmann equation with discrete velocities, a new LBGK scheme with general propagation step is proposed to address this problem. In this model, two free parameters are introduced into the propagation step, which can be adjusted to obtain a small kinematic viscosity and improved numerical stability as well. Numerical simulations of the two-dimensional Taylor vortex and the unsteady Womersley flow are carried out to test the kinematic viscosity, numerical diffusion, and numerical stability of the proposed scheme.     

Numerical diffusive terms in weakly-compressible SPH schemes

A discussion on the use of numerical diffusive terms in SPH models is proposed. Such terms are, generally, added in the continuity equation, in order to reduce the spurious numerical noise that affects the density and pressure fields in weakly-compressible SPH schemes. Specific focus has been given to the theoretical analysis of the diffusive term structure, highlighting the main benefits and drawbacks of the most widespread formulations. Finally, specific test cases have been used to compare such formulations and to confirm the theoretical findings.     

Real-time Simulation of Gas Based on Smoothed Particle Hydrodynamics

This paper extends the SPH method to gas simulation. The SPH(Smoothed Particles Hydrodynamics) method is the most popular method of flow simulation, which is widely used in large-scale liquid simulation. However, it is not found to apply to gas simulation, since those methods based on SPH can’t be used in real-time simulation due to their enormous particles and huge computation. This paper proposes a method for gas simulation based on SPH with a small number of particles. Firstly, the method computes the position and density of each particle in each point-in-time, and outlines the shape of the simulated gas based on those particles. Secondly the method uses the grid technique to refine the shape with the diffusion of particle’s density under the control of grid, and get more lifelike simulation result. Each grid will be assigned density according to the particles in it. The density determines the final appearance of the grid. For ensuring the natural transition of the color between adjacent grids, we give a diffuse process of density between these grids and assign appropriate values to vertexes of these grids. The experimental results show that the proposed method can give better gas simulation and meet the request of real-time.     

On the SPH tensile instability in forming viscous liquid drops

Smoothed Particle Hydrodynamics (SPH) simulations of elastic solids and viscous fluids may suffer from unphysical clustering of particles due to the tensile instability. Recent work has shown that in simulations of elastic or brittle solids the instability can be removed by an artificial stress whose form is derived from a linear perturbation analysis of the full set of governing SPH equations. While a linear analysis cannot be used to derive the corresponding form of the artificial stress for a viscous fluid, here we show that the same construction which applies to elastic solids may also work for viscous fluids provided that the constant parameter ? entering in the definition of the artificial stress is properly chosen. As a suitable test case, we model the formation of a circular van der Waals liquid drop and show that the tensile instability is removed when an artificial viscous force and energy generation term are added to the standard SPH equations of motion and energy, respectively. The optimal value of the constant ? is constrained by the ability of the model simulation to reproduce both a sufficiently smoothed density profile and the van der Waals phase diagram.     

A two-steps time discretization scheme using the SPH method for shock wave propagation

The Smoothed Particle Hydrodynamics (SPH) is a meshfree method which has been applied to a wide range of problems. In the present work, a new time integration algorithm using a corrected SPH spatial discretization for small deformations is applied to solve the propagation of shock waves in viscoplastic continua. In the method presented herein the equations are formulated in terms of stress and velocity. A corrected Lagrangian kernel is employed and two different sets of particles are used for the time discretization. Numerical instabilities are not present when using this new SPH formulation. The method proposed here has been proved to be efficient and it provides solutions of good accuracy.     

基于粒子分解的SPH并行算法研究与应用

作为一种典型的拉格朗日型无网格数值方法,光滑粒子流体动力学(SPH)方法在模拟自由表面流问题时具有天然优势。但是,该方法计算量大、耗时长,为此提出了一种基于粒子分解的SPH并行算法。该算法将所有粒子平均分配到各个进程进行计算,每个时间步通信仅调用一次发送、接收和广播函数,因此易于实现且可扩展性较好。应用该并行算法对二维溃坝流和三维液滴冲击液膜问题进行数值模拟,结果表明：该并行算法能显著减少模拟所消耗的计算时间,有利于进行三维大规模计算问题的数值模拟;当粒子数大于百万时,最大加速比可达30以上。     

Navier-Stokes calculations with a coupled strongly implicit method—I: Finite-difference solutions

Stone's unconditionally stable, strongly implicit numerical method is extended to the 2 x 2 coupled vorticity-stream function form of the Navier-Stokes equations. The solution algorithm allows for complete coupling of the boundary conditions. Solution for arbitrary large time steps, and for cell Reynolds numbers much greater than two have been obtained. The method converges quite rapidly without adding artificial viscosity or the necessity for under-relaxation. This technique is used here to solve for a variety of internal and external flow problems. Moderate to large Reynolds numbers are considered for both separated and unseparated flows. The procedure is extended to higher-order splines in Part 2 of this study.     

Fourth-order central compact scheme for the numerical solution of incompressible Navier–Stokes equations

This paper provides an implicit central compact scheme for the numerical solution of incompressible Navier–Stokes equations. The solution procedure is based on the artificial compressibility method that transforms the governing equations into a hyperbolic-parabolic form. A fourth-order central compact scheme with a sixth-order numerical filtering is used for the discretization of convective terms and fourth-order central compact scheme for the viscous terms. Dual-time stepping approach is applied to time discretization with backward Euler difference scheme to the pseudo-time derivative, and three point second-order backward difference scheme to the physical time derivative. An approximate factorization-based alternating direction implicit scheme is used to solve the resulting block tridiagonal system of equations. The accuracy and efficiency of the proposed numerical method is verified by simulating several two-dimensional steady and unsteady benchmark problems.