Stability and accuracy analysis for viscous flow simulation by the moving particle semi-implicit method

The viscosity term, which is the explicit part of the MPS method, is analyzed here. In this section, the term 'viscosity' and 'diffusion' can replace each other freely because they have similar physical mechanisms. To simplify the analysis, the one-dimensional diffusion or viscosity problem is considered as follows:

$$\begin{equation} \frac{{\partial \phi }}{{\partial t}} = v\frac{{\partial ^{\rm{2}} \phi }}{{\partial x^{\rm{2}} }}\quad \left( { - \infty < x < + \infty ,t > {\rm{0}}} \right), \end{equation} \tag{ 12 }$$

where ϕ is also an arbitrary scalar variable and ν is the diffusion coefficient. Given an initial scalar field $\left. \phi \right|_{t = {\rm{0}}} = \phi _0 \left( x \right)\,\left( { - \infty < x < + \infty } \right)$, ϕ will vary with time according to equation (12). The theoretical solution of equation (12) is obtained from the Fourier transformation as follows:

$$\begin{equation} {\phi \left( {x,t} \right) = \int_{ - \infty }^{ + \infty } G \left( {x,t;\xi ,{\rm{0}}} \right)\phi _0 \left( \xi \right)\,{\rm{d}}\xi }, \end{equation} \tag{ 13 }$$

where the kernel function G(x,t;ξ,0) is a Gaussian function whose symmetric center is ξ. So the Gaussian function is a theoretical weight function in the Laplacian model. The specific expression of the Gaussian function is

$$\begin{equation} {G\left( {x,\Delta t;\xi ,{\rm{0}}} \right) = G\left( {x - \xi ,\Delta t - {\rm{0}}} \right) = \left( {\frac{{\rm{1}}}{{\sqrt {4\pi \,{\rm{d}}v\Delta t} }}} \right)^d \,{\rm{e}}^{\frac{{ - \left( {x - \xi } \right)^2 }}{{4\,{\rm{d}}v\Delta t}}} }, \end{equation} \tag{ 14 }$$

where Δt is the diffusion time. When d = 1, equation (14) is the solution for the one-dimensional problem described by equation (12). The diffusion time, Δt, is a key parameter affecting the shape of the Gaussian function. As depicted in figure 1, with the increase of diffusion time Δt, the Gaussian function is dumpier, which means the effective range of diffusion is longer.

Zoom In Zoom Out Reset image size

The effective range of diffusion is defined quantitatively as below. σ, which is the root-mean-square deviation of the Gaussian function, is defined as follows:

$$\begin{equation} \sigma ^2 (\Delta t) = \int_{ - \infty }^{ + \infty } {x^2 G\left( {x,\Delta t;0,{\rm{0}}} \right){\rm{d}}x = 2\,{\rm{d}}v\Delta t} . \end{equation} \tag{ 15 }$$

Here dx is the distance, area, volume of a point in one dimension, two dimensions, three dimensions, respectively. Take the heat diffusion problem as an example. In the diffusion time Δt, 99.7% of the heat diffused from the heat source at the symmetric center is still in the range of three times root-mean-square deviation (3σ) from the symmetric center as shown in figure 2. So the diffusion or viscosity is a local effect and 3σ is the effective range of diffusion. Moreover, 3σ is the effective range of the Gaussian function and 3σ should be the interaction radius in the Laplacian model.

Zoom In Zoom Out Reset image size

We can rewrite the Gaussian function as follows and deduce a Laplacian model based on the Gaussian function:

$$\begin{equation} {G\left( {\sigma ,r} \right) = \left( {\frac{{\rm{1}}}{{\sqrt {2\pi \sigma ^2 } }}} \right)^d {\rm{e}}^{\frac{{ - r^2 }}{{2\sigma ^2 }}} }. \end{equation} \tag{ 16 }$$

The increase of ϕ_i at particle i can be expressed as follows:

