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')
end
if MassResidual <= 10^-3
alphap = 0.1
disp('alphap is now 0.3')
end

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

I am also solving the equations using Jacobi:
kk=1;
ip=0;
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);
%end
%end
kk=kk+1 ;
mx = max(max(max(abs(pcheck-pdash))));
ip=ip+1;
if (kk>50&&mx<10^-6), break, end
end

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

I'd appreciate any assistance!

Best Regards,
Michael
0565113@gmail.com

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

This is the momentum equation that i solve for, including under-relaxation:
kk=1;
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);
end
end
kk=kk+1; % Counter
end

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);