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

1、二维数值传热学C+程序代码#include#include#include#include/*/从此开始至下一个/*/标记带处之间的代码可以写在头文件里。/计算参数宏定义 开始/#define GRID_ROW 82#define GRID_COLUMN 102 /网格横纵数目，边界处多出一列处理#define DX 0.05#define DY 0.05/网格单元尺寸，m#define DEN_COPPER 8900#define DEN_STEEL 7850#define DEN_CONCRETE 2430 /材料的密度，kg/(m3)#define C_COPPER 390#define

2、 C_STEEL 460#define C_CONCRETE 970/材料比热容，J/(kg*K)#define COND_COPPER 400#define COND_STEEL 20#define COND_CONCRETE 1.63/材料导热系数，W/(m*K)#define CONV 5/空气与各种材料的对流换热系数，W/(m2*K)#define INITL_TMPRT 298/初始温度，K#define T_C2K(x) (x)+273)/C to K#define T_K2C(x) (x)-273)/K to C/摄氏温度与绝对温度相互转换#define T_AIR 298;/空

3、气温度,K#define t_END 3600/计算时长，s#define t_STEP 1/时间步长，s#define COND_COP_CON 3.2467#define COND_COP_STL 38.0952#define COND_CON_STL 3.0143/不同材料的两个相邻结点，导热系数取调和平均数，W/(m*K)/计算参数宏定义 结束/enum Property copper, steel, concrete, iso834, adiabat, air ;/结点属性的枚举类型声明struct Grid double mNodeTGRID_ROWGRID_COLUMN;/结点温