$$\begin{eqnarray} &&\fl \Delta \phi _i = \frac{1}{{n_G^0 }}\sum\limits_j {\left( {\Delta \phi _{j \to i} - \Delta \phi _{i \to j} } \right) = } \frac{1}{{n_G^0 }}\sum\limits_j {\left( {\phi _j G\left( {\sigma ,\left| {{\bf{r}}_j - {\bf{r}}_i } \right|} \right) - \phi _i G\left( {\sigma ,\left| {{\bf{r}}_j - {\bf{r}}_i } \right|} \right)} \right)} \nonumber\\ &&= \frac{1}{{n_G^0 }}\sum\limits_j {\left( {\left( {\phi _j - \phi _i } \right)G\left( {\sigma ,\left| {{\bf{r}}_j - {\bf{r}}_i } \right|} \right)} \right)} . \end{eqnarray} \tag{ 17 }$$

Note that $n_G^0 = \sum\nolimits_{j,i} {G(\sigma ,| {{\bf{r}}_j - {\bf{r}}_i } |)}$, where label i is one of the summation index j. So the Laplacian model can be transformed to

$$\begin{equation} \left\langle {\nabla ^2 \phi } \right\rangle _i = \frac{{\Delta \phi _i }}{{v\Delta t}} = \frac{1}{{v\Delta t\,n_G^0 }}\sum\limits_{j \ne i} {\left( {\left( {\phi _j - \phi _i } \right)G\left( {\sigma ,\left| {{\bf{r}}_j - {\bf{r}}_i } \right|} \right)} \right)} . \end{equation} \tag{ 18 }$$

The term νΔt above should be substituted by the interaction radius r_e according to the relationship $r_e = 3\sigma = 3\sqrt {2d\,v\,\Delta t}$ as discussed above. So the Laplacian model based on the Gaussian function is

$$\begin{equation} \left\langle {\nabla ^2 \phi } \right\rangle _i = \frac{{18d}}{{r_e ^2 \;n_G^0 }}\sum\limits_{j \ne i} {\left( {\left( {\phi _j - \phi _i } \right)G\left( {\left| {{\bf{r}}_j - {\bf{r}}_i } \right|} \right)} \right)} . \end{equation} \tag{ 19 }$$

The widely used Laplacian model, equation (9), can be obtained by substituting the Gaussian function and (3σ)² (or r_e²) in equation (19) by a hyperbolic function and λ (equation (10)) respectively because both σ² and λ represent the variance of the respective weight function and have the same units. As shown in figure 3, the shape of the Gaussian function is similar to that of the hyperbolic function when their effective radius is the same, which guarantees the accuracy of the conventional Laplacian model with the λ parameter.

Zoom In Zoom Out Reset image size

Here we discuss how the configurations of time step Δt and space step l₀ affect the stability and accuracy of the viscosity term in the MPS method. Firstly, two conditions for stable and accurate simulation of the viscosity term are proposed and then the reasons why the accuracy condition provides the highest accuracy of the viscosity term and why numerical viscosity in the stability condition of the viscosity term is larger than fluid viscosity are explained.

When considering the stability and accuracy issues, one can compare the realistic effective range of viscosity with interaction radius in the Laplacian model. The realistic effective range of fluid viscosity during time step Δt is ${{\rm{3}}\sigma {\rm{ = 3}}\sqrt {{\rm{2}}dv\Delta t} }$ where Δt is the time step and ν is the fluid viscosity. This realistic effective range is shown by the solid line circle B in figure 4. On the other hand, each radius of circles A–C in figure 4 can be selected as the interaction radius of the Laplacian model for the viscosity term. When the realistic effective range (3σ) is selected as the interaction radius (r_e) for the Laplacian model, the viscosity term will be simulated realistically and accurately. So the relationship

$$\begin{equation} {{\rm{3}}\sqrt {{\rm{2}}dv\Delta t} = r_e } \end{equation} \tag{ 20 }$$

