计算传热学课程设计报告正六面体的非稳态导热模拟

1、计算传热学课程设计报告论文题目：正六面体的非稳态导热模拟专业写全称，仿宋GB-2312，三号字体。专业写全称，仿宋GB-2312，三号字体。 2014年 7 月 5 日目 录一、设计题目3二、控制方程4三、数学模型5四、数值方法及程序流程：54.1隐式格式下离散过程54.2显式格式下离散过程64.3程序流程：8五、集总参数法的应用8六、导热方程的无量纲化10七、结果的分析与处理117.1 显式格式处理结果117.2 隐式格式处理结果15参考文献18附录A：组员分工19附录B：计算环境及源程序19 正六面体的非稳态导热模拟摘要 本报告对正六面体的的非稳态导热问题进行了模拟，运用了数值模拟的方法在

2、三位直角坐标系中对模型的导热方程进行了离散, 根据离散方程进行了编程，最终模拟出了温度随时间变化的关系，得到了相关的图像。关键词：三维；数值模拟；非稳态导热 Regular hexahedron unsteady heatconduction simulationAbstractThis report for regular hexahedron of unsteady heat conduction problem is simulated, using the method of numerical simulation in the three rectangular coordinat

3、e system to the discrete model of heat conduction equation based on the discrete equations programming, finally simulate the relationship between temperature change over time, obtained the relevant images.Keywords：Three dimension；Numerical simulation；Unsteady heat conduction一、 设计题目正六面体的非稳态导热模拟。如图所示为

4、一个正六面体，该正六面体由不锈钢材料制成，三个边的边长分别为a=20mm，b=10mm，c=30mm，初始时刻温度为0，突然将该六面体置于温度为40的热空气中，外表面的对流换热系数为20W/（m²·K），不锈钢的物性参数为16W/（m·K），=7200kg/m³，c=460J/（kg·K）。请完成以下题目：建立该问题的数学方程，并给出定解条件；采用集总参数法对该问题进行分析，获得温度随时间变化的函数关系，并计算该物体的时间常数和达到稳态过程所需要的时间；选择合适的参数将该问题无量纲化，得到无量纲导热方程及边界条件；分别用隐式和显式格式对该问题进

5、行数学模拟，得到三维温度随时间的变化关系；将数值结果与集总参数法结果进行对比，看看二者差距如何？二、控制方程由导热微分方程可得控制方程为：（1）物理特性：直角坐标系，三维，非稳态，无内热源，导热（2）材料特性：常物性（3）边界条件：（六个面两两对称）=0， =0, =0, 三、数学模型四、 数值方法及程序流程：4.1隐式格式下离散过程用内节点法将所求的区域在均分网格的情况下离散化（1） 非稳态项积分：（2扩散项积分：整理可得=+ + 对隐式格式进行整理得：其中：4.2显式格式下离散过程=+ (1)非稳态项积分：(2)扩散项积分：整理得:=+ 对显式格式进行整理得：其中：4.3程序流程：五、集总

6、参数法的应用（1）根据能量平衡，在时间内物体内能的增量等于物体与环境的对流换热量，即：式中，为在时间段内物体的内能的增量，式中，V为物体体积，和c分别是物体的密度和比热，而T是物体在任意时刻时的温度。（2） 物体在时间内的对流换热量为，其中， ，A=整理得：引入过余温度：初始条件为：故，温度随时间变化的关系式：（3）时间常数：， 带入数据，（3） 取稳态温度T=39，则1666s六、导热方程的无量纲化（1）控制方程：（2）引入下述无量纲量特征长度：无量纲时间t：无量纲温度：， ， 则有T=，（3）将上述无量纲数带入控制方程中，经整理化简得：（4）无量纲化的边界条件：由无量纲温度：，故，七、结果

7、的分析与处理7.1 显式格式处理结果取中间某一截面绘图4s时温度分布8s时温度分布12s时温度分布16s时温度分布20s时温度分布24s时温度分布256秒时的温度分布7.2 隐式格式处理结果选取i=5截面进行绘图4s时的温度场分布8s时的温度场分布12s时的温度场分布16s时的温度场分布20s时的温度场分布24s时的温度场分布148s时的温度场分布由上面不同时刻的截图可以得出如下结论：（1）在该模型的导热过程中热量由四周向中心传递（2）沿同一坐标方向上，同一截面上不同温度的分布区域呈椭圆形（3）随着时间的增加，截面上低温区域的面积越来越小（4）中心温度越来越高，最终各不同温度区域的温差会越来越

8、小（5）整个区域的温度也会越来越接近环境温度。(6)显示格式和隐式格式下的温度变化情况大致相同。(7)随着时间的推移，温度升高的速度趋于平缓。(8)整个矩形截面的对称轴也是温度的对称轴。将数值结果与集总参数法的结果进行对比，集总参数法的时间大于数值方法的时间，造成这种结果的原因是，由于集总参数法应用的前提是0.1，而题目中取，所以误差较大。参考文献1 黄善波 刘中良编 计算传热学基础.青岛：中国石油大学（华东）出版社，2009. 2杨世铭  陶文铨编著  传热学（第四版） 高等教育出版社，2006 3陶文

9、铨编著.数值传热学（第二版）.西安：西安交通大学出版社，2001. 4 袁恩熙主编工程流体力学（第一版）.石油工业出版社，19865 魏东平 朱连章 于广斌编 C程序设计语言.北京：电子工业出版社,2009附录A：组员分工组 长：杨伟茂 编写程序 方程离散组 员：蔡佳鑫 设计报告 方程离散组 员：李培强 制作ppt 建立模型附录B：计算环境及源程序计算环境：Microsoft Visual C+ 6.0隐式格式主程序：#include<stdio.h>#include<math.h># define NX 15# define NY 15# define NZ

