MAC Projection in Projection method?

Question

My question concerns the following paper: A Conservative Adaptive Projection Method for the Variable Density Incompressible Navier–Stokes Equations (http://www.sciencedirect.com/science/article/pii/S0021999198958909)

It describes a second order accurate projection method to solve incompressible Navier Stokes equations. Usually projection methods are based on implicit schemes, but as far as I can figure out, this is based on explicit methods. My motive is to solve incompressible NS equations on a uniform grid using this method.

To solve for the non-linear convection term of momentum equation, the authors first extrapolate the cell-centered velocities to cell faces using predictor-corrector upwind schemes. This is an extrapolation in space (from cell center to cell face) and also an extrapolation in time (from time level $n$ to level $n+\frac{1}{2}$). After they calculate normal velocities on cell faces, they enforce the divergence constraint at level $n+\frac{1}{2}$. For this, they are using the MAC projection in which they write:

where D implies divergence, G implies gradient. This equation is used to solve for $\phi^{MAC}$ But I do not understand how to get $\phi^{MAC}$ from this equation. The term on the right contains two terms in two dimensional case (corresponding to u and v velocity each). How should I write the term on LHS in order to get $\phi^{MAC}$ from this equation?

Answer 1

Without entering much into the details of your specific problem, let's write your problem in a simplified manner:

$\nabla \cdot \left(\frac{\nabla \phi}{\rho}\right) = F $

Note that $F(U,\delta t)$ is already given by the predictor step. The important thing to remember is that $U$ is given on the faces of the cell, thus some interpolation has to be made in order to have all coefficients "in place". Assuming that the coefficients are sufficiently smooth, the Gauss-Seidel iteration scheme in one of the directions (x, say) for the problem in $\Omega=[0,1]^2$ and $u=\phi/\rho $:

$d^2u/dx^2 = f$

$u(0)=u(1)$ -Periodic B.C.

is, using a second order centered discretization, written as:

$u_i^{(k+1)} = 2u_i^{(k)} - \frac{1}{2}(h^2F - u_{i-1}^{(k+1)}-u_{i+1}^{(k)})$

Where $h$ is the grid spacing and $u^{(k)}$ is the value of the grid function in the kth iteration. From here, is up to you and the problem at hand which solver to use.

A good, general reference is the book by Saab (2003). The Author has put this version for free in his webpage. On multigrid, the classic paper of Achi Brandtl (1977) is pretty good on the subject.

APPENDIX

For a linear system $A x = b$, where $A = (a_{ij})$; a fixed point iteration scheme, known as Gauss-Seidel, can be written as:

$ x^{(k+1)}_i = x^{(k)}_i + \frac{1}{a_{ii}}\left(b-\sum_{i=1}^{j-1}a_{ij}x^{(k+1)}_i -\sum_{i=1}^{j}a_{ij}x^{(k)}_i \right) $

Note that, in fact, no iteration over $i$ is needed since the resulting system is explicit in nature (i.e.: the $k+1$ value on the RHS is known from the previous iteration).