is termed the accuracy condition of the viscosity term. In this situation, the accuracy of the viscosity term will be highest, which will be explained in the next paragraph. When circle A in figure 4 is adopted in the Laplacian model for the viscous term, the particles in a larger area than the realistic diffusion area will be selected to interact with a central particle. The simulation tends to be stable when the numerical diffusion area is larger than the realistic diffusion area (Anderson 1995). Also, considering that diffusion can average velocity field, the velocity will be over-averaged in a larger area when circle A is adopted, which helps to stabilize the simulation. So the simulation will be stable when the interaction radius for the viscosity term is larger than the realistic effective range. The relationship

$$\begin{equation} {r_e \geqslant {\rm{3}}\sqrt {{\rm{2}}dv\Delta t} },\quad {\rm{or}}\quad {\frac{{\nu \Delta t}}{{l_0 ^2 }} \leqslant \frac{{k^2 }}{{18d}}\quad \left( {r_e = k\,l_0 } \right)} \end{equation} \tag{ 21 }$$

is termed as the stability condition of the viscosity term. Comparing condition (21) with conditions (1)–(3), one can find that the parameter at the right-hand side of condition (21) relates closely with the interaction radius of the Laplacian model. This parameter increases with the interaction radius, r_e, which is also reasonable. This stability condition is especially for the particle method, so there is no need to borrow the stability condition from the grid-based method.

Zoom In Zoom Out Reset image size

Now let us turn to the situation where r_e is adopted as the interaction radius of the viscosity term in the MPS method. One can calculate numerical diffusion time Δt_d from the interaction radius (r_e) of the Laplacian model by the condition ${r_e {\rm{ = 3}}\sqrt {{\rm{2}}dv\Delta t_{\rm{d}} } }$ where ν is regarded as fluid viscosity. From the comparison between the calculated numerical diffusion time Δt_d and the time step Δt, the reason why the accuracy condition of the viscosity term produces the highest accuracy of the viscosity term can be explained. From equation (18), one can conclude that the Laplacian model in MPS is the average change rate of a variable during the diffusion time (Δt_d) caused by the viscosity term when considering ν as a constant coefficient to adjust units in Newtonian fluid simulations. On the other hand, when calculating the velocity variation Δu caused by the viscosity term during a time step Δt based on the following equation,

$$\begin{equation} \Delta {\bf{u}} = {\bar{\bf a}} \cdot \Delta t, \end{equation} \tag{ 22 }$$

one just wants ${\bar{\bf a}}$, the average change rate of velocity during the time step Δt. When the average change rate during a longer or shorter time interval than time step Δt is adopted as ${\bar{\bf a}}$ in calculating velocity variation Δu during time step Δt according to equation (22), the accuracy of Δu will decrease so the average change rate during the time step Δt can exactly provide the highest accuracy of velocity variation. When the accuracy condition (20) is adopted to set up Δt and l₀ (further setting up r_e because r_e = kl₀), the calculated numerical diffusion time Δt_d is exactly the time step Δt, which means the calculated Laplacian model (namely, the average change rate of velocity during numerical diffusion time Δt_d) is exactly the wanted average change rate during the time step Δt. So in this situation, the viscosity term is simulated most accurately and realistically. Otherwise, irrespective of whether circles A or B in figure 4 are adopted, the calculated diffusion time Δt_d is either longer (circle A) or shorter (circle C) than the time step Δt causing that the calculated Laplacian model is the average change rate during either longer (circle A) or shorter (circle C) diffusion time Δt_d, so that the accuracy of the viscosity term will decrease when this calculated Laplacian model is adopted. Therefore, the accuracy condition of the viscosity term produced the highest accuracy for the viscosity term just because the Laplacian model calculated from r_e in this situation is exactly the wanted average change rate during time step Δt. This is the reason why condition (20) is addressed as the accuracy condition of the viscosity term.

When stability condition (21) or circle A is adopted, the calculated Laplacian model is over-averaged during a longer time interval (Δt_d) than the time step (Δt) so that the simulation tends to be stable. On the other hand, since the circle B in figure 4 provides the highest accuracy of viscosity, the accuracy will decrease a little when circle C is adopted just a little smaller than circle B to setting up Δt and l₀. In this situation, the simulation may be stable because the slightly decreased accuracy is not enough to break down the simulation. However, when circle C in figure 4 is much smaller than circle B, the simulation tends to be unstable because of the insufficient diffusion. Therefore, the stability condition of the viscosity term is just a sufficient rather than necessary condition of stable calculation of the viscosity term.

