Implementing the pressure correction method using finite elements

Question

Ok so I am nearing the completion of my finite element Navier-Stokes solver that uses the $\theta$-method for time stepping and the pressure correction method for the pressure. I am following the notes http://ta.twi.tudelft.nl/users/vuik/burgers/fem_notes.pdf by Ir. A. Segal (page 74-75). So far everything seems to work except there appears to be a bug in how I implement the pressure projection. Allow me to explain using some pictures...

Consider flow going over a step. When I run my simulation I get:

At first glance everything appears to be ok. However you may notice that in the bottom left corner the pressure is slightly larger (I know its hard to see). This bottom corner is exactly where I specified my one Dirichlet boundary condition on the Poisson equation necessary for the system to be non-singular. If we look at the corresponding u-velocity the error becomes much more clear:

Clearly the velocity should not be peaking at the bottom left-hand corner. I have tried moving the Dirichlet condition away from the corner, but this does not fix the problem. For example moving it up the left-hand side boundary results in the pressure and u-velocity looking like:

I think the issue is with how I am solving the Poisson equation that arises from the pressure correction method. In particular with how I am implementing the boundary conditions. Ok lets look at the math now (hopefully someone can shed some light on what I am doing wrong). From the above cited paper, we need to solve the equation:

$\nabla^{2}(p^{n+1}-p^{n}) = \frac{\nabla{}\cdot{\vec{u}^{*}}}{\theta{}dt}$

Doing this with finite elements we multiply by a test function $\phi_{j}$ and integrate over the domain (I replace $\frac{1}{\theta{}dt}$ with the constant c):

$\int\int_{\Omega}{\nabla^{2}(p^{n+1}-p^{n})\phi_{j}}dA = c\int\int_{\Omega}{\nabla{}\cdot{\vec{u}^{*}}\phi_{j}}dA$

Using Greens theorem,

$\int\int_{\Omega}{\nabla{}\cdot{\vec{F}}}dA = \int_{\partial\Omega}{\vec{F}\cdot{\vec{n}}}dS$

We get:

$\int\int_{\Omega}{\nabla^{2}(p^{n+1}-p^{n})\phi_{j}}dA = -\int\int_{\Omega}{\nabla{(p^{n+1}-p^{n})}\cdot{\nabla{\phi_{j}}}}dA + \int_{\partial\Omega}{(\nabla{(p^{n+1}-p^{n})\phi_{j}})\cdot{\vec{n}}}dS$

Now we prescribe Neumann boundary conditions for pressure on all boundaries except for one point on the boundary for which we prescribe a Dirichlet condition. This is necessary since otherwise the resulting Poisson matrix system would be singular. Now my theory as to what is going wrong is currently in the above equation I just get rid of the last term (because I assume a Neumann boundary condition on all boundaries for pressure, i.e. $\nabla{(p^{n+1}-p^{n})}\cdot{n} = 0$. However I am now beginning to think that I must account for that one point on the boundary. So for that one Dirichlet point I need to retain:

$\int_{\partial\Omega}{(\nabla{(p^{n+1}-p^{n})\phi_{j}})\cdot{\vec{n}}}dS$ ?

Does anyone know if this is right? Has anyone come accross similar problems when implementing the pressure correction method? I have been racking my brain trying to figure out what could be going wrong. In the equation:

$\nabla^{2}(p^{n+1}-p^{n}) = \frac{\nabla{}\cdot{\vec{u}^{*}}}{\theta{}dt}$

what, if any, boundary conditions do you apply to $\frac{\nabla{}\cdot{\vec{u}^{*}}}{\theta{}dt}$? Should it be zero on the boundary? How does one enforce the Dirichlet condition at one point for pressure in the global stiffness matrix? Any help is really appreciated. I have run out of ideas!

Edit

Ok so after reading the comments and running so more tests I am more convinced my problem is in how I implement the one Dirichlet boundary condition on the pressure. So far I have been following http://math.mit.edu/~gs/cse/codes/mit18086_navierstokes.pdf (pg 12-13) in prescribing the one Dirichlet condition on pressure. Basically if Kp is my global stiffness matrix for the Poisson equation then I simply add a 1 to the first entry, i.e. Kp[0][0] = Kp[0][0] + 1. Is this correct? I dont do anything special on the r.h.s of the Poisson equation, i.e. $\frac{\nabla{}\cdot{\vec{u}^{*}}}{\theta{}dt}$.

Answer 1

It depends on how you set your Dirichlet conditions, but they are often equivalent to setting $\phi_j=0$ on those nodes, so you shouldn't have to worry about that equation as long as you set the Dirichlet condition correctly.

Now, that being said, I've never been a huge fan of pinning one pressure node in traditional mixed formulations for Navier-Stokes using Taylor-Hood elements. You have other choices. If your linear solver is a Krylov method, you can usually get away with doing nothing. That is, form the matrix with the Neumann/natural boundary condition everywhere, and just feed it to the solver. As long as your preconditioner is non-singular, you will usually iterate to a good solution with the Krylov method "picking" your pressure level for you. Your other choice is to impose an average pressure zero condition across the whole domain $$ \int_\Omega p \; dx =0 $$ which, when you plug in your pressure representation will give you a constraint equation for your system.

Edited to add based on the question edit:

Suppose you have the $j^{\rm th}$ node where you want to apply the boundary condition $p_j=c$, then you can replace row $j$ of your matrix with all zeros and a 1 on the diagonal. Then you put $c$ in the RHS vector on row $j$ as well (don't add, replace). Now on whatever other rows $p_j$ appeared in your LHS, i.e. non-zeros row entries in original column $j$, you need to subtract $c \times K_p[i][j]$ from the RHS vector on row $i$. You then set $K_p[i][j] = 0$ in the original matrix. That will enforce the Dirichlet condition.

You can get away with not modifying anything other than row $j$, but you will destroy the symmetry of your matrix. If your solver required symmetry (like Conjugate Gradient or Cholesky), then you have to pick a non-symmetric solver. If you can get the bookkeeping right, it's probably better to modify the matrix.

Answer 2

You should use the discrete counterpart of the pressure-correction method so that no boundary conditions are needed for pressure (like in the original problem). Then this problem is avoided.