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}$.
do nothing
conditions. $\endgroup$