In an alternative perspective, one can calculate a numerical viscosity coefficient ν_n from the interaction radius (r_e) of the Laplacian model by the condition ${r_e {\rm{ = 3}}\sqrt {{\rm{2}}dv_{\rm{n}} \Delta t} }$ where Δt is regarded as time step. When the accuracy condition of the viscosity term (${{\rm{3}}\sqrt {{\rm{2}}dv\Delta t} = r_e }$) is adopted to set up Δt and l₀ (further setting up r_e because r_e = kl₀), the numerical viscosity ν_n is exactly the fluid viscosity ν so the viscosity term is realistically and completely modeled without extra numerical viscosity. While the stability condition of the viscosity term (circle C) is adopted for setting up Δt and l₀, the calculated numerical viscosity coefficient ν_n will be larger than the fluid viscosity so extra numerical viscosity is added, which also indicates why the stability condition of the viscosity term provides stable simulation. The over-averaged velocity field in which the stability condition of the viscosity term is adopted indicates extra numerical viscosity. When circle C in figure 4 is adopted in the Laplacian model for the viscous term, the numerical viscosity is even smaller than the fluid viscosity so the simulation may be unstable. The two proposed conditions mainly apply to the flow where viscous forces are dominant or comparable with other forces.

The CFL to stabilize the simulation is

$$\begin{equation} {u_{{\rm{max}}} {\rm{\Delta }}t = C_{\rm{0}} l_{\rm{0}} },\quad C_0 \leqslant C_{{\rm{max}}} , \end{equation} \tag{ 23 }$$

where u_max is the maximal velocity, C₀ is the Courant number and C_max (typically 0.25) is a maximal Courant number in particle methods. Both CFL and the proposed accuracy or stability condition of the viscosity term are necessary to set up time step Δt and space step l₀. The CFL condition which is a straight line and the accuracy condition of the viscosity term which is a parabolic curve are plotted in figure 5. The areas A and B are stable areas guaranteed by CFL. The areas B and C are stable areas guaranteed by the stability condition of the viscosity term. So area B is the common and final stable area for viscous flow. The accuracy and stability conditions of the viscosity term are very similar to each other, but they differ as shown in figure 5. The accuracy condition of the viscosity term states that the combinations of l₀ and Δt on the curve ${ {\mathord{\buildrel{\lower3pt\hbox{$\frown$}} \over {\rm ode} } } }$ in figure 5 produce the most accurate simulation of the viscosity term as discussed above. When the accuracy condition of the viscosity term is adopted to realistically and accurately simulate the viscosity term causing that the point must be selected on the curve ${ {\mathord{\buildrel{\lower3pt\hbox{$\frown$}} \over {\rm ode} } } }$, point (e) in figure 5 can be selected to save computational load because Δt and l₀ at point (e) are the largest ones on the curve ${ {\mathord{\buildrel{\lower3pt\hbox{$\frown$}} \over {\rm ode} } } }$ in the stable area. The computational load is still very large even though point (e) is selected, which will be discussed in the next subsection. So the accuracy condition of the viscosity term is a very stringent condition because the viscosity term is realistically and accurately simulated. On the other hand, the stability condition of the viscosity term states that the combinations of l₀ and Δt in the area B in figure 5 produce stable simulation for viscous flow. Area B in figure 5 is a very large area so that much larger Δt and l₀ than Δt and l₀ at point (e) can be selected to produce stable simulation. So the stability condition of the viscosity term can save computational load obviously but extra numerical simulation is introduced as discussed above. For the high Reynolds number flows, the combinations of l₀ and Δt in the area B are usually selected to save computational load.

Zoom In Zoom Out Reset image size

