📄 tempf3d.h
字号:
if(*(TBorderCon+loop*(TElem_node+3)+10)!=4 && *(TBorderCon+loop*(TElem_node+3)+10)!=5)
{
Bord6_b=fabs(*(Node_Coor+3*Enode3+2)-*(Node_Coor+3*Enode1+2))/2.0;
Bord6_a=fabs(*(Node_Coor+3*Enode2+1)-*(Node_Coor+3*Enode1+1))/2.0;
if(fabs(*(Node_Coor+3*Enode2)-*(Node_Coor+3*Enode1))/2.0>Bord6_a)
{
Bord6_a=fabs(*(Node_Coor+3*Enode2)-*(Node_Coor+3*Enode1))/2.0;
}
}
//受火面上单元温度取4个节点的平均温度,且将温度化为绝对温度
TempAverage=(*(Temp1+Enode1)+*(Temp1+Enode2)+*(Temp1+Enode3)+*(Temp1+Enode4))/4.0+273.15;
if(*(TBorderCon+loop*(TElem_node+3)+9)==2) //常温边界条件的处理
{
ChangHeatCT=hc_Nort+hr_Nort*Ste_Bol*(TempAverage*TempAverage+(iniT0+273.0)*(iniT0+273.0))*(TempAverage+(iniT0+273.0));
}
if(*(TBorderCon+loop*(TElem_node+3)+9)==3) //受火边界条件的处理;过镇海书.P93,(6-38b),综述P8,式(7)
{
ChangHeatCT=hr_Heat*Ste_Bol*(TempAverage*TempAverage+(Envtemper+273.15)*Defer*(Envtemper+273.15)*Defer)*(TempAverage+(Envtemper+273.15)*Defer);
}
//单元传热矩阵有4×4个元素;只有4个代表性的元素
CKBe[0][0]=4.0*Bord6_a*Bord6_b*ChangHeatCT/9.0;
CKBe[0][1]=4.0*Bord6_a*Bord6_b*ChangHeatCT/18.0;
CKBe[0][2]=4.0*Bord6_a*Bord6_b*ChangHeatCT/18.0;
CKBe[0][3]=4.0*Bord6_a*Bord6_b*ChangHeatCT/36.0;
//第二行
CKBe[1][0]=CKBe[0][1];
CKBe[1][1]=CKBe[0][0];
CKBe[1][2]=CKBe[0][3];
CKBe[1][3]=CKBe[0][1];
//第二行
CKBe[2][0]=CKBe[0][2];
CKBe[2][1]=CKBe[1][2];
CKBe[2][2]=CKBe[0][0];
CKBe[2][3]=CKBe[0][1];
//第二行
CKBe[3][0]=CKBe[0][3];
CKBe[3][1]=CKBe[1][3];
CKBe[3][2]=CKBe[2][3];
CKBe[3][3]=CKBe[0][0];
//六面体单元的传热荷载向量
if(*(TBorderCon+loop*(TElem_node+3)+9)==2)//常温边界条件的处理
{
PBe[0]=Bord6_a*Bord6_b*ChangHeatCT*iniT0;
PBe[1]=Bord6_a*Bord6_b*ChangHeatCT*iniT0;
PBe[2]=Bord6_a*Bord6_b*ChangHeatCT*iniT0;
PBe[3]=Bord6_a*Bord6_b*ChangHeatCT*iniT0;
}
if(*(TBorderCon+loop*(TElem_node+3)+9)==3) //受火边界条件的处理;过镇海书.P93,(6-38b),综述P8,式(7)
{
PBe[0]=Bord6_a*Bord6_b*hc_Heat*Envtemper+Bord6_a*Bord6_b*ChangHeatCT*Envtemper*Defer; //?
PBe[1]=PBe[0];
PBe[2]=PBe[0];
PBe[3]=PBe[0];
}
}
for(loop1=0; loop1<TElem_node; loop1++) //组装边界刚度矩阵和传热荷载向量
{
iRow=*(TCElemNode_Info+TElem_node*loop+loop1);
for(loop2=0; loop2<=loop1; loop2++)
{
iCol=*(TCElemNode_Info+TElem_node*loop+loop2);
if(iRow>=iCol)
ibuf=DiagAddress[iRow]-(iRow-iCol);
else
ibuf=DiagAddress[iCol]-(iCol-iRow);
KBG[ibuf]=KBG[ibuf]+CKBe[loop1][loop2];
}
PBG[iRow]=PBG[iRow]+PBe[loop1];
}
}//!!!边界单元的和传热刚度矩阵和传热荷载向量组装完毕!!!!!
if(SteelExit==1) //三维温度场分析时考虑钢筋的影响；
{
;
}
for(loop=0; loop<TotalLen; loop++)
KG[loop]=KBG[loop]+KTG[loop];
//备份t0时刻第一次迭代的刚度矩阵,质量热容矩阵和热荷载向量
if(ComputeMethod==1 && IterCount==1)//1表示权函数按Galerkin取值时得到的隐格式
{
for(loop=0;loop<TotalLen;loop++)
{
*(KG0+loop)=*(KG+loop); //t0时刻的总刚度矩阵；
*(CcG0+loop)=*(CcG+loop); //t0时刻的质量热容矩阵
}
for(loop=0;loop<TNode_num;loop++)
{
*(PBG0+loop)=*(PBG+loop); //t0时刻的总热荷载向量
}
}
if(ComputeMethod==1) //1表示权函数按Galerkin取值时得到的隐格式
{
if(IterCount==1)
{
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]+PBG0[loop];
Gauss_1D(A,B,Temp); //高斯法解方程组Ax=B
}
if(IterCount!=1)
{
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] +=PBG0[loop]/3.0+PBG[loop]*2.0/3.0;
Gauss_1D(A,B,Temp); //高斯法解方程组Ax=B
}
}
if(ComputeMethod==2) //1表示采用两点后差分格式
{
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]+PBG0[loop];
Gauss_1D(A,B,Linshi); //高斯法解方程组Ax=B
for(loop=0; loop<TNode_num; loop++)
Temp[loop]=Temp0[loop]+Linshi[loop];
}
cout<<" 第 "<<NumTime<<" 次时间步,第 "<<IterCount<<" 子步计算完成\n"<<endl;
WResult<<" 第 "<<NumTime<<" 次时间步,第 "<<IterCount<<" 子步计算完成\n"<<endl;
if(IterCount!=1) //下面判断迭代是否收敛
{
IterStop=0; error=0.0;
for(loop=0;loop<TNode_num;loop++)
{
if(fabs(*(Temp+loop)-*(Temp1+loop))>FNMaxErr)
{
IterStop=1;
}
if(fabs(*(Temp1+loop)-*(Temp+loop))>error)
{
error=fabs(*(Temp1+loop)-*(Temp+loop));
}
}
}
cout<<" 第 "<<NumTime<<" 次时间步,第 "<<IterCount<<" 子步的最大收敛误差为："<<error<<",容许误差为："<<FNMaxErr<<endl;
WResult<<" 第 "<<NumTime<<" 次时间步,第 "<<IterCount<<" 子步的最大收敛误差为："<<error<<",容许误差为："<<FNMaxErr<<endl;
//超过最大允许迭代次数
if(IterCount>IterMax)
{
IterStop=0;
}
for(loop=0;loop<TNode_num;loop++) //计算t+delt时刻的温度
{
*(Temp1+loop)=*(Temp+loop);
}
}while(IterStop==1);
if(IterCount<=IterMax)
{
cout<<"********第 "<<NumTime<<" 次时间步计算成功收敛! ********\n"<<endl;
WResult<<"********第 "<<NumTime<<" 次时间步计算成功收敛! ********\n"<<endl;
}
if(IterCount>IterMax)
{
cout<<"********第 "<<NumTime<<" 次时间步超过最大迭代次数 "<<IterMax<<" ********\n"<<endl;
WResult<<"********第 "<<NumTime<<" 次时间步超过最大迭代次数 "<<IterMax<<" ********\n"<<endl;
}
for(loop=0;loop<TNode_num;loop++) //计算t+delt时刻的温度
{
if(TemperAlsub==1)//子空间过滤
{
if(*(Temp1+loop)>=*(Temp0+loop))//最小值过滤
{
*(Temp0+loop)=*(Temp1+loop);
}
else
*(Temp1+loop)=*(Temp0+loop);
if(*(Temp1+loop)>=Envtemper)//最大值过滤
{
*(Temp0+loop)=Envtemper;
*(Temp1+loop)=Envtemper;
}
}
}
//********************输出结果***节点温度****************
OutResults<<"\n"<<"时间="<<Time/60.0<<"温度="<<Envtemper<<"\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);
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";
}
}
}
}
}
//***********钢筋的体积比热随温度的变化**************//
double CalcuSSpecilHeat(int SSpecilHeat_Type,double Temp)
{
double SSpecilHeat;
if(SSpecilHeat_Type==1) //按Lie的建议公式取值
{
if(Temp<=650.0)
SSpecilHeat=(0.004*Temp+3.3)*1.0e6;
else if(Temp>=650.0 && Temp<=725.0)
SSpecilHeat=(0.0680*Temp-38.3)*1.0e6;
else if(Temp>=725.0 && Temp<=800.0)
SSpecilHeat=(-0.0860*Temp+73.35)*1.0e6;
else if(Temp>=800.0)
SSpecilHeat=4.55*1.0e6;
}
else if(SSpecilHeat_Type==2) //2表示按同济大学陆州导建议公式取
{
if(Temp<=500.0)
SSpecilHeat=7850.0*(4.77*Temp/10000+0.48)*1.0e3;
else if(Temp>=5⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -