دینامیک سیالات محاسباتی (CFD)

clc
clear all
dx=1;dy=1;
L=1;roh=1;%%%%%%%    
a=10;b=2;
nx=25;ny=25;%%%%%  number of nod in x and y direction
u=1;v=2;
FI=zeros(ny,nx);
fi0=0;fi1=100;
FI(1,:)=fi1;FI(:,1)=fi1;FI(ny,:)=fi0;FI(:,nx)=fi0;
%%%%%%%%%%%%%%%%%5   diffution term
ae_d=L*dy/dx;
aw_d=L*dy/dx;
an_d=L*dx/dy;
as_d=L*dx/dy;
%%%%%%%%%%%%%%%%%  convection term by Central
ae_c=-roh*u*dy/2;
aw_c=roh*u*dy/2;
an_c=-roh*v*dx/2;
as_c=roh*v*dx/2;
%%%%%%%%%%%%%%%%%%  ae=convection + diffution
ae=ae_d+ae_c;
as=as_c+as_d;
aw=aw_c+aw_d;
an=an_c+an_d;
%%%%%%%%%%%%%%%%
s=1;     %%%%%%%%%%%  parameter for account 
for j=1:ny
    for i=1:nx
        if (i==1||j==1)
            A(s,s)=1;
            B(s)=fi1;
        elseif(i==nx||j==ny)
            A(s,s)=1;
            B(s)=fi0;
        else
            ap=-(ae+as+aw+an)-b*dx*dy;
            A(s,[s-nx  s-1 s s+1 s+nx])=[as aw ap ae an];
            B(s)=-a*dx*dy;
        end
        s=s+1;
    end
end
sol=solution_matrix(A,B);
s=1;
for j=1:ny
    for i=1:nx
        FI(j,i)=sol(s);
        s=s+1;
    end
end
xlswrite('Central.xls',FI)
%%%%%   plot
x=dx:dx:nx*dx;
y=dy:dy:ny*dy;
[xx,yy]=meshgrid(x,y);
surf(xx,yy,FI) 
title(' Function \phi By Central')
xlabel('X')
ylabel('Y')
view(45,45)
colorbar