Now, the computational cost when the accuracy condition of the viscosity term is adopted is discussed. To ensure stability of particles' motion and the pressure term, CFL is still necessary. Assuming u_max = βu₀, r_e = kl₀ and Reynolds number Re = Du₀/ν (D is the characteristic length, u₀ is the characteristic velocity), l₀ and Δt can be set up by coupling the CFL (23) and the accuracy condition of the viscosity term (20). The result is rearranged as follows:

$$\begin{equation} {\left\{ {\begin{array}{@{}c@{}} \frac{D}{{l_{\rm{0}} }} = \frac{{k^{\rm{2}} \beta }}{{{\rm{18}}C_{\rm{0}} d}}Re, \\ \frac{D}{{u_{\rm{0}} {\rm{\Delta }}t}} = \frac{{k^{\rm{2}} \beta ^{\rm{2}} }}{{{\rm{18}}C_{\rm{0}}^{\rm{2}} d}}Re. \\ \end{array}} \right\}} \end{equation} \tag{ 24 }$$

The above results are stated as follows:

• (1)  
  D/l₀, which represents the number of particles used to discretize characteristic length D, is proportional to Reynolds number, Re.
• (2)  
  D/(u₀Δt), which represents the number of iterations from the initial state to the fully developed state, is also proportional to Reynolds number, Re.

The computational cost shown above is very large because the molecular viscosity is realistically simulated. Assuming that u_max is proportional to fluctuating velocity u' and that D is the energetic scale, Reynolds number is redefined as Re_l = u'D/v. One can also conclude that D/l₀ is proportional to Re_l in the particle method. So the computational cost of the accuracy condition of the viscosity term in the MPS method is even larger than that of the theoretical direct numerical simulation (DNS) method, which states that the grid number in one dimension is proportional to (Re_l)^3/4 (Wilcox et al 1998).

As discussed above, the accuracy condition of the viscosity term actually provides a realistic and complete simulation of the molecular viscosity term without extra numerical viscosity. So when the accuracy condition of the viscosity term is used to set up l₀ and Δt for the steady laminar flow simulation, the results will be the most accurate because time step and space step match each other. When the accuracy condition of the viscosity term is adopted for the unsteady turbulent flow simulation, no turbulence model is needed. These simulations can be regarded as DNS in the MPS method. But the authors fail to provide an unsteady turbulent flow simulation because of the unacceptable computational load.

The Poiseuille flow is simulated by the MPS method to evaluate the accuracy and stability condition of the viscosity term because the viscous force is dominant in Poiseuille flow. The initial particle configuration is shown is figure 6. The top and bottom wall is set as the wall boundary condition (Koshizuka and Oka 1996). The inflow boundary condition is similar to Shakibaeinia's setting (Shakibaeinia and Jin 2010) and the outflow boundary is set as a free surface boundary (Koshizuka and Oka 1996). In the simulation, gravity is neglected. The fluid property and some parameters are shown in table 1, where D is the distance between the top and bottom walls and L is the length from inlet to outlet as shown in figure 6. Reynolds number Re is small, so the velocity at 6D's cross-section is fully developed. We compare the simulated velocity profile with the theoretical one from 6.5D to 7.5D in the fully developed area as shown in figure 6. The criterion of accuracy is defined by the mean relative error as follows:

$$\begin{equation} {\rm{Err}} = \frac{1}{n}\sum\limits_{i = 1}^n {\frac{{u - u_{{\rm{theory}}} }}{{u_{{\rm{theory}}} }}} , \end{equation} \tag{ 25 }$$

where i is the label of particles which located between 6.5D and 7.5D at the fully developed time layers and n is the number of all particles in the specific area above. Two simulation cases (simulation case 1 and simulation case 2) of Poiseuille flow are simulated to evaluate the proposed conditions of the viscosity term. In simulation case 1, the uniform mean velocity profile is configured for the inlet particles which are used to push the fluid static particles between the two solid walls until the flow is fully developed. There is always about 5% relative error which cannot be eliminated effectively even though we decrease the initial particle distance (l₀). Then all the initial fluid particles are set up with the theoretical fully developed parabolic velocity to test whether the parabolic profile can maintain or not, which is simulation case 2.

