2.1 Governing Equations
The governing equations of the present flow solver are the three-dimensional incompressible Navier–Stokes equations and the equation of continuity. No averaging or filtering process is involved and the flows are solved without any turbulence models:
$\begin{aligned} & \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial x}+v\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2}\right) \\ & \frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial y}+v\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}\right) \\ & \frac{\partial w}{\partial t}+u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial z}+v\left(\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}+\frac{\partial^2 w}{\partial z^2}\right)\end{aligned}$
(1)
$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$
(2)
where u, v, w, p, ρ and ν are the fluid velocity, pressure, density, kinematic viscosity, respectively. The fractional step method is applied for time marching. The grid is defined by an equally spaced Cartesian mesh. The convection term is evaluated by the second-order skew–symmetric scheme [2]. The pressure and diffusion term are calculated by the second-order finite-difference method. The Poisson equation of the pressure is calculated by the successive over-relaxation (SOR) method.
2.2 Immersed Boundary Method
Our flow solver represents objects by the level set and ghost-cell (GC) methods [3]. Each cell is classified as a fluid cell (FC), a GC, or an object cell (OC) by the classification manner of the level set function $\varphi$ expressed as eqx. (3).
$\begin{aligned} & 0<\phi_{F C}, \\ & -2 \sqrt{3} \Delta x \leq \phi_{G C} \leq 0, \\ & \phi_{O C}<-2 \sqrt{3} \Delta x .\end{aligned}$
(3)
The GCs are assigned in two layers as shown in Fig. 1. The flow equations in a GC are determined by the flow in an image point within the region occupied by fluid cells. To avoid recursive references, the probe length is set to 1.75 times the mesh size Δx. The primitive variables on the image point are interpolated of the values in the surrounding fluid. Finally, the value of a GC is defined by the value at the image point. To determine the velocity components of the GCs, non-slip boundary conditions are assumed at the object surface. The density and pressure are subjected to Neumann boundary conditions as follows
$\begin{aligned} & u_{G C}=u_{I P}-\frac{d_{I P}+\left|\phi_{G C}\right|}{d_{I P}}\left(u_{I P}-u_{I B}\right), \\ & v_{G C}=v_{I P}-\frac{d_{I P}+\left|\phi_{G C}\right|}{d_{I P}}\left(v_{I P}-v_{I B}\right), \\ & w_{G C}=w_{I P}-\frac{d_{I P}+\left|\phi_{G C}\right|}{d_{I P}}\left(w_{I P}-w_{I B}\right), \\ & P_{G C}=P_{I P} .\end{aligned}$
(4)
Figure 1. Domain around an object
2.3 The Motion of the Particle
The particles movements obey Newton’s equations for linear and transportation of a rigid body. Note that this simulation ignores the rotation of the particle. The velocities of two particles (labels 1 and 2) after their collision are, respectively, given by
$\begin{aligned} & \mathbf{u}_{p, 1}^{\prime}=\frac{m_{p, 2}}{m_{p, 1}+m_{p, 2}}(1+e)\left[\left(\mathbf{u}_{p, 2}-\mathbf{u}_{p, 1}\right) \cdot \mathbf{c}\right] \mathbf{c}+\mathbf{u}_{p, 1}, \\ & \mathbf{u}_{p, 2}^{\prime}=\frac{m_{p, 1}}{m_{p, 1}+m_{p, 2}}(1+e)\left[\left(\mathbf{u}_{p, 1}-\mathbf{u}_{p, 2}\right) \cdot \mathbf{c}\right] \mathbf{c}+\mathbf{u}_{p, 2},\end{aligned}$
(5)
where mp, up = (up, vp, wp) and e are the mass, velocity of the particle coefficient restitution respectively, and c is a standardization vector [4, 5]. The collision detection of particles 1 and 2 is determined by Pythagoras’ theorem as follows
$\sqrt{\left(x_2-x_1\right)^2+\left(y_2-y_1\right)^2+\left(z_2-z_1\right)^2} \leq \sqrt{\left(r_1+r_2\right)^2}$
(6)
where r is the particle radius.
After the particle–wall collisions, the particle velocity becomes
$\mathbf{u}_p^{\prime}=(1+e)\left[\left(-\mathbf{u}_p\right) \cdot \mathbf{c}\right] \mathbf{c}+\mathbf{u}_p$
(7)
To detect collision between a particle and the wall, the central position of the particle and the normal vector are obtained from the level set function. When the image point is appeared in the GC of the other object, the velocity and pressure components are defined by those of the nearest object.