📄 tempf2d.h
字号:
void TempF2D();
//*****钢筋混凝土矩形截面瞬态温度场非线性有限元程序*****//
//***华南理工大学建筑学院土木工程系结构工程2003级硕士研究生唐贵和编制***//
//***说明:***
//***1.本程序可用来计算在国际ISO834、美国升温曲线ASTE E119、自定义升温曲线(每一分//
//钟输入一个环境温度)下的钢筋混凝土截面二维温度场;//
//***2.程序中温度单元采用四节点矩形单元(过镇海);//
//***3.采用定时间步长求解策略和变时间步长两种求解策略//
//***4.可以计算构件在单面、相邻两面、相对两面、三面、四面受火形式下的温度场;//
//***5.可以计算普通混凝土构件温度场,近似计算高强混凝土构件的温度场;//
//***6.为简化分析过程,考虑如下基本假定:忽略钢筋表面与混凝土之间的热阻,钢筋的温度
//直接采用钢筋质心处的混凝土温度;混凝土中裂缝对温度场的影响忽略不计
//***7.(a).不考虑含水量
//(b).按照Lie-韩林海方法。
//(c).按照欧洲规范(EC94)近似考虑含湿率对温度场的影响。
//(d).徐玉野方法
//***8.(a).Galerkin法,参照过镇海
//(b).两点向后差分法(implicit Euler method)
long *BandWidth; //一维变带变带宽存储法的带宽
long *DiagAddress; //一维变带变带宽存储法的对角地址
long TotalLen=0; //系数矩阵在一维数组存储的数量
//**********函数声明***************************//
double CalcuCSpecialHeat(double,int);
double CalcuKc(double,int);
void CalcuBt(double &,double);
void Form_KTG(double *,int,double *);
void Form_CcG(double *,int,double *);
void Form_KBG_PBG(double *,double *,double *,double *);
void Form_krG_PrG(double,double *,double *,double *,double *,double *);
void Multiply_1D(double *,double *,double *);
void Gauss_1D(double *,double *,double *);
void Cal_Temp_TiQu(int,double *,double *,double *,double *);
void Cal_Spall_Ind(double,double,int &);
double Cal_TimeStep(double,double);
void TempF2D()
{
ifstream Tempfin("Temp2D.txt");
ifstream BoundTemp("边界温度.xls");
ifstream CircumTemp("环境温度.xls");
ofstream OutResults("Struc_Tin.txt");
ofstream TimeOut("TimeHistory.txt");
ofstream OutTemp("温度提取.csv");
ofstream Temperp("Temperp.csv");
ofstream YunTu1("云图提取1.csv");
ofstream YunTu2("云图提取2.csv");
ofstream YunTu("云图.csv");
ofstream WResult("WResult.txt");
if(Tempfin.fail())
{
cerr<<"error opening the Temp2D.txt!\n";
exit(0);
}
WResult<<"\n********!!!钢筋混凝土矩形截面瞬态温度场非线性有限元计算程序!!!*********\n"<<endl;
WResult<<" **************************************************************\n"<<endl;
WResult<<"\n***!!!华南理工大学建筑学院土木工程系结构工程2003级硕士研究生唐贵和编制!!!***\n"<<endl;
WResult<<"***!!!下面进行瞬态温度场分析!!!***\n"<<endl;
cout<<"\n********!!!钢筋混凝土矩形截面瞬态温度场非线性有限元计算程序!!!*********\n"<<endl;
cout<<" **************************************************************"<<endl;
cout<<"\n***!!!华南理工大学建筑学院土木工程系结构工程2003级硕士研究生唐贵和编制!!!***\n"<<endl;
cout<<"******!!!!!下面进行瞬态温度场分析!!!!!******\n"<<endl;
TNode_num=(Tdiv1+divProt*2+1)*(Tdiv2+divProt*2+1); //温度场分析时节点总数目
TConElem_num=(Tdiv1+divProt*2)*(Tdiv2+divProt*2); //温度场分析时单元总数目
double AfSectwid,AfSectlen; //爆裂后的构件截面边长
double BoundTemper; //不受火边的温度(注意:随时间变化)
int Bound; //
double TotalTime,SpallTime; //控制柱子火烧总时间,爆裂时间 时间以秒为单位
double TimeStep=0.0; //时间步长
int CCondu_Type,CSpeciheat_Type; //混凝土导热系数和热容的选择
double *Hc,*Hr; //对流换热系数,辐射换热系数
double Defer; //炉对构件表面辐射延迟系数
int Bord_heatType=0; //边界的处理方法,1-分别考虑Hc和Hr,2-综合考虑
int ComputeMethod; //计算方法,1-Galerkin法,2-两点向后差分法
int Spall_Ind=0; //判断点是否已经处于明火中,已爆
int NumTime=0; //累计时间步的数量
int NumDieDai=0; //累计时间步内迭代的次数
int IterMax; //每一时间步内迭代的最大次数
double Time=0.0; //用来累积时间
double Temper; //用来累积温度(炉子内部的温度)
double FNMaxErr; //前后两次迭代容差
double error=0; //前后两次迭代最大温差
int Num_TiQu=0; //需要提取温度点的个数
int Points=0; //节点号提取法的点数
double *CorX_TiQu,*CorY_TiQu; //需要提取温度点的坐标
double *Temp_TiQu; //需要提取温度点的温度
double *PointsTemp; //由节点号提取点的温度
int *PointsNum; //提取点的节点号
int ContourNum; //提取云图的时刻个数
double *Contour; //提取云图的时刻
double *AfNode_Coor; //爆裂后单元节点坐标(相对左下角点)
int loop,loop1; //循坏变量
int ibuf=0;
double *KTG,*CcG,*KG,*KBG,*KrG; //刚度矩阵
double *PBG,*PrG,*PG;
double *A,*B,*C; //解方程用来简化表达的数组
double *CcG0,*KG0,*KTG0,*KBG0,*KrG0;
double *PBG0,*PrG0,*PG0;
double *EleAverTemp; //单元平均温度
double *AireTemp; //柱四面的空气温度
double *Temp0,*Temp1,*Temp,*TempLin;
double *Linshi; //两点后差法时的温度增量TNode_num*1
char NoUse[200]; //为输入设置的无用的变量
Hc=(double *)calloc(4,8);
Hr=(double *)calloc(4,8);
//**********从文件读取初始数据**********
Tempfin.getline(NoUse,sizeof(NoUse)); //受火形式
Tempfin>>Fire0>>Fire1>>Fire2>>Fire3;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Bound;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>TotalTime>>SpallTime>>TimeStep;
TotalTime=TotalTime*60.0; //把时间变成秒单位制
SpallTime=SpallTime*60.0;
TimeStep=TimeStep*60.0;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>IterMax>>FNMaxErr;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>CCondu_Type>>CSpeciheat_Type;
Tempfin.getline(NoUse,sizeof(NoUse));//含水量热工参数
Tempfin>>CWaterContent;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Qw>>Cw>>PWater;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>CWaterPeakT;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Water_Method;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>PConcrete;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Bord_heatType;
Tempfin.getline(NoUse,sizeof(NoUse));//对流换热系数
Tempfin>>Hc[0]>>Hc[1]>>Hc[2]>>Hc[3];
Tempfin.getline(NoUse,sizeof(NoUse));//辐射换热系数
Tempfin>>Hr[0]>>Hr[1]>>Hr[2]>>Hr[3];
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Defer;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>ComputeMethod;
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>ContourNum;
if(ContourNum>=1)
{
Contour=(double *)calloc(ContourNum,8);
Tempfin.getline(NoUse,sizeof(NoUse));
for(loop=0; loop<ContourNum; loop++)
Tempfin>>Contour[loop];
}
if(ContourNum<1) //不需要提取云图
{
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin.getline(NoUse,sizeof(NoUse));
}
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Num_TiQu;
if(Num_TiQu>=1)
{
CorX_TiQu=(double *)calloc(Num_TiQu,8);
CorY_TiQu=(double *)calloc(Num_TiQu,8);
Temp_TiQu=(double *)calloc(Num_TiQu,8);
Tempfin.getline(NoUse,sizeof(NoUse));
for(loop=0; loop<Num_TiQu; loop++)
Tempfin>>CorX_TiQu[loop];
Tempfin.getline(NoUse,sizeof(NoUse));
for(loop=0; loop<Num_TiQu; loop++)
Tempfin>>CorY_TiQu[loop];
}
if(Num_TiQu<1) //没有依靠坐标提取点
{
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin.getline(NoUse,sizeof(NoUse));
}
Tempfin.getline(NoUse,sizeof(NoUse));
Tempfin>>Points;
if(Points>=1)
{
PointsTemp=(double *)calloc(Points,8);
PointsNum=(int *)calloc(Points,2);
Tempfin.getline(NoUse,sizeof(NoUse));
for(loop=0;loop<Points;loop++)
Tempfin>>PointsNum[loop];
}
WResult<<"***************************************************************"<<endl;
WResult<<"***!!!程序状态: 输入的信息已经读入!!!***"<<endl;
if(ComputeMethod==1)
{
WResult<<"***!!!采用过镇海推导的Galerkin法迭代式,ComputeMethod= "<<ComputeMethod<<endl;
cout<<"***!!!采用过镇海推导的Galerkin法迭代式,ComputeMethod= "<<ComputeMethod<<endl;
}
else if(ComputeMethod==2)
{
WResult<<"****!!!采用两点向后差分法而推导的迭代式,ComputeMethod= "<<ComputeMethod<<endl;
cout<<"***!!!采用两点向后差分法而推导的迭代式,ComputeMethod= "<<ComputeMethod<<endl;
}
else
{
WResult<<"****!!!计算方法选择有错,ComputeMethod= "<<ComputeMethod<<endl;
cout<<"***!!!计算方法选择有错,ComputeMethod= "<<ComputeMethod<<endl;
exit(0);
}
if(VarTstep==1)
{
WResult<<"***!!!采用定时间步长的求解策略,VarTstep= "<<VarTstep<<endl;
cout<<"***!!!采用定时间步长的求解策略,VarTstep= "<<VarTstep<<endl;
}
else if(VarTstep==2)
{
WResult<<"***!!!采用变时间步长的求解策略,VarTstep= "<<VarTstep<<endl;
cout<<"***!!!采用变时间步长的求解策略,VarTstep= "<<VarTstep<<endl;
}
else
{
WResult<<"****!!!求解策略选择有错,VarTstep= "<<VarTstep<<endl;
cout<<"***!!!求解策略选择有错,VarTstep= "<<VarTstep<<endl;
exit(0);
}
if(Water_Method==0)
{
WResult<<"***!!!不考虑钢筋混凝土含湿量对温度场的影响,Water_Method= "<<Water_Method<<endl;
cout<<"***!!!不考虑钢筋混凝土含湿量对温度场的影响,Water_Method= "<<Water_Method<<endl;
}
else if(Water_Method==1 || Water_Method==2)
{
WResult<<"***!!!考虑钢筋混凝土含湿量对温度场的影响****"<<endl;
WResult<<"***!!!:含湿量!!!***: "<<CWaterContent*100<<"%"<<endl;
WResult<<"***!!!:含湿量采用的计算方法是!!!***: ";
cout<<"***!!!考虑钢筋混凝土含湿量对温度场的影响*****"<<endl;
cout<<"***!!!:含湿量!!!***: "<<CWaterContent*100<<"%"<<endl;
cout<<"***!!!:含湿量采用的计算方法是!!!***: ";
if(Water_Method==1)
{
WResult<<" 参考Lie-韩林海的方法!!!"<<endl;
cout<<" 参考Lie-韩林海的方法!!!"<<endl;
}
else if(Water_Method==2)
{
WResult<<" 参考欧洲规范(EC94)的方法!!!"<<endl;
cout<<" 参考欧洲规范(EC94)的方法!!!"<<endl;
}
else if(Water_Method==3){
WResult<<" 徐玉野建议的方法!!!"<<endl;
cout<<" 徐玉野建议的方法!!!"<<endl;
}
}
else
{
WResult<<"***!!!含水量方法选择有错(Water_Method)!!! "<<Water_Method<<endl;
cout<<"***!!!含水量方法选择有错(Water_Method)!!! "<<Water_Method<<endl;
exit(0);
}
if(Bord_heatType==1)
{
WResult<<"***!!!分开考虑辐射和对流换热******Bord_heatType= "<<Bord_heatType<<endl;
cout<<"***!!!分开考虑辐射和对流换热******Bord_heatType= "<<Bord_heatType<<endl;
}
else if(Bord_heatType==2)
{
WResult<<"***!!!采用综合换热系数(段文玺建议的1985)***Bord_heatType= "<<Bord_heatType<<endl;
cout<<"***!!!采用综合换热系数(段文玺建议的1985)****Bord_heatType= "<<Bord_heatType<<endl;
}
else
{
WResult<<"***!!!换热系数方法选择有错(Bord_heatType)!!! "<<Bord_heatType<<endl;
cout<<"***!!!换热系数方法选择有错(Bord_heatType)!!! "<<Bord_heatType<<endl;
exit(0);
}
AfNode_Coor=(double *)calloc(TConElem_num*2,8);
Node_Coor=(double *)calloc(TConElem_num*2,8);
EleAverTemp=(double *)calloc(TConElem_num,8);
AireTemp =(double *)calloc(4,8);
Temp0=(double *)calloc(TNode_num,8);
Temp1=(double *)calloc(TNode_num,8);
Temp =(double *)calloc(TNode_num,8);
TempLin=(double *)calloc(TNode_num,8);
Linshi=(double *)calloc(TNode_num,8);
for(loop=0; loop<TConElem_num; loop++)
EleAverTemp[loop]=iniT0;
for(loop=0; loop<TNode_num; loop++)
{
Temp1[loop]=iniT0; Temp[loop] =iniT0;
Temp0[loop]=iniT0; TempLin[loop]=iniT0;
Linshi[loop]=0.0;
}
//输出提取点温度
OutTemp<<setiosflags(ios::fixed);
OutTemp<<setprecision(5)<<0.0<<",";
for(loop=0; loop<Num_TiQu; loop++)
{
OutTemp<<setiosflags(ios::fixed);
OutTemp<<setprecision(5)<<iniT0<<",";
}
OutTemp<<"\n";
Temperp<<setiosflags(ios::fixed);
Temperp<<setprecision(5)<<0.0<<",";
for(loop=0; loop<Points; loop++)
{
Temperp<<setiosflags(ios::fixed);
Temperp<<setprecision(5)<<iniT0<<",";
}
Temperp<<"\n";
TCElemNode_Info=(long *)calloc(TConElem_num*4,4);//单元结点整体编码
Node_Coor=(double *)calloc(TNode_num*2,8);
BandWidth=(long *)calloc(TNode_num,4);
DiagAddress=(long *)calloc(TNode_num,4);
//单元节点编号信息
/*for(loop=0;loop<2*divProt+Tdiv1;loop++)
{
for(loop1=0;loop1<2*divProt+Tdiv2;loop1++)
{
*(TCElemNode_Info+4*(loop*(2*divProt+Tdiv2)+loop1)+0)=loop*(2*divProt+Tdiv2+1)+loop1;
*(TCElemNode_Info+4*(loop*(2*divProt+Tdiv2)+loop1)+1)=loop*(2*divProt+Tdiv2+1)+loop1+1;
*(TCElemNode_Info+4*(loop*(2*divProt+Tdiv2)+loop1)+2)=(loop+1)*(2*divProt+Tdiv2+1)+loop1;
*(TCElemNode_Info+4*(loop*(2*divProt+Tdiv2)+loop1)+3)=(loop+1)*(2*divProt+Tdiv2+1)+loop1+1;
}
}*/
for(loop=0; loop<TConElem_num; loop++)
{
*(TCElemNode_Info+4*loop+0)=loop+int(loop/(Tdiv1+divProt*2));
*(TCElemNode_Info+4*loop+1)=loop+int(loop/(Tdiv1+divProt*2))+1;
*(TCElemNode_Info+4*loop+2)=loop+int(loop/(Tdiv1+divProt*2))+(Tdiv1+divProt*2)+1;
*(TCElemNode_Info+4*loop+3)=loop+int(loop/(Tdiv1+divProt*2))+(Tdiv1+divProt*2)+2;
}
//计算带宽
for(loop=0; loop<TConElem_num; loop++)//带宽的计算
{
ibuf=*(TCElemNode_Info+4*loop+0); //选出单元中节点编号最小的节点
if(ibuf>*(TCElemNode_Info+4*loop+1)) ibuf=*(TCElemNode_Info+4*loop+1);
if(ibuf>*(TCElemNode_Info+4*loop+2)) ibuf=*(TCElemNode_Info+4*loop+2);
if(ibuf>*(TCElemNode_Info+4*loop+3)) ibuf=*(TCElemNode_Info+4*loop+3);
if(BandWidth[*(TCElemNode_Info+4*loop+0)]<=(*(TCElemNode_Info+4*loop+0)-ibuf))
BandWidth[*(TCElemNode_Info+4*loop+0)]=*(TCElemNode_Info+4*loop+0)-ibuf+1;
if(BandWidth[*(TCElemNode_Info+4*loop+1)]<=(*(TCElemNode_Info+4*loop+1)-ibuf))
BandWidth[*(TCElemNode_Info+4*loop+1)]=*(TCElemNode_Info+4*loop+1)-ibuf+1;
if(BandWidth[*(TCElemNode_Info+4*loop+2)]<=(*(TCElemNode_Info+4*loop+2)-ibuf))
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -