清华大学传热大作业-无限大平壁简单数值计算报告

传热大作业报告2011010***热动*****一、大作业题目一厚度为0.1m的无限大平壁，两侧均为对流换热边界条件，初始时两侧流体温度与壁内温度一致，tf1=tf2=t0=5℃；已知两侧对流换热系数分别为h1=11W/m2K、h2=23W/m2K,平壁材料的导热系数=0.43W/mK，导温系数a=0.3437×10-6m2/s。如果一侧的环境温度tf1突然升高为50℃并维持不变，计算在其它参数不变的条件下，平壁内温度分布及两侧壁面热流密度随时间的变化规律（用图形表示）。要求：将全部计算内容（包括网格的划分、节点方程组、计算框图、程序及计算结果）用A4纸打印。二、网格划分如图，将无限大平板作为一维处理，本题为一维非稳态导热问题，对流换热边界条件。空间网格划分：平板总厚度为delta=0.1m，定义空间步长为dx=0.005m，则距离份数为N=delta/dx=20份。定义x{n}为以0为首项，以dx为公差的等差数列，尾项为delta=0.1m，共有N+1项，则x{n}中的每一项即表示一个沿平板厚度方向中的划分点。时间网格划分：设总时间长度为T=100000s，定义时间步长为dtao=20s，则时间份数为M=5000份。定义tao{m}是以0为首项，以dtao为公差的等差数列，尾项为T=100000s，共有M+1项,则tao{m}中每一项即表示一个时刻。三、计算框图程序中的各个变量的名称及意义：1.题设中各个常数lambda=0.43导热系数；a=0.3437e-6热扩散率；h1=11边界对流换热系数；h2=23边界对流换热系数2；t0=5初始温度；tf1=50初始流体温度；tf2=5初始流体温度2；delta=0.1总距离长度（无限大平板厚度）；2.网格划分所设的变量T=100000总时间长度（在T时间内考虑本问题）；dtao=20定义时间步长；dx=0.005定义距离步长；M=floor(T/dtao)时间份数=总时间/时间步长（向下取整）；N=floor(delta/dx距离份数=总厚度/距离步长(向下取整）；tao=0:dtao:T定义时间划分单元（以0为首项，以dtao为公差的等差数列，尾项为T），共有M+1项；x=0:dx:delta定义距离划分单元（以0为首项，以dx为公差的等差数列，尾项为delta），共有N+1项；3.判定稳定性的准则数Bi1=h1*dx/lambda边界节点网格毕渥数；Bi2=h2*dx/lambda边界节点网格毕渥数2；Fo=a*dtao/dx^2傅里叶数；程序计算框图NoYesNoYes四、程序代码本程序在MATLABR2008a中运行通过，以下是源代码（%后为注释）：lambda=0.43;%导热系数a=0.3437e-6;%热扩散率h1=11;%边界对流换热系数h2=23;%边界对流换热系数2输入delta,T,dtao,dx,M,N,tao,xFo1/(2*Bi1+2)&&Fo1/(2*Bi2+2)开始建立t(M+1,N+1)温度矩阵，令t(1,:)=t0,令m=1m=m+1t(m,1)=2*Fo*(t(m-1,2)+Bi1*tf1)+(1-2*Bi1*Fo-2*Fo)*t(m-1,1)t(m,N+1)=2*Fo*(t(m-1,N)+Bi2*tf2)+(1-2*Bi2*Fo-2*Fo)*t(m-1,N+1)t(m,n)=Fo*(t(m-1,n-1)+t(m-1,n+1))+(1-2*Fo)*t(m-1,n)mM+1输出温度矩阵t(M+1,N+1)和相应图象停机打印“不稳定”t0=5;%初始温度tf1=50;%初始流体温度tf2=5;%初始流体温度2delta=0.1;%总距离长度（无限大平板厚度）T=100000;%总时间长度（在T时间内考虑本问题）dtao=20;%定义时间步长dx=0.005;%定义距离步长M=floor(T/dtao);%时间份数=总时间/时间步长（向下取整）N=floor(delta/dx);%距离份数=总厚度/距离步长(向下取整）tao=0:dtao:T;%定义时间划分单元（以0为首项，以dtao为公差的等差数列，尾项为T），共有M+1项x=0:dx:delta;%定义距离划分单元（以0为首项，以dx为公差的等差数列，尾项为delta），共有N+1项Bi1=h1*dx/lambda;%边界节点网格毕渥数Bi2=h2*dx/lambda;%边界节点网格毕渥数2Fo=a*dtao/dx^2;%傅里叶数ifFo1/(2*Bi1+2)&&Fo1/(2*Bi2+2)%判断稳定性，不稳定则显示毕渥数、傅里叶数disp('不稳定');disp(Bi1);disp(Bi2);disp(Fo);disp(1/(2*Bi1+2));disp(1/(2*Bi2+2));else%若稳定，则进行迭代计算t=zeros(M+1,N+1);%建立一个(M+1)*(N+1)的温度矩阵，M+1为时间节点个数，N+1为空间节点个数，以便进行迭代计算q1=zeros(M+1,1);%根据题目要求算两壁面处热流密度q2=zeros(M+1,1);t(1,:)=t0;%初始温度均为t0=5℃form=2:M+1%m=1时是初值上一行已计算出，则从m=2一直计算到m=M+1，m对应的时刻是tao=（m-1)dtaot(m,1)=2*Fo*(t(m-1,2)+Bi1*tf1)+(1-2*Bi1*Fo-2*Fo)*t(m-1,1);%首先计算一边界这个时刻温度t(m,N+1)=2*Fo*(t(m-1,N)+Bi2*tf2)+(1-2*Bi2*Fo-2*Fo)*t(m-1,N+1);%再计算另一边界这个时刻的温度q1(m)=h1*(tf1-t(m,1));q2(m)=h2*(t(m,N+1)-tf2);forn=2:N%然后计算内部，n=1和n=N+1时是边界节点温度，上面两行已经计算出，n对应的坐标是x=(n-1)*dxt(m,n)=Fo*(t(m-1,n-1)+t(m-1,n+1))+(1-2*Fo)*t(m-1,n);endend%以下是画图figureplot(x,t(1,:),x,t(11,:),x,t(21,:),x,t(51,:),x,t(101,:),x,t(1001,:),x,t(5001,:));legend('t=0s','t=200s','t=400s','t=1000s','t=2000s','t=20000s','t=100000',0);title('一定时间下温度随距离的分布','fontsize',12,'fontweight','bold','fontname','楷体');axis([0,0.1,0,40]);figureplot(tao,t(:,1),tao,t(:,6),tao,t(:,11),tao,t(:,16),tao,t(:,21));legend('x=0','x=0.025','x=0.05','x=0.075','0.1',0);title('一定位置处温度随时间的分布','fontsize',12,'fontweight','bold','fontname','楷体');axis([0,100000,0,40]);figuremesh(x,tao,t);title('温度随时间和空间的分布','fontsize',12,'fontweight','bold','fontname','楷体');figureplot(tao,q1,tao,q2);legend('q1','q2');title('两壁面热流密度随时间变化曲线','fontsize',12,'fontweight','bold','fontname','楷体');end五、计算结果及图表最终t（M+1,N+1）矩阵数据因为太庞大，详见“传热大作业数据.xls”。两壁面热流密度随时间变化曲线（横轴时间，纵轴热流密度）不同位置处温度随时间分布曲线(横轴时间，纵轴温度)不同时间下温度随距离分布曲线（横轴距离，纵轴温度）温度随距离、时间分布的三维网格图