二维数值传热学C 程序代码 (2)

精心整理精心整理精心整理#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<<"X,Y,TEMPERATURE\n";for(inti=GRID_ROW-2;i>0;--i){for(intj=1;j<GRID_COLUMN-1;++j){if(y>100)y=1;fout<<x<<","<<y<<","<<A.mNodeT[i][j]<<"\n";++y;}++x;}}voidSetNodePara(Grid&A){for(inti=0;i<GRID_ROW;++i)for(intj=0;j<GRID_COLUMN;++j)switch(A.mNodeProperty[i][j]){casecopper:A.mNodeC[i][j]=C_COPPER;A.mNodeDen[i][j]=DEN_COPPER;break;casesteel:A.mNodeC[i][j]=C_STEEL;A.mNodeDen[i][j]=DEN_STEEL;break;caseconcrete:A.mNodeC[i][j]=C_CONCRETE;A.mNodeDen[i][j]=DEN_CONCRETE;break;default:A.mNodeC[i][j]=1;A.mNodeDen[i][j]=1;}}/////////////////////////////////////////////////////////////////////////////////******************************************************************************///////////////////////////////////////////////////////////////////////////////////主函数开始intmain(){GridP;SetProperty(P);InitlTMPRT(P);SetNodePara(P);for(intt=1;t<=t_END;++t){for(inti=41;i<=60;++i)P.mNodeT[i][0]=T_C2K(345*log10(8*t/60+1)+26);CalculateTMPRT(P);if(!(t%100))write(t,P);}}