Problem solving pressure-correction equation
 

Hi everybody,
I'm trying to simulate an incompressible unsteady flow. I assume inlet BC (given velocities), wall, symetry plane and outlet (Neumann condition for velocities). So I use zero normal derivative condition for pressure correction on whole boundary, but I'm not sure that this is right. The pressure-correction equation is to be solved by Gauss-Seidel or SOR solver. Solver behaves quite a curious way. It does not converge neither diverge. More accurately residuum = A fi - q is falling for a while and then it stops on a value (too high). It seems problem is in a source term but I can't resolve it. It's needed for a poisson equation with Neumann condition source term to be zero in average. I tried to fix it by correcting computed mass fluxes on outlet to giver the same total mass flux as on inlet in accordance with sugestion in Ferziger, Peric. This treatment ensures that sum of source terms is about 1.e-13. I'm not sure whether it isn't too much and doesn't cause the bad convergence. Or could be the cause of the problem large gradient of the source term near outlet boundary caused by this correction?
What do you think?

Thank you for your suggestions

yes you are right. pressure correction equation with a insulation bc

Hi Alex,

I would suggest that you enforce Dirichlet conditions on some parts (at least one point) of your whole boundary set. Then, the problem is well-posed and you are not bound by the compatibility integral constraint.

If the convergence does improve (and goes down to machine-accuracy), then it is probably this integral condition which is not fullfilled. In this case, make sure that your discrete scheme is fully conservative... that might help.

If it does not improve, then the problem is somewhere else.

Hope this helps...

Franck, thank you! Introducing dirichlet bc in only single point really improved accuracy to machine precision. What did you mean by conservation?

Hi Alex.

It seems that the integral compatibility constraint is not fulfilled then...

What I mean by conservative: since the laplacian operator is formally the divergence of a gradient, you can express it as flux differences (of the gradient component normal to the cell faces) between the cell faces of your mesh cells. Conservative means that the flux leaving one cell should be identical (in your discrete scheme) to the flux entering the neighbouring one. If your equation reads like this:
div(grad(fi))=q
adding this over all cells of the entire mesh will result in:
sum of q over all cells = sum of normal flux over all boundary faces
which is the integral constraints which should be fulfilled. If the scheme is not conservative, then you cannot be sure that this constraints holds...

As you wrote before, you should as well ensure that your Neumann boundary conditions are set so that the constraint is fulfilled. You can do that by correcting them as you mentionned in your previous post.

Does it make sense?