[michaelmoor.aero]

Originally Posted by houkensjtu

I read through your description briefly and I think you should apply the boundary condition of Wall3(Which as you said, ensure the continuity condition) EVERY TIME before you solve the pressure correction equation. Coefficient of outlet velocity is cut as you said but they will enter the source term of pressure correction equation.
Also, I think BICGSTAB is ok but not necessary, because your grid number is certainly not so huge(maybe several hundreds or thousands). Usually backslash is far more quick for small problem. Matlab's backslash is very capable.
good luck!

Thanks again for the help! so you are saying that after solving the momentum equations yielding u* and v*, that i should apply the outlet boundary condition to these guessed velocity fields u* and v*? That makes sense, as that velocity cv (at i=NI), is the location of the east face of the pressure control volume at I=NI-1, i.e. the last pressure cv in the domain.

On one quick note about the initial guesses, which initialise the solution... at the moment, I have a guessed pressure field p* which is a linear distribution from west to east, i.e. 10-8-6-4-2-0, and I have also set the intial velocity field u* equal to the inlet vlaue, and v*=0. Anderson suggests putting in a velocity spike in a v control volume to create a 2D flow also.
If i set the u velocity field to zero (i=3 onward to i=NI-1), my outlet b.c returns NaN, and the solution fails due to this piece of code:

u(nX,2:nYp-1)=u(nX-1,2:nYp-1); % First extrapolate values from the domain

u(nX,2:nYp-1)=u(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1)); % This multiplier is the mass in divided by the mass out, which ensures continuity

So I guess the question is, do i assign velocity values that correspond with p*, or is it just sufficient that the u-velocities be non-zero?

This is all much appreciated and I am making more progress than I have in a long time! Best Regards, Michael

[houkensjtu]

Originally Posted by michaelmoor.aero

Thanks again for the help! so you are saying that after solving the momentum equations yielding u* and v*, that i should apply the outlet boundary condition to these guessed velocity fields u* and v*? That makes sense, as that velocity cv (at i=NI), is the location of the east face of the pressure control volume at I=NI-1, i.e. the last pressure cv in the domain.

On one quick note about the initial guesses, which initialise the solution... at the moment, I have a guessed pressure field p* which is a linear distribution from west to east, i.e. 10-8-6-4-2-0, and I have also set the intial velocity field u* equal to the inlet vlaue, and v*=0. Anderson suggests putting in a velocity spike in a v control volume to create a 2D flow also.
If i set the u velocity field to zero (i=3 onward to i=NI-1), my outlet b.c returns NaN, and the solution fails due to this piece of code:

u(nX,2:nYp-1)=u(nX-1,2:nYp-1); % First extrapolate values from the domain

u(nX,2:nYp-1)=u(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1)); % This multiplier is the mass in divided by the mass out, which ensures continuity

So I guess the question is, do i assign velocity values that correspond with p*, or is it just sufficient that the u-velocities be non-zero?

This is all much appreciated and I am making more progress than I have in a long time! Best Regards, Michael

Yeah you should ensure continuity on outlet every time before solving pressure correction.
As for NaN problem, I guess, because you applied 0 to all u velocity as initial guess, so when you apply

Originally Posted by michaelmoor.aero

u(nX,2:nYp-1)=u(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1)); % This multiplier is the mass in divided by the mass out, which ensures continuity

actually you are dividing by zero,or sth very near zero.
As far as I know...pressure initial guess is not so crucial problem.
good luck

[michaelmoor.aero]

I assume then that this is what you mean:

%% Apply Boundary Conditions
boundary_conditions();
%% Solve for u*
[ustar, diJ]=u_momentum();

%% Check continuity again
% u(nX,2:nYp-1)=ustar(nX-1,2:nYp-1); % First extrapolate values from the domain
% u(nX,2:nYp-1)=ustar(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1));
%% Solve for v*
[vstar, dIj]=v_momentum();

%% Solve for p'
[pdash] = pressure_correction(diJ, dIj);


Where now I have used the values of u* along i=NI-1, whereas in the boundary conditions, I had used the values of the initial guess.

[Mos]

Thank you very much for your notes!
My cavity code wasn't working till I used your suggested relaxation factors! And now it is working like a charm!!

[shivam vakhariya]

Originally Posted by houkensjtu

Hey guys!
I am a Master student and i'm programming SIMPLE method on myself.(duct flow)
I have been confused for a very very long time,and now I eventually get familiar with SIMPLE and I'd like to share some of experience to all of you.

1.Use Poiseuille flow to check your momentum equation.Apply the exact solution of velocity and pressure drop as boundary and initial condition.If your momentum equation is correct, it should never change the velocity field.Especially pay attention to the wall boundary condition.

2.Check mass balance over all flow field.This is VERY important and have NOT been stressed in most books,like Patankar's and HK.Versteeg's text book.
For example,u assumed a pressure field to be even as an initial condition.Then u solve the momentum equation,the result is that velocity field will definitely NOT fulfill the overall mass balance,i.e.,the mass flux on outlet will not equal flux you set on the inlet.
Terrible thing is,If u don't correct the velocity on outlet to make overall mass balance,the pressure correction equation will BLOW UP!You can never get a pressure correction solution, and the pressure correction will grow up and never stop.

By the way,you may wonder, that why all the sample codes are written for a cavity flow, one important reason I guess, cavity flow is INTERNAL flow and u never worry about overall mass balance,so there is no necessary to correct. However, if u modify them to a duct flow you will definitely FAIL,because of the blow up of pressure correction equation.

Hope this can help someone,good luck!

iam also facing same issue would you please share what action you taken for this or if it possible share your code