4、度 Property mNodePropertyGRID_ROWGRID_COLUMN;/结点材料属性 double mNodeCGRID_ROWGRID_COLUMN;/结点比热容 double mNodeDenGRID_ROWGRID_COLUMN;/结点密度;/结构体Grid定义，此为计算的主要承载数据结构，将在计算中存储温度数据/函数全局声明开始void SetProperty(Grid &);/属性设置函数void InitlTMPRT(Grid &);/温度初始化函数void SetNodePara(Grid &);/结点物性参数设置函数void CalculateTMPRT(Gr

5、id &);/温度计算函数，几乎全部数值计算的承担者void write(int,Grid &);/写文件函数，每隔100秒输出温度分布文件到当前目录/函数全局说明结束void SetProperty(Grid &A) for (int j = 0; j GRID_COLUMN;+j) A.mNodeProperty0j = air; /第一行为空气边界，计算区域外 for (int j = 0; j GRID_COLUMN; +j) A.mNodePropertyGRID_ROW-1j = adiabat; /最后一行为绝热边界，计算区域外 for (int i = 0; i = 0 & i

6、 = 61 & i GRID_ROW) A.mNodePropertyi0 = air; else A.mNodePropertyi0 = iso834; /第一列的空气边界和ISO834边界，计算区域外 for (int i = 0; i GRID_ROW; +i) A.mNodePropertyiGRID_COLUMN - 1 = air; /最后一列，空气边界，计算区域外 for (int i = 1; i GRID_ROW - 1; +i) for (int j = 1; j = 1 & i = 1 & j = 61 & i = 1 & j = (100 - (80 - i)*(49

7、/ 19) A.mNodePropertyij = concrete; /混凝土区域 else if (i - 20)*(i - 20) + (j - 25)*(j - 25) = 41 & i = 1 & j = 20 * 20) A.mNodePropertyij = copper; /金属铜区域的一部分 else if (i = 1 & i = 21 & j = 90) A.mNodePropertyij = copper; /金属铜区域的另一部分 else A.mNodePropertyij = steel; /剩下的即为金属铁的区域 /属性设置函数，将Grid结构体A的所有网格赋予材

8、料属性（枚举类型 Property）void InitlTMPRT(Grid &A) for (int i = 0; i GRID_ROW ; +i) for (int j = 0; j GRID_COLUMN ; +j) A.mNodeTij = INITL_TMPRT;/温度初始化函数，将Grid结构体A的全部结点温度设置为初始温度INITL_TMPRTvoid CalculateTMPRT(Grid &A) double QE=0, QW=0, QN=0, QS=0;/四个方向传来的热量 for (int i = 1; i GRID_ROW - 1; +i) for (int j = 1

9、; j GRID_COLUMN - 1; +j) switch (A.mNodePropertyij) switch (A.mNodePropertyij + 1)/判断右侧属性，并计算QE case copper: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_COPPER; break; case air: QE = (A.mNodeTij + 1 - A.mNodeTij)*CONV*DY; break; case steel: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_COP_STL; break; case co

10、ncrete: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_COP_CON; break; case iso834: QE = 2 * (A.mNodeTij + 1 - A.mNodeTij)*COND_COPPER; break; case adiabat: QE = 0; break; default: std:cout i , j + 1 结点属性未指定n; getchar(); / std:cout QE n; switch (A.mNodePropertyij - 1)/判断左侧属性，并计算QW case copper: QW = (A.mNod

11、eTij - 1 - A.mNodeTij)*COND_COPPER; break; case air: QW = (A.mNodeTij - 1 - A.mNodeTij)*CONV*DY; break; case steel: QW = (A.mNodeTij - 1 - A.mNodeTij)*COND_COP_STL; break; case concrete: QW = (A.mNodeTij - 1 - A.mNodeTij)*COND_COP_CON; break; case iso834: QW = 2 * (A.mNodeTij - 1 - A.mNodeTij)*COND_

12、COPPER; break; case adiabat: QW = 0; break; default: std:cout i , j - 1 结点属性未指定n; getchar(); switch (A.mNodePropertyi - 1j)/判断上侧属性，并计算QN case copper: QN = (A.mNodeTi - 1j - A.mNodeTij)*COND_COPPER; break; case air: QN = (A.mNodeTi - 1j - A.mNodeTij)*CONV*DX; break; case steel: QN = (A.mNodeTi - 1j -

13、 A.mNodeTij)*COND_COP_STL; break; case concrete: QN = (A.mNodeTi - 1j - A.mNodeTij)*COND_COP_CON; break; case iso834: QN = 2 * (A.mNodeTi - 1j - A.mNodeTij)*COND_COPPER; break; case adiabat: QN = 0; break; default: std:cout i - 1 , j 结点属性未指定n; getchar(); switch (A.mNodePropertyi + 1j)/判断下侧属性，并计算QS c

14、ase copper: QS = (A.mNodeTi + 1j - A.mNodeTij)*COND_COPPER; break; case air: QS = (A.mNodeTi + 1j - A.mNodeTij)*CONV*DX; break; case steel: QS = (A.mNodeTi + 1j - A.mNodeTij)*COND_COP_STL; break; case concrete: QS = (A.mNodeTi + 1j - A.mNodeTij)*COND_COP_CON; break; case iso834: QS = 2 * (A.mNodeTi

15、+ 1j - A.mNodeTij)*COND_COPPER; break; case adiabat: QS = 0; break; default: std:cout i + 1 , j 结点属性未指定n; getchar(); break; switch (A.mNodePropertyij + 1)/判断右侧属性，并计算QE case copper: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_COP_STL; break; case air: QE = (A.mNodeTij + 1 - A.mNodeTij)*CONV*DY; break; ca

16、se steel: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_STEEL; break; case concrete: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_CON_STL; break; case iso834: QE = 2 * (A.mNodeTij + 1 - A.mNodeTij)*COND_STEEL; break; case adiabat: QE = 0; break; default: std:cout i , j + 1 结点属性未指定n; getchar(); switch (A.mNodeP

17、ropertyij - 1)/判断左侧属性，并计算QW case copper: QW = (A.mNodeTij - 1 - A.mNodeTij)*COND_COP_STL; break; case air: QW = (A.mNodeTij - 1 - A.mNodeTij)*CONV*DY; break; case steel: QW = (A.mNodeTij - 1 - A.mNodeTij)*COND_STEEL; break; case concrete: QW = (A.mNodeTij - 1 - A.mNodeTij)*COND_CON_STL; break; case

18、iso834: QW = 2 * (A.mNodeTij - 1 - A.mNodeTij)*COND_STEEL; break; case adiabat: QW = 0; break; default: std:cout i , j - 1 结点属性未指定n; getchar(); switch (A.mNodePropertyi - 1j)/判断上侧属性，并计算QN case copper: QN = (A.mNodeTi - 1j - A.mNodeTij)*COND_COP_STL; break; case air: QN = (A.mNodeTi - 1j - A.mNodeTij

19、)*CONV*DX; break; case steel: QN = (A.mNodeTi - 1j - A.mNodeTij)*COND_STEEL; break; case concrete: QN = (A.mNodeTi - 1j - A.mNodeTij)*COND_CON_STL; break; case iso834: QN = 2 * (A.mNodeTi - 1j - A.mNodeTij)*COND_STEEL; break; case adiabat: QN = 0; break; default: std:cout i - 1 , j 结点属性未指定n; getchar

20、(); switch (A.mNodePropertyi + 1j)/判断下侧属性，并计算QS case copper: QS = (A.mNodeTi + 1j - A.mNodeTij)*COND_COP_STL; break; case air: QS = (A.mNodeTi + 1j - A.mNodeTij)*CONV*DX; break; case steel: QS = (A.mNodeTi + 1j - A.mNodeTij)*COND_STEEL; break; case concrete: QS = (A.mNodeTi + 1j - A.mNodeTij)*COND_C

21、ON_STL; break; case iso834: QS = 2 * (A.mNodeTi + 1j - A.mNodeTij)*COND_STEEL; break; case adiabat: QS = 0; break; default: std:cout i + 1 , j 结点属性未指定n; getchar(); break; switch (A.mNodePropertyij + 1)/判断右侧属性，并计算QE case copper: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_COP_CON; break; case air: QE = (

22、A.mNodeTij + 1 - A.mNodeTij)*CONV*DY; break; case steel: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_CON_STL; break; case concrete: QE = (A.mNodeTij + 1 - A.mNodeTij)*COND_CONCRETE; break; case iso834: break; case adiabat: QE = 0; break; default: std:cout i , j + 1 结点属性未指定n; getchar(); switch (A.mNodePr

23、opertyij - 1)/判断左侧属性，并计算QW case copper: break; case air: QW = (A.mNodeTij - 1 - A.mNodeTij)*CONV*DY; break; case steel: break; case concrete: QW = (A.mNodeTij - 1 - A.mNodeTij)*COND_CONCRETE; break; case iso834: break; case adiabat: QW = 0; break; default: std:cout i , j - 1 结点属性未指定n; getchar(); switch (A.mNodePr