一维非稳态导热CRANK.docx

一维非稳态导热CRANK-NICOLSON解法题目：数值计算一维非稳态导热，长度1米的不锈钢棒原来温度都是0度，一端温度突然变为300度，并保存不变，采用CRANK-NICOLSON 方法数值计算不锈钢内温度分布随时间的变化解法：一维导热微分方程IGU 「©% 小 n 气亍广沃，身四）边界条件为 u（0，t）=0； u（a0，t）=300初值 u（X，0）=0；主程序clc clearuX=1; %不锈钢长1米uT=2000; %时长 2000 秒M=10;%空间轴等分区间数N=1000; %时间轴等分区间数rou=8030; %不锈钢密度cp=502.48; %不锈钢热容kk=16.27; %不锈钢导热率D=kk/rou/cp;% 扩散系数phi=inline（'0'）; % 初值psi1=inline（'0'）; % 左边界psi2=inline（'300'）; % 右边界%计算步长dx=uX/M;%x 的步长dt=uT/N;%t 的步长x=(0:M)*dx;t=(0:N)*dt;r=D*dt/dx/dx;% 步长比Diag=zeros(1,M-1);%矩阵的对角线元素Low=zeros(1/M-2);%矩阵的下对角线元素Up=zeros(1,M-2);%矩阵的上对角线元素for i=1:M-2Diag(i)=1+r;Low(i)=-r/2;Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));endB=zeros(M-1,M-1);for i=1:M-2B(i,i)=1-r;B(i,i+1)=r/2;B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐层求解，需要使用追赶法(调用函数EqtsForwardAndBackward)for j=1:Nb1=zeros(M-1,1);b1(1)=r*(U(1,j+1)+U(1,j))/2;b1(M-1)=r*(U(M+1,j+1)+U(M+1,j))/2;b=B*U(2:M,j)+b1;U(2:M,j+1)=zhuiganfa(Low,Diag,Up,b);endU=U';%作出图形 mesh(x,t,U);xlabel('空间变量x')ylabel('时间变量t')shading interp程序用到了追赶法子程序，代码如下function x=zhuiganfa(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%检查参数的输入是否正确n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 ~= 1 || n-n2 ~= 1 || n ~= mdisp('输入参数有误！')x='';return;end%追的过程for i=2:nL(i-1)=L(i-1)/D(i-1);D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:nx(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);for i=n-1:-1:1x(i)=(x(i)-U(i)*x(i+1))/D(i);endreturn;运行主程序，最终得到如图所示结果。