CFD Online Discussion Forums

CFD Online Discussion Forums (
-   Main CFD Forum (
-   -   mesh size & urf SIMPLE FVM MATLAB ( August 14, 2012 08:52

mesh size & urf SIMPLE FVM MATLAB
Hello everyone,
I have developed a laminar 2D solver for Poiseulle flow using FVM and SIMPLE, and the results match nicely to the Hagen-Poiseulle velocity profile, but if i increase the mesh size, the solution blows up... My results are coming out nicely, but i cannot seem to be able to increase the mesh size beyond 36*36... I have to decrease the pressure under-relaxation to 0.01 to get it to iterate, but then it just oscillates at a continuity residual of around e-2... do you know the reason for this and why it happens? or how to fix it?

I am trying to maybe implement a method such that:
if MassResidual <= 10^-2
alphap = 0.05
disp('alphap is now 0.1')
if MassResidual <= 10^-3
alphap = 0.1
disp('alphap is now 0.3')

in an effort to try and march the process along? is this a correct approach?

I am also solving the equations using Jacobi:
while kk<=200;

pcheck = pdash;
%for i=2:nXp-1;
%for j=2:nYp-1;
%pdash(i,j)=(bdash(i,j)+saE(i,j)*pdash(i+1,j)+saW( i,j)*pdash(i-1,j)+saN(i,j)*pdash(i,j+1)+saS(i,j)*pdash(i,j-1))/saP(i,j);
kk=kk+1 ;
mx = max(max(max(abs(pcheck-pdash))));
if (kk>50&&mx<10^-6), break, end

and for residuals:
for i=2:nXp-1;
for j=2:nYp-1;
MassResidual=max([MassResidual abs(bdash(i,j))]);

I'd appreciate any assistance!

Best Regards,
Michael August 14, 2012 10:12

Apologies, I did not add all of the relevant information...

This is the momentum equation that i solve for, including under-relaxation:
while kk<=100;
%Solve iteratively using the Jacobi Method to yield u, which will now
%be the most up to date value, i.e. u*
for i=3:nX-1;
for j=2:nYp-1;
u(i,j)=((biJ(i,j)+uaE(i,j)*u(i+1,j)+uaW(i,j)*u(i-1,j)+uaN(i,j)*u(i,j+1)+uaS(i,j)*u(i,j-1) + (1-alphau)*uaP(i,j)*u(i,j)/alphau))*alphau/uaP(i,j);
kk=kk+1; % Counter

and i correct the pressure with under-relaxation:
%% Correct the Pressure, including under-relaxation

pnew(2:nXp,2:nYp-1) = pstar(2:nXp,2:nYp-1) + alphap*pdash(2:nXp,2:nYp-1);

All times are GMT -4. The time now is 11:51.