一、简介
二维对流扩散温度场
二、源代码
function A=secdim2
% 本程序可实现两种可选边界条件的对流扩散实时动态模拟
% 采用乘方格式
% 初始温度场为10,速度可选
% 时间差分采用全隐式格式
clc
dt=0.001;
bc=input('南北边界条件:1-第一类,2-第二类: ');
q=0;
ddt=input('迭代次数: ');
L1=0.02;r=1000;k=0.017;sp=0;L2=0.01;
u=input('横向速度: ');
if bc==2
v=0;
else
v=input('纵向速度 : ');
end
n=100;m=100;
dx=L1/n;dy=L2/m;F1=r*u*dy;F2=r*v*dx;
ap0=r*dx*dy/dt;
D1=k*dy/dx;pe1=F1/D1;D2=k*dx/dy;pe2=F2/D2;
tp0=10*ones(n*m,1);%初场温度
tw=10;te=10; %边界条件
A=sparse(n*m);
b=ap0*tp0;
%二维影响系数采用乘方格式
ae=D1*max(0,(1-0.1*abs(pe1^5)))+max(-F1,0);
aw=D1*max(0,(1-0.1*abs(pe1^5)))+max(0,F1);
as=D2*max(0,(1-0.1*abs(pe2^5)))+max(0,F2);
an=D2*max(0,(1-0.1*abs(pe2^5)))+max(0,-F2);
ap=aw+ae+as+an+ap0-sp;
spw=-1*2*k*dy/dx;
spe=-1*2*k*dy/dx;
%边界条件源项处理
if bc==2
q=input('热流为: ');
spn=0;sps=0;
tn=0;ts=0;te=0;spe=0;
else
tn=input('北侧温度: ');
ts=input('南侧温度: ');
spn=-1*2*k*dx/dy;
sps=-1*2*k*dx/dy;
end
sun=zeros(n,1);
sus=zeros(n,1);
sun(5:64)=2*k*tn*dx/dy+q*dx;
sus(5:64)=2*k*ts*dx/dy+q*dx;
suw=2*k*tw*dy/dx;
sue=2*k*te*dy/dx;
%系数矩阵,常数向量赋值
for i=1:n*m
%南s边界点赋值
if(i<n&&i>1)
A(i,i)=ae+an+aw+ap0-sps;
A(i,i+1)=-ae; A(i,i-1)=-aw;A(i,i+n)=-an;
b(i)=b(i)+sus(i);
%内节点inner赋值
else if(i>n&&rem(i,n)~=1&&rem(i,n)~=0&&i<n*(m-1)+1)
A(i,i)=ap;
A(i,i+n)=-an; A(i,i-n)=-as; A(i,i+1)=-ae; A(i,i-1)=-aw;
%北n边界点赋值
else if(i>n*m-n+1&&i<n*m)
A(i,i)=ae+aw+as+ap0-spn;
A(i,i+1)=-ae; A(i,i-1)=-aw; A(i,i-n)=-as;
b(i)=b(i)+sun(i-n*m+n);
%西w边界点赋值
else if(i>n&&i<n*m-n+1&&rem(i,n)==1)
A(i,i)=ae+an+as+ap0-spw;
A(i,i+n)=-an; A(i,i-n)=-as; A(i,i+1)=-ae;
b(i)=b(i)+suw;
%东e边界点赋值
else if(i>n&&i<n*m&&rem(i,n)==0)
A(i,i)=aw+an+as+ap0-spe;
A(i,i+n)=-an; A(i,i-n)=-as; A(i,i-1)=-aw;
b(i)=b(i)+sue;
else
end
end
end
end
end
end
%四个角点赋值
%西南角点
A(1,1)=1;
%东南角点
A(n,n)=1;
%西北角点
A(n*m-n+1,n*m-n+1)=1;
%东北角点
A(n*m,n*m)=1;
%
end
三、备注
完整代码或者代写添加QQ912100926
往期回顾>>>>>>
【物理应用】杨氏双缝干涉实验【Matlab 310期】
【物理应用】井筒多相流【Matlab 311期】