Oscillating drop under a central force field
The test case of oscillating drop under a central force field is a widely considered benchmark test case for validation of stability, accuracy, (energy and volume) conservation properties and convergence of particle methods. This test was originally proposed by Monaghan and Rafiee (2013) and has been considered in many single-phase (e.g., Antuono et al., 2015; Khayyer et al., 2017a,b) as well as multi-phase SPH simulations (e.g., Lind et al., 2016; Khayyer et al., 2019). Due to presence of simple boundary conditions (e.g., absence of wall boundaries) and existence of an exact theoretical solution, this classical test case is recommended as a basic benchmark for validation of SPH codes, especially for long-duration validations.
Single-Phase or Two-Phase Flow
The geometry is a 2D initially circular drop of radius R that comprises concentric regions of fluids with different densities, where the lighter fluid surrounds the heavier one. The heavier fluid has an initial radius of r with R = z r.
No specific boundary condition needs to be imposed explicitly. Kinematic free-surface boundary condition is automatically satisfied by Lagrangian meshfree nature of SPH. Dynamic free-surface boundary condition should also be satisfied if a proper SPH scheme is implemented so that the density at the free-surface remains invariant, resulting in zero gauge pressure at the free-surface.
The initial setup is presented in Fig. 1 for a multi-phase (two-fluid) case. A drop of an outer radius R comprises concentric elliptical regions of fluids with different densities, where the lighter fluid surrounds the heavier one. The fluids are assumed as incompressible and inviscid. The drop is subjected to the following initial velocity field:
u = σ x; v = - σ y
where (x, y) and (u, v) are the components of the position vector r = (x, y) and velocity vector u = (u, v), respectively, and s is a constant value. The drop is also subjected to a conservative force field corresponding to a potential field of ½(Ωr)², where Ω is a dimensional parameter and r is the position vector.
Fig. 1 Schematic sketch of the benchmark test – multi-phase case
Initial pressure distribution
In simulations by projection methods, imposition of initial pressure field would not be necessary because a global Poisson Pressure Equation (PPE) would be solved along with explicit imposition of dynamic free-surface boundary condition. For weakly compressible SPH or explicit SPH simulations, however, it is recommended to set the initial pressure field with respect to analytical solution to ensure minimization of errors associated with the initially obtained numerical pressure field (as highlighted by Hammani et al., 2020). According to Monaghan and Rafiee (2013), the pressure in the outer and inner fluids, Pₒ and Pᵢ, must be a quadratic function of the coordinates as follows:
where Qₒ, Qᵢ and C are constants. From the fact that the pressure at the boundary of inner and outer fluids must be the same, we have the following relation:
Also, for both outer and inner fluids, we have the following relations of σ and Ω:
In case of the initial state (t = 0 s), we have the relations of A = B = R and a = b = r, and therefore:
Therefore, for the outer fluid, we have the initial pressure field of:
and for the inner fluid, the initial pressure field would be expressed as:
If a single-phase calculation is considered, then the initial pressure field would be:
Outer and inner elliptical surfaces are expressed as:
where a, A, b and B are inner semi-major, outer semi-major, inner semi-minor and outer semi-minor axes, respectively. For these variables, we have the following relation:
Time evolutions of a, A, b and B are expressed as:
where t is time; the variable σ is evolved in time as:
Theoretical solutions for semi-major and semi-minor axes can be obtained by simultaneously evolving the above equations (Eqs. 1 and 2) in time.
The theoretical solutions for kinetic and potential energies (EK and EP) are expressed as follows:
where stands for the ratio between A to a or B to b (i.e., A = ζ a ; B = ζ b ; ζ is usually set to 2 for multi-phase case); and are densities of inner and outer fluids. The time evolutions of kinetic and potential energies can be obtained by inserting the semi-major and semi-minor axes, A and B, into Eq. (3) in time.
For a single-phase calculation (Fig. 2), densities of both inner and outer regions can be set the same resulting in a single-fluid drop. Accordingly, time evolutions of A, B and σ are obtained by solving the following Eqs. (4) and (5):
The time variations of kinetic and potential energies for single-phase flow can be obtained from Eq. (3) with ρᵢ and ρₒ begin equal to ρ (ρᵢ = ρₒ = ρ), i.e.,
Fig. 2 Schematic sketch of the benchmark test – single-phase case
For a drop of radius R = 1.0 m, an initial setup of particle spacing can be considered to be dₒ = 0.01 m. Certainly, as previously stated, this test serves as an excellent benchmark to investigate convergence property of the model and thus, simulations are suggested to be conducted by considering finer and coarser initial particle spacings.
Validations are suggested to be made with respect to theoretical solutions corresponding to time variations of semi-major and semi-minor axes as well as time histories of energy components and total energy. Certainly, to assess the local volume conservation of an SPH scheme, investigations can be made regarding spatial distribution of velocity divergence. Global volume conservation can be verified by simultaneous consideration of evolution of semi-major and semi-minor axes.
An excel file containing the theoretical solutions in terms of time evolution of semi-minor axis, semi-major axis and energy components for the single-phase case is included as a supplementary material of this test. The corresponding Fortran file to generate it is also provided. Another Fortran file for generation of theoretical solution corresponding to multi-phase case is also provided.
In this section, some typical numerical results corresponding to this benchmark test are provided for both single-phase and multi-phase configurations. Fig. 3 presents typical results for the single-phase version by improved MPS and ISPH discussed in the paper by Khayyer et al. (2017b), in specific, snapshots of particles and pressure field by improved MPS (Fig. 3a) as well as time history of semi-major axis A by improved ISPH (Fig. 3b). For the simulations performed in this paper, s and W were set as 0.4 and 1.2, respectively. Simulations were conducted by considering the CFL stability condition as well as a maximum allowable time step set as Dtₘₐₓ = 0.001 s. The diameter of particles or particle spacing was considered as dₒ = 0.005 m.
Fig. 4 presents typical results for the multi-phase case by improved multi-phase MPS incorporating the Optimized Particle Shifting (OPS) scheme (improved MPS + OPS, Khayyer et al. 2019). In particular, Fig. 4b shows a typical snapshot of particles and density field by improved MPS+OPS. Fig. 4b portrays time history of outer semi-major axis A by the same method. For this simulation, s and Ω were set equal to 0.5 and 1.0, respectively. A density ratio of 1 to 1000 is considered (ρₕ/ρₗ = 1000). The particles are 1 cm in diameter (dₒ = 0.01 m). The calculation time step is set according to the CFL stability condition and a maximum allowable time step of Dtₘₐₓ = 1.0E-4 s.
Fig. 3 Typical results for single-phase case by improved MPS and ISPH (Khayyer et al., 2017b)
Fig. 4 Typical results for multi-phase case by improved multi-phase MPS incorporating Optimised Particle Shifting Scheme, referred to as Improved MPS+OPS (Khayyer et al., 2019)
You can download the full test case below:
Monaghan, J.J., Rafiee, A., 2013. A simple SPH algorithm for multi-fluid flow with high density ratios. International Journal for Numerical Methods in Fluids 71, 537–561. [Link]
Antuono, M., Marrone, S., Colagrossi, A., & Bouscasse, B., 2015. Energy balance in the δ-SPH scheme. Computer Methods in Applied Mechanics and Engineering, 289, 209-226. [Link]
Khayyer, A., Gotoh, H., Shimizu, Y., 2017a. Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context. Journal of Computational Physics 332, 236–256. [Link]
Khayyer, A., Gotoh, H., Shimizu, Y., Gotoh, K., 2017b. On enhancement of energy conservation properties of projection-based particle methods. European Journal of Mechanics - B/Fluids 66, 20–37. [Link]
Khayyer, A., Gotoh, H., Shimizu, Y., 2019. A projection-based particle method with optimized particle shifting for multiphase flows with large density ratios and discontinuous density fields. Computers & Fluids 179, 356–371. [
Khayyer, A., Gotoh, H., Shimizu, Y., 2019. A projection-based particle method with optimized particle shifting for multiphase flows with large density ratios and discontinuous density fields. Computers & Fluids 179, 356–371. [Link]
Hammani, I., Marrone, S., Colagrossi, A., Oger, G., Le Touzé, D., 2020. Detailed study on the extension of the δ-SPH model to multi-phase flow. Computer Methods in Applied Mechanics and Engineering 368, 113189. [Link]
Other publications considering this benchmark test:
You, Y., Khayyer, A., Zheng, X., Gotoh, H., Ma, Q., 2021. Enhancement of δ-SPH for ocean engineering applications through incorporation of a background mesh scheme. Applied Ocean Research 110, 102508. [Link]
Sun, P.N., Colagrossi, A., Marrone, S., Zhang, A.M., 2017. The δplus-SPH model: Simple procedures for a further improvement of the SPH scheme. Computer Methods in Applied Mechanics and Engineering 315, 25–49. [Link]
Wang, P.-P., Meng, Z.-F., Zhang, A.-M., Ming, F.-R., Sun, P.-N., 2019. Improved particle shifting technology and optimized free-surface detection method for free-surface flows in smoothed particle hydrodynamics. Computer Methods in Applied Mechanics and Engineering 357, 112580. [Link]
Lyu, H.-G., Sun, P.-N., 2022. Further enhancement of the particle shifting technique: Towards better volume conservation and particle distribution in SPH simulations of violent free-surface flows. Applied Mathematical Modelling 101, 214–238. [Link]
Meringolo, D.D., Colagrossi, A., Marrone, S., Aristodemo, F., 2017. On the filtering of acoustic components in weakly-compressible SPH simulations. Journal of Fluids and Structures 70, 1–23. [Link]
Krimi, A., Jandaghian, M., Shakibaeinia, A., 2020. A WCSPH Particle Shifting Strategy for Simulating Violent Free Surface Flows. Water 12, 3189. [Link]
Jandaghian, M., Siaben, H.M., Shakibaeinia, A., 2022. Stability and accuracy of the weakly compressible SPH with particle regularization techniques. European Journal of Mechanics - B/Fluids. [Link]
Rezavand, M., Taeibi-Rahni, M., Rauch, W., 2018. An ISPH scheme for numerical simulation of multiphase flows with complex interfaces and high density ratios. Computers & Mathematics with Applications 75, 2658–2677. [Link]
Shimizu, Y., Gotoh, H., Khayyer, A., 2018. An MPS-based particle method for simulation of multiphase flows characterized by high density ratios by incorporation of space potential particle concept. Computers & Mathematics with Applications 76, 1108–1129. [Link]
Meng, Z.-F., Wang, P.-P., Zhang, A.-M., Ming, F.-R., Sun, P.-N., 2020. A multiphase SPH model based on Roe’s approximate Riemann solver for hydraulic flows with complex interface. Computer Methods in Applied Mechanics and Engineering 365, 112999. [Link]