📄 tempf2d.h
字号:
OutTemp<<setprecision(6)<<*(Temp_TiQu+loop)<<",";
}
OutTemp<<"\n";
for(loop=0; loop<Points; loop++)
{
for(loop1=0; loop1<TNode_num; loop1++)
{
if(*(PointsNum+loop)==loop1)
PointsTemp[loop]=Temp0[loop1];
}
}
Temperp<<setiosflags(ios::fixed)<<setprecision(5)<<Time/60<<",";
for(loop=0; loop<Points; loop++)
{
Temperp<<setiosflags(ios::fixed);
Temperp<<setprecision(5)<<*(PointsTemp+loop)<<",";
}
Temperp<<"\n";
//*******************指定时刻的云图***************
for(loop=0; loop<ContourNum; loop++)
{
if(Contour[loop]==Time)
{
YunTu<<"t= "<<Time<<" 时刻的云图 "<<endl;
for(loop=0; loop<Tdiv2+divProt*2+1; loop++)
{
for(loop1=0; loop1<Tdiv1+divProt*2+1; loop1++)
YunTu<<setw(8)<<Temp0[loop*(Tdiv1+divProt*2+1)+loop1]<<" ";
YunTu<<"\n";
}
}
}
}//******!!!!!!!!!时间推进已经结束!!!!!!!!!*****
//******************温度场计算结束时画云图准备***************
for(loop=0; loop<Tdiv2+divProt*2+1; loop++)
{
for(loop1=0; loop1<Tdiv1+divProt*2+1; loop1++)
YunTu1<<setw(8)<<Temp0[loop*(Tdiv1+divProt*2+1)+loop1]<<" ";
YunTu1<<"\n";
}
cout<<"\n"<<"*************成功**!!!!爆裂前温度场计算**成功!!!!!***************"<<endl;
cout<<"\n"<<"*************成功**!!!!爆裂前温度场计算**成功!!!!!***************"<<endl;
WResult<<"\n"<<"*************成功**!!!!爆裂前温度场计算**成功!!!!!***************"<<endl;
WResult<<"\n"<<"*************成功**!!!!爆裂前温度场计算**成功!!!!!***************"<<endl;
if(SpallTime==TotalTime)
{
cout<<"\n"<<" *******!!!!本次温度场计算不考虑高强混凝土的爆裂效!!!!!******"<<endl;
cout<<"\n"<<" *******!!!!本次温度场计算不考虑高强混凝土的爆裂效!!!!!******"<<endl;
WResult<<"\n"<<" *******!!!!本次温度场计算不考虑高强混凝土的爆裂效!!!!!******"<<endl;
WResult<<"\n"<<" *******!!!!本次温度场计算不考虑高强混凝土的爆裂效!!!!!******"<<endl;
}
else
{
cout<<"\n"<<" ******!!!!下面计算考虑高强混凝土的爆裂效应的温度场!!!!!*****"<<endl;
cout<<"\n"<<" ******!!!!下面计算考虑高强混凝土的爆裂效应的温度场!!!!!*****"<<endl;
WResult<<"\n"<<" ******!!!!下面计算考虑高强混凝土的爆裂效应的温度场!!!!!*****"<<endl;
WResult<<"\n"<<" ******!!!!下面计算考虑高强混凝土的爆裂效应的温度场!!!!!*****"<<endl;
}
//***********************************************************************
//************************爆裂后温度场的计算*****************************
//***********************************************************************
//****爆裂后划分单元数不变(Tdiv1+divProt*2,Tdiv2+divProt*2),单元编号和边界节点编号都没有变化,不必再
//****重新带宽和对角地址的计算********
//****温度场结构输出难度加大,爆裂后的温度场节点已经不是爆裂前的节点,要换算
while(Time<TotalTime)
{ //为爆裂后温度场的计算做准备
Time=Time+TotalTime-SpallTime;
for(loop=0; loop<TNode_num; loop++)
{
Temp[loop] =0; //重新清空温度数组
Temp1[loop]=0;
}
if(Fire0==0 && Fire1==1 && Fire2==0 && Fire3==0)//一面受火 右
{
AfSectwid=Sectwid-Spall_thick;
AfSectlen=Sectlen;
}
else if(Fire0==0 && Fire1==1 && Fire2==1 && Fire3==0)//相邻两面受火 右-上
{
AfSectwid=Sectwid-Spall_thick;
AfSectlen=Sectlen-Spall_thick;
}
else if(Fire0==1 && Fire1==1 && Fire2==1 && Fire3==0)//三面受火 右-上-下
{
AfSectwid=Sectwid-Spall_thick;
AfSectlen=Sectlen-Spall_thick*2;
}
else if(Fire0==1 && Fire1==1 && Fire2==1 && Fire3==1)//四面受火 右-上-下-左
{
AfSectwid=Sectwid-Spall_thick*2;
AfSectlen=Sectlen-Spall_thick*2;
}
else if(Fire0==0 && Fire1==1 && Fire2==0 && Fire3==1)//相对受火 右-左
{
AfSectwid=Sectwid-Spall_thick*2;
}
else
{
cout<<"受火形式(Fire0,Fire1,Fire2,Fire3)输入有错,请检查"<<endl;
exit(0);
}
for(loop=0; loop<(Tdiv1+divProt*2+1)*(Tdiv2+divProt*2+1); loop++)
{
PBG[loop]=0.0; B[loop]=0.0;
PBG0[loop]=0.0; PrG0[loop]=0.0;
PrG[loop]=0.0; PG0[loop]=0.0;
PG[loop]=0.0;
}
for(loop=0; loop<TotalLen; loop++)
{
KTG[loop]=0.0; KG[loop]=0.0;
KBG[loop]=0.0; CcG[loop]=0.0;
A[loop]=0.0; C[loop]=0.0;
CcG0[loop]=0.0; KG0[loop]=0.0;
KTG0[loop]=0.0; KBG0[loop]=0.0;
KrG[loop]=0.0; KrG0[loop]=0.0;
}
//***************提取爆裂后节点初始温度***********************
//计算爆裂后的相对节点坐标
for(loop=0; loop<Tdiv2+2*divProt+1; loop++) //行号
{
for(loop1=0; loop1<Tdiv1+2*divProt+1; loop1++)//
{
//节点z坐标
if(loop1<divProt+1)
AfNode_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1]=Prot_thick*loop1/divProt;
if(loop1>divProt && loop1<Tdiv1+divProt+1)
AfNode_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1]=Prot_thick+(AfSectwid-Prot_thick*2)*(loop1-divProt)/Tdiv1;
if(loop1>Tdiv1+divProt)
AfNode_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1]=AfSectwid-Prot_thick+(loop1-divProt-Tdiv1)*Prot_thick/divProt;
//节点y坐标
if(loop<divProt+1)
AfNode_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1+1]=Prot_thick*loop/divProt;
if(loop>divProt && loop<Tdiv2+divProt+1)
AfNode_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1+1]=Prot_thick+(AfSectlen-Prot_thick*2)*(loop-divProt)/Tdiv2;
if(loop>Tdiv2+divProt)
AfNode_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1+1]=AfSectlen-Prot_thick+(loop-divProt-Tdiv2)*Prot_thick/divProt;
}
}
for(loop=0; loop<TNode_num; loop++)
{
Temp0[loop]=Temp[loop]; Temp[loop]=0.0;
}
NumDieDai=0;
cout<<"\n"<<"****!!!爆裂后x方向划分单元总数!!!***: "<<Tdiv1+divProt*2<<endl;
cout<<"\n"<<"****!!!爆裂后y方向划分单元总数!!!***: "<<Tdiv2+divProt*2<<endl;
cout<<"\n"<<"****!!!爆裂后二维温度场单元总数!!!***: "<<TNode_num<<endl;
cout<<"\n"<<"****!!!爆裂后二维温度场节点总数!!!***: "<<TNode_num<<endl;
cout<<"\n"<<"****!1!爆裂后一维半带宽存储的元素个数为!!!***: "<<TotalLen<<"\n"<<endl;
WResult<<"\n"<<"****!!!爆裂后x方向划分单元总数!!!***: "<<Tdiv1+divProt*2<<endl;
WResult<<"\n"<<"****!!!爆裂后y方向划分单元总数!!!***: "<<Tdiv2+divProt*2<<endl;
WResult<<"\n"<<"****!!!爆裂后二维温度场单元总数!!!***: "<<TNode_num<<endl;
WResult<<"\n"<<"****!!!爆裂后二维温度场节点总数!!!***: "<<TNode_num<<endl;
WResult<<"\n"<<"****!1!爆裂后一维半带宽存储的元素个数为!!!***: "<<TotalLen<<"\n"<<endl;
cout<<"\n"<<"****************!!!!爆裂后温度场计算!!!!!****************"<<endl;
cout<<"\n"<<"****************!!!!爆裂后温度场计算!!!!!****************"<<endl;
WResult<<"\n"<<"****************!!!!爆裂后温度场计算!!!!!****************"<<endl;
WResult<<"\n"<<"****************!!!!爆裂后温度场计算!!!!!****************"<<endl;
}
Time=Time-(TotalTime-SpallTime);
//******************时间推进************主循环******************************
while(Time<TotalTime)
{
NumTime++;
if(VarTstep==1) //定时间步长
{
Time +=TimeStep;
TimeHis[loop]=Time; //记录时间步信息
if(Env_heatType==1) //炉内空气温度,即受火面的环境温度
Temper=iniT0+345.0*log10(8.0*Time/60.0+1.0);
if(Env_heatType==2)
Temper=iniT0+750.0*(1-exp(-3.79553*sqrt(Time/3600)))+170.41*sqrt(Time/3600);
if(Env_heatType==3)
CircumTemp>>Temper;
}
if(VarTstep==2) //采用变时间步长的求解策略
{
if(Env_heatType==1)
TimeStep=(8.0*Time/60.0+1.0)*(pow(10.0,DelTemp/345)-1)*60.0/8.0;
if(Env_heatType==2)
TimeStep=Cal_TimeStep(Time,DelTemp);
if(TimeStep>5*60.0){
TimeStep=5*60.0; //不大于5分钟
Time +=TimeStep;
TimeHis[NumTime-1]=Time; //记录时间步信息
if(Env_heatType==1)
Temper=iniT0+345.0*log10(8.0*Time/60.0+1.0);
if(Env_heatType==2)
Temper=iniT0+750.0*(1.0-exp(-3.79553*sqrt(Time/3600)))+170.41*sqrt(Time/3600.0);
}
if(TimeStep<=5*60.0){
Time +=TimeStep;
TimeHis[NumTime-1]=Time; //记录时间步信息
Temper +=DelTemp;
}
}
if(Fire0==1) //确定受火和不受火边界的坏境温度
AireTemp[0]=Temper;
else
{
if(Bound==1)
AireTemp[0]=BoundTemper;
else
AireTemp[0]=iniT0;
}
if(Fire1==1)
AireTemp[1]=Temper;
else
{
if(Bound==1)
AireTemp[1]=BoundTemper;
else
AireTemp[1]=iniT0;
}
if(Fire2==1)
AireTemp[2]=Temper;
else
{
if(Bound==1)
AireTemp[2]=BoundTemper;
else
AireTemp[2]=iniT0;
}
if(Fire3==1)
AireTemp[3]=Temper;
else
{
if(Bound==1)
AireTemp[3]=BoundTemper;
else
AireTemp[3]=iniT0;
}
if(Bound==1) //读入下一时间的边界温度
{
BoundTemp>>BoundTemper;
}
cout<<"***!爆裂后!***第"<<NumTime<<"个时间步 "<<" 时间:"<<Time/60<<"受火的环境温度:"<<Temper<<"摄氏度"<<endl;
WResult<<"***!爆裂后!***第"<<NumTime<<"个时间步 "<<" 时间:"<<Time/60<<"受火的环境温度:"<<Temper<<"摄氏度"<<endl;
if(Bord_heatType==2)//求综合换热系数
{
for(loop=0; loop<4; loop++)
{
Hc[loop]=0.0;
CalcuBt(Hc[loop],AireTemp[loop]);
}
}
//***叠加总刚度矩阵***
for(loop=0; loop<TConElem_num; loop++)
{
EleAverTemp[loop]=(Temp0[*(TCElemNode_Info+4*loop)]+Temp0[*(TCElemNode_Info+4*loop+1)]+Temp0[*(TCElemNode_Info+4*loop+2)]+Temp0[*(TCElemNode_Info+4*loop+3)])/4.0;
}
Form_KTG(EleAverTemp,CCondu_Type,KTG0);
Form_CcG(EleAverTemp,CSpeciheat_Type,CcG0);
if(Bord_heatType==1)
{
Form_KBG_PBG(Hc,AireTemp,KBG0,PBG0);
Form_krG_PrG(Defer,AireTemp,Temp0,Hr,KrG0,PrG0);
}
if(Bord_heatType==2)
{
Form_KBG_PBG(Hc,AireTemp,KBG0,PBG0);
}
for(loop=0; loop<TotalLen; loop++)
KG0[loop]=KBG0[loop]+KTG0[loop]+KrG0[loop];
for(loop=0; loop<TNode_num; loop++)
PG0[loop]=PBG0[loop]+PrG0[loop];
if(ComputeMethod==1) //Galerkin法
{
for(loop=0; loop<TotalLen; loop++)
{
A[loop]=CcG0[loop]/TimeStep+KG0[loop]*2/3;
C[loop]=CcG0[loop]/TimeStep-KG0[loop]/3;
}
Multiply_1D(B,C,Temp0);
for(loop=0; loop<TNode_num; loop++)
B[loop]=B[loop]+PG0[loop];
Gauss_1D(A,B,Temp1);
}
if(ComputeMethod==2) //两点向后差分法
{
for(loop=0; loop<TotalLen; loop++)
{
A[loop]=CcG0[loop]/TimeStep+KG0[loop];
}
Multiply_1D(B,KG0,Temp0);
for(loop=0; loop<TNode_num; loop++)
B[loop]=-B[loop]+PG0[loop];
Gauss_1D(A,B,Linshi); //Ax=B
for(loop=0; loop<TNode_num; loop++)
Temp1[loop]=Temp0[loop]+Linshi[loop];
}
NumDieDai=0; error=0.0;
//************************时*间*步*内*迭*代************************
do
{
error=0.0;
NumDieDai++;
for(loop=0; loop<(Tdiv1+divProt*2)*(Tdiv2+divProt*2); loop++)
{
EleAverTemp[loop]=(Temp1[*(TCElemNode_Info+4*loop)]+Temp1[*(TCElemNode_Info+4*loop+1)]+Temp1[*(TCElemNode_Info+4*loop+2)]+Temp1[*(TCElemNode_Info+4*loop+3)])/4.0;
}
//***叠加总刚度矩阵***
Form_KTG(EleAverTemp,CCondu_Type,KTG);
Form_CcG(EleAverTemp,CSpeciheat_Type,CcG);
if(Bord_heatType==1)
{
Form_KBG_PBG(Hc,AireTemp,KBG,PBG);
Form_krG_PrG(Defer,AireTemp,Temp1,Hr,KrG,PrG);
}
if(Bord_heatType==2)
{
Form_KBG_PBG(Hc,AireTemp,KBG,PBG);
}
for(loop=0; loop<TotalLen; loop++)
KG[loop]=KBG[loop]+KTG[loop]+KrG[loop];
for(loop=0; loop<(Tdiv1+divProt*2+1)*(Tdiv2+divProt*2+1); loop++)
PG[loop]=PBG[loop]+PrG[loop];
if(ComputeMethod==1) //Galerkin法
{
for(loop=0; loop<TotalLen; loop++)
{
A[loop]=(CcG0[loop]/3+CcG[loop]*2/3)/TimeStep+KG0[loop]/6+KG[loop]/2;
C[loop]=(CcG0[loop]/3+CcG[loop]*2/3)/TimeStep-KG0[loop]/6-KG[loop]/6;
}
Multiply_1D(B,C,Temp0);
for(loop=0; loop<TNode_num; loop++)
B[loop] +=PG0[loop]/3+PG[loop]*2/3;
Gauss_1D(A,B,Temp);
}
if(ComputeMethod==2) //两点向后差分法
{
for(loop=0; loop<TotalLen; loop++)
{
A[loop]=CcG[loop]/TimeStep+KG[loop];
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -