📄 tempf2d.h
字号:
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++)
Temp0[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(4)<<error<<" 容许误差为:"<<FNMaxErr<<endl;
WResult<<" 第 "<<NumTime<<"时间步 "<<"第"<<NumDieDai<<"迭代步 "<<"误差error为:"<<setiosflags(ios::fixed)<<setprecision(4)<<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];
//******提取爆裂前节点温度,根据前后坐标的对应关系********
for(loop=0; loop<TNode_num; loop++)
{
Spall_Ind=0;
Cal_Spall_Ind(Node_Coor[loop*2],Node_Coor[loop*2+1],Spall_Ind);
if(Spall_Ind==1)//判断提取点是否已经爆裂
{ //已经爆裂处于明火中的单元
if(Env_heatType==1)
TempLin[loop]=iniT0+345.0*log10(8.0*Time/60.0+1.0);
if(Env_heatType==2)
TempLin[loop]=iniT0+750.0*(1-exp(-3.79553*sqrt(Time/3600)))+170.41*sqrt(Time/3600);
if(Env_heatType==3)
TempLin[loop]=Temper;
}
else //尚未爆裂的单元
{
if(Fire0==1) //爆裂后相对原点的位置 把原坐标换成后坐标
Node_Coor[loop*2+1]-=Spall_thick;
if(Fire3==1)
Node_Coor[loop*2]-=Spall_thick;
Cal_Temp_TiQu(1,Node_Coor+loop*2,Node_Coor+loop*2+1,Temp0,TempLin+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<<TempLin[loop]<<endl;
if((ibuf%(Tdiv1+divProt*2+1))==0) OutResults<<"\n";
}
//********************提取提取点的温度****************
for(loop=0; loop<Num_TiQu; loop++)
{
Spall_Ind=0;
Cal_Spall_Ind(CorX_TiQu[loop],CorY_TiQu[loop],Spall_Ind);
if(Spall_Ind==1)//判断提取点是否已经爆裂
{ //已经爆裂处于明火中的单元
if(Env_heatType==1)
Temp_TiQu[loop]=iniT0+345.0*log10(8.0*Time/60.0+1.0);
if(Env_heatType==2)
Temp_TiQu[loop]=iniT0+750.0*(1-exp(-3.79553*sqrt(Time/3600)))+170.41*sqrt(Time/3600);
if(Env_heatType==3)
Temp_TiQu[loop]=Temper;
}
else //尚未爆裂的单元
{
if(Fire0==1) //爆裂后相对原点的位置 把原坐标换成后坐标
CorY_TiQu[loop]-=Spall_thick;
if(Fire3==1)
CorX_TiQu[loop]-=Spall_thick;
Cal_Temp_TiQu(1,CorX_TiQu+loop,CorY_TiQu+loop,Temp0,Temp_TiQu+loop);
if(Fire0==1) //为下一步计算做准备
CorY_TiQu[loop]+=Spall_thick;
if(Fire3==1)
CorX_TiQu[loop]+=Spall_thick;
}
}
//**********************输出***提取点的温度*********************
OutTemp<<setiosflags(ios::fixed)<<setprecision(3)<<Time/60<<",";
for(loop=0; loop<Num_TiQu; loop++)
{
OutTemp<<setiosflags(ios::fixed);
OutTemp<<setprecision(3)<<Temp_TiQu[loop]<<",";
}
OutTemp<<"\n";
}//******时间推进已经结束!!!!!*****
//*********************画云图准备*****************
for(loop=0; loop<Tdiv2+divProt*2+1; loop++)
{
for(loop1=0; loop1<Tdiv1+divProt*2+1; loop1++)
YunTu2<<setw(8)<<Temp0[loop*(Tdiv1+divProt*2+1)+loop1]<<" ";
YunTu2<<"\n";
}
//为结构计算做准备
Step=NumTime;
for(loop=0; loop<Step; loop++)
{
TimeOut<<TimeHis[loop]<<endl;
}
cout<<"\n"<<" ***********!!!!本次温度场计算全部结束!!!!!**********"<<endl;
cout<<"\n"<<" ***********!!!!本次温度场计算全部结束!!!!!**********"<<endl;
WResult<<"\n"<<" ***********!!!!本次温度场计算全部结束!!!!!**********"<<endl;
WResult<<"\n"<<" ***********!!!!本次温度场计算全部结束!!!!!**********"<<endl;
}
//******************将单元刚度矩阵**KTe**合成整体刚度矩阵***************//
void Form_KTG(double *EleAverTemp,int CCondu_Type,double *KTG)
{
int loop,loop1,loop2;
int iRow,iCol,ibuf;
double KTe[4][4];
double a=0.0,b=0.0;
double Kc=0.0;
for(loop=0; loop<TotalLen; loop++)
KTG[loop]=0.0;
for(loop=0; loop<TConElem_num; loop++)
{
a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;
Kc=CalcuKc(EleAverTemp[loop],CCondu_Type);
KTe[0][0]=KTe[1][1]=KTe[2][2]=KTe[3][3]=(2.0*b/a+2.0*a/b)*Kc/6.0;
KTe[0][1]=KTe[1][0]=KTe[3][2]=KTe[2][3]=(-2.0*b/a+a/b)*Kc/6.0;
KTe[0][2]=KTe[2][0]=KTe[1][3]=KTe[3][1]=(b/a-2.0*a/b)*Kc/6.0;
KTe[3][0]=KTe[2][1]=KTe[1][2]=KTe[0][3]=(-b/a-a/b)*Kc/6.0;
for(loop1=0; loop1<4; loop1++)
{
iRow=*(TCElemNode_Info+4*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+4*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
KTG[ibuf] +=KTe[loop1][loop2];
}
}
}
}
//******************将单元刚度矩阵**Cce**合成整体刚度矩阵***************//
void Form_CcG(double *EleAverTemp,int CSpeciheat_Type,double *CcG)
{
int loop,loop1,loop2;
int iRow,iCol,ibuf;
double SpecialHeat=0.0;
double Cce[4][4];
double a=0.0,b=0.0;
for(loop=0; loop<TotalLen; loop++)
CcG[loop]=0.0;
for(loop=0; loop<TConElem_num; loop++)
{
SpecialHeat=CalcuCSpecialHeat(EleAverTemp[loop],CSpeciheat_Type);
a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;
Cce[0][0]=Cce[1][1]=Cce[2][2]=Cce[3][3]=4.0*a*b*SpecialHeat/9.0;
Cce[1][0]=Cce[2][0]=Cce[3][1]=Cce[3][2]=2.0*a*b*SpecialHeat/9.0;
Cce[0][1]=Cce[0][2]=Cce[1][3]=Cce[2][3]=2.0*a*b*SpecialHeat/9.0;
Cce[2][1]=Cce[3][0]=Cce[0][3]=Cce[1][2]=1.0*a*b*SpecialHeat/9.0;
for(loop1=0; loop1<4; loop1++)
{
iRow=*(TCElemNode_Info+4*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+4*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
CcG[ibuf] +=Cce[loop1][loop2];
}
}
}
}
//***************将对流边界KBe,PBe矩阵合成整体刚度矩阵*********************//
void Form_KBG_PBG(double *Hc,double *AireTemp,double *KBG,double *PBG)
{
int loop,loop1,loop2;
int iRow=0,iCol=0,ibuf=0;
double a=0.0,b=0.0;
double KBe[4][4]; //单元传热矩阵
double PBe[4]; //单元传热荷载向量
for(loop=0; loop<4; loop++) //初始化
PBe[loop]=0.0;
for(loop=0; loop<4; loop++)
for(loop1=0; loop1<4; loop1++)
KBe[loop][loop1]=0.0;
for(loop=0; loop<TotalLen; loop++)
KBG[loop]=0.0;
for(loop=0; loop<TNode_num; loop++)
PBG[loop]=0.0;
//底面受火
for(loop=0; loop<Tdiv1+divProt*2; loop++)
{
a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;
KBe[0][0]=KBe[1][1]=2.0*Hc[0]*a/3.0; //
KBe[0][1]=KBe[1][0]=1.0*Hc[0]*a/3.0;
KBe[0][2]=KBe[0][3]=KBe[1][2]=KBe[1][3]=0.0;
KBe[2][0]=KBe[2][1]=KBe[2][2]=KBe[2][3]=0.0;
KBe[3][0]=KBe[3][1]=KBe[3][2]=KBe[3][3]=0.0;
PBe[0]=PBe[1]=Hc[0]*AireTemp[0]*a; //
PBe[2]=PBe[3]=0.0;
for(loop1=0; loop1<4; loop1++)
{
iRow=*(TCElemNode_Info+4*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+4*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
KBG[ibuf] +=KBe[loop1][loop2];
}
PBG[iRow]+=PBe[loop1];
}
}
//计算试验柱的第二个受火面 叠加柱子右边受火时的刚度矩阵
for(loop=Tdiv1+divProt*2-1; loop<(Tdiv1+divProt*2)*(Tdiv2+divProt*2); )
{
a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;
KBe[1][1]=KBe[3][3]=2.0*Hc[1]*b/3.0; //
KBe[1][3]=KBe[3][1]=1.0*Hc[1]*b/3.0;
KBe[0][0]=KBe[0][1]=KBe[0][2]=KBe[0][3]=0.0;
KBe[2][0]=KBe[2][1]=KBe[2][2]=KBe[2][3]=0.0;
KBe[3][0]=KBe[3][2]=KBe[1][0]=KBe[1][2]=0.0;
PBe[1]=PBe[3]=Hc[1]*AireTemp[1]*b; //
PBe[0]=PBe[2]=0.0;
for(loop1=0; loop1<4; loop1++)
{
iRow=*(TCElemNode_Info+4*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+4*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
KBG[ibuf] +=KBe[loop1][loop2];
}
PBG[iRow] +=PBe[loop1];
}
loop +=Tdiv1+divProt*2;
}
//计算试验柱的第三个受火面 叠加柱子上边受火时的刚度矩阵
for(loop=(Tdiv1+divProt*2)*(Tdiv2+divProt*2-1); loop<(Tdiv1+divProt*2)*(Tdiv2+divProt*2); loop++)
{
a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;
KBe[2][2]=KBe[3][3]=2.0*Hc[2]*a/3.0;
KBe[2][3]=KBe[3][2]=1.0*Hc[2]*a/3.0;
KBe[0][0]=KBe[0][1]=KBe[0][2]=KBe[0][3]=0.0;
KBe[1][0]=KBe[1][1]=KBe[1][2]=KBe[1][3]=0.0;
KBe[2][0]=KBe[2][1]=KBe[3][0]=KBe[3][1]=0.0;
PBe[2]=PBe[3]=Hc[2]*AireTemp[2]*a; //
PBe[0]=PBe[1]=0.0;
for(loop1=0; loop1<4; loop1++)
{
iRow=*(TCElemNode_Info+4*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+4*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
KBG[ibuf] +=KBe[loop1][loop2];
}
PBG[iRow] +=PBe[loop1];
}
}
//计算试验柱的第四个受火面 叠加柱子左边受火时的刚度矩阵
for(loop=0; loop<=(Tdiv1+divProt*2)*(Tdiv2+divProt*2-1); )
{
a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;
KBe[0][0]=KBe[2][2]=2.0*Hc[3]*b/3.0;
KBe[0][2]=KBe[2][0]=1.0*Hc[3]*b/3.0;
KBe[0][1]=KBe[0][3]=KBe[2][1]=KBe[2][3]=0.0;
KBe[1][0]=KBe[1][1]=KBe[1][2]=KBe[1][3]=0.0;
KBe[3][0]=KBe[3][1]=KBe[3][2]=KBe[3][3]=0.0;
PBe[0]=PBe[2]=Hc[3]*AireTemp[3]*b; //
PBe[1]=PBe[3]=0.0;
for(loop1=0; loop1<4; loop1++)
{
iRow=*(TCElemNode_Info+4*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+4*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
KBG[ibuf] +=KBe[loop1][loop2];
}
PBG[iRow] +=PBe[loop1];
}
loop +=Tdiv1+divProt*2;
}
}
//****************叠加辐射边界条件到[KrG]矩阵************//
void Form_krG_PrG(double Defer,double *AireTemp,double *NTemp,double *Hr,double *KrG,double *PrG)
{
int loop,loop1,loop2;
int iRow=0,iCol=0,ibuf=0;
double Ratio=0.0,BJTemp=0.0,AireT[4];
double Kre[4][4],Pre[4]; //单元辐射传热矩阵,辐射传热荷载向量
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -