Hi all,

I'm trying to find a bug in my code. My code seems to work for a range of cases (airfoil in a channel) but messes up on the most simple! I have a uniform 2D Cartesian grid with U=10 and V=0 on the left boundary. All other boundaries are set to zero gradient. Obviously the converged solution should be U=10, V=0, P=const. everywhere (I did say it was a simple case!). I initialise my solution to U=10 etc. (so the initialised solution is actually the correct result!) and set it going. After 10 iterations or so, the two right corners of the domain have a large pressure. This pressure alternates between +ve and -ve from one iteration to the next and keeps increasing in magnitude until NAN (after 10,000 iterations of so).

I've been trying to find whats up with this for some time now with no luck so any tests or suggestions you can think of to try and isolate the problem would be much appreciated.

FYI: the code is steady, 2D, finite volume, SIMPLE, with Rhie-Chow.

Thanks in advance

Have you tried debugging the code using gdb or any other debugger ?

Sounds like a boundary condition problem. What boundary conditions are being applied at these corners? Is everything being extrapolated or is something anchored?

Hi, thanks for the replies. I do use gdb but have not found it to be helpful in this case as I don't really know what to be looking for.

As for boundary conditions, I have simply set U(i,j)=U(i-1,j) on the right boundary, U(i,j)=U(i,j-1) on the top boundary and similarly for the bottom boundary and similarly for the V BC's. On the left boundary, U=10, V=0 is used. Pressure is extrapolated on to all boundaries.

Upon further investigation, I've found that writing all the x,y values of each node to a file and then re-reading this information from the file strangely works SOMETIMES. The only difference this step makes is in the floating point values being slightly different e.g. 1.000000000 becomes 0.9999999999. It seems that some tiny computational round-off error is causing a pressure gradient that should not be there and for some reason this grows over several thousand iterations to become problematic. As I mentioned, the pressure goes from +ve to -ve from one iteration to the next so it seems the pressure correction step may be over correcting (despite it being under-relaxed). Humm, this one has really got me!

Thanks for any help again.

"he only difference this step makes is in the floating point values being slightly different e.g. 1.000000000 becomes 0.9999999999"

I can assure you this is not the cause of your very ligh values at the corners.

It seems you are doing some very silly mistake somewhere. (which is difficult for me to point out).

Check your code again very carefully, you might be skipping something thinking that it is obviously correct.

(once it took me 2 weeks to figure out what was wrong with my code, and the mistake was instead of 'i' i typed 'ii' in the loop (ii was another int defined). I kept overlooking it thinking this part is correct) and kept wondering where the bug is).

Check again. Pressure is very sensitive to mass imbalance and AP values. (I feel you are getting your AP values wrong).

You should not extrapolate the pressure on ALL boundaries. Either you should arbitrarily set the pressure to be zero at some arbitrary point of the domain, or else if you want only Neumann boundary conditions on the pressure, you should make sure that the pressure satisfies satisfies an integral condition. Check the literature for this latter condition.

Also, try setting v=5 boundary conditions at the inflow, and initialize v=5 everywhere accordingly.