An Anti-Clustering Model for Stability Enhancement of a 3D Moving Particle Semi-Implicit Method and Two-Phase Coupling between MPS and Euler Grids

Abstract

As a Lagrangian gridless particle method, the MPS (Moving Particle Semi-implicit) method has a wide engineering application. However, for complex 3D flows, unphysical pressure oscillations often occur and result in the failure of simulations. This paper compares the stability enhancement methods proposed by different researchers to develop a 3D, stable MPS method. The results indicate that the proposed methods are incapable of eliminating the particle clustering that leads to instability as the main source in coarser particle spacing cases. An anti-clustering model, referring to the SPH (Smoothed Particle Hydrodynamics) artificial viscosity model, is proposed to further reduce instability. Combining various proposed methods and models, several typical examples are simulated comparatively. The results are compared with those of the VOF (Volume of Fluid) model using commercial software to validate the accuracy and stability of the combination of the proposed methods and models. It is concluded that (1) 3D cases that adopt a high-order Laplacian model and high-order source terms in PPE are more accurate than those adopting the low-order operators; (2) the proposed anti-clustering model can produce a tuned interparticle force to prevent particle clustering and introduce no additional viscosity effects in the flow of the normal state, which plays a very positive role for further stability enhancement of MPS; (3) particle resolution significantly maintains simulation accuracy given the stable algorithms by the combination of stability enhancement methods. The 3D MPS method is coupled with the Euler grid (FLUENT V17 software, ANSYS, Pittsburgh, PA, USA) in two phases. In particular, the 3D MPS algorithm is used to calculate the liquid-phase change from the continuous to the dispersed, and the finite volume method based on the Euler grid is adopted to measure the corresponding gas-phase motion. The atomization of the liquid jet under static air flow is calculated and compared with the results of the VOF method, which can capture the continuous interface.

1. Introduction

The MPS (Moving Particle Semi-implicit) method is a kind of meshless method resolving incompressible fluid dynamics, firstly proposed by Koshizuka et al. [1] in 1995, by which the governing N–S (Navier–Stokes) equations are described in a fully Lagrange frame, and the fluid is discretized into particles and then directly tracked in the Lagrange frame. By the virtue of this method, no numerical diffusion can be produced in the solution since the advection term has disappeared in the governing equations and no complicated mesh generation is needed due to its meshless nature. Therefore, the MPS method has been applied in a wide range of engineering fields in the past years, especially in dealing with the fluid problem with violent free-surface deformation [2,3,4].

However, with the application of the MPS method, some built-in problems of this method came out gradually, especially its instability. The instability of the MPS method mainly lies in the spurious pressure fluctuations in the solution of the pressure Poisson equations, leading to an unphysical velocity and divergence of the solution. In the literature, six kinds of techniques to improve this instability have been proposed.

The first is just to average the pressure fluctuations in space and time simply, which was proposed by Sueyoshi and Naito [5]. Although it produces some effects in view of the result, it is not an ideal method to tackle the problem. The second technique involves the modification of the kernel function. In fact, Ashtiani and Farhadi [6] has studied the effects of kernel functions on stability by comparing a variety of kernel functions; they suggest one from SPH (Smooth Particle Hydrodynamics), which has some effects to improve the instability. Similar works were done by Pan et al. [7]. The third technique focuses on the solution of the PPE (Poisson Pressure Equation) [8,9,10]. Masayuki et al. [11] introduced a Quasi-Compressibility term in PPE to stabilize the pressure solution. Kondo and Koshizuka [12] derived a more generic formulation of the source term for PPE, including one main part and two error-compensating parts. It was proved that with proper selection of the coefficients for the error-compensating parts, the unphysical pressure oscillation can be suppressed. With more attention on the PPE solution, Khayyer and Gotoh [13,14,15] proposed a higher-order Laplacian model and a higher-order source term for PPE recently. They discussed the instability of the MPS-based simulations in presence of attractive interparticle forces, which is similar to the so-called tensile instability in the SPH method, and proposed that the modified PPE and a corrective matrix for the pressure gradient model are helpful to stabilize and enhance the performance of the MPS method. The fourth technique is to modify the pressure gradient model [16,17]. The first modification of the pressure gradient model was proposed by Koshizuka et al. [18,19], in which the pressure of the target particle in the pressure gradient model is replaced by the minimum pressure value in the influence circle of the particle. This modification is tested to be able to improve the numerical stability. Based on this modification, Khayyer and Gotoh [15] further proposed a new pressure gradient model with full anti-symmetry, which conserves the linear momentum exactly. Besides, based on the 2nd-order Taylor series expansion of pressure, a corrective matrix for the pressure gradient model was also proposed by Khayyer and Gotoh [15] to enhance the exactness and stability of the pressure solution. The fifth technique presents a scheme for Dynamic Stabilization [20,21], which meticulously reproduces adequate repulsive forces to attenuate the interparticle penetration and thus stabilizes the calculations. The sixth techniques is based on the regularization techniques (such as particle shifting and collision) [22,23,24,25,26]. On the other hand, the stability algorithm and application of SPH [27,28,29], which is also a meshless method, attracted the authors’ attention, such as the high accuracy pressure model proposed by Zheng et al. [30] and artificial viscosity model used by Monaghan et al. and Daniel et al. [31,32,33].

As addressed previously, although many techniques have been developed to enhance the stability of MPS, there still lacks a complete method to essentially resolve this problem. In the present investigation, several selected stabilization enhancement methods are tested and extended to the 3D case, and their adaptability in 2D and 3D in relatively coarse particle spacing cases are tested firstly. However, we found that the instability also occurs due to particle clustering. Thus, an anti-clustering model referring to the SPH artificial viscosity model is proposed to further enhance the stability in the 3D case. Encouraging results are obtained.

The atomization mechanism of a liquid jet under supersonic flow [34] is very complicated. Direct numerical observations of the spray process under this condition contribute to the establishment of a general theoretical model. This paper develops an integrated calculation method to include the atomization process. The advantage of the MPS method lies in its treatment on fluid interfaces, which simulates the entire process from the continuous to the discrete interfaces. Therefore, the two-phase coupling between the MPS and the Euler grid (FLUENT V17 software) is adopted to calculate the liquid water jetting in the boundary layer of the static gas flow and explore the atomization of the liquid jet. This calculation is compared to the result of the VOF method that captures the continuous interface. The method provides a feasible path for the integrated simulation of liquid jet atomization in supersonic flow.

2. Standard MPS Method and Stabilization Methods Proposed Previously

2.1. Governing Equation

The governing equations are incompressible Navier–Stokes equations, listed as follows:

$$\frac{D\rho }{Dt}=0$$

(1)

$$\frac{Du}{Dt}=-\frac{\nabla p}{\rho }+\upsilon {\nabla }^{2}u+F$$

(2)

where ρ is the fluid density, u is velocity, p is pressure, υ is dynamic viscosity, t is time, and F is the sum of the forces, including the volume forces and surface tension, etc.

2.2. Kernel Function and Particle Interaction Models

The original kernel function of the MPS method [19] is

$$w\left(r\right)=\begin{array}{ll}\frac{r_e}{r^{}}-1& 0\le r<r_e\\ 0& r>r_e\end{array}$$

(3)

where r_e is the radius of the influence region of each particle, r = |r_j − r_i|.

The gradient operator and the Laplacian model are listed as follows:

$${〈\nabla \phi 〉}_i=\frac{d^{’}}{n_0}{\displaystyle \sum _{j\ne i}\frac{{\phi }_j-{\phi }_i{}^{}}{{\left|r_j-r_i\right|}^2}}\left(r_j-r_i\right)w\left(\left|r_j-r_i\right|,r_e\right)$$

(4)

$${〈{\nabla }^2\phi 〉}_i=\frac{2d^{\prime }}{n_0\lambda }{\displaystyle \sum _{j\ne i}\left({\phi }_j-{\phi }_i\right)w\left(\left|r_j-r_i\right|,r_e\right)}$$

(5)

where $\lambda =\frac{{\displaystyle \sum _{j\ne i}w\left(\left|r_j-r_i\right|,r_e\right){\left|r_j-r_i\right|}^2}}{{\displaystyle \sum _{j\ne i}w\left(\left|r_j-r_i\right|,r_e\right)}}$, and φ can be either a scalar or a vector physical quantity, ${\phi }_{ij}={\phi }_j-{\phi }_i$. d’ is the number of space dimensions, and n₀ is the constant particle number density $n_0={\displaystyle \sum _{j\ne i}{w}_{ij}}$. The gradient operator (Equation (4)) is related to the non-exact conservation of momentum as the resulting interparticle pressure forces are not anti-symmetric (equal in magnitude, opposite in direction). However, the numerical errors arising from a non-exact conservation of momentum might be quite less than those corresponding to the approximation of pressure gradient, especially in the presence of negative pressure fields [15].

2.3. Several Methods for Stability Enhancement Proposed Previously

Method 1: The Laplacian operator of pressure in PPE is adopted as Equation (6), while the source term is adopted as Equation (7); the pressure gradient model is adopted as Equation (8).

$${〈{\nabla }^2\phi 〉}_i=\frac{2d^{\prime }}{n_0}{\displaystyle \sum _{j\ne i}\frac{{\phi }_j-{\phi }_i}{{\left|r_j-r_i\right|}^2}}w(\left|r_j-r_i\right|,r_e)$$

(6)

$$\frac{Dn}{Dt}=\frac{1}{n_0}\frac{n_{*}-n_0}{\mathsf{\Delta }t}$$

(7)

$${〈\nabla p〉}_i^{}=\frac{d^{\prime }}{n_0}{\displaystyle \sum _{j\ne i}\frac{p_j-p_{\min }{}^{}}{{\left|r_j-r_i\right|}^2}\left(r_j-r_i\right)w\left(\left|r_j-r_i\right|,r_e\right)}$$

(8)

where $p_{\min }=\min \left(p_j\right),for\left\{j\left|w\left(\left|r_j-r_i\right|\right)\ne 0\right.\right\}$.

For the convenience of addressing, the MPS case using Equation (6) in PPE is referred to as MPS-IL (Initial Laplacian model), while the one using Equation (8) [35] as the gradient model is referred to as MPS-GM (Gradient model with Minimum pressure) in this paper.

Method 2: The Laplacian operator of pressure in PPE is adopted as Equation (10), while the source term is adopted as Equation (9) or Equation (11), and the pressure gradient model is adopted as Equation (12). Also for convenience of addressing, the MPS case using Equation (10) in PPE is referred to as MPS-HL (Higher-order Laplacion model) [15], while one using Equation (9) as the source term of PPE is referred to as MPS-HS (Higher source term); the one using Equation (11) as the source term is referred to as MPS-HS-ECS (Error Compensating parts in the source term) [15], and one using Equation (12) as the source term is referred to as MPS-GC (Gradient Correction) [15]. For the 2D case [15], the formulas are listed as follows:

MPS-HS:

