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

精心整理#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/(m^3)#defineC_COPPER390#defineC_STEEL460#defineC_CONCRETE970材料比热容，J/(kg*K)#defineCOND_COPPER400#defineCOND_STEEL20#defineCOND_CONCRETE1.63材料导热系数，W/(m*K)#defineCONV5空气与各种材料的对流换热系数，W/(m^2*K)#defineINITL_TMPRT298初始温度，K#defineT_C2K(x)((x)+273)//CtoK#defineT_K2C(x)((x)-273)//KtoC摄氏温度与绝对温度相互转换#defineT_AIR298;空气温度,K#definet_END3600计算时长，s#definet_STEP1时间步长，s#defineCOND_COP_CON3.2467#defineCOND_COP_STL38.0952#defineCOND_CON_STL3.0143精心整理//不同材料的两个相邻结点，导热系数取调和平均数，

W/(m*K)////////////////////////////////////////////////////////////////////////////////

计算参数宏定义 结束////////////////////////////////////////////////////////////////////////////////enumProperty{copper,steel,concrete,iso834,adiabat,air};结点属性的枚举类型声明structGrid{doublemNodeT[GRID_ROW][GRID_COLUMN];//结点温度PropertymNodeProperty[GRID_ROW][GRID_COLUMN];//结点材料属性doublemNodeC[GRID_ROW][GRID_COLUMN];//结点比热容doublemNodeDen[GRID_ROW][GRID_COLUMN];//结点密度};结构体Grid定义，此为计算的主要承载数据结构，将在计算中存储温度数据函数全局声明开始voidSetProperty(Grid&);//属性设置函数voidInitlTMPRT(Grid&);//温度初始化函数voidSetNodePara(Grid&);//结点物性参数设置函数voidCalculateTMPRT(Grid&);//温度计算函数，几乎全部数值计算的承担者voidwrite(int,Grid&);//写文件函数，每隔 100秒输出温度分布文件到当前目录函数全局说明结束voidSetProperty(Grid&A){for(intj=0;j<GRID_COLUMN;++j)A.mNodeProperty[0][j]=air;第一行为空气边界，计算区域外for(intj=0;j<GRID_COLUMN;++j)A.mNodeProperty[GRID_ROW-1][j]=adiabat;最后一行为绝热边界，计算区域外for(inti=0;i<GRID_ROW;++i)if((i>=0&&i<=40)||(i>=61&&i<GRID_ROW))A.mNodeProperty[i][0]=air;elseA.mNodeProperty[i][0]=iso834;第一列的空气边界和ISO834边界，计算区域外for(inti=0;i<GRID_ROW;++i)A.mNodeProperty[i][GRID_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>=61&&i<GRID_ROW-2&&j>=1&&j<=(100-(80-i)*(49/19))))A.mNodeProperty[i][j]=concrete;混凝土区域elseif((i-20)*(i-20)+(j-25)*(j-25)<=20*20)A.mNodeProperty[i][j]=air;圆孔空气区域精心整理elseif(i>=41&&i<=60&&j>=1&&j<=90&&((i-20)*(i-20)+(j-25)*(j-25)>=20*20))A.mNodeProperty[i][j]=copper;金属铜区域的一部分elseif(i>=1&&i<+40&&j>=21&&j<=90)A.mNodeProperty[i][j]=copper;金属铜区域的另一部分elseA.mNodeProperty[i][j]=steel;剩下的即为金属铁的区域}}//属性设置函数，将Grid结构体A的所有网格赋予材料属性（枚举类型Property）voidInitlTMPRT(Grid&A){for(inti=0;i<GRID_ROW;++i)for(intj=0;j<GRID_COLUMN;++j)A.mNodeT[i][j]=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.mNodeProperty[i][j]){{switch(A.mNodeProperty[i][j+1])//判断右侧属性，并计算 QE{casecopper:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_COPPER;break;caseair:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*CONV*DY;break;casesteel:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_COP_STL;break;caseconcrete:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_COP_CON;break;caseiso834:QE=2*(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_COPPER;break;caseadiabat:QE=0;精心整理break;default:std::cout<<i<<","<<j+1<<结"点属性未指定\n";getchar();}std::cout<<QE<<"\n";switch(A.mNodeProperty[i][j-1])//判断左侧属性，并计算 QW{casecopper:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_COPPER;break;caseair:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*CONV*DY;break;casesteel:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_COP_STL;break;caseconcrete:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_COP_CON;break;caseiso834:QW=2*(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_COPPER;break;caseadiabat:QW=0;break;default:{std::cout<<i<<","<<j-1<<结"点属性未指定\n";getchar();}}switch(A.mNodeProperty[i-1][j])//判断上侧属性，并计算 QN{casecopper:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_COPPER;break;caseair:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*CONV*DX;break;casesteel:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_COP_STL;break;caseconcrete:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_COP_CON;break;caseiso834:QN=2*(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_COPPER;break;精心整理caseadiabat:QN=0;break;default:std::cout<<i-1<<","<<j<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i+1][j])//判断下侧属性，并计算 QS{casecopper:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_COPPER;break;caseair:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*CONV*DX;break;casesteel:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_COP_STL;break;caseconcrete:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_COP_CON;break;caseiso834:QS=2*(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_COPPER;break;caseadiabat:QS=0;break;default:std::cout<<i+1<<","<<j<<结"点属性未指定\n";getchar();}break;}{switch(A.mNodeProperty[i][j+1])//判断右侧属性，并计算 QE{casecopper:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_COP_STL;break;caseair:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*CONV*DY;break;casesteel:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_STEEL;break;caseconcrete:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_CON_STL;break;caseiso834:QE=2*(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_STEEL;精心整理break;caseadiabat:QE=0;break;default:std::cout<<i<<","<<j+1<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i][j-1])//判断左侧属性，并计算 QW{casecopper:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_COP_STL;break;caseair:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*CONV*DY;break;casesteel:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_STEEL;break;caseconcrete:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_CON_STL;break;caseiso834:QW=2*(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_STEEL;break;caseadiabat:QW=0;break;default:std::cout<<i<<","<<j-1<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i-1][j])//判断上侧属性，并计算 QN{casecopper:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_COP_STL;break;caseair:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*CONV*DX;break;casesteel:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_STEEL;break;caseconcrete:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_CON_STL;break;caseiso834:精心整理QN=2*(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_STEEL;break;caseadiabat:QN=0;break;default:std::cout<<i-1<<","<<j<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i+1][j])//判断下侧属性，并计算 QS{casecopper:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_COP_STL;break;caseair:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*CONV*DX;break;casesteel:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_STEEL;break;caseconcrete:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_CON_STL;break;caseiso834:QS=2*(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_STEEL;break;caseadiabat:QS=0;break;default:std::cout<<i+1<<","<<j<<结"点属性未指定\n";getchar();}break;}{switch(A.mNodeProperty[i][j+1])//判断右侧属性，并计算 QE{casecopper:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_COP_CON;break;caseair:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*CONV*DY;break;casesteel:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_CON_STL;break;caseconcrete:QE=(A.mNodeT[i][j+1]-A.mNodeT[i][j])*COND_CONCRETE;精心整理break;caseiso834:break;caseadiabat:QE=0;break;default:std::cout<<i<<","<<j+1<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i][j-1])//判断左侧属性，并计算 QW{casecopper:break;caseair:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*CONV*DY;break;casesteel:break;caseconcrete:QW=(A.mNodeT[i][j-1]-A.mNodeT[i][j])*COND_CONCRETE;break;caseiso834:break;caseadiabat:QW=0;break;default:std::cout<<i<<","<<j-1<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i-1][j])//判断上侧属性，并计算 QN{casecopper:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_COP_CON;break;caseair:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*CONV*DX;break;casesteel:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_CON_STL;break;caseconcrete:QN=(A.mNodeT[i-1][j]-A.mNodeT[i][j])*COND_CONCRETE;break;caseiso834:break;精心整理caseadiabat:QN=0;break;default:std::cout<<i-1<<","<<j<<结"点属性未指定\n";getchar();}switch(A.mNodeProperty[i+1][j])//判断下侧属性，并计算 QS{casecopper:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_COP_CON;break;caseair:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*CONV*DX;break;casesteel:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_CON_STL;break;caseconcrete:QS=(A.mNodeT[i+1][j]-A.mNodeT[i][j])*COND_CONCRETE;break;caseiso834:break;caseadiabat:QS=0;break;default:std::cout<<i+1<<","<<j<<结"点属性未指定\n";getchar();}break;}{A.mNodeT[i][j]=T_AIR;break;}default:break;}A.mNodeT[i][j]=(QE+QN+QS+QW)/(A.mNodeC[i][j]*A.mNodeDen[i][j]*DX*DY);}}函数会判断每个结点相邻4个结点的属性，然后将四个结点传来的热量加和,最后根据能量守恒计算出中间结点的温度值（时间全显式）voidwrite(intt,Grid&A){charfilename[40];intx=1,y=1;sprintf_s(filename,"temperature%d.csv",t);精心整理std::ofstreamfout(filename);fout<<"