10、 15double TNXNYNZ,bNXNYNZ,T0NXNYNZ,TTNXNYNZ;double aP,aE,aW,aS,aN,aF,aB;double dx,dy,dz,dt,XL,YL,ZL;int N,M,W,p,c,lamd,Tw,Tf,h;double Tao;void main()FILE *fp;int i,j,k;XL=0.02;/三边距离YL=0.03; ZL=0.01;N=12;/节点数M=12;W=12;dt=0.5;/时间步长Tw=273;/物体开始温度p=7200;/物体密度c=460;Tf=313;/环境温度h=20;/对流传热系数lamd=16;/导热系数Ta

11、o=0;/时间 for (i=1;i<N-1;i+)for (j=1;j<M-1;j+)for(k=1;k<W-1;k+)T0ijk=Tw;Tijk=Tw;/步长 dx=1.0*XL/(N-1); dy=1.0*YL/(M-1);dz=1.0*ZL/(W-1);/赋初始值for (j=0;j<M;j+)for(k=0;k<W;k+) T00jk=Tf;T0jk=Tf; T0N-1jk=Tf;TN-1jk=Tf;for (i=0;i<N;i+)for(k=0;k<W;k+) T0i0k=Tf;Ti0k=Tf;T0iM-1k=Tf;TiM-1k=Tf;fo

12、r(i=0;i<N;i+)for (j=0;j<M;j+) T0ij0=Tf;Tij0=Tf;T0ijW-1=Tf;TijW-1=Tf;fp=fopen("result.txt","w");do aE=aW=lamd*dz*dy/dx;aN=aS=lamd*dx*dz/dy;aF=aB=lamd*dx*dy/dz;aP=p*c*dx*dy*dz/dt+2*aE+2*aN+2*aF;/结果loop:for(i=1;i<N-1;i+)for(j=1;j<M-1;j+)for(k=1;k<W-1;k+) bijk=p*c*dx*d

13、y*dz*T0ijk/dt;Tijk=aE*Ti+1jk+aW*Ti-1jk+aN*Tij+1k+aS*Tij-1k+aF*Tijk+1+aB*Tijk-1;Tijk=Tijk+bijk;Tijk=Tijk/aP;for(i=0;i<N;i+)for(j=0;j<W;j+) for(k=0;k<M;k+) T0ijk=Tijk;Tao=Tao+dt;if(fmod(Tao,4)!=0) goto loop;/输出 if(fmod(Tao,4)=0) fprintf(fp,"%lf时刻的温度场为n",Tao); for(i=0;i<N;i+)for(j

14、=0;j<W;j+)for(k=0;k<M;k+)fprintf(fp,"T%2d%2d%2d=%lf %lf %lf %lfn",i,j,k,i*dx,j*dy,k*dz,Tijk);for(i=5;i<6;i+)for(j=0;j<W;j+)for(k=0;k<M;k+)fprintf(fp,"%lf %lf %lfn",j*dy,k*dz,Tijk); while(Tao<30);显式格式主程序：#include<stdio.h>#include<math.h># define NX 15

15、# define NY 15# define NZ 15double TNXNYNZ,bNXNYNZ,T0NXNYNZ,TTNXNYNZ;double aP,aE,aW,aS,aN,aF,aB;double dx,dy,dz,dt,XL,YL,ZL;int N,M,W,p,c,lamd,Tw,Tf,h;double Tao;void main()FILE *fp;int i,j,k;XL=0.02;/三边距离YL=0.03; ZL=0.01;N=12;/节点数M=12;W=12;dt=0.5;/时间步长Tw=273;/物体开始温度p=7200;/物体密度c=460;Tf=313;/环境温度h=

16、20;/对流传热系数lamd=16;/导热系数Tao=0;/时间 for (i=1;i<N-1;i+)for (j=1;j<M-1;j+)for(k=1;k<W-1;k+)T0ijk=Tw;Tijk=Tw;/步长 dx=1.0*XL/(N-1); dy=1.0*YL/(M-1);dz=1.0*ZL/(W-1);/赋初始值for (j=0;j<M;j+)for(k=0;k<W;k+) T00jk=Tf;T0jk=Tf; T0N-1jk=Tf;TN-1jk=Tf;for (i=0;i<N;i+)for(k=0;k<W;k+) T0i0k=Tf;Ti0k=T

17、f;T0iM-1k=Tf;TiM-1k=Tf;for(i=0;i<N;i+)for (j=0;j<M;j+) T0ij0=Tf;Tij0=Tf;T0ijW-1=Tf;TijW-1=Tf;fp=fopen("result.txt","w");do aE=aW=lamd*dz*dy/dx;aN=aS=lamd*dx*dz/dy;aF=aB=lamd*dx*dy/dz;aP=p*c*dx*dy*dz/dt+2*aE+2*aN+2*aF;/结果loop:for(i=1;i<N-1;i+)for(j=1;j<M-1;j+)for(k=1;k<W-1;k+) bijk=p*c*dx*dy*dz*T0ijk/dt;Tijk=aE*Ti+1jk+aW*Ti-1jk+aN*Tij+1k+aS*Tij-1k+aF*Tijk+1+aB*Tijk-1;Tijk=Tijk+bijk;Tijk=Tijk/aP;for(i=0;i<N;i+)for(j=0;j<W;j+) for(k=0;k<M;k+) T0ijk=Tijk;Tao=Tao+dt;if(fmod(Tao,4)!=0) goto loop;/输出 if(fmod(Tao,4)=0) fprintf(fp,&quo