$$\frac{Dn}{Dt}={\displaystyle \sum _{i\ne j}\frac{D\mathrm{w}(|r_j-r_i|)}{Dt}={\displaystyle \sum _{i\ne j}\frac{D{w}_{i\mathrm{j}}}{Dt}}}={\displaystyle \sum _{i\ne j}\frac{\partial w_{ij}}{\partial r_{ij}}\left(\frac{x_{ij}}{r_{ij}}{u}^{*}{}_{ij}+\frac{y_{ij}}{r_{ij}}{v}^{*}{}_{ij}\right)}$$

(9)

MPS-HL:

$$\nabla {〈\nabla {\phi }_{ij}〉}_i=\nabla {〈{\phi }_{ij}.\nabla w_{ij}〉}_i=\frac{1}{n_0}{\displaystyle \sum _{i\ne j}\left({\phi }_{ij}.\frac{{\partial }^2{w}_{ij}}{\partial r_{ij}^2}-\frac{{\phi }_{ij}}{r_{ij}}\frac{\partial w_{ij}}{\partial r_{ij}}\right)}$$

(10)

MPS-HS-ECS:

$$\{\begin{array}{cc}\frac{\mathsf{\Delta }t}{\rho }({\nabla }^{2}Pk+1\hfill & )i=\frac{1}{n_0}{\left(\frac{Dn}{Dt}\right)}_i^{*}+\left|\left(\frac{n^k-n_0}{n_0}\right)\right|\left[\frac{1}{n_0}{\left(\frac{Dn}{Dt}\right)}_i^k\right]\hfill \\ & +\left|\left(\frac{\mathsf{\Delta }t}{n_0}{\left(\frac{Dn}{Dt}\right)}_i^k\right)\right|\left[\frac{1}{\mathsf{\Delta }t}\left(\frac{n^k-n_0}{n_0}\right)\right] {\mathrm{n}}^{\mathrm{k}}> \gamma {\mathrm{n}}_0\hfill \\ \frac{\mathsf{\Delta }t}{\rho }({\nabla }^{2}Pk+1\hfill & )i=\frac{1}{n_0}{\left(\frac{Dn}{Dt}\right)}_i^{*}+\left|\left(\frac{n^k-{\left(n_i\right)}_s}{{\left(n_i\right)}_s}\right)\right|\left[\frac{1}{n_0}{\left(\frac{Dn}{Dt}\right)}_i^k\right]\hfill \\ & +\left|\left(\frac{\mathsf{\Delta }t}{n_0}{\left(\frac{Dn}{Dt}\right)}_i^k\right)\right|\left[\frac{1}{\mathsf{\Delta }t}\left(\frac{n^k-{\left(n_i\right)}_s}{{\left(n_i\right)}_s}\right)\right] {\mathrm{n}}^{\mathrm{k}}\le \gamma {\mathrm{n}}_0\hfill \end{array}$$

(11)

MPS-GC:

$${〈\nabla P〉}_i=\frac{d^{’}}{n_0}{\displaystyle \sum _{}\frac{p_j-p_i}{{\left|r_j-r_i\right|}^2}\left(r_j-r_i\right)C_{ij}w\left(\left|r_j-r_i\right|\right)}$$

(12)

3. The 3D Extending of the Proposed Stability Enhancement Methods and an Anti-Clustering Model

3.1. 3D Higher-Order Laplacian Model for PPE

Based on the 2D higher-order Laplacian model proposed by Khayyer and Gotoh [13,15], a 3D higher-order Laplacian model is derived as

$$\begin{array}{l}\nabla {〈\nabla {\phi }_{ij}〉}_i=\nabla {〈{\phi }_{ij}.\nabla w_{ij}〉}_i\\ =\frac{1}{n_0}{\displaystyle \sum _{i\ne j}\left(\nabla {\phi }_{ij}.\nabla w_{ij}+{\phi }_{ij}.{\nabla }^2{w}_{ij}\right)}\\ =\frac{1}{n_0}{\displaystyle \sum _{i\ne j}\left(\frac{2{\phi }_{ji}}{r_{ij}}\frac{\partial w_{ij}}{\partial r_{ij}}+{\phi }_{ij}.\frac{{\partial }^2{w}_{ij}}{\partial r_{ij}^2}+\frac{2{\phi }_{ij}}{r_{ij}}\frac{\partial w_{ij}}{\partial r_{ij}}\right)}\\ =\frac{1}{n_0}{\displaystyle \sum _{i\ne j}\left({\phi }_{ij}.\frac{{\partial }^2{w}_{ij}}{\partial r_{ij}^2}\right)}\end{array}$$

(13)

3.2. 3D Higher-Order Source Term for PPE

The 3D higher-order source term for PPE was also derived as

$$\frac{Dn}{Dt}={\displaystyle \sum _{i\ne j}\frac{D\mathrm{w}(|r_j-r_i|)}{Dt}={\displaystyle \sum _{i\ne j}\frac{D{w}_{i\mathrm{j}}}{Dt}}}={\displaystyle \sum _{i\ne j}\frac{\partial w_{ij}}{\partial r_{ij}}\left(\frac{x_{ij}}{r_{ij}}{u}^{*}{}_{ij}+\frac{y_{ij}}{r_{ij}}{v}^{*}{}_{ij}+\frac{z_{ij}}{r_{ij}}{w}^{*}{}_{ij}\right)}$$

(14)

3.3. 3D Correction Matrix for the Pressure Gradient Model

As discussed by Khayyer and Gotoh [15], an accurate pressure gradient lead to the more accurate motions of particles, therefore reducing the disorder and confusion of the particles movements. They proposed a matrix C_ij in the 2D case for gradient correction. In the present investigation, a 3D correction matrix for the pressure gradient model is derived as

$$C_{ij}=\frac{1}{d^{’}}{\left(\begin{array}{lll}{\displaystyle \sum V_{ij}\frac{w_{ij}{x}_{ij}^2}{r_{ij}^2}}& {\displaystyle \sum V_{ij}\frac{w_{ij}{x}_{ij}{y}_{ij}}{r_{ij}^2}}& {\displaystyle \sum V_{ij}\frac{w_{ij}{x}_{ij}{z}_{ij}}{r_{ij}^2}}\\ {\displaystyle \sum V_{ij}\frac{w_{ij}{x}_{ij}{y}_{ij}}{r_{ij}^2}}& {\displaystyle \sum V_{ij}\frac{w_{ij}{y}_{ij}^2}{r_{ij}^2}}& {\displaystyle \sum V_{ij}\frac{w_{ij}{y}_{ij}{z}_{ij}}{r_{ij}^2}}\\ {\displaystyle \sum V_{ij}\frac{w_{ij}{x}_{ij}{z}_{ij}}{r_{ij}^2}}& {\displaystyle \sum V_{ij}\frac{w_{ij}{y}_{ij}{z}_{ij}}{r_{ij}^2}}& {\displaystyle \sum V_{ij}\frac{w_{ij}{z}_{ij}^2}{r_{ij}^2}}\end{array}\right)}^{-1}$$

(15)

where $V_{ij}=\frac{1}{{\displaystyle {\sum }_{i\ne j}{w}_{ij}}}.$

3.4. B-Spline Kernel Function

As discussed by Ashtiani and Farhadi [6], the original kernel function and its derivatives approach an unlimited value when the two particles tend to be clustered, which may produce spurious pressure fluctuations, since an unphysical pressure gradient and velocities may be produced in that case. They suggest an M₄ (cubic) B-spline kernel function, as listed in Equation (16) (the MPS using Equation (16) as kernel function is referred to as MPS-BK (B-spline Kernel function) while using Equation (3) is MPS-OK (Old Kernel function).

$$w(r,r_e)=\left\{\begin{array}{l}\frac{8}{\pi r_e{}^3}\left(1-6{\left(\frac{r}{r_e}\right)}^2+6{\left(\frac{r}{r_e}\right)}^3\right) 0\le r<0.5r_e\\ \frac{2}{\pi r_e{}^3}{\left(2-2\frac{r}{r_e}\right)}^3 0.5r_e<r<r_e\\ 0 r>r_e\end{array}\right.$$

(16)

3.5. Anti-Clustering Model for the MPS Method

Different from the 2D case, the particle movements in the 3D case are more complex. Due to particle disorder and clustering, the unphysical pressure oscillations still occur in the 3D MPS simulation even with the most elaborated stability enhancement method addressed previously. The reason for the particle clustering lies in that there lacks tuned forces to prevent the particles approaching violent 3D flows according to their local states. To overcome such a disadvantage, an anti-clustering model referred to the SPH artificial viscosity concept is introduced in the present investigation.

In fact, artificial viscosity models are commonly used in the SPH method for discontinuity capturing, as discussed by Daniel et al. [31,32,33], as well as used in wave breaking by Edmond [36]. Specifically, the proposed anti-clustering model for the MPS method is adopted as follows:

$$\begin{array}{l}{\left(\frac{d{v}_i}{dt}\right)}_{diss}={\displaystyle \sum _j{m}_j}\frac{\alpha {\nu }_{sig}(v_i-v_j)\cdot {\widehat{r}}_{ij}}{{\overline{\rho }}_{ij}}{\widehat{r}}_{ij}\overline{F_{ij}}\\ ={\displaystyle \sum _j\frac{\alpha {\nu }_{sig}(u_{ij}*x_{ij}+v_{ij}*y_{ij}+w_{ij}*z_{ij})}{n_0{\left|r_i-r_j\right|}^2}\cdot }\frac{\partial w_{ij}}{\partial r_{ij}}\left(x_{ij}i+y_{ij}j+z_{ij}k\right)\end{array}$$

(17)

where ${\widehat{r}}_{ij}\overline{F_{ij}}=\overline{{\nabla }_{ij}{w}_{ij}}$,$\overline{{\nabla }_{ij}{w}_{ij}}=\frac{1}{2}\left[\nabla w_{ij}(r_{ij},r_e)+\nabla w_{ij}(r_{ij},r_e)\right]$, ${\widehat{r}}_{ij}=\raisebox{1ex}{$\left(r_i-r_j\right)$}\!\left/ \!\raisebox{-1ex}{$\left|r_i-r_j\right|$}\right.$, and ν_sig refers to the maximum (averaged) signal speed between a particle pair, which is calculated as

$$\begin{array}{l}{\nu }_{sig}=\left\{\begin{array}{cc}\frac{1}{2}\left[\left(c_{s,a}+c_{s,b}\right)-\beta v_{ab}\cdot {\widehat{r}}_{ab}\right]\cdot l_0^4/r_{ab}{}^4;& v_{ab}\cdot {\widehat{r}}_{ab}<0,and,r_{ab}<0.99l_0\\ 0;& v_{ab}\cdot {\widehat{r}}_{ab}\ge 0\end{array}\right.\\ =\left\{\begin{array}{cc}\frac{1}{2}\left[2c_s-\beta v_{ab}\cdot \raisebox{1ex}{$\left(r_a-r_b\right)$}\!\left/ \!\raisebox{-1ex}{$\left|r_a-r_b\right|$}\right.\right]\cdot l_0^4/r_{ab}{}^4;& v_{ab}\cdot {\widehat{r}}_{ab}<0,and,r_{ab}<0.99l_0\\ 0;& v_{ab}\cdot {\widehat{r}}_{ab}\ge 0\end{array}\right.\end{array}$$

(18)

The “standard” formulation of artificial viscosity given by [32] conserves the total linear and angular momentum. The artificial viscosity is equivalent to a physical viscosity [32]. Equation (17) provides a standard artificial viscosity term similar to [32] and it conserves density, specific momentum, and energy [33]. It is noted that ν_sig is worked as a switch from the anti-clustering term through the sign $v_{ab}\cdot {\widehat{r}}_{ab}$ ($v_{ab}\cdot {\widehat{r}}_{ab}$ < 0 means the approaching of two particles) and the distance of two particles. Here, 0.99l₀ is chosen as a threshold, which is in accordance with incompressibility of the governing equations (l₀ is the initial spacing of particles), i.e., the anti-clustering model only starts working when the distance of the two approaching particles is less than l₀, violating the incompressibility situation. Besides, the magnitude of anti-clustering forces is adjusted by the terms of the relative velocities and $l_0^4/r_{ab}{}^4$, which means a tuned increase in the interparticle forces with a decrease in particle distance. As results, the clustering of particles will be prevented in a smoother way, and so a local pressure extremum can be avoided.

As indicated by Equation (17), the anti-clustering model is added in the particle velocity updating step under the condition of incompressibility violation (local abnormal state), which will introduce no viscosity effects in the fluid flow in normal states, although the model utilizes the concept of artificial viscosity.

Finally, the MPS case using Equations (17) and (18) as the anti-clustering model is referred to as MPS-AC (Artificial viscosity) in this paper.

4. Results and Discussions

The evolution of a swirling 2D square patch of fluid has often been used as a benchmark test for evaluation of SPH [37,38] and MPS [15] numerical schemes in suppression of tensile instability. In the present investigation, this case was also chosen as the benchmark example for different stability enhancement methods applications. As a comparison, the VOF (volume of fluid) simulation under similar conditions with the available commercial software was also conducted. The square patch (cubic patch for 3D) of the fluid was 1.0 m in length (L = 1.0 m) and the material was water. The initial swirling conditions were set as

$$\left\{\begin{array}{c}u=\omega y,\\ \begin{array}{l}v=-\omega x,\\ w=0.\end{array}\end{array}\right.$$

(19)

where ω = 10.0 s⁻¹ and u, ν, and w are the x, y, and z velocity components.

Table 1. Cases of the selected stability enhancement methods’ combination.

Case	Abbreviation	Equations Used in Modified MPS Methods
		Laplacian Operator in PPE	Source Term of PPE	$\nabla \mathit{p}$	Kernel Function	Anti-Clustering Model
1	MPS-HS-HL-ECS-GC-OK	(13)	(14)	(12)	(3)	×
2	MPS-HS-HL-ECS-GC-BK	(13)	(14)	(12)	(16)	×
3	MPS-IL-GM-BK	(6)	(7)	(8)	(16)	×
4	MPS-HS-HL-ECS-GC-BK-AC	(13)	(14)	(12)	(16)	(17)(18)
5	MPS-HS-HL-ECS-GC-OK-AC	(13)	(14)	(12)	(3)	(17)(18)

MPS-HS-HL-ECS-GC-OK: MPS with a Higher order Source term-Higher order Laplacian of pressure-Error Compensating parts in the Source term of PPE-Gradient Correction-Old Kernal function. MPS-HS-HL-ECS-GC-BK: MPS with a Higher order Source term-Higher order Laplacian of pressure-Error Compensating parts in the Source term of PPE-Gradient Correction-B-spline Kernal function. MPS-IL-GM-BK: MPS with an Initial Laplacian of pressure-Gradient operator with a Minimum pressure-Bspline Kernal function. HS-HL-ECS-GC-BK-AC: MPS with a Higher order Source term-Higher order Laplacian of pressure-Error Compensating parts in the Source term of PPE-Gradient Correction-B-spline Kernal function-Artificial viscosity. HS-HL-ECS-GC-OK-AC: MPS with a Higher order Source term-Higher order Laplacian of pressure-Error Compensating parts in the Source term of PPE-Gradient Correction-Old Kernal function-Artificial viscosity.

lists the five cases of the selected stability enhancement methods’ combination. According to Table 1, the stability enhancement methods adopted by Case 1 are the same as that of Khayyer and Gotoh [15], who proved its effects in the 2D case. For Case 2, only the kernel function is changed to the B-spline type based on Case 1. While for Case 4 and Case 5, the anti-clustering model is added based on Case 2 and Case 1, respectively. As for Case 3, two early proposed stability enhancement methods were applied.

4.1. Simulation of a Swirling 2D Square Patch of Fluid

4.1.1. Comparison of Selected Stability Enhancement Methods in the 2D Case

In this section, the initial distance between the particles is 0.025 m, and 1681 particles were generated for the present cases.

Figure 1 shows deforming states for four cases for t = 0.2 s of computation. By comparison between cases, it is clear that:

(1)

For Case 1 and Case 2, although the stability enhancement methods are adopted the same as that by Khayyer and Gotoh [15], no stable deforming pattern as shown in Reference [9] is observed in the present cases. The reason may lie in the coarse initial particle spacing (0.025 m, while 0.002 m is adopted in Reference [9]) for the present cases.

(2)

For Case 4 and Case 5, the anti-clustering model is adopted on the basis of Case 1 and Case 2, respectively. The results show an impressive improvement of the stable deforming pattern, especially for Case 5. Note that the only difference between Case 4 and Case 5 is the kernel function. For the present cases, adopting the kernel function in Equation (3) shows better effects than adopting the kernel function in Equation (16).

(3)

The results of the four cases preliminarily validate the effects of the anti-clustering model proposed in the present investigation on stability.

Figure 2 shows the time histories of the pressure and particle number density of a particle initially located at the center of the square fluid patch in the 2D MPS simulations, where the unit of tw is 1. As a comparison, the time history of the pressure at the same location of the VOF simulation and BEM simulation result [39] is also plotted in Figure 2a. It is clear that:

(1)

Pressure oscillations are observed for Case 1 and Case 2, especially after t = 0.15 s, while for Case 4 and Case 5, a relatively smooth pressure change in the whole computation process is observed, indicating the anti-clustering model played a key role in the stability of the simulating process.

(2)

The pressure change patterns observed for the four cases are also in accordance with the change patterns in Figure 2b, in which the change in particle number density vary smoothly around zero for Case 4 and Case 5 while rapidly deviate from zero after 0.2 s for Case 1 and Case 2.

(3)

Comparing the pressure histories between that of the MPS cases and VOF, the results calculated by Case 4 and Case 5 are observed in relatively good agreement with that calculated by VOF generally, although small pressure oscillations are observed for Case 4 and Case 5, indicating no negative effects on accuracy for the presently proposed anti-clustering model.

4.1.2. The Effects of Particle Spacing on the Accuracy of the Proposed Stability Method

To check the accuracy of proposed stability method in Case 5, the deformation process of the swirling square fluid patch was computed with three particle spacings: 50 mm, 25 mm, and 10 mm. The results are compared with those of VOF at the same conditions, as shown in Figure 3. The deformation shape calculated by VOF at t = 0.2 s is represented by iso-surface lines of the water volume fraction, with a = 0.001~0.5 (the interface of the VOF is not clear due to numerical diffusion, while the deformation shape of the fluid calculated by the MPS method is represented directly by the particles). According to Figure 3, the evolution of a square fluid patch calculated by the MPS method is in good agreement with the results calculated by VOF generally, especially for the case of d = 10 mm, indicating the good accuracy and reliability of the proposed stability MPS simulation. It is noted that although a relatively coarse particle spacing is adopted for the present simulations (only d = 2 mm is adopted by Khayyer and Gotoh [15]), stable and relatively accurate results are obtained, which should be attributed to the anti-clustering model proposed in the present investigation. Besides, with the reducing of particle spacing in Figure 3, a clearer interface and a smoother pressure field are also observed.

Figure 4 shows the time histories of the pressure and particle number density of a particle initially located at the center of the square fluid patch with a different particle spacing in Case 5 of MPS method; the pressure histories of the VOF at the same conditions are also plotted in Figure 4a for comparison. As shown in Figure 4a, it is clear that the results of all the calculated cases of the MPS are generally in agreement with that of VOF, except that some fluctuations are observed for the MPS cases, which should be attributed to the sensitivity of the particle number density with particle spacing. In fact, as shown in Figure 4b, large magnitude fluctuations of the particle number density are observed for the case with d = 50 mm and d = 25 mm, while most of the smoother particle number density history is observed for the case with d = 10 mm, which is in accordance with smoother pressure history and in better agreement with that of the VOF for case with d = 10 mm in Figure 4a. Therefore, it is worthy to be pointed out that choosing a proper particle spacing is very important to obtain an accurate MPS simulation, which should be balanced between accuracy and cost.

4.1.3. Calculation Verification of Hydrostatic Pressure

In order to verify the correctness of the MPS algorithm in this paper, the two-dimensional hydrostatic pressure problem has been calculated. The calculation domain is a rectangular water tank with a liquid level of 0.2 m. A pressure measuring point is arranged at the center of the bottom of the water tank, and the range of bottom pressure is recorded by the point in the calculation process. The initial particle spacing is 5 mm, the density of water is 1000 kg/m³, and the acceleration of gravity is 9.8 m/s². The calculation results are shown in Figure 5 and Figure 6.

Figure 5 shows the snapshots of the hydrostatic pressure at different times. It can be seen that as the liquid depth is increasing the pressure gradually increases. Figure 6 shows the pressure at the bottom pressure measurement point during this calculation. It is obvious that the curve is consistent with the theoretical value. Because the particles are constantly exercising during the calculation, the pressure still has certain oscillations.

4.2. Comparison of Selected Stability Enhancement Methods in the 3D Case

The computational models for the MPS as well as VOF in 3D are shown in Figure 7. For the MPS case, the initial distance between the particles was 0.05 m, and 9261 particles were generated for present case. For VOF case, the computational domain was a 7 m cubic air region with 1 m cubic water zone at the center. A total of 2,744,000 mesh cells were generated, with a uniform mesh size of 0.05 m.

4.2.1. Simulation of a Swirling 3D Cubic Patch of Fluid

Figure 8 shows deforming states for four cases at the end of computation. For Case 1, the divergence occurs at t = 0.134 s in the simulation; for Case 2, the divergence occurs at t = 0.190 s; for Case 3, the divergence does not occur, but several particles are flying out of deforming region at t = 0.47 s with large velocities, which should be caused by an unphysical pressure extremum—besides, the deforming shape at the free surface of Case 3 is also incorrect (should collapse to the center); and for Case 4, no divergence and unphysical deforming shape occur, indicating a stable and correct 3D computation.

Theoretically, the deforming of a swirling cubic fluid patch should evolve as follows: in the process of swirling, the cubic fluid patch was extended as a whole within a plane perpendicularly to the rotating axis by centrifugal force, while the negative pressure produced by centrifugal movement push the particles back to center. The movements of particles are determined by the combination of both forces addressed above. For particles far from the axis, a larger centrifugal force will be exerted, while for particles closer to axis, a smaller centrifugal force will be exerted. Therefore, for particles on a free surface cubic fluid patch, a deforming shape of collapsing to the center will be observed in the process of outward extending. Based on such a theoretical analysis, it is obvious to observe that only the deforming evolution of Case 1, Case 2, and Case 4 are theoretically correct.

Further comparison between cases reveals that:

(1)

For Case 1 and Case 2, only the kernel function is different, while the results show that the stability of Case 2 is more positive than that of Case 1. So, unlike in the 2D case, the kernel function (16) is more suitable for the 3D case.

(2)

For Case 2 and Case 3, different PPE and pressure gradient models were adopted, while the results show a better stability but worse deforming evolution for Case 3, verifying the exactness of the high-order PPE, as discussed by Khayyer and Gotoh [9].

(3)

For Case 2 and Case 4, the difference lies in the anti-clustering model adopted by Case 4, while the results show only Case 4 obtain both a stable simulation and correct deforming evolution of cubic fluid significantly, verifying the positive effect of the anti-clustering model to stability in the 3D case.

In general, as discussed previously, there exists a more complex situation in the 3D case for stability and accuracy. Some stability enhancement methods tested in the 2D case cannot fully meet the requirements in the 3D case. According to the comparison of the four cases, it is concluded preliminarily that adopting a comprehensive combination of techniques, such as a high-order PPE, pressure gradient model, kernel function, as well as the anti-clustering model proposed in present investigation, has both good effects on the stability and accuracy in the 3D case.

Figure 9 shows the time histories of the pressure and particle number density of a particle initially located at the center of the cubic fluid patch in the MPS simulations. As a comparison, the time history of the pressure at the same location of the VOF simulation is also plotted in Figure 9a. It is clear that:

(1)

Pressure peaks are observed for Case 1 and Case 2 at dimensionless time 1.2, which corresponds to the time of divergence. For Case 3 and Case 4, a relatively smooth pressure change in the whole computation process is observed.

(2)

Comparing the pressure histories between that of the MPS cases and VOF, a relatively good agreement between that of Case 4 and VOF is observed generally, although a little over-predicted negative pressure at the beginning and a small pressure oscillation at dimensionless time 1.4~3.0 are also observed. Here, it is worthy to note that for MPS computation, all stability enhancement methods tend to reduce but not get rid of pressure fluctuation completely due to the fact that the calculated pressure values of the particles are in charge of correcting all errors produced in a previous explicit step. Hence, for violent particle movements, with the disorder of the particles’ position, especially in the 3D case, pressure fluctuation is unavoidable. As observed for Case 4 at dimensionless time 1.4~3.0, it corresponds to the time zone of collapse of the upper and down free surfaces. Particles tend to be disordered at this time. Obviously, the anti-clustering model plays a positive effect on the stability for Case 4 as compared with that of Case 2, and therefore computation is going on for Case 4 while divergence occurs for Case 2.

(3)

Further comparing the pressure histories of Case 2 and Case 4 before the divergence time of Case 2, they are the same with each other, verifying the fact that the addition of an anti-clustering model proposed in the present investigation will not affect the accuracy of the MPS method but enhance the stability of the solution.

In order to evaluate the accuracy of the MPS method in the 3D case, the deforming process of a swirling cubic fluid patch calculated by MPS was compared with that of VOF, as shown in Figure 10. The deformation shape calculated by the VOF is represented by the iso-surface of the water volume fraction between 0.001 and 0.5, as done in the 2D case, while the deformation shape of the fluid calculated by MPS is represented directly by the particles. According to Figure 10, it can be observed that the evolution of a cubic fluid patch calculated by MPS is in good agreement with the results calculated by VOF, verifying the accuracy and reliability of the MPS simulation. On the other hand, owing to the limitation of the interface-capturing accuracy and the effect of the numerical diffusion of VOF, no exact deforming interface shape of the VOF can be observed, while the deforming interface shape of the MPS method can be observed accurately.

4.2.2. Simulation of 3D Fluid Extruding

With the combined method of Case 4, we simulated another case of a 3D cubic fluid patch extruding, as shown in Figure 11, in which the half patch of fluid where z > 0 is firstly initialized with a velocity of 0.1 m/s along a negative z axis while the left patch is kept stationary. In the fluid extruding case, the particle is forced to be approaching, which more easily cause particle disorder and clustering, leading to instability. Figure 11 shows the whole evolution of a fluid patch, as indicated by the red particles. The evolution pattern simulated by the VOF was also shown by the iso-surfaces of the water volume fraction between 0.01 and 0.5. It can be observed that the deformation of the cubic patch of fluid calculated by the MPS method is relatively smooth, and the free surface evolution is in good agreement with the results calculated by the VOF. No spurious instability occurs in the simulation process, verifying the accuracy and reliability of the MPS simulation by the proposed methods’ combination.

4.3. Two-Phase Coupling between the MPS Liquid Phase and the Euler-Grid Gas-Phase Compressible Flow

Liquid and gas phases were coupled in two phases by the 3D MPS and the FLUENT V17 software (Euler grid), respectively. Liquid water jetting in the boundary layer of the static gas flow was calculated and the atomization of the liquid jet was investigated. The results were compared to those of the VOF method, which captures the continuous interface, under the same conditions.

Figure 12 demonstrates the two-phase coupling diagram that may be created in the spray calculation. This can be divided into 3 areas: (1) the large interface area, which is mainly located in the primary crushing zone. In this area, the liquid phase occupies several gas-phase calculation grids with a clear gas–liquid interface; (2) small interface area, mainly in the secondary crushing area, in which the liquid phase is not enough to occupy a single gas-phase calculation grid; (3) free particle area, primarily formed by particles separated from the interface as the smallest atomized ones, whose represented actual droplet size is determined by the initial particle distribution. These three situations require to be handled separately, and the coupling methods are proposed as follows:

The large interface area is treated by the Immersed Boundary Method (IBM) [40]. The small interface and the free particle areas use particle dynamics to establish the two-phase coupling conditions (including momentum, mass, and energy).

The initial mesh and boundary conditions are provided in Figure 13. The square jet had a size of 0.5 mm × 0.5 mm, with a particle spacing of 0.05 mm. The liquid jet velocity was 20 m/s, and the initial pressure of the gas phase field was 50,000 Pa.

Figure 14 illustrates a comparison of the liquid jet calculations in the static airflow between the MPS integration and the VOF methods. The jet simulated by the VOF indicates the iso-surfaces of the water volume fraction to be 0.5, which demonstrates a shape of the unbroken continuous interface, as shown in the gray translucent surface (see Figure 14). The calculation using the MPS integration method is indicated by particles in the figure, in which red particles represent free particles (B = 0), while blue particles refer to internal particles (B = −1), or wall virtual particles in particular if at the nozzle (B = −2). Figure 14 demonstrates that the position of the continuous fluid interface calculated by the MPS integration method is substantially coincident with that of the interface by the VOF method. In the early stage of jetting, the VOF captures a clear interface since the liquid interface has not been significantly broken. However, as the jetting continues to penetrate due to the shearing effect between the liquid jet and the surrounding air, some free surface particles gradually fall off from the continuous fluid interface to form atomized droplets. The VOF fails to describe this process in detail. In contrast, the MPS integration method captures the gas–liquid interface and the movement of atomized particles, forming a clear simulation picture of the atomization process. This reflects the characteristics of the MPS integration method. Besides, a central continuous jet column is formed in the late stage of jetting, which is consistent with the results by both the MPS and VOF methods. This further verifies the reliability of the method adopted in this paper.

Figure 15a,b, respectively, show the pressure maps and velocity maps at a liquid jet velocity of 20 m/s in the stationary airflow. The figure indicates that large droplets are created by the surface tension when the liquid fuel is injected into the gas-phase flow field at the time of 1.5 × 10⁻⁵ s and 3.2 × 10⁻⁵ s. With further liquid injection, the liquid jet column presents a dynamic process, such as crushing and peeling. Simultaneously, the velocity field near the nozzle changes locally, and the local vortex caused by the jetting can be clearly identified. Driven by the surrounding airflow, the particles on the free interface are gradually falling off, forming atomized droplet, and moving downstream at a certain angle to create a distinct spray angle.

5. Conclusions

Aiming at developing a 3D, stable MPS method, the selected stability enhancement methods proposed by different researchers was extended to the 3D case, and then an anti-clustering model was proposed to further reduce the instability in the MPS method. Meanwhile, the effect of the anti-clustering model is validated, and the effect of the particle spacing is also examined. With different combinations of the proposed methods and models, the deforming process of a swirling cubic fluid patch was simulated comparatively. The results in 2D and 3D were also compared with that of the VOF model. It is concluded that:

(1)

For the 3D case, only extending and applying previous stability enhancement method is still dissatisfactory, while good effects in accuracy and stability can be obtained by incorporating the anti-clustering model proposed in the present paper.

(2)

In view of accuracy, adopting a high-order Laplacian model and high-order source term in PPE are more accurate than that adopting a low-order operator.

(3)

No matter if in the 2D or 3D case, the anti-clustering model presented in this paper apparently plays a very positive role for further stability enhancement of MPS. The results show that the model can produce a tuned interparticle force to prevent particle clustering and introduce no additional viscosity effects in the flow of the normal state.

(4)

Finer particle spacing is very important to obtain an accurate MPS simulation, although the anti-clustering model helps to stabilize the computation of cases with a coarser particle spacing.

In summary, a more complex situation exists in the 3D case for stability and accuracy. Some stability enhancement methods tested in the 2D case cannot fully meet the requirements in the 3D case. According to the results of the present investigation, it is encouraging to summarize that adopting a comprehensive combination of techniques, such as a high-order PPE, pressure gradient model, kernel function, as well as the anti-clustering model proposed by authors, has prominent effects on stability and accuracy in complex 3D flow situations, which will be helpful to many engineering applications.

Liquid and gas phases were coupled in two phases by the 3D MPS method and the FLUENT V17 software (Euler grid), respectively. Liquid water jetting in the boundary layer of the static gas flow was calculated by this integrated method and compared to the VOF method. The results show that both the MPS integration and the VOF methods capture a clear interface at the beginning. As the jet penetrates due to the shearing effect between the liquid jet and the surrounding air, the VOF cannot describe the process in detail, while the MPS integration method still can capture the gas–liquid interface and the movement of the atomized particles, creating a clear simulation of the atomization process.