Help with lid-driven cavity problem
 

Hi, im using a collocated non uniform grid for solving the lid-driven cavity problem and i dont know very much about CFD (i started this week learning!), i hope if someone can help me with my code (Matlab):

function rhs=Advection(q,jcb,csi_x,eta_y,im,jm)

for i=1:im
for j=1:jm
U1(i,j)=q(i,j,2)*csi_x(i,j);
F(i,j,1)=U1(i,j)/jcb(i,j);
F(i,j,2)=(q(i,j,2)*U1(i,j)+q(i,j,1)*csi_x(i,j))/jcb(i,j);
F(i,j,3)=q(i,j,3)*U1(i,j)/jcb(i,j);

U2(i,j)=q(i,j,3)*eta_y(i,j);
G(i,j,1)=U2(i,j)/jcb(i,j);
G(i,j,2)=q(i,j,2)*U2(i,j)/jcb(i,j);
G(i,j,3)=(q(i,j,3)*U2(i,j)+q(i,j,1)*eta_y(i,j))/jcb(i,j);
end,end

for i=2:im-1
for j=2:jm-1
rhs(i,j,1:3)=0.5*(F(i+1,j,1:3)-F(i-1,j,1:3)+G(i,j+1,1:3)-G(i,j-1,1:3));
end,end


function rhs=Viscous(q,Re,jcb,csi_x,eta_y,im,jm,rhs)

for j=2:jm-1
for i=2:im-1

g11=(0.5*(csi_x(i+1,j)+csi_x(i,j)))^2;
g22=(0.5*(eta_y(i,j+1)+eta_y(i,j)))^2;
J=0.5*(jcb(i+1,j)+jcb(i,j));

FV(i,j,1)=0;
FV(i,j,2)=1/Re*g11/J*(q(i+1,j,2)-q(i,j,2));
FV(i,j,3)=1/Re*g11/J*(q(i+1,j,3)-q(i,j,3));

GV(i,j,1)=0;
GV(i,j,2)=1/Re*g22/J*(q(i,j+1,2)-q(i,j,2));
GV(i,j,3)=1/Re*g22/J*(q(i,j+1,3)-q(i,j,3));
end,end

for j=2:jm-1
for i=2:im-1

rhs(i,j,1:3)=rhs(i,j,1:3)-(FV(i,j,1:3)-FV(i-1,j,1:3))...
-(GV(i,j,1:3)-GV(i,j-1,1:3));
end,end


MAIN CODE:


clear all
clc
close all

%%%% Navier-Stokes

Re=50;

% MALLA

N=200;
L=1;
r=1.1; % 1 <= r <= 1.3

d=L/2*(r-1)/(r^((N-1)/2)-1); %%Geometric progression
for i=2:round((N-1)/2+1)
x(i,1)=d*(r^(i-1)-1)/(r-1);
x(N+2-i,1)=L-x(i-1);
end
if r==1 x=linspace(0,L,N);x=x';end
y=x;

for i=2:N-1
xc(i)=0.5*(x(i+1)-x(i-1));
end
xc(1)=x(2)-x(1);
xc(N)=x(N)-x(N-1);
ye=xc;

for i=1:N
for j=1:N
jcb(i,j)=1/(xc(i)*ye(j)); %%Jacobian
csi(i,j)=1/xc(i);
eta(i,j)=1/ye(j);
end,end
[X,Y]=meshgrid(x,y);

% boundary conditions
q(1:N,1:N,3)=0;
q(1:N,1:N,2)=0;
q(1,1:N,2)=1;
q(1,1,1)=0; % refrence preassure
pfix=q(1,1,1);

% N-S

itmax=5;
CFL=0.9;

for it=1:itmax

qn=q;

%Runge-Kutta

for l=1:4
rhs=Advective(q,jcb,csi,eta,N,N);
rhs=Viscous(q,Re,jcb,csi,eta,N,N,rhs);

%avance del paso R-K
for i=2:N-1
for j=2:N-1

aux=dt*jcb(i,j)*rhs(i,j,1:3)/(5-l);

q(i,j,1:3)=qn(i,j,1:3)-aux;
q(i,j,1)=q(i,j,1)-pfix;

end,end

q(:,1,1)=(q(:,3,1)-q(:,2,1))/(X(:,3)-X(:,2))*(X(:,1)-X(:,2))+q(:,2,1);
q(:,end,1)=(q(:,end-1,1)-q(:,end-2,1))/...
(X(:,end-1)-X(:,end-2))*(X(:,end)-X(:,end-2))+q(:,end-2,1);
q(1,:,1)=(q(3,:,1)-q(2,:,1))/(Y(3,:)-Y(2,:))*(Y(1,:)-Y(2,:))+q(2,:,1);
q(end,:,1)=(q(end-1,:,1)-q(end-2,:,1))/...
(Y(end-1,:)-Y(end-2,:))*(Y(end,:)-Y(end-2,:))+q(end-2,:,1);
end

end

figure(1)
hold on
u=(flipud(q(:,:,2)));
v=(flipud(q(:,:,3)));
quiver(X,Y,u,v)
rectangle('Position',[0,0,L,L])

axis([-0.2 1.2 -0.2 1.2])


i have trouble defining my CFL time step also :(