Improvement of Calculation Stability for Slow Fluid Flow Analysis Using Particle Method
2017 Volume 58 Issue 3 Pages 479-484
Details
2017 Volume 58 Issue 3 Pages 479-484
MPS (Moving particle semi-implicit) method is one of the popular particle methods. It is well known that calculation instability occurs especially for a slow fluid flow analysis, and the instability appears as unnatural pressure and velocity oscillation. This is because the MPS method adopts predictor-corrector method, and tentative position of particles are not corrected enough at the corrector phase. Therefore an improved multiple relaxation method is proposed in this study to improve the calculation stability. Pressure and velocity correction is conducted multiply in the correction phase, and also adjusting coefficient is adopted in the source term in the Poisson equation for the pressure. As a result, the proposed method improved the pressure and velocity distribution instability. The method was applied to complex shaped casting with a shoulder. The conventional method showed unnatural high pressure and the fluid hopped, whereas smooth and natural pressure field could be obtained using the proposed method.
This Paper was Originally Published in Japanese in J. JFS. 87 (2015) 702–708.
Recently, particle methods, which are based on a fully Lagrangian method, have been developing to calculate complex phenomena occurring in the casting processes1–4). Nowadays, the typical particle methods for computational continuum dynamics are usually categorized into two methods, SPH (Smoothed Particle Hydrodynamics)2) method and MPS (Moving Particle Semi-implicit)3,4) method. Discrete objects are used as calculation elements in the particle methods, and they can be placed and moved freely in the space. This feature allows the particle methods to simulate the heat and mass transfer phenomena that are observed in the casting process more easily and directly than other methods that use the calculation lattice.
However, the particle methods have still some problems. Calculation stability is one of the significant problems in the particle methods. Casting process involves a variety of phenomena which occur simultaneously from the pouring to the solidification. Because it is difficult to expect what kind of phenomena will occur at where and at which timing, higher calculation stability is required for the integrated simulation using particle methods. It is known that the oscillation of velocity and pressure occurs in the flow simulation, and the oscillation decreases the calculation stability especially for a slow flow5). Hirata et al. have proposed a stable and rapid calculation method for the slow flow by ignoring inertia force and adjusting the gravitational force after the fluid flow can be assumed to be almost stopped5). However, the procedure requires precise estimation of the flow stop time. The original flow calculation by the particle method3) is unstable for the slow fluid flow. Therefore the flow calculation without any improvement will decrease the calculation stability when the flow becomes slow during the casting processes. Recently, Hirata et al. reported the multiple relaxation method which improves calculation stability and speed6) in the case of slow fluid flow. However, the oscillation of pressure and velocity are still not negligible. Therefore further improvement in the accuracy is expected. In this study, we tried further improvement of the calculation stability and accuracy of the particle method in the case of slow fluid flow by improving the multiple relaxation method.
We adopted MPS method to combine flow and solidification simulations. In the MPS method, the calculations such as gradient operation are evaluated by calculating an interaction among the surrounding particles. The interaction is calculated with the particles by weighted averaging within a certain distance, which is referred to as the kernel size. A summary of the proposed calculation algorithms is shown in this section.
2.1 Weight function and particle interaction modelsA weight function w is used in the MPS method to calculate the interaction of particles by weighted averaging. Equation (1) was used in this study because of its higher stability for the flow calculation7).
| \[\begin{split} & w(r,{c_k}{r_0}) \\ & = \left\{ \begin{array}{cc} \displaystyle\frac{40}{{7 \pi c_{k}}^{2}} \left( 1 - \frac{6r^{2}}{{c_{k}}^{2} {r_{0}}^{2}} + \frac{6r^{3}}{{c_{k}}^{3} {r_{0}}^{3}} \right) & (0 \le r < 0.5c_{k}r_{0}) \\[2mm] \displaystyle\frac{10}{{7\pi c_{k}}^{2}} \left( 2 - \frac{2r}{c_{k} r_{0}} \right)^3 & (0.5c_{k} r_{0} \le r < c_{k} r_{0}) \\[2mm] 0 & (c_{k} r_{0} < r) \end{array} \right. \end{split} \] | (1) |
Here, r is a distance between particles i and j. ckr0 is the kernel size which represents interaction calculation range. The MPS method uses particle size r0 which has the same significance as the lattice size in the Eulerian methods such as FDM (finite difference method). ck is the kernel size coefficient and usually varies between 2 and 4. ck = 2.1 was used in this study as a commonly used value3,7).
Particle number density ni of the particle i is used in the particle method to calculate the interaction between particle i and the surrounding particles. ni is the sum of the weight function of the particles surrounding the particle i, and is expressed as follows.
| \[ n_{i} = \sum_{i \ne j} w \left( \left| \boldsymbol{r}_{j} - \boldsymbol{r}_{i} \right|, c_{k} r_{0} \right) \] | (2) |
Here, ri and rj are the position vectors of particle i and j, respectively. Assuming that the fluid is incompressible and considering that the particle number density is directly related to the fluid density, ni is equal to n0 in the case that a particle has no surface particle around it. And we can use n0 for the mass conservation condition in the incompressible fluid flow analysis using the MPS method.
Particle interaction models are used to describe the differential operators in the MPS method. If ϕi and ϕj are arbitrary scalars at positions ri and rj, then the particle interaction models for the differential operators at the particle i can be expressed as follows.
| \[ \nabla \phi_{i} = \frac{d}{n_{i}} \sum_{i \ne j} \frac{(\phi_{i} - \phi_{j})}{\left| \boldsymbol{r}_{j} - \boldsymbol{r}_{i} \right|} \cdot \frac{(\boldsymbol{r}_{j} - \boldsymbol{r}_{i})}{\left| \boldsymbol{r}_{j} - \boldsymbol{r}_{i} \right|} w \left( \left| \boldsymbol{r}_{j} - \boldsymbol{r}_{i} \right|, c_{k} r_{0} \right) \] | (3) |
| \[ \nabla^{2} \phi_{i} = \frac{2d}{n_{i}} \sum_{i \ne j} \frac{(\phi_{i} - \phi_{j})}{\left| \boldsymbol{r}_{j} - \boldsymbol{r}_{i} \right|^2} w \left( \left| \boldsymbol{r}_{j} - \boldsymbol{r}_{i} \right|, c_{k} r_{0} \right) \] | (4) |
Equation (3) is a gradient model, and eq. (4) is a Laplacian model. Suffixes i and j represent the assigned numbers of particles. The d is the number of space dimensions.
2.2 Flow simulationThe governing equations for incompressible flow calculation by the MPS method consist of the continuity equation and the Navier-Stokes equation as follows.
| \[ \frac{D\rho}{Dt} = 0, \quad \frac{D\boldsymbol{u}}{Dt} = - \frac{1}{\rho} \nabla p + \nu \nabla^{2} \boldsymbol{u} + \boldsymbol{f} \] | (5) |
Here, ρ indicates density (kg/m3); u, velocity vector (m/s); p, pressure (Pa); ν, kinematic viscosity (m2/s); and f (m/s2), the body force acceleration vector including gravity. The left-hand side of the second equation, D/Dt, is the Lagrangian derivative, having the same significance as temporal differentiation.
Flow calculation by the MPS method is based on the predictor-corrector method similar to that in the case of the FDM or other Eulerian methods. The major differences between the MPS method and the FDM are the formulation of the Poisson equation of pressure and the calculation of particle position. In the MPS method, the tentative particle position is calculated using tentative velocity in the prediction phase. Therefore, the correction in position is calculated along with the correction in velocity in the correction phase. The Poisson equation for pressure in the MPS method is described using the tentative particle number density n* and n0 as follows1).
| \[ \nabla^{2} p^{k + 1} = - \frac{\rho}{\Delta t^{2}} \frac{n^{*} - n_{0}}{n_{0}} \] | (6) |
Here, n* is the tentative particle number density calculated using the tentative particle positions r* and eq. (2).
In the MPS method, the pressure is calculated so as to make the tentative particle number density n* to the constant value n0, whereas making the divergence of velocity to zero in the FDM.
In the flow simulation with a free surface, the pressure at the free surface is assumed to be zero. The criteria to determine the free surface particles is described as following equation because the particle number density decrease at the free surface.
| \[ \frac{n_{i}}{n_{0}} < \beta \] | (7) |
β = 0.97 is usually used in the MPS method with a different weight function to the eq. (1)3). In this study, β = 0.985 was used for eq. (1), which shows the most natural free surface particle criteria.
The time increment Δt is calculated by using Courant number Cn using the following equation.
| \[ \Delta t = C_{n} \frac{r_{0}}{\left| \boldsymbol{u}_{\max} \right|} \] | (8) |
Here, umax is the maximum velocity of the particles. It is known that the use of small Δt sometimes worsens the calculation stability. This is because the source term of the Poisson equation includes Δt2 in its denominator, and too small Δt produces unnatural high pressure. In this study, constant value Cn = 0.1 was used.
The authors have reported that the multiple relaxation method improves calculation stability and speed in the previous report6). However unnatural pressure oscillation and unnatural movement of the surface particle are still at issues. Pressure oscillation should be avoided as far as possible to predict a feeding behavior to the solidification shrinkage precisely. Therefore an improved multiple relaxation method was proposed in this study to avoid unnatural pressure oscillation.
3.1 Improved multiple relaxation methodThe Poisson equation of the pressure in the MPS method uses the source term based on the particle number density ni, and the tentative value n* is expected to be corrected to the constant value n0. Therefore, if an excess concentration of particles occurs in the prediction phase, excess positional correction occurs because of excess high pressure. Such excess correction is one of the reasons of velocity and pressure oscillation occurred in the particle method calculation. The multiple relaxation method improves the stability and accuracy of flow simulation using the particle method by calculating the Poisson equation more than once, so as to the particle number density ni become closer to the constant value n06). However, a possibility of excess correction of the position is still remain even by using multiple relaxations; sometimes the excess correction will let particles move a distance more than expected, and the oscillation of velocity and position will occur during the multiple relaxation procedures. Therefore, an artificial adjustment coefficient of the source term in the Poisson equation cpoi was introduced to reduce the amount of positional correction in one relaxation process to settle the particles at the expected position without excess movement. The flowchart for the flow simulation with the improved multiple relaxation method is shown in the Fig. 1. The Poisson equation is solved more than once in the multiple relaxation method. Accordingly the velocity and position are also corrected more than once during one time step. The improved method uses the adjustment coefficient cpoi in the source term in the Poisson equation for the pressure.
Flow chart for the flow simulation.
Two-dimensional flow calculation is conducted for still fluid as shown in the Fig. 2. The width of casting (fluid) is 40 mm and the height is 80 mm. Particle size is 2 mm; thus, we use 800 particles for the casting, 384 particles for the mold. The calculation is conducted for 2 seconds. The maximum time increment Δtmax, the number of multiple relaxations (MR) and the adjustment coefficient in the source term of the Poisson equation cpoi, and the pressure value at three points were measured. The measuring points were at the center in the horizontal directions and the depth from the casting surface were h = 30, 50 70 mm. The physical properties and the calculation conditions are shown in Table 1 and 2. The relationship between Δtmax and the number of multiple relaxation is determined to optimize the stability and the efficiency of calculation referring the past research6). Four cpoi values (cpoi = 1.0, 0.5, 0.2, 0.1) were used for each MR and Δtmax conditions.
Calculation model.
| Casting | |
|---|---|
| Materials | Pure Al |
| Density, kg/m3 | 2500 |
| Kinematic Viscosity, m2/s | 10−5 |
Figure 3 shows a pressure distribution at 0.1 s. From now on, the calculation conditions are expressed using the combination of the number of multiple relaxation and cpoi; for example, in the case of two times of multiple relaxation and cpoi = 0.5, the condition is expressed as MR2-0.5. An ideal hydrostatic pressure at the bottom of the casting is 1960 Pa. Figure 3(a) is the result in the case of the basic MPS method without multiple relaxations and cpoi = 1.0. The pressure distribution was non-uniform, and the maximum value exceeded 5000 Pa. Such non-uniform distribution was not settled during 2 s calculation, and the particles continued to flow in whole unnaturally. Figure 3(b) is the result in the case of two times multiple relaxations with cpoi = 1.0. The pressure at the bottom was around 2000 Pa which is similar to the ideal value, however remains non-uniform distribution. Figure 3(c) is the results in the case of five times multiple relaxations and cpoi = 0.5. In this case, significant improvement was found compared with the former two results. We could obtain almost similar results under the conditions shown in the Table 2, except for MR1-1.0 and MR2-1.0 as shown in the Fig. 3(a) and (b).
Pressure distribution at 0.1 s.
| Δtmax (s) | MR (Multiple relaxation) | cpoi |
|---|---|---|
| 0.0001 | 1 | 1.0 |
| 0.0005 | 2 | 0.5 |
| 0.001 | 5 | 0.2 |
| 0.002 | 10 | 0.1 |
Next, an average pressure value during 2 s at h = 70 mm and a standard deviation of the pressure are shown in the Fig. 4. Black dots are the average value and the error bars correspond to the standard deviation. The ideal pressure value at h = 70 mm is 1715 Pa. The average values indicated by the black dots agreed well with the ideal value in each condition. On the other hand, the standard deviation of MR1-1.0 and MR2-1.0 were larger than the other conditions, showing a similar tendency in the Fig. 3. As a result, the ideal pressure distribution could be achieved in the case that the standard deviation was less than 200 Pa in this study.
Pressure at h = 0.07 m.
For further discussion, the pressure change in time under typical six conditions are shown in Fig. 5. The figures show the pressure oscillation at h = 30, 50, 70 mm. Figure 5(a) to (c) show the results in the case of cpoi = 1.0. The results show that the multiple relaxation reduced the pressure oscillation, and the effect was more obvious if the number of multiple relaxations is higher. Figure 5(d) to (f) are the case of cpoi = 0.5. Figure 5(d) shows that the pressure oscillation could not be reduced enough without multiple relaxations. From the Fig. 5(e) and (f), the combination of the multiple relaxations and adjustment coefficient cpoi significantly reduces the pressure oscillation.
Pressure oscillation.
A relationship between calculation conditions and the fluid behavior was summarized in Fig. 6. The circle indicates the conditions showing the excellent stability of still fluid with less pressure oscillation. Cross mark means poor stability conditions which show an unnatural rapid increase of pressure, and the calculation resulted in failure with a divergence of the velocity field. Diamond shows the conditions which do not result in failure, however, unnatural behavior such as hopping of fluid occurs. Overall, a use of multiple relaxations and less cpoi conditions improved calculation stability. Whereas in the case of 2 times multiple relaxations, too small cpoi could reduce the stability. This is because 2 times multiple relaxations is not enough for pressure, velocity and positional correction in the correction phase if cpoi is too small.
Effect of multiple relaxation and source term adjustment on the pressure stability of the flow simulation for gently placed fluid. Circle indicates excellent stability, diamond shows good stability but large oscillation of pressure occurs, cross means poor stability.
Figure 7 shows a variation of the maximum velocity of particles. Figure 7(a) and (b) are the results in the case of MR1 and MR5, respectively. The results in the case of MR2 and MR10 were similar to the case of MR5. In the case without multiple relaxations shown in the Fig. 7(a), the use of smaller cpoi decrease the maximum velocity, whereas the values were larger than any case with multiple relaxations shown in the Fig. 7(b). However, the maximum velocity was over 0.05 m/s and the average velocity was over 0.02 m/s even using the multiple relaxations. Therefore, if the maximum or average velocity is expected to be under such values, some artificial operation will be useful; such as ignoring the inertia force of fluid particles after pouring as reported by Hirata et al.5)
Maximum velocity of particles.
Finally, we discuss the pressure stability in the complex shaped casting. In this section, three-dimensional flow calculation was conducted for a cylindrical casting with a shoulder. The height of the casting was 180 mm, and the upper half of the casting was 25 mm and the lower half was 40 mm in diameter, respectively. The physical properties used in the calculation are shown in the Table 1, and the fluid was poured from the upper side with a rate of 100 g/s and pouring finish within 4 s. The results are shown in the Fig. 8. Figure 8(a) is the result of MR1-1.0 at 4.5 s, showing unnatural pressure distribution at the bottom of the casting, and hopping of the casting particles. At the last moment, the fluid flow was very slow. However the pressure in the casting was unnaturally high under the shoulder. In such shapes with shoulders, positional correction in the correction phase tends to be insufficient, and the pressure under the shoulder will be accumulated unnaturally. As a result, unnatural hopping occurs, and the calculation will result in failure. In the Fig. 8(b), we can find that the proposed method improves the calculation stability. The calculated pressure value at the bottom agreed well with the theoretical value, 4410 Pa. It is possible to improve calculation stability by using smaller Δt. However the pressure or velocity oscillation problem will remain, and also, longer calculation time will be required because of higher maximum velocity.
Pressure distribution in the casting with shoulder.
The series of results show that the multiple relaxation method and the use of adjustment coefficient cpoi are effective to reduce unnatural oscillation of the pressure and the velocity, and to improve calculation stability for slow fluid flow, even for the complex shaped casting with the shoulder. Improved multiple relaxation method can improve the unnatural velocity decay and calculation stability for faster fluid flow. The results of the faster flow calculation are discussed in the next report8).
The improved multiple relaxation method was proposed to improve calculation stability for slow fluid flow based on the MPS method; one of the particle methods. The multiple relaxation method is combined with the adjustment coefficient of the source term in the Poisson equation of pressure. As a result, the pressure oscillation was significantly reduced by the proposed method, and the calculation stability was also improved. The proposed method can calculate stably even for the complex shaped casting with the shoulder.
The authors wish to thank the young researchers support program 2012 of JFS for its support.
