# MAC Projection in Projection method?

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?

• it seems you need an iterative solver to compute the numerical "pressure" from this equation which, in principle, is the laplace equation. BTW, what is the E-->C superscript? Sep 21 '15 at 13:14
• @SantiagoLópezCastaño I absolutely have no idea. The paper writes "we apply the MAC projection to the face-based velocity field before construction of the conservative updates" So, C might stand for conservation? But I can not be very sure. Sep 21 '15 at 13:18
• It seems to me that the superscript E-->C refers to an interpolation from the faces to the cell center of the grid, and C-->E viceversa. Note that the contravariant (auxiliary,etc) )fluxes (U) are on the cell faces, whereas the comp. Pressure is usually located on the cell center, along with the velocity terms. In any case, the Poisson condition needs of an iterative method to be solved. Sep 21 '15 at 13:25
• As far as I know, pressure in this case is specified at cell corners and is also staggered in time while velocity, density are defined at cell centers and is defined at integer times (not staggered). Sep 21 '15 at 13:40
• A classic iterative solver for the Poisson equation (I said laplace earlier, I take that back) is a v- (or w-) cycle multigrid using Gauss-Seidel as smoother, and some sort of Krylov Solver as accelerator (to speed up the convergence). Note that for simple cases, (2D, uniform orthogonal grids) using such solvers can be a bit of an overshoot; instead, a simple SSOR, SOR, or Gauss-Seidel will suffice. Sep 21 '15 at 13:44

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).

• Thanks for the detailed analysis. Though I have some queries: 1) Why did you let $u = \frac{\phi}{\rho}$ and does this u corresponds to x-velocity or is it just any variable? 2) Since F is a scalar (divergence of some vector). Therefore when I use this scheme in two dimension, say x and y, how do I get the value of $\phi$ at any point? Do I use a second order discretization to laplacian on LHS and calculat $u_{i,j}$ from that? Sep 22 '15 at 3:18
• @TanmayAgrawal Responding to your queries: (1) I just used another name for $\phi/\rho$; (2) Indeed, F is a scalar, but $U$ is a vector so the divergence has to be taken in both directions. Note that you will get $\phi/\rho$ from the scheme above. In the simple scheme shown above, you just have to move first in the i-direction and then in the j-direction. Sep 22 '15 at 7:36
• @TanmayAgrawal Another thing: discretize the 2D Laplacian, and calculate the coefficient matrix $A$ as shown in the answer. Sep 22 '15 at 7:59