二维数值传热学C++程序代码

1、#include<stdio.h>#include<iostream>#include<math.h>#include<fstream>/*/从此开始至下一个/*/标记带处之间的代码可以写在头文件里。/计算参数宏定义开始/#defineGRID_ROW82#defineGRID_COLUMN102网格横纵数目，边界处多出一列处理#defineDX0.05#defineDY0.05网格单元尺寸，m#defineDEN_COPPER8900#defineDEN_STEEL7850#defineDEN_CONCRETE2430材料的密度，kg/(mA3

2、)#defineC_COPPER390#defineC_STEEL460#defineC_CONCRETE970材料比热容，J/(kg*K)#defineCOND_COPPER400#defineCOND_STEEL20#defineCOND_CONCRETE1.63/材料导热系正W/(m*K)W/(mA2*K)#defineCONV5/空气与各种材料的对流换热系数，#defineINITL_TMPRT298初始温度，K#defineT_C2K(x)(x)+273)/CtoK#defineT_K2c(x)(x)-273)/KtoC/摄氏温度与绝对温度相互转换#defineT_AIR298;空气

3、温度,K#definet_END3600计算时k,s#definet_STEP1时间步k,s#defineCOND_COP_CON3.2467#defineCOND_COP_STL38.0952#defineCONDCONSTL3.0143精心整理/不同材料的两个相邻结点，导热系数取调和平均数，W/(m*K)/计算参数宏定义结束/enumPropertycopper,steel,concrete,iso834,adiabat,air;结点属性的枚举类型声明structGriddoublemNodeTGRID_ROWGRID_COLUMN缩点温度PropertymNodePropertyGRID

4、_ROWGRID_COLUMN站点材料属性doublemNodeCGRID_ROWGRID_COLUMN缁点比热容doublemNodeDenGRID_ROWGRID_COLUMN结点密度;结构体Grid定义，此为计算的主要承载数据结构，将在计算中存储温度数据函数全局声明开始voidSetProperty(Grid&);/属性设置函数voidInitlTMPRT(Grid&);/温度初始化函数voidSetNodePara(Grid&);/结点物性参数设置函数voidCalculateTMPRT(Grid&);/温度计算函数，几乎全部数值计算的承担者voidwr

5、ite(int,Grid&);/写文件函数，每隔100秒输出温度分布文件到当前目录/函数全局说明结束voidSetProperty(Grid&A)for(intj=0;j<GRID_COLUMN;+j)A.mNodeProperty0j=air;第一行为空气边界，计算区域外for(intj=0;j<GRID_COLUMN;+j)A.mNodePropertyGRID_ROW-1j=adiabat;最后一行为绝热边界，计K域外for(inti=0;i<GRID_ROW;+i)if(i>=0&&i<=40)|(i>=61&

6、&i<GRID_ROW)A.mNodePropertyi0=air;elseA.mNodePropertyi0=iso834;第一列的空气边界和ISO834边界，计算区域外for(inti=0;i<GRID_ROW;+i)A.mNodePropertyiGRID_COLUMN-1=air;最后一列，空气边界，计算由域外for(inti=1;i<GRID_ROW-1;+i)for(intj=1;j<GRID_COLUMN-1;+j)if(i>=1&&i<=44&&j>=1&&j<=20)|(i

7、>=61&&i<GRID_ROW-2&&j>=1&&j<=(100-(80-i)*(49/19)A.mNodeProperty皿=concrete;混凝土区域elseif(i-20)*(i-20)+(j-25)*(j-25)<=20*20)A.mNodeProperty皿=air;/圆孔空气区域elseif(i>=41&&i<=60&&j>=1&&j<=90&&(i-20)*(i-20)+(j-25)*(j-25)>=20*

8、20)A.mNodePropertyij=copper;金属铜区域的一部分elseif(i>=1&&i<+40&&j>=21&&j<=90)A.mNodeProperty皿=copper;金属铜区域的另一部分elseA.mNodeProperty皿=steel;剩下的即为金属铁的区域属性设置函数，将Grid结构体A的所有网格赋予材料属性(枚举类型Property)voidInitlTMPRT(Grid&A)for(inti=0;i<GRID_ROW;+i)for(intj=0;j<GRID_COLUMN

9、;+j)A.mNodeTij=INITL_TMPRT;一温度初始化函数，将Grid结构体A的全部结点温度设置为初始温度INITL_TMPRTvoidCalculateTMPRT(Grid&A)doubleQE=0,QW=0,QN=0,QS=0;网个方向传来的热量for(inti=1;i<GRID_ROW-1;+i)for(intj=1;j<GRID_COLUMN-1;+j)switch(A.mNodePropertyij)switch(A.mNodePropertyij+1)判断右侧属性，并计算QEcasecopper:QE=(A.mNodeTij+1-A.mNodeTij

10、)*COND_COPPER;break;caseair:QE=(A.mNodeTij+1-A.mNodeTij)*CONV*DY;break;casesteel:QE=(A.mNodeTij+1-A.mNodeTij)*COND_COP_STL;break;caseconcrete:QE=(A.mNodeTij+1-A.mNodeTij)*COND_COP_CON;break;caseiso834:QE=2*(A.mNodeTij+1-A.mNodeTij)*COND_COPPER;break;caseadiabat:QE=0;break;default:std:cout<<i&l

11、t;<","<<j+1<<结点属性未指定n”;getchar();/std:cout<<QE<<"n"switch(A.mNodePropertyij-1)/判断左侧属性，并计算QWcasecopper:QW=(A.mNodeTij-1-A.mNodeTij)*COND_COPPER;break;caseair:QW=(A.mNodeTij-1-A.mNodeTij)*CONV*DY;break;casesteel:QW=(A.mNodeTij-1-A.mNodeTij)*COND_COP_STL;b

12、reak;.-caseconcrete:QW=(A.mNodeTij-1-A.mNodeTij)*COND_COP_CON;break;caseiso834:QW=2*(A.mNodeTij-1-A.mNodeTij)*COND_COPPER;break;caseadiabat:QW=0;break;default:std:cout<<i<<","<<j-1<<结点属性未指定n"getchar();switch(A.mNodePropertyi-1j)/判断上侧属性，并计算QNcasecopper:QN=(A.mNo

13、deTi-1j-A.mNodeTij)*COND_COPPER;break;caseair:QN=(A.mNodeTi-1j-A.mNodeTij)*CONV*DX;break;casesteel:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_COP_STL;break;caseconcrete:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_COP_CON;break;caseiso834:QN=2*(A.mNodeTi-1j-A.mNodeTij)*COND_COPPER;break;caseadiabat:QN=0;break;default:s

14、td:cout<<i-1<<","<<j<<结点属性未指定n"getchar();switch(A.mNodePropertyi+1j)判断下侧属性，并计算QScasecopper:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_COPPER;break;caseair:QS=(A.mNodeTi+1j-A.mNodeTij)*CONV*DX;break;casesteel:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_COP_STL;Ibreak;caseconcrete:

15、QS=(A.mNodeTi+1j-A.mNodeTij)*COND_COP_CON;break;Icaseiso834:QS=2*(A.mNodeTi+1j-A.mNodeTij)*COND_COPPER;break;caseadiabat:QS=0;break;default:std:cout<<i+1<<","<<j<<结点属性未指定n”;getchar();break;switch(A.mNodePropertyij+1)/判断右侧属性，并计算QEcasecopper:QE=(A.mNodeTij+1-A.mNodeTi

16、j)*COND_COP_STL;break;caseair:QE=(A.mNodeTij+1-A.mNodeTij)*CONV*DY;break;casesteel:QE=(A.mNodeTij+1-A.mNodeTij)*COND_STEEL;break;caseconcrete:QE=(A.mNodeTij+1-A.mNodeTij)*COND_CON_STL;break;caseiso834:QE=2*(A.mNodeTij+1-A.mNodeTij)*COND_STEEL;break;caseadiabat:QE=0;break;default:std:cout<<i<

