[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!!