The numerical method used in this study is the full implicit moving particle simulation (MPS) method, solves the Navier-Stokes equations in a Lagrangian mesh free framework. The fluid domain is represented by particles, and the motion of each particle is calculated based on the interactions with the neighboring particles by means of a kernel function. The viscosity and pressure calculation of Navier Stokes equation are solved by implicit method. Besides that, the calculation of Kernel function and satisfaction of incompressible condition are done similarly as MPS method, which was originally proposed by Koshizuka and Oka (1996) for incompressible flow simulation. In particular, this implicit particle method used in high viscous flow simulation makes the calculation stable at high kinematic viscosity (1~10 m2/sec) and can solve in large time steps (0.0001~0.00001 sec). The MPS method is similar to SPH (Smooth particle hydrodynamics) method developed by Gingold and Monaghan (1977). In last decade particle based method like MPS/SPH method has been applied in wide range of Nuclear Engineering, Ocean Engineering, Automotive and Chemical engineering etc.
The governing equations for incompressible viscous flows are the continuity and Navier–Stokes equations as follows:
The symbol is the density, t the time, the velocity vector, P the pressure, and the kinematic viscosity.
Particle interaction models
The governing equations written with partial differential equation are transformed to the equation of particle interaction. During interaction the discretized particle domain is smoothed by the smoothing function called Kernel function. In this study, the following function is employed:
Fig.1. Kernel Function
Linked-list Search AlgorithmTo carry out search of nearest neighboring particles, we had used the Linked-list search algorithm, for which a temporal square mesh is overlaid on the computational domain and the mesh spacing is defined by the kernel size. Then, for a particle in a cell, its neighboring particles within the kernel size are considered in the same cell; otherwise, they belong to the surrounding cells. This search method is very efficient and robust, especially for the simulation with a large number of particles.
A gradient vector between two particles i and j possessing scalar quantities fi and fj at coordines ri and rj is represented as (fi – fj) (ri-rj)/( ri-rj)2 as shown in Fig. 2. The gradient vector of each particle i is given by the weighted average of these gradient vectors:
The diffusiono at particle i is described as
Fig.3. Diffusion Model
Where, Lambda(?) is the parameter by which the varience increase is equal to that of the analytical solution. The diffusion can be modeled by distribution of a physical value from a particle to its neighboring particles by use of the kernel function. As shown in Fig. 3, the model is conservative since the quantity lost by particle i is obtained by the neighboring particle j.
Incompressibility is represented by keeping the particle number density constant. After calculation of viscosity and gravity it is assumed that the particle number density n* deviates from n0. The particle number density is to be modified to n0 by n`:
Thus, the continuity equation (1) is fulfilled with fixing the particle number density through the simulation. And the modification of the particle number density is related to the modification of velicity v`.
In present study, the kinematic and dynamic boundary conditions are imposed in free surface particles. Free surface particles are judged by particle number density by satisfying condition n*< βno. Where β is chossen as 0.98 in this study. The kinematic condition can be directly satisfied by moving particles on the free-surface and dynamic condition is satisfied by making atmospheric Pressure (Patm = zero) at free surface particles.
Fig.4 Free surface particles
Algorithm for Full Implicit MPS
The algorithm of incompressibility for the MPS method is similar to that of the SMAC method in a grid system. In which viscosity, external forces and convetion terms are explicitly calculated. But while solving viscosity by explicit process seems unstable and we have to choose extremely small time step which makes calculation, time consuming and extremely unstable to calculate. So, in this paper we had choose to compute viscosity by Implicit method.
In each time there are two stages: in the first stage, the temporal velocity components and coordinates of particle i are obtained using diffusion, external forces, and convection terms, which are implicitly calculated with the values in the nth timestep.
Thus, the temporal coordinates ri of particle i can be written using the temporal velocity ui as follows:
The left hand side of Eq(11) is discretized by diffusion model Eq(6). Since, Velocity is a vector quantity, it has 3 direction (x,y,z). So we need to calculate implicitly in every direction. For simplicity, expanding above equation in x direction only we found:
Implicit derivation of above equation becomes
Rearranging above equation
Matrix formulation of A x = b is represented below
This formation of system of linear equation is solved by matrix iterative method called Conjugate Gradient Method. Thus, conjugate gradient method is solved for 3 times for three direction of velocity vector.
Due to the movement of particles in above implicit mehod, the particle number densities might be changed i.e. n* deviates from no. Fig. 5 shows the algorithm procedure of the present method.
In second stage, the temporal particle number densities are calculated from the temporal coordinates. The poisson equation for pressure is again calculated implicitly by conjugate gradient method (Koshizuka and Oka, 1996).
Rearranging above two equation
Matrix formulation of Ax = B is represented below
The right-handside of Eq.(18) is represented by the deviation of the temporal particle number density from the constant. Which maintains the particle number densities constant during the simulation. The left-handside of Eq. (18) is discretized by the diffusion model(Eq.(6)). Thus, we have simultaneous equations expressed by a linear symmetric matrix and they are solved by the Congjugate Gradient iteration method.
After calculating the pressure field the velocity is updated by the gradient of pressure which is governed by following Eq. (20).