17、;<","<<j+1<<结点属性未指定n”;getchar();switch(A.mNodePropertyij-1)/判断左侧属性，并计算QWcasecopper:QW=(A.mNodeTij-1-A.mNodeTij)*COND_COP_STL;break;caseair:QW=(A.mNodeTij-1-A.mNodeTij)*CONV*DY;break;casesteel:QW=(A.mNodeTij-1-A.mNodeTij)*COND_STEEL;break;caseconcrete:QW=(A.mNodeTij-1-A.mNode

18、Tij)*COND_CON_STL;break;caseiso834:QW=2*(A.mNodeTij-1-A.mNodeTij)*COND_STEEL;break;caseadiabat:QW=0;break;default:std:cout<<i<<","<<j-1<<结点属性未指定n"getchar();:二，二:Iswitch(A.mNodePropertyi-1j)/判断上侧属性，并计算QNcasecopper:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_COP_STL;break;

19、caseair:QN=(A.mNodeTi-1j-A.mNodeTij)*CONV*DX;break;casesteel:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_STEEL;break;caseconcrete:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_CON_STL;break;caseiso834:QN=2*(A.mNodeTi-1j-A.mNodeTij)*COND_STEEL;break;caseadiabat:QN=0;break;default:std:cout<<i-1<<","&l

20、t;<j<<结点属性未指定n"getchar();switch(A.mNodePropertyi+1j)/判断下侧属性，并计算QScasecopper:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_COP_STL;break;caseair:QS=(A.mNodeTi+1j-A.mNodeTij)*CONV*DX;zbreak;.°.Icasesteel:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_STEEL;break;caseconcrete:QS=(A.mNodeTi+1j-A.mNodeTij)*CON

21、D_CON_STL;break;caseiso834:QS=2*(A.mNodeTi+1j-A.mNodeTij)*COND_STEEL;break;caseadiabat:QS=0;break;default:std:cout<<i+1<<","<<j<<结点属性未指定n”;getchar();break;switch(A.mNodePropertyij+1)/判断右侧属性，并计算QEcasecopper:QE=(A.mNodeTij+1-A.mNodeTij)*COND_COP_CON;break;caseair:QE=(

22、A.mNodeTij+1-A.mNodeTij)*CONV*DY;break;casesteel:QE=(A.mNodeTij+1-A.mNodeTij)*COND_CON_STL;break;caseconcrete:QE=(A.mNodeTij+1-A.mNodeTij)*COND_CONCRETE;break;caseiso834:break;caseadiabat:QE=0;break;default:std:cout<<i<<","<<j+1<<结点属性未指定n”;getchar();switch(A.mNodePr

23、opertyij-1)/判断左侧属性，并计算QWcasecopper:break;caseair:QW=(A.mNodeTij-1-A.mNodeTij)*CONV*DY;break;.°.Icasesteel:break;caseconcrete:QW=(A.mNodeTij-1-A.mNodeTij)*COND_CONCRETE;break;caseiso834:break;caseadiabat:QW=0;break;default:std:cout<<i<<","<<j-1<<结点属性未指定n"ge

24、tchar();switch(A.mNodePropertyi-1j)/判断上侧属性，并计算QNcasecopper:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_COP_CON;break;caseair:QN=(A.mNodeTi-1j-A.mNodeTij)*CONV*DX;break;casesteel:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_CON_STL;break;caseconcrete:QN=(A.mNodeTi-1j-A.mNodeTij)*COND_CONCRETE;break;caseiso834:break;casea

25、diabat:QN=0;break;default:std:cout<<i-1<<","<<j<<结点属性未指定n"getchar();switch(A.mNodePropertyi+1j)判断下侧属性，并计算QScasecopper:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_COP_CON;break;caseair:QS=(A.mNodeTi+1j-A.mNodeTij)*CONV*DX;break;casesteel:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_

26、CON_STL;Ibreak;caseconcrete:QS=(A.mNodeTi+1j-A.mNodeTij)*COND_CONCRETE;break;Icaseiso834:break;caseadiabat:QS=0;break;default:std:cout<<i+1<<","<<j<<结点属性未指定n"getchar();break;A.mNodeT皿=T_AIR;break;default:break;A.mNodeTij=(QE+QN+QS+QW)/(A.mNodeC皿*A.mNodeDenij*DX*DY);函数会判断每个结点相邻4个结点的属性，然后将四个结点传来的热量加和，/最后根据能量守恒计算出中间结点的温度值(时间全显式)voidwrite(intt,Grid&A)c