二维瞬态导热时间显式格式MATLAB程序

最新推荐文章于 2024-04-14 13:36:19 发布

心同流水

最新推荐文章于 2024-04-14 13:36:19 发布

阅读量2.1k

点赞数 2

分类专栏：数值传热学文章标签： matlab

本文链接：https://blog.csdn.net/dashen_daifei/article/details/110229606

版权

数值传热学专栏收录该内容

7 篇文章 26 订阅

订阅专栏

二维瞬态导热时间显式格式MATLAB程序

代码改编自计算传热学–11-二维瞬态导热时间显式格式C程序-数值传热学
下面是代码

%
%%
%
%%
clc
clear
LengthX = 3;
LengthY = 2;
Tleft = 3;
Tright = 5;
Tbottom = 2;
Tup = 4;
den = 100;
s = 10;
k=10;
cv = 1000;

dt = 200;
N = 5;%X
M = 5;%Y
maxstep = 200;
%%
dx =LengthX/N;
dy =LengthY/M;
ae0 = zeros(M+2,N+2);
aw0 = zeros(M+2,N+2);
an0 = zeros(M+2,N+2);
as0 = zeros(M+2,N+2);
ap0 = zeros(M+2,N+2);
ap1 = zeros(M+2,N+2);
b= zeros(M+2,N+2);
T0 =zeros(M+2,N+2);
T1 = zeros(M+2,N+2);
Tall = zeros(N+2,maxstep+1);
%% 边界条件
T0(2:M+1,2:N+1) = 3;
T0(1:M+2,1) = Tleft;
T0(1:M+2,N+2) = Tright;
T0(1,1:N+2) = Tup;
T0(M+2,1:N+2) = Tbottom;

T1(1:M+2,1) = Tleft;
T1(1:M+2,N+2) = Tright;
T1(1,1:N+2) = Tup;
T1(M+2,1:N+2) = Tbottom;
%% 网格系数计算
%内部网格计算
ae0(3:M,3:N) = k*dy/dx;
aw0(3:M,3:N) = k*dy/dx;
an0(3:M,3:N) = k*dx/dy;
as0(3:M,3:N) = k*dx/dy;
ap0(3:M,3:N) = den*cv*dx*dy/dt-aw0(3:M,3:N)-ae0(3:M,3:N)-an0(3:M,3:N)-as0(3:M,3:N);
ap1(3:M,3:N) = ap0(3:M,3:N)+ae0(3:M,3:N)+aw0(3:M,3:N)+an0(3:M,3:N)+as0(3:M,3:N);
b(3:M,3:N) = s*dx*dy;
%边界网格计算
%第一个
i=2;
for j=3:N
ae0(i,j) = k*dy/dx;
aw0(i,j) = k*dy/dx;
as0(i,j) = k*dx/dy;
an0(i,j) = k*dx/(dy/2);
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;
end

i=M+1;
for j=3:N
ae0(i,j) = k*dy/dx;
aw0(i,j) = k*dy/dx;
as0(i,j) = k*dx/(dy/2);
an0(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;
end
j=2;
for i=3:M
ae0(i,j) = k*dy/dx;
aw0(i,j) = k*dy/(dx/2);
an0(i,j) = k*dx/dy;
as0(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;
end
j=N+1;
for i=3:M
ae0(i,j) = k*dy/(dx/2);
aw0(i,j) = k*dy/dx;
an0(i,j) = k*dx/dy;
as0(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;
end

j=2;
i=2;
ae0(i,j) = k*dy/dx;
aw0(i,j) = k*dy/(dx/2);
as0(i,j) = k*dx/dy;
an0(i,j) = k*dx/(dy/2);
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;

j=2;
i=M+1;
ae0(i,j) = k*dy/dx;
aw0(i,j) = k*dy/(dx/2);
as0(i,j) = k*dx/(dy/2);
an0(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;

i=2;
j=N+1;
ae0(i,j) = k*dy/(dx/2);
aw0(i,j) = k*dy/dx;
as0(i,j) = k*dx/dy;
an0(i,j) = k*dx/(dy/2);
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;

i=M+1;
j=N+1;
ae0(i,j) = k*dy/(dx/2);
aw0(i,j) = k*dy/dx;
as0(i,j) = k*dx/(dy/2);
an0(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt-aw0(i,j)-ae0(i,j)-an0(i,j)-as0(i,j);
ap1(i,j) = ap0(i,j)+ae0(i,j)+aw0(i,j)+an0(i,j)+as0(i,j);
b(i,j) = s*dx*dy;


Tall(1:M+2,1:N+2) = T0;
for h = 1:maxstep
T1(2:M+1,2:N+1) = (ae0(2:M+1,2:N+1).*T0(2:M+1,3:N+2)+aw0(2:M+1,2:N+1).*T0(2:M+1,1:N)...
    +an0(2:M+1,2:N+1).*T0(1:M,2:N+1)+as0(2:M+1,2:N+1).*T0(3:M+2,2:N+1)...
    +ap0(2:M+1,2:N+1).*T0(2:M+1,2:N+1)+b(2:M+1,2:N+1))./ap1(2:M+1,2:N+1);
T0=T1;
Tall((h)*(M+2)+1:(h+1)*(M+2),1:N+2) = T0;
end
T11=zeros(maxstep,1);
for i =1:maxstep+1
    T11(i) = Tall((i-1)*(M+2)+6,2);
end
plot(dt:dt:dt*maxstep,T11(2:maxstep+1));
hold on