📄 tempf2d.h
字号:
BandWidth[*(TCElemNode_Info+4*loop+2)]=*(TCElemNode_Info+4*loop+2)-ibuf+1;
if(BandWidth[*(TCElemNode_Info+4*loop+3)]<=(*(TCElemNode_Info+4*loop+3)-ibuf))
BandWidth[*(TCElemNode_Info+4*loop+3)]=*(TCElemNode_Info+4*loop+3)-ibuf+1;
}
//对角地址的计算按一维压缩数组贮存
for(loop=0; loop<TNode_num; loop++)
DiagAddress[loop]=1;
for(loop=1; loop<TNode_num; loop++)
{
DiagAddress[loop]=0;
for(loop1=0; loop1<=loop; loop1++){
DiagAddress[loop]+=BandWidth[loop1];
}
}
for(loop=0; loop<TNode_num; loop++)
DiagAddress[loop]--;
TotalLen=DiagAddress[TNode_num-1]+1;
cout<<"***!!!x方向划分单元总数!!!***: "<<Tdiv1+divProt*2<<endl;
cout<<"***!!!y方向划分单元总数!!!***: "<<Tdiv2+divProt*2<<endl;
cout<<"***!!!二维温度场单元总数!!!***: "<<TConElem_num<<endl;
cout<<"***!!!二维温度场节点总数!!!***: "<<TNode_num<<endl;
cout<<"***!!!变带宽下三角一维存储总刚度矩阵存储的总刚元素的个数为!!!***: "<<TotalLen<<endl;
WResult<<"***!!!x方向划分单元总数!!!***: "<<Tdiv1+divProt*2<<endl;
WResult<<"***!!!y方向划分单元总数!!!***: "<<Tdiv2+divProt*2<<endl;
WResult<<"***!!!二维温度场单元总数!!!***: "<<TConElem_num<<endl;
WResult<<"***!!!二维温度场节点总数!!!***: "<<TNode_num<<endl;
WResult<<"***!!!变带宽下三角一维存储总刚度矩阵存储的总刚元素的个数为!!!***: "<<TotalLen<<endl;
//计算相对节点坐标
for(loop=0; loop<Tdiv2+2*divProt+1; loop++) //行号
{
for(loop1=0; loop1<Tdiv1+2*divProt+1; loop1++)//
{
//节点z坐标
if(loop1<divProt+1)
Node_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1]=Prot_thick*loop1/divProt;
if(loop1>divProt && loop1<Tdiv1+divProt+1)
Node_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1]=Prot_thick+(Sectwid-Prot_thick*2)*(loop1-divProt)/Tdiv1;
if(loop1>Tdiv1+divProt)
Node_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1]=Sectwid-Prot_thick+(loop1-divProt-Tdiv1)*Prot_thick/divProt;
//节点y坐标
if(loop<divProt+1)
Node_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1+1]=Prot_thick*loop/divProt;
if(loop>divProt && loop<Tdiv2+divProt+1)
Node_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1+1]=Prot_thick+(Sectlen-Prot_thick*2)*(loop-divProt)/Tdiv2;
if(loop>Tdiv2+divProt)
Node_Coor[2*loop*(divProt*2+Tdiv1+1)+2*loop1+1]=Sectlen-Prot_thick+(loop-divProt-Tdiv2)*Prot_thick/divProt;
}
}
//给各刚度矩阵分配内存
KTG=(double *)calloc(TotalLen,8);
KG=(double *)calloc(TotalLen,8);
KBG=(double *)calloc(TotalLen,8);
KrG=(double *)calloc(TotalLen,8);
CcG=(double *)calloc(TotalLen,8);
PBG=(double *)calloc(TNode_num,8);
PrG=(double *)calloc(TNode_num,8);
PG =(double *)calloc(TNode_num,8);
A=(double *)calloc(TotalLen,8);
B=(double *)calloc(TNode_num,8);
C=(double *)calloc(TotalLen,8);
CcG0=(double *)calloc(TotalLen,8);
KG0 =(double *)calloc(TotalLen,8);
KTG0=(double *)calloc(TotalLen,8);
KBG0=(double *)calloc(TotalLen,8);
KrG0=(double *)calloc(TotalLen,8);
PBG0=(double *)calloc(TNode_num,8);
PrG0=(double *)calloc(TNode_num,8);
PG0=(double *)calloc(TNode_num,8);
for(loop=0; loop<Num_TiQu; loop++)
Temp_TiQu[loop]=0.0;
for(loop=0; loop<Points; loop++)
PointsTemp[loop]=0.0;
for(loop=0; loop<4; loop++)
AireTemp[loop]=iniT0;
if(Bound==1)
{
BoundTemper=iniT0;
}
TimeHis=(double *)calloc(Step,8);
Temper=iniT0;
Step=0;
//************************爆裂前温度场的计算****************************
cout<<"\n****************!!!!爆裂前温度场计算!!!!!****************\n"<<endl;
WResult<<"\n****************!!!!爆裂前温度场计算!!!!!****************\n"<<endl;
//**********************时间推进******主循环****************************
while(Time < SpallTime)
{
NumTime++;
if(VarTstep==1) //采用定时间步长的求解策略
{
Time +=TimeStep;
TimeHis[NumTime-1]=Time; //记录时间步信息
if(Env_heatType==1) //炉内空气温度,即受火面的环境温度
Temper=iniT0+345.0*log10(8.0*Time/60.0+1.0);
else if(Env_heatType==2)
Temper=iniT0+750.0*(1.0-exp(-3.79553*sqrt(Time/3600.0)))+170.41*sqrt(Time/3600.0);
else if(Env_heatType==3)
CircumTemp>>Temper;
else
{
cout<<"没有这样的升温曲线(Env_heatType)"<<endl;
exit(0);
}
}
if(VarTstep==2) //采用变时间步长的求解策略
{
if(Env_heatType==1)
TimeStep=(8.0*Time/60.0+1.0)*(pow(10.0,DelTemp/345.0)-1)*60.0/8.0;
if(Env_heatType==2)
TimeStep=Cal_TimeStep(Time,DelTemp);
if(TimeStep>5.0*60.0){
TimeStep=5.0*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.0)))+170.41*sqrt(Time/3600.0);
}
if(TimeStep<=5.0*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.0<<" 受火的环境温度:"<<Temper<<"摄氏度"<<endl;
WResult<<"第"<<NumTime<<"个时间步"<<" 时间:"<<Time/60.0<<" 受火的环境温度:"<<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);
}
else if(Bord_heatType==2)
{
Form_KBG_PBG(Hc,AireTemp,KBG0,PBG0);
}
else
{
cout<<"换热方法选择有错(Bord_heatType)!!!"<<endl;
cout<<"Bord_heatType= "<<Bord_heatType<<endl;
exit(0);
}
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.0/3.0;
C[loop]=CcG0[loop]/TimeStep-KG0[loop]/3.0;
}
Multiply_1D(B,C,Temp0);
for(loop=0; loop<TNode_num; loop++)
B[loop]=B[loop]+PG0[loop];
Gauss_1D(A,B,Temp1); //Ax=B
}
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<TConElem_num; 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<TNode_num; 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.0/3.0)/TimeStep+KG0[loop]/6+KG[loop]/2.0;
C[loop]=(CcG0[loop]/3+CcG[loop]*2.0/3.0)/TimeStep-KG0[loop]/6-KG[loop]/6.0;
}
Multiply_1D(B,C,Temp0);
for(loop=0; loop<TNode_num; loop++)
B[loop] +=PG0[loop]/3.0+PG[loop]*2.0/3.0;
Gauss_1D(A,B,Temp);
}
if(ComputeMethod==2) //两点向后差分法
{
for(loop=0; loop<TotalLen; loop++)
{
A[loop]=CcG[loop]/TimeStep+KG[loop];
}
Multiply_1D(B,KG,Temp0);
for(loop=0; loop<TNode_num; loop++)
B[loop]=-B[loop]+PG[loop];
Gauss_1D(A,B,Linshi); //Ax=B
for(loop=0; loop<TNode_num; loop++)
Temp[loop]=Temp0[loop]+Linshi[loop];
}
for(loop=0; loop<TNode_num; loop++) //选出前后两次迭代节点最大温差
{
if(error<fabs(Temp[loop]-Temp1[loop]))
error=fabs(Temp[loop]-Temp1[loop]);
if(Temp[loop]==Temp1[loop]) //前后前后两次迭代基本上相等
break;
}
for(loop=0; loop<TNode_num; loop++)
Temp1[loop]=Temp[loop];
for(loop=0; loop<TNode_num; loop++) //子空间法去掉不合理解
{
if(Temp1[loop]<iniT0) //不低于初始温度
Temp1[loop]=iniT0;
if(Temp1[loop]>Temper) //不高于介质温度
Temp1[loop]=Temper;
}
cout<<" 第 "<<NumTime<<"时间步 "<<"第"<<NumDieDai<<"迭代步 "<<"误差error为:"<<setiosflags(ios::fixed)<<setprecision(5)<<error<<" 容许误差为:"<<FNMaxErr<<endl;
WResult<<" 第 "<<NumTime<<"时间步 "<<"第"<<NumDieDai<<"迭代步 "<<"误差error为:"<<setiosflags(ios::fixed)<<setprecision(5)<<error<<" 容许误差为:"<<FNMaxErr<<endl;
if(NumDieDai==IterMax)
{
cout<<"!!!!!!!!!迭代了100次,不收敛!!!!!!!!!"<<endl;
WResult<<"!!!!!!!!!迭代了100次,不收敛!!!!!!!!!"<<endl;
}
if(error<FNMaxErr && NumDieDai<IterMax)
{
cout<<"!!!!!第"<<NumTime<<"时间步成功收敛!!!!!!!!!"<<endl;
WResult<<"!!!!!第"<<NumTime<<"时间步成功收敛!!!!!!!!!"<<endl;
}
}while(error>FNMaxErr && NumDieDai<IterMax); //控制迭代循坏
for(loop=0; loop<TNode_num; loop++)
Temp0[loop]=Temp1[loop];
//********************输出结果***节点温度****************
// OutResults<<"\n"<<"时间="<<Time/60<<"温度="<<Temper<<"\n";
for(loop=0; loop<TNode_num; loop++)
{
// ibuf=loop+1;
// OutResults<<setiosflags(ios::fixed);
// OutResults<<setprecision(3)<<Temp0[loop]<<" ";
OutResults<<Temp0[loop]<<endl;
// if((ibuf%(Tdiv1+divProt*2+1))==0) OutResults<<"\n";
}
//**********************输出***提取点的温度*********************
OutTemp<<setiosflags(ios::fixed)<<setprecision(5)<<Time/60<<",";
Cal_Temp_TiQu(Num_TiQu,CorX_TiQu,CorY_TiQu,Temp0,Temp_TiQu);
for(loop=0; loop<Num_TiQu; loop++)
{
OutTemp<<setiosflags(ios::fixed);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -