Fluid Time Step

In the absence of reactions, we assume that the fluid density is unchanged.

We compute the fluid volume fraction directly from the particle locations.

Thus here we focus on the discretization of the momentum equation

In the predictor

• Define $U^{MAC,n}$, the face-centered (staggered) MAC velocity which is used for advection, using $U_g^n$

• Define an approximation to the new-time state, $({\epsilon }_g{\rho }_g{U}_g{)}^{\ast }$ by setting

$$\begin{array}{rl}({\epsilon }_g{\rho }_g{U}_g{)}^{\ast }& =({\epsilon }_g{\rho }_g{U}_g{)}^n-\mathrm{\Delta }t(\mathrm{\nabla }\cdot ({\epsilon }_g{\rho }_g{U}^{MAC}{U}_g)+{\epsilon }_g\mathrm{\nabla }{p_g}^{n-1/2})\\ & +\mathrm{\Delta }t(\mathrm{\nabla }\cdot {\tau }^n+\sum _p{\beta }_p(V_p-{U_g}^{\ast })+{\rho }_g{\epsilon }_{g}g)\end{array}$$

• Project $U_g^{\ast }$ by solving

$$\mathrm{\nabla }\cdot \frac{{\epsilon }_g}{{\rho }_g}\mathrm{\nabla }\varphi =\mathrm{\nabla }\cdot (\frac{1}{\mathrm{\Delta }t}({\epsilon }_g{U}_g{)}^{\ast }+\frac{{\epsilon }_g}{{\rho }_g}\mathrm{\nabla }{p_g}^{n-1/2})$$

then defining

$$U_g^{\ast \ast }=U_g^{\ast }-\frac{\mathrm{\Delta }t}{{\rho }_g}\mathrm{\nabla }\varphi$$

and

$${p_g}^{n+1/2,\ast }=\varphi$$

In the corrector

• Define $U^{MAC,\ast \ast }$ at the “new” time using $U_g^{\ast \ast }$

• Define a new approximation to the new-time state, $({\epsilon }_g{\rho }_g{U}_g{)}^{\ast \ast \ast }$ by setting

$$\begin{array}{rl}({\epsilon }_g{\rho }_g{U}_g{)}^{\ast \ast \ast }& =({\epsilon }_g{\rho }_g{U}_g{)}^n-\frac{\mathrm{\Delta }t}{2}(\mathrm{\nabla }\cdot ({\epsilon }_g{\rho }_g{U}^{MAC}{U}_g{)}^n+\mathrm{\nabla }\cdot ({\epsilon }_g{\rho }_g{U}^{MAC}{U}_g{)}^{\ast \ast })+\\ & +\frac{\mathrm{\Delta }t}{2}(\mathrm{\nabla }\cdot {\tau }^n+\mathrm{\nabla }\cdot {\tau }^{\ast \ast \ast })+\mathrm{\Delta }t(-{\epsilon }_g\mathrm{\nabla }{p_g}^{n+1/2,\ast }+\sum _p{\beta }_p(V_p-{U_g}^{\ast \ast \ast })+{\epsilon }_g{\rho }_{g}g)\end{array}$$

• Project $U_g^{\ast \ast \ast }$ by solving

$$\mathrm{\nabla }\cdot \frac{{\epsilon }_g}{{\rho }_g}\mathrm{\nabla }\varphi =\mathrm{\nabla }\cdot (\frac{1}{\mathrm{\Delta }t}({\epsilon }_g{U}_g{)}^{\ast \ast \ast }+\frac{{\epsilon }_g}{{\rho }_g}\mathrm{\nabla }{p_g}^{n+1/2,\ast })$$

then defining

$$U_g^{n+1}=U_g^{\ast \ast \ast }-\frac{\mathrm{\Delta }t}{{\rho }_g}\mathrm{\nabla }\varphi$$

and

$${p_g}^{n+1/2}=\varphi$$