Zoom In Zoom Out Reset image size

The time step Δt is fixed in the simulation so that the effect of the proposed stability and accuracy condition of the viscosity term on stability and accuracy is more obvious. Different setups of initial particle distance l₀ and time step Δt are employed to evaluate the stability and accuracy condition of the viscosity term as shown in figure 7. As discussed above in section 3.2, when C₀ ⩽ C_max (C_max = 0.25) in CFL (equation (23)), the CFL is guaranteed. As shown in figure 7, here we choose CFL with smaller C₀ = 0.068 to couple with the accuracy condition to solve the basic Δt and l₀ (point (c)) so that when Δt is enlarged (points (d) and (e)) or l₀ is reduced (points (f) and (g)), the CFL (C₀ ⩽ C_max) is always guaranteed. In figure 7, only point (c) is on the curve of the accuracy condition of the viscosity term. So the accuracy of the simulation with l₀ and Δt at point (c) is expected to be the highest. As shown by the vertical line in figure 7, keeping the space step l₀ still and changing the time step can produce four cases: much shorter Δt (point (a)), shorter Δt (point (b)), longer Δt (point (d)), much longer Δt (point (e)), which are expressed in table 2. As shown by the horizontal line in figure 7, keeping the time step Δt still and changing the space step can produce another four cases: much shorter l₀ (point (f)), shorter l₀ (point (g)), longer l₀ (point (h)), much longer l₀ (point (i)), which are expressed in table 3. Also we compare the accuracy of the Laplacian model based on the Gaussian function (equation (19)) and the conventional Laplacian model (equation (9)). We combine the two Laplacian models and the two simulation cases to produce four schemes: G1, G2, λ1 and λ2, where G represents the Laplacian model based on the Gaussian function (equation (19)), λ represents the conventional Laplacian model with the λ parameter (equation (9)), 1 is short for simulation case 1 and 2 is short for simulation case 2. So G1 means the viscosity term is modeled by the Laplacian model based on the Gaussian function and that the type of the simulation is simulation case 1. These symbols are also explained in detail as the footnote of tables 2 and 3. In G1 and G2, the Laplacian model based on the Gaussian function is only used to calculate the viscosity term. The simulated results are shown in tables 2 and 3, where points (a)–(i) in figure 7 are also shown.

Zoom In Zoom Out Reset image size

Firstly, the relative error of the simulation with l₀ and Δt at point (c) is always the smallest as shown by the bold numbers in tables 2 and 3, because point (c) is on the curve of the accuracy condition of the viscosity term. So the accuracy condition of the viscosity term provides the highest accuracy in the MPS method.

Secondly, comparing figures 5 and 7, one can find that points (a)–(c), (h), (i) in figure 7 are in area B in figure 5, which is the stable area in viscous flow. The simulations with Δt and l₀ at points (a)–(c), (h), (i) are all stable which prove that the stability condition of the viscosity term proposed in this paper is effective to provide stable simulation. Points (d)–(g) in figure 7 are in area A in figure 5, where the CFL (C₀ ⩽ C_max) is guaranteed and the stability condition of viscosity is violated. So the simulations at points (d)–(g) in figure 7 are mainly divergent, which is consistent with most results in tables 2 and 3. However, the simulation at point (d) in the 'λ1' and 'λ2' columns in table 2 and simulations at point (g) in the 'λ' and 'λ2' columns in table 3 are stable, which support that the stability condition of the viscous term proposed in this paper is just a sufficient condition rather than a necessary condition.

Thirdly, one can find that in tables 2 and 3, the relative error of simulations at point (c) in the 'G1' and 'G2' columns is smaller than that in the 'λ1' and 'λ2' columns. So the Laplacian model based on the Gaussian function is more accurate than the conventional Laplacian model with the λ parameter. Yet, the former is more unstable than the latter by comparing the number of divergent cases in the 'G1' and 'G2' columns with those in the 'λ1' and 'λ2' columns.