Example The Simulation of a 1D diffusion case using Runge-Kutta for time stepping
clc
clear
% Constants Input
L=1;
N=32;
ALPHA=32*(10)^(-6);
DT=1;
DX=L/N;
umax=100;
% velocity Profile
p=0;
for pp=-N/2:N/2;
p=p+1;
X(p)=pp*DX;
u(p)=umax*(1-(X(p))^2);
end
% checking paramters
% plot(u)
% x'
% u'
pp=0;
j=1;
for p=-N/2:N/2;
pp=pp+1;
if (p==-N/2)
pp;
RHS(pp,j)=ALPHA*((-2*u(pp)+u(pp+1))/(DX)^2);
elseif ((p>-N/2) & (p<N/2) )
pp;
RHS(pp,j)=ALPHA*((-2*u(pp)+u(pp+1)+u(pp-1))/(DX)^2);
elseif (p==N/2)
pp;
RHS(pp,j)=ALPHA*((-2*u(pp)+u(pp-1))/(DX)^2);
end
end
% checking paramters
% plot(RHS)
% RHS
% p'
% pp'
pp=0
j=1
for p=-N/2:N/2;
pp=pp+1
U(pp,j)=DT*RHS(pp,j)+u(pp)
end
%plot(X,U)
%grid on
%Time Stepping Section
pp=0;
for j=1:10;
pp=0;
for p=-N/2:N/2;
pp=pp+1;
if (p==-N/2)
pp;
RHS(pp,j)=0;
elseif ((p>-N/2) & (p<N/2) )
pp;
RHS(pp,j)=ALPHA*((-2*U(pp,j)+U(pp+1)+U(pp-1,j))/(DX)^2);
elseif (p==N/2)
pp;
RHS(pp,j)=0;
end
end
pp=0;
for p=-N/2:N/2;
pp=pp+1;
k1(pp,j)=DT*RHS(pp,j);
k2(pp,j)=DT*(0.5*k1(pp,j)+DT*RHS(pp,j));
k3(pp,j)=DT*(0.5*k2(pp,j)+DT*RHS(pp,j));
k4(pp,j)=DT*(k3(pp,j)+DT*RHS(pp,j));
U(pp,j+1)=U(pp,j)+(1/6)*k1(pp,j)+(1/3)*(k2(pp,j)+k3(pp,j))+(1/6)*k4(pp,j)
end
plot(X,U)
grid on
hold on
end
clear
% Constants Input
L=1;
N=32;
ALPHA=32*(10)^(-6);
DT=1;
DX=L/N;
umax=100;
% velocity Profile
p=0;
for pp=-N/2:N/2;
p=p+1;
X(p)=pp*DX;
u(p)=umax*(1-(X(p))^2);
end
% checking paramters
% plot(u)
% x'
% u'
pp=0;
j=1;
for p=-N/2:N/2;
pp=pp+1;
if (p==-N/2)
pp;
RHS(pp,j)=ALPHA*((-2*u(pp)+u(pp+1))/(DX)^2);
elseif ((p>-N/2) & (p<N/2) )
pp;
RHS(pp,j)=ALPHA*((-2*u(pp)+u(pp+1)+u(pp-1))/(DX)^2);
elseif (p==N/2)
pp;
RHS(pp,j)=ALPHA*((-2*u(pp)+u(pp-1))/(DX)^2);
end
end
% checking paramters
% plot(RHS)
% RHS
% p'
% pp'
pp=0
j=1
for p=-N/2:N/2;
pp=pp+1
U(pp,j)=DT*RHS(pp,j)+u(pp)
end
%plot(X,U)
%grid on
%Time Stepping Section
pp=0;
for j=1:10;
pp=0;
for p=-N/2:N/2;
pp=pp+1;
if (p==-N/2)
pp;
RHS(pp,j)=0;
elseif ((p>-N/2) & (p<N/2) )
pp;
RHS(pp,j)=ALPHA*((-2*U(pp,j)+U(pp+1)+U(pp-1,j))/(DX)^2);
elseif (p==N/2)
pp;
RHS(pp,j)=0;
end
end
pp=0;
for p=-N/2:N/2;
pp=pp+1;
k1(pp,j)=DT*RHS(pp,j);
k2(pp,j)=DT*(0.5*k1(pp,j)+DT*RHS(pp,j));
k3(pp,j)=DT*(0.5*k2(pp,j)+DT*RHS(pp,j));
k4(pp,j)=DT*(k3(pp,j)+DT*RHS(pp,j));
U(pp,j+1)=U(pp,j)+(1/6)*k1(pp,j)+(1/3)*(k2(pp,j)+k3(pp,j))+(1/6)*k4(pp,j)
end
plot(X,U)
grid on
hold on
end
Example The Simulation of a 2D diffusion case using the Crank Nicolson Method for time stepping and TDMA Solver
clc
clear
DX=0.1;
DY=0.1;
DZ=0.01;
QW(1:4)=500000;
AN=DY*DZ;
AS=DY*DZ;
AW=DX*DZ;
AE=DX*DZ;
K=1000;
TN=100;
I=1;
AAN(I)=(K*AN)/DY;
AAS(I)=0;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(1,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=2;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(2,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=3;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(3,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=4;
AAN(I)=0;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=((2*K*AN*TN)/DY);
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=(-2*K*AN)/DY;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(4,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=5;
AAN(I)=(K*AN)/DY;
AAS(I)=0;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(5,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=6;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DX;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(6,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=7;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DX;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(7,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=8;
AAN(I)=0;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=((2*K*AN*TN)/DY);
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=(-2*K*AN)/DY;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(8,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=9;
AAN(I)=(K*AN)/DY;
AAS(I)=0;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(9,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=10;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(10,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=11;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(11,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=12;
AAN(I)=0;
AAS(I)=(K*AS)/DY;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=((2*K*AN*TN)/DY);
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=(-2*K*AN)/DY;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(12,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
for I=1:4;
ALPHA(I)=AAN(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==1)
AJ(I)=ALPHA(I)/DJ(I);
else if (I>1)
AJ(I)=ALPHA(I)/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=1:4;
CJ(I)=SU(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==1)
CJP(I)=CJ(I)/DJ(I);
else if (I>1)
CJP(I)=(BETA(I)*CJP(I-1)+CJ(I))/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=1:4;
VECTOR2(I,1:6)=[ BETA(I) DJ(I) ALPHA(I) CJ(I) AJ(I) CJP(I) ];
end
I=5;
for P=1:4;
I=I-1;
if (I==4)
T(I)=CJP(I);
else if (I<4)
T(I)=AJ(I)*T(I+1)+CJP(I);
end
end
end
for I=5:8;
ALPHA(I)=AAN(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==5)
AJ(I)=ALPHA(I)/DJ(I);
else if (I>5)
AJ(I)=ALPHA(I)/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
P=0;
for I=5:8;
P=P+1;
CJ(I)=AAW(I)*T(P)+SU(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==5)
CJP(I)=CJ(I)/DJ(I);
else if (I>5)
CJP(I)=(BETA(I)*CJP(I-1)+CJ(I))/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=5:8;
VECTOR2(I,1:4)=[ BETA(I) DJ(I) ALPHA(I) CJ(I) ];
end
I=9;
for P=1:4;
I=I-1;
if (I==8)
T(I)=CJP(I);
else if (I<8)
T(I)=AJ(I)*T(I+1)+CJP(I);
end
end
end
for I=9:12;
ALPHA(I)=AAN(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==9)
AJ(I)=ALPHA(I)/DJ(I);
else if (I>9)
AJ(I)=ALPHA(I)/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
P=4;
for I=9:12;
P=P+1;
CJ(I)=AAW(I)*T(P)+SU(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==9)
CJP(I)=CJ(I)/DJ(I);
else if (I>9)
CJP(I)=(BETA(I)*CJP(I-1)+CJ(I))/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=9:12;
VECTOR2(I,1:4)=[ BETA(I) DJ(I) ALPHA(I) CJ(I) ];
end
I=13;
for P=1:4;
I=I-1;
if (I==12)
T(I)=CJP(I);
else if (I<12)
T(I)=AJ(I)*T(I+1)+CJP(I);
end
end
end
%OUTPUT
T'
clear
DX=0.1;
DY=0.1;
DZ=0.01;
QW(1:4)=500000;
AN=DY*DZ;
AS=DY*DZ;
AW=DX*DZ;
AE=DX*DZ;
K=1000;
TN=100;
I=1;
AAN(I)=(K*AN)/DY;
AAS(I)=0;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(1,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=2;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(2,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=3;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(3,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=4;
AAN(I)=0;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=0;
BN(I)=((2*K*AN*TN)/DY);
BS(I)=0;
BE(I)=0;
BW(I)=QW(I)*AW;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=(-2*K*AN)/DY;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(4,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=5;
AAN(I)=(K*AN)/DY;
AAS(I)=0;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(5,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=6;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DX;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(6,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=7;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DX;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(7,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=8;
AAN(I)=0;
AAS(I)=(K*AS)/DY;
AAE(I)=(K*AE)/DX;
AAW(I)=(K*AW)/DX;
BN(I)=((2*K*AN*TN)/DY);
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=(-2*K*AN)/DY;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(8,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=9;
AAN(I)=(K*AN)/DY;
AAS(I)=0;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(9,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=10;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(10,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=11;
AAN(I)=(K*AN)/DY;
AAS(I)=(K*AS)/DY;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=0;
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=0;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(11,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
I=12;
AAN(I)=0;
AAS(I)=(K*AS)/DY;
AAE(I)=0;
AAW(I)=(K*AW)/DX;
BN(I)=((2*K*AN*TN)/DY);
BS(I)=0;
BE(I)=0;
BW(I)=0;
SU(I)=BN(I)+BS(I)+BE(I)+BW(I);
SPN(I)=(-2*K*AN)/DY;
SPS(I)=0;
SPE(I)=0;
SPW(I)=0;
SP(I)=SPN(I)+SPS(I)+SPE(I)+SPW(I);
AAP(I)=AAW(I)+AAE(I)+AAN(I)+AAS(I)-SP(I);
VECTOR(12,1:6)=[ AAN(I) AAS(I) AAW(I) AAE(I) AAP(I) SU(I)];
for I=1:4;
ALPHA(I)=AAN(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==1)
AJ(I)=ALPHA(I)/DJ(I);
else if (I>1)
AJ(I)=ALPHA(I)/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=1:4;
CJ(I)=SU(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==1)
CJP(I)=CJ(I)/DJ(I);
else if (I>1)
CJP(I)=(BETA(I)*CJP(I-1)+CJ(I))/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=1:4;
VECTOR2(I,1:6)=[ BETA(I) DJ(I) ALPHA(I) CJ(I) AJ(I) CJP(I) ];
end
I=5;
for P=1:4;
I=I-1;
if (I==4)
T(I)=CJP(I);
else if (I<4)
T(I)=AJ(I)*T(I+1)+CJP(I);
end
end
end
for I=5:8;
ALPHA(I)=AAN(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==5)
AJ(I)=ALPHA(I)/DJ(I);
else if (I>5)
AJ(I)=ALPHA(I)/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
P=0;
for I=5:8;
P=P+1;
CJ(I)=AAW(I)*T(P)+SU(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==5)
CJP(I)=CJ(I)/DJ(I);
else if (I>5)
CJP(I)=(BETA(I)*CJP(I-1)+CJ(I))/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=5:8;
VECTOR2(I,1:4)=[ BETA(I) DJ(I) ALPHA(I) CJ(I) ];
end
I=9;
for P=1:4;
I=I-1;
if (I==8)
T(I)=CJP(I);
else if (I<8)
T(I)=AJ(I)*T(I+1)+CJP(I);
end
end
end
for I=9:12;
ALPHA(I)=AAN(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==9)
AJ(I)=ALPHA(I)/DJ(I);
else if (I>9)
AJ(I)=ALPHA(I)/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
P=4;
for I=9:12;
P=P+1;
CJ(I)=AAW(I)*T(P)+SU(I);
BETA(I)=AAS(I);
DJ(I)=AAP(I);
if (I==9)
CJP(I)=CJ(I)/DJ(I);
else if (I>9)
CJP(I)=(BETA(I)*CJP(I-1)+CJ(I))/(DJ(I)-BETA(I)*AJ(I-1));
end
end
end
for I=9:12;
VECTOR2(I,1:4)=[ BETA(I) DJ(I) ALPHA(I) CJ(I) ];
end
I=13;
for P=1:4;
I=I-1;
if (I==12)
T(I)=CJP(I);
else if (I<12)
T(I)=AJ(I)*T(I+1)+CJP(I);
end
end
end
%OUTPUT
T'
Unless otherwise noted, all content on this site is @Copyright by Ahmed Al Makky 2012-2013 - http://cfd2012.com