[vaibhavkhanna]

Hey, I am quite new to CFD programming and was writing my own Matlab code for flow between parallel plates to get more familiar with the simple algorithm which I will be using for my thesis. I am trying to simulate the entrance development region which results in the parabolic velocity profile.
So I am getting NaN values after about 50 iterations of my solver. I have given constant velocity at inlet and u=0 at the walls along with zero pressure correction at the outlet. My guess is that it is related to some fundamental mistakes regarding the coefficients in the momentum equations since I am getting a centreline velocity greater than inlet before solving the pressure correction (continuity) equation but have been unsuccessful in finding that mistake.

[tas38]

Asking for others to sift through your codes for bugs, be they theoretical or in numerical implementation, is asking too much. I suggest you narrow your question or provide some sample results (e.g. contour plots before divergence) to better solicit assistance.

[vaibhavkhanna]

Can you provide some insight on how to debug my code? I have been trying to find the problem but there might be a mistake in the formulation? I am not asking you to debug my code. I am asking how I should proceed with debugging because I don't even know yet what is causing the error

[manle0312]

Your code is full of missing errors and I am not even digging into the main SIMPLE algorithm. Please dont name your matlab file as .txt. I am not going to do the whole thing for you but here are a few pointers:

1. You are calling the wrong function name from salgo1.m: it is gauss_seidel1, not gauss_seidel.

2. Your gauss_seidel1 function has an error on line 24, it is A(i,j), not a(i,j)

3. I can confirm that your coefficients Feu, Feu_r, Fev....etc. are NaN. My suggestion is to run it on a smaller grid with nicer number: like Lx = 1 and Nx = 5; so we have nicer numbers to work with.

4. Print stuff out. I see you only print out gauss iteration and count.. Your first big for loop to solve for the X momentum involves coefficients Feu, Feu_r....etc. Print them out along with iteration number of the for loop so you can see when things become NaN. Even better, make an if statement that if Feu is NaN then break then print out other important variables: iteration number, dx, dy...etc.


5. It looks like you have 4 BCs for top, bottom, left and right wall. For your first implementation, keep them simple: perhaps only left and right wall are constants. I dont know, I dont work with parallel pipe but make keep them simple. The best thing is to get your SIMPLE algorithm works first before tweaking.

6. Make sure your Gauss Seidel scheme works perfectly first before plugging it in the SIMPLE algorithm.

TL,DR: just print stuff out, and keep track of when stuff becomes NaN. When in doubt, try to do a hand calculation of the expression that becomes NaN and backtrack.
A side note on gauss_seidel1.m, not sure why you do: PHI = phi, use a different variable name would help I think. Matlab is case sensitive and you might have written a wrong variable somewhere.

[vaibhavkhanna]

Thanks Man Le for your advice. I found several bugs and now it finally seems to converge. However, the issue I am facing with my simulations now is that the residuals (L2 norm of mass imbalance) for different Re and different fluids all seem to oscillate about a fixed value instead of decreasing as the iterations proceed. The termination criteria for my simulation is never met. What are your thoughts on this?
Here are some of the residual plots from my simulation. The velocity values are oscillating about the right solution (the centreline velocity is 1.5 times the inlet velocity) but the residuals do not decrease with iterations.

https://drive.google.com/file/d/16fb...ew?usp=sharing