📄 tempf3d.h
字号:
}
if(BorderFz2>=Colu_h1+Colu_h2+Colu_h3+Colu_h4 && BorderFz2<Colu_high)
{
BFendceng=zdiv1+zdiv22+zdiv33+zdiv44+(int)((BorderFz2-Colu_h1-Colu_h2-Colu_h3-Colu_h4)/(Colu_h5/zdiv5)+0.00001);
}
//考虑2种边界条件,1-表示第一类边界条件温度为给定值,2-空气保持常温的边界,3-受火边界条件
Bord_num=0;
for(loop=0;loop<2*divProt+Tdiv2;loop++) //柱底面和柱顶面的边界的编码
{
for(loop1=0;loop1<2*divProt+Tdiv1;loop1++)//柱底面的边界
{
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+0)=loop*(2*divProt+Tdiv1)+loop1;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+1)=loop*(2*divProt+Tdiv1+1)+loop1;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+2)=loop*(2*divProt+Tdiv1+1)+loop1+1;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+3)=(loop+1)*(2*divProt+Tdiv1+1)+loop1;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+4)=(loop+1)*(2*divProt+Tdiv1+1)+loop1+1;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+5)=0;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+6)=1;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+7)=2;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+8)=3;
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+9)=*(BodreArea+2*4+1);
*(TBorderCon+(loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+10)=*(BodreArea+2*4+0);
}
for(loop1=0;loop1<2*divProt+Tdiv1;loop1++)//柱顶面的边界
{
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+0)=loop*(2*divProt+Tdiv1)+loop1+TStoryElem_num*(zdiv-1);
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+1)=loop*(2*divProt+Tdiv1+1)+loop1+TStoryNode_num*zdiv;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+2)=loop*(2*divProt+Tdiv1+1)+loop1+1+TStoryNode_num*zdiv;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+3)=(loop+1)*(2*divProt+Tdiv1+1)+loop1+TStoryNode_num*zdiv;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+4)=(loop+1)*(2*divProt+Tdiv1+1)+loop1+1+TStoryNode_num*zdiv;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+5)=4;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+6)=5;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+7)=6;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+8)=7;
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+9)=*(BodreArea+2*5+1);
*(TBorderCon+(TStoryElem_num+loop*(2*divProt+Tdiv1)+loop1)*(TElem_node+3)+10)=*(BodreArea+2*5+0);
}
}
Bord_num=TStoryElem_num*2; //柱底和柱顶的边界单元数
for(loop=0;loop<zdiv;loop++) //柱的边界单元的编码
{
for(loop1=0;loop1<2*divProt+Tdiv1;loop1++) //0面 侧面的底边面
{
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+0)=loop1+loop*TStoryElem_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+1)=loop1+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+2)=loop1+1+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+3)=loop1+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+4)=loop1+1+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+5)=0;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+6)=1;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+7)=4;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+8)=5;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+9)=2;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+10)=*(BodreArea+2*0+0);
}
Bord_num=Bord_num+2*divProt+Tdiv1;
}
//并修正边界信息 !!!!将常温边界2修正成受火边界3!!!
if(*(BodreArea+2*0+1)==3) //修改柱0面的边界条件
{
for(loop=BFbegceng;loop<BFendceng;loop++)
{
for(loop1=0;loop1<2*divProt+Tdiv1;loop1++)
{
*(TBorderCon+(TStoryElem_num*2+(2*divProt+Tdiv1)*loop+loop1)*(TElem_node+3)+9)=3;
}
}
}
for(loop=0;loop<zdiv;loop++) //柱的边界单元的编码
{
for(loop1=0;loop1<2*divProt+Tdiv2;loop1++) //1面 侧面的右边面
{
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+0)=(loop1+1)*(2*divProt+Tdiv1)-1+loop*TStoryElem_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+1)=(loop1+1)*(2*divProt+Tdiv1+1)-1+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+2)=(loop1+2)*(2*divProt+Tdiv1+1)-1+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+3)=(loop1+1)*(2*divProt+Tdiv1+1)-1+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+4)=(loop1+2)*(2*divProt+Tdiv1+1)-1+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+5)=1;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+6)=3;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+7)=5;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+8)=7;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+9)=2;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+10)=*(BodreArea+2*1+0);
}
Bord_num=Bord_num+2*divProt+Tdiv2;
}
//并修正边界信息 !!!!将常温边界2修正成受火边界3!!!
if(*(BodreArea+2*1+1)==3) //修改柱1面的边界条件
{
for(loop=BFbegceng;loop<BFendceng;loop++)
{
for(loop1=0;loop1<2*divProt+Tdiv2;loop1++)
{
*(TBorderCon+(TStoryElem_num*2+(2*divProt+Tdiv1)*zdiv+(2*divProt+Tdiv2)*loop+loop1)*(TElem_node+3)+9)=3;
}
}
}
for(loop=0;loop<zdiv;loop++) //柱的边界单元的编码
{
for(loop1=0;loop1<2*divProt+Tdiv1;loop1++) //2面 侧面的上边面
{
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+0)=loop1+(2*divProt+Tdiv1)*(2*divProt+Tdiv2-1)+loop*TStoryElem_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+1)=loop1+(2*divProt+Tdiv2)*(2*divProt+Tdiv1+1)+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+2)=loop1+(2*divProt+Tdiv2)*(2*divProt+Tdiv1+1)+1+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+3)=loop1+(2*divProt+Tdiv2)*(2*divProt+Tdiv1+1)+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+4)=loop1+(2*divProt+Tdiv2)*(2*divProt+Tdiv1+1)+1+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+5)=2;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+6)=3;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+7)=6;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+8)=7;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+9)=2;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+10)=*(BodreArea+2*2+0);
}
Bord_num=Bord_num+2*divProt+Tdiv1;
}
//并修正边界信息 !!!!将常温边界2修正成受火边界3!!!
if(*(BodreArea+2*2+1)==3) //修改柱2面的边界条件
{
for(loop=BFbegceng;loop<BFendceng;loop++)
{
for(loop1=0;loop1<2*divProt+Tdiv1;loop1++)
{
*(TBorderCon+(TStoryElem_num*2+(4*divProt+Tdiv1+Tdiv2)*zdiv+(2*divProt+Tdiv1)*loop+loop1)*(TElem_node+3)+9)=3;
}
}
}
for(loop=0;loop<zdiv;loop++) //柱的边界单元的编码
{
for(loop1=0;loop1<2*divProt+Tdiv2;loop1++) //3面 侧面的左边面
{
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+0)=loop1*(2*divProt+Tdiv1)+loop*TStoryElem_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+1)=loop1*(2*divProt+Tdiv1+1)+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+2)=(loop1+1)*(2*divProt+Tdiv1+1)+loop*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+3)=loop1*(2*divProt+Tdiv1+1)+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+4)=(loop1+1)*(2*divProt+Tdiv1+1)+(loop+1)*TStoryNode_num;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+5)=0;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+6)=2;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+7)=4;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+8)=6;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+9)=2;
*(TBorderCon+(Bord_num+loop1)*(TElem_node+3)+10)=*(BodreArea+2*3+0);
}
Bord_num=Bord_num+2*divProt+Tdiv2;
}
//并修正边界信息 !!!!将常温边界2修正成受火边界3!!!
if(*(BodreArea+2*3+1)==3) //修改柱3面的边界条件
{
for(loop=BFbegceng;loop<BFendceng;loop++)
{
for(loop1=0;loop1<2*divProt+Tdiv2;loop1++)
{
*(TBorderCon+(TStoryElem_num*2+(6*divProt+Tdiv1*2+Tdiv2)*zdiv+(2*divProt+Tdiv2)*loop+loop1)*(TElem_node+3)+9)=3;
}
}
}
//给各刚度矩阵分配内存
KTG=(double *)calloc(TotalLen,8);
KG=(double *)calloc(TotalLen,8);
KBG=(double *)calloc(TotalLen,8);
CcG=(double *)calloc(TotalLen,8);
PBG=(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);
PBG0=(double *)calloc(TNode_num,8);
Temp0=(double *)calloc(TNode_num,8);
Temp1=(double *)calloc(TNode_num,8);
Temp =(double *)calloc(TNode_num,8);
Linshi=(double *)calloc(TNode_num,8);
for(loop=0; loop<TNode_num; loop++)
{
Temp1[loop]=iniT0; Temp[loop]=iniT0;
Temp0[loop]=iniT0; Linshi[loop]=0.0;
}
//给导热矩阵,传热矩阵,总温度荷载向量赋初值
for(loop=0; loop<TotalLen; loop++)
{
KG[loop]=0.0; KG0[loop]=0.0;
KBG[loop]=0.0; KBG0[loop]=0.0;
KTG[loop]=0.0; KTG0[loop]=0.0;
CcG[loop]=0.0; CcG0[loop]=0.0;
A[loop]=0.0; C[loop]=0.0;
}
for(loop=0; loop<TNode_num; loop++)
{
PBG[loop]=0.0; PBG0[loop]=0.0;
B[loop]=0.0;
}
Envtemper=iniT0;
Step=int(TotalTime/TimeStep+0.00001);
TimeHis=(double *)calloc(Step,8);
Time=0.0;
while(Time<TotalTime)
{
NumTime++;
IterCount=0; //初始迭代次数1
IterStop=1; //设置初始值,1表示继续迭代
if(VarTstep==1) //采用定时间步长的求解策略
{
Time +=TimeStep;
TimeHis[NumTime-1]=Time; //记录时间步信息
if(Env_heatType==1) //炉内空气温度,即受火面的环境温度
Envtemper=iniT0+345.0*log10(8.0*Time/60.0+1.0);
else if(Env_heatType==2)
Envtemper=iniT0+750.0*(1.0-exp(-3.79553*sqrt(Time/3600.0)))+170.41*sqrt(Time/3600.0);
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)
Envtemper=iniT0+345.0*log10(8.0*Time/60.0+1.0);
if(Env_heatType==2)
Envtemper=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; //记录时间步信息
Envtemper +=DelTemp;
}
}
cout<<"第 "<<NumTime<<" 次时间步,升温时间为: "<<Time/60.0<<" 分钟,受火的环境温度为: "<<Envtemper<<" 摄氏度"<<endl;
WResult<<"第 "<<NumTime<<" 次时间步,升温时间为: "<<Time/60.0<<" 分钟,环境温度为："<<Envtemper<<" 摄氏度"<<endl;
for(loop=0; loop<TNode_num; loop++)
{
Temp1[loop]=iniT0; Temp[loop]=iniT0;
Temp0[loop]=iniT0; Linshi[loop]=0.0;
}
//************************时*间*步*内*迭*代************************
do
{
IterCount++; //迭代计数器
//给刚度矩阵和温度荷载向量赋初值
for(loop=0; loop<TotalLen; loop++)
{
KTG[loop]=0.0; KG[loop]=0.0;
CcG[loop]=0.0;
if(ComputeMethod==1)
{
C[loop]=0.0;
}
if(ComputeMethod==1 && IterCount==1)
{
KG0[loop]=0.0;
CcG0[loop]=0.0;
}
}
for(loop=0; loop<TNode_num; loop++)
{
PBG[loop]=0.0; B[loop]=0.0;
Temp[loop]=0.0;
if(ComputeMethod==1 && IterCount==1)
{
PBG0[loop]=0.0;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -