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
 

注：完整代码或者代写添加QQ2449341593
往期回顾>>>>>>