Hello all,

I have still not got around this problem.

Arjun - could you please advise me on how you would go about fixing the value of pressure at one node? No matter how I do this it always leads to the same problem of increased pressure gradients and therefore increased velocities and continuity is no longer obeyed in the cell at which I fix the pressure.

I have tested the different iterative solvers available in Matlab. They give the same results. I have also used the approach of Allan from another thread which is:

'setting a Dirichlet condition at an interior node and after that solve another Poisson equation in order to remove the errors around that node at which the Dirichlet condition was imposed.'

Although this does remove the sharp pressure gradients, the field is still not correct as it induces false velocities of around 10E-5 over the entire field and continuity is no longer obeyed, particualrly at the boundaries.

I have checked my discretization of the poisson equation many times and I am sure it is correct. I am not using a software package for this. I have written all the code myself apart from the iterative solver which is the biconjugate gradient method taken from Numerical Recipes in FORTRAN. I have no pressure correction equation. When the Navier-Stokes equations were first temporally discretised the pressure terms were treated implicitly which leads to the implicit pressure poisson equation with Neumann boundary conditions which is solved first, then the explicit velocity and density updates are carried out. This completes a time step.

I am trying to recreate the code from the paper http://www-math.mit.edu/~bush/fub.pdf

@Tom, thanks for posting the paper, now it would be easy to understand what you are doing.

I will try to write in more detail as i get time but here is a question.
Are you solving page 79 A4 for pressure?

If you are then it seems the paper also mentions fixing pressure at one point.

"The boundary conditions for (A4) are of Neumann type, provided in terms of
^ n
rp+ by the velocity boundary conditions, and p+ = 0 is set at some point
for uniqueness."


This is exactly what I do with my FVM code too. But I fix the value at one point on any one of the boundary.

Important: I noticed that if I fix value at flow out then it usually works better than just randomly picking point and fixing it.

I will write more later.

At a different point in the paper is states that p=0 should be set at p(1,1). However, this is not an outflow, there are no outflows as the system is closed, it describes the axisymmetrical spin-up from rest in a cylindrical container. The fluid is bound between 0<=r<=R and 0<=z<=2H. However, only the bottom right corner, 0<=r<=R and 0<=z<=H,of the container is modelled due to the symmetry of the problem.

Therefore at r=0 and z=H there are symmetry conditions, which gives pressure boundary conditions of dp/dr=0 at r=0 and dp/dz=0 at z=H. The remaining boundary conditions are found, as you have qouted, by inserting the velocity boundary conditions into equation (A1). This gives dp/dr=gamma*phi*r at r=R and dp/dz=-Froude^-2 *phi at z=H, these come from the no slip condition of velocity being zero at the walls.

Therefore, this full Neumman boundary condition equation leads to the singular matrix I am trying to solve. As the paper states, p+=0 is set at p(1,1) for uniqueness which allows the iterative biconjugate gradient solver to find a solution. However, setting p+=0 at this or any other point leads to the problem I have already described whereby I have large, incorrect pressure gradients concentrated around the static pressure point. I have gone through my code countless times to check that everything is correct, I have even rewritten my code completely from the start and it has led to the same issues.

Thank you for reading and replying arjun. I am desperate to sort this problem out. I am only a lowly Civil Engineering graduate now 2 years into my PhD. Before starting this I had never done any numerical modelling before and also my knowledge of matrices is not the best as maths is not a real priority on an engineering course - its more a case of 'here is the equation you need and this is how to use it'.

A quick addition to the above post.

You said that you fix the value at flow out and this works well. Wherever I fix the value it creates the effect of a sink or outflow with fluid being pulled towards the point. Maybe this is why you have not had problems as fixing the value of pressure at your flow out location would be a correct boundary condition for your physical problem?

I tried to read the pdf you posted and here is my initial comment.

Do you have to do this thing exactly the way this paper mentions or are you free to do the same in any other way too.

why I am asking this is that by first look it feels that you are supposed to solve navier stokes with a separate equation for density function.

The paper you are using as base write navier stokes in axis-symetric form and creates a procedure for solving it. (Please correct me if I am wrong).


Further by looking at the procedure, I expect the solution process to get into trouble. Why?? The reason is it seems in his procedure pressure is solved and then it is used the way he gets it.

If I say the same thing for SIMPLE algorithm then it would mean using pressure correction without any under-relaxation factor. which makes solution process very difficult to converge.

Or in other words if this procedure could be converted to a pressure correction form it would give you a lot more stability and lot less headache with better convergence.

More later on.

You are correct on all these points.

It would be good if i could match the results of this paper using the exact same method but I am free to solve it another way. I just need to make sure my numerical results are in good agreement with experimental ones.

I have just started reading up on how to use the SIMPLE algorithm.

Don't bother to read about it as yet. If you can write down in equation forms with Navier stokes, I can add it onto iNavier and can solve that. If results match then you are free to write your own code or just use iNavier. With iNavier you won't be confined to cartesian mesh and will be able to use any grid.

Do you mean you want me to give you my update equations for the pressure, velocity and density? If so, in what form? Differential or spatially discretized?

Also, use of iNavier is not an option for me.

Sorry for not writing, just got a bit busy with things.

If you want to do with your code then my advice would be convert that into pressure correction type equation and correct pressure slowly. Off course it also means solving pressure correction on each time steps few times, that is may be 4-5 times before moving to next time level.

S.o.r
 

Just incase anyone is having similair problems I have now found my way around this issue.

I spent so long trying to get the Biconjugate Gradient Method to work as I was trying to recreate the code of someone else. However, I gave up trying to with BiCG and decided to try some alternative iterative solvers. The first one I tried was Successive Over Relaxation as it is a very simple method, easily programmed and tested.

SOR has worked brilliantly for me. There is no need to fix the pressure at an abitrary point in order to obtain a solution and as such I am not experiencing the problem of false fluid flow which is described earlier in this thread. Also, I am finding that SOR is running quicker than BiCG with a relaxation parameter of 1.25 which was just my first test. Therefore it may run even faster once I have had a chance to play around with the relaxation parameter.

Basically, SOR has solved my problems and I am now able to get proper results from my code and as such I can get on with my life.

Thanks again to Arjun who is a very willing helper, and also thanks to anyone else that has posted their advice. I hope this thread can help others in the future.

Tom.