Projection Method: Boundary condition on intermediate velocity field

Question

I'm trying to solve variable density and viscosity Navier-Stokes equation using lagged pressure projection method. I'm solving for cavity problem as a test case now (once I get projection right, I will be solving a multiphase problem). The problem is I get oscillations developed in the solution as I march in time. .

I'm using collocated grid for solving this problem.

I'm solving variable density, viscosity N-S equation: $$ \dfrac{u^* - u^n}{\Delta t} + \left(u_i\frac{\partial u_i}{\partial x_j} \right)^n = \left(-\frac{1}{\rho}\frac{\partial p}{\partial x_i} \right)^n + \frac{1}{\rho} \left[\frac{\partial}{\partial x_j}\left( \mu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right) \right]^n $$

and on boundary, $U^* = U^{n+1}$.

I'm approximating convective terms with 1st order godunov as of now, diffusion term and pressure gradient with central difference.

Then I'm solving the projection step,

$$ \frac{\nabla.u^*}{\Delta t} = \nabla.(\frac{1}{\rho}\nabla\phi) $$

using boundary condition $$ \frac{\partial \phi}{\partial n} = 0 $$

I'm discretizing $\nabla.u^*$ as $\dfrac{u^*_{i+1,j} - u^*_{i-1,j}}{2\Delta x} + \dfrac{v^*_{i,j+1} - v^*_{i,j-1}}{2\Delta y}$ and applying boundary condition using ghost node. and for pressure, I'm discretizing as $\dfrac{\left(\frac{1}{\rho}\frac{\partial p}{\partial x} \right)_{i+1/2,j} - \left(\frac{1}{\rho}\frac{\partial p}{\partial x} \right)_{i-1/2,j}}{h} + \dfrac{\left(\frac{1}{\rho}\frac{\partial p}{\partial y} \right)_{i,j+1/2} - \left(\frac{1}{\rho}\frac{\partial p}{\partial y} \right)_{i,j-1/2}}{h}$

For solving the singular poisson equation, I'm removing the nullspace from the RHS ($\nabla.u^*$), operator matrix ($\nabla.(1/\rho$)). I'm also removing the nullspace from the solution vector.

I'm using petsc to solve the poisson equation and using GMRES ksp solver.

I,m updating velocity using $U^{n+1} = U^* - \Delta t*(\nabla\phi)$ and $P^{n+1} = P^* + \phi$. I'm fine with first order accurate pressure.

However, I get the oscillations in my solution. And then the solver diverges if kept iterating. I strongly feel there is something wrong with the boundary conditions, especially on $U^*$ as the solution is very sensitive to what I specify $U^*$ on boundaries. I also tried with $U^* = 0$ but no improvement.

Is there anything wrong in the method I'm using here? It'll great help if you can point it out.

Thank you.

See results here:

Found the solution:

I'm aware of the fact that the staggered grid is the most straight forward solution to this checkerboard problem but I wanted to stick to the collocated grid as it is easier to code and less prone to error, plus having velocities at cell centers will greatly reduce the level set formulation efforts compared to staggered grid.

Now, instead of lagged-pressure projection method, I've changed it to the (Kim and Moin's) pressure-free projection method which is working great as of now. However, I've coded it for constant density and viscosity Navier-Stokes I'm hoping it to work for the variable density and viscosity as well.

While implementing, now there is no pressure gradient term in the momentum equation and I used following boundary conditions on intermediate velocity field: $$ U^*.\hat{n} = U^{n+1}.\hat{n} $$ $$ U^*.\hat{t} = (U^{n+1} + \Delta t*\nabla\phi).\hat{t} $$ where $\hat{n}$ and $\hat{t}$ are normal and tangent vectors. And poisson equation has the neumann B.C. as before.

As I have to get going with my project, I have not thoroughly investigated the error in lagged-pressure projection method But I will do so as soon as I get some time.

Thank you.

Here's the solution (at t = 10sec with Re400) with pressure-free projection method:

Answer 1

What you are observing is odd-even decoupling between the velocity and pressure. This is a well-known problem on collocated grids. At any point on the grid, the gradient of the pressure is calculated using $p_{i+1}$ and $p_{i-1}$:

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} + (u_i^n.\nabla)u_i^n = -\frac{p_{i+1}-p_{i-1}}{2\Delta x}+\nu\nabla^2 u_i^n$$

Hence the velocity $u_i$ only depends on the pressure at the adjacent grid points, and not on $p_i$ itself. This can cause the numerical pressure to behave strangely. As you can see in your results, a "checkerboard" pattern is observed, where only the pressures at alternate grid points influence each other.

This problem can be solved by using a staggered grid, where the velocities and pressures are calculated at different points in the domain. The most commonly used staggered grid on a Cartesian mesh is the one introduced by Harlow and Welch [1], but it will make your code more complex. If you want to make minimal changes to your code, then use the same collocated grid for the velocities, but compute the pressures at the center of each cell. Take care to set up your pressure equation carefully, and interpolate the velocities as required when you calculate the right-hand side $\nabla\cdot u^*$.

Here is another page that discusses this phenomenon: http://www.imm.dtu.dk/~pcha/Projekter/Shock/main/node17.html

Regarding the boundary conditions for the intermediate velocity field, rather than $u^*=0$, you should be using $u^*=u^{n+1}$. In the problem you have considered, this makes a difference on the top surface where the lid is moving. For more information about the boundary conditions using a projection method, read the paper by Perot [2] and follow the references.

[1] Harlow, Francis H., and J. Eddie Welch. "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface." Physics of fluids 8.12 (1965): 2182.

[2] Perot, J. Blair. "An analysis of the fractional step method." Journal of Computational Physics 108.1 (1993): 51-58.

Answer 2

I think your problem lies in the boundary conditions for $\mathbf{u^*}$. Since your grid is collocated, $\mathbf{u^*}$ is stored at cell centers; and the boundary lies half a cell distance from the place where $\mathbf{u^*}$ is actually computed. Thus a higher-order BC for $\mathbf{u^*}$ should do the trick. Also, I would like to mention that the Kim-Moin projection method does not yield a high order of accuracy in pressure [1]. [1] will also help with many other order of accuracy issues you might face.

[1] Brown D. L., Cortex R. and Minion M. L. "Accurate projection methods for the incompressible Navier-Stokes equations"