📄 tempf3d.h
字号:
}
if(ComputeMethod==2)
{
Linshi[loop]=0.0;
}
}
//将内部单元导热矩阵和质量热容矩阵组装总刚度矩阵
for(loop=0; loop<TConElem_num; loop++)
{
TempAverage=0.0;
for(loop1=0; loop1<TElem_node; loop1++)//求单元平均温度
{
TempAverage +=Temp1[*(TCElemNode_Info+TElem_node*loop+loop1)]/double(TElem_node);
}
//下面取出六面体单元的边长的一半 x y z
Solid_a=fabs(*(Node_Coor+*(TCElemNode_Info+loop*TElem_node+1)*3+0)-*(Node_Coor+*(TCElemNode_Info+loop*TElem_node+0)*3+0))/2.0;
Solid_b=fabs(*(Node_Coor+*(TCElemNode_Info+loop*TElem_node+2)*3+1)-*(Node_Coor+*(TCElemNode_Info+loop*TElem_node+0)*3+1))/2.0;
Solid_c=fabs(*(Node_Coor+*(TCElemNode_Info+loop*TElem_node+4)*3+2)-*(Node_Coor+*(TCElemNode_Info+loop*TElem_node+0)*3+2))/2.0;
CK=CalcuKc(TempAverage,CCondu_Type); //当前温度的导热系数
//刚度矩阵的第一行
CKTe[0][0]=CK*2.0/9.0*(Solid_a*Solid_c/Solid_b+Solid_b*Solid_c/Solid_a+Solid_a*Solid_b/Solid_c);
CKTe[0][1]=CK*(Solid_a*Solid_c/Solid_b/9.0-2.0/9.0*Solid_b*Solid_c/Solid_a+Solid_a*Solid_b/Solid_c/9.0);
CKTe[0][2]=CK*(Solid_b*Solid_c/Solid_a/9.0-2.0/9.0*Solid_a*Solid_c/Solid_b+Solid_a*Solid_b/Solid_c/9.0);
CKTe[0][3]=CK*(-Solid_b*Solid_c/Solid_a/9.0-Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/18.0);
CKTe[0][4]=CK*(Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0-2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[0][5]=CK*(-Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/9.0);
CKTe[0][6]=CK*(Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/9.0-Solid_a*Solid_b/Solid_c/9.0);
CKTe[0][7]=CK*(-Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/18.0);
//第二行
CKTe[1][0]=CKTe[0][1];
CKTe[1][1]=CK*(2.0*Solid_a*Solid_c/Solid_b/9.0+2.0/9.0*Solid_b*Solid_c/Solid_a+2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[1][2]=CK*(-Solid_b*Solid_c/Solid_a/9.0-Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/18.0);
CKTe[1][3]=CK*(Solid_b*Solid_c/Solid_a/9.0-2.0*Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/9.0);
CKTe[1][4]=CK*(-Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/9.0);
CKTe[1][5]=CK*(Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0-2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[1][6]=CK*(-Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/18.0);
CKTe[1][7]=CK*(Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/9.0-Solid_a*Solid_b/Solid_c/9.0);
//第三行
CKTe[2][0]=CKTe[0][2];
CKTe[2][1]=CKTe[1][2];
CKTe[2][2]=CK*(2.0*Solid_b*Solid_c/Solid_a/9.0+2.0/9.0*Solid_a*Solid_c/Solid_b+Solid_a*Solid_b/Solid_c/9.0);
CKTe[2][3]=CK*(-2.0*Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/9.0);
CKTe[2][4]=CK*(Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/9.0-Solid_a*Solid_b/Solid_c/9.0);
CKTe[2][5]=CK*(-Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/18.0);
CKTe[2][6]=CK*(Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0-2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[2][7]=CK*(-Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/9.0);
//第四行
CKTe[3][0]=CKTe[0][3];
CKTe[3][1]=CKTe[1][3];
CKTe[3][2]=CKTe[2][3];
CKTe[3][3]=CK*(2.0*Solid_b*Solid_c/Solid_a/9.0+2.0*Solid_a*Solid_c/Solid_b/9.0+2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[3][4]=CK*(-Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/18.0);
CKTe[3][5]=CK*(Solid_b*Solid_c/Solid_a/18.0-Solid_a*Solid_c/Solid_b/9.0-Solid_a*Solid_b/Solid_c/9.0);
CKTe[3][6]=CK*(-Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/18.0-Solid_a*Solid_b/Solid_c/9.0);
CKTe[3][7]=CK*(Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0-2.0*Solid_a*Solid_b/Solid_c/9.0);
//第五行
CKTe[4][0]=CKTe[0][4];
CKTe[4][1]=CKTe[1][4];
CKTe[4][2]=CKTe[2][4];
CKTe[4][3]=CKTe[3][4];
CKTe[4][4]=CK*(2.0*Solid_b*Solid_c/Solid_a/9.0+2.0*Solid_a*Solid_c/Solid_b/9.0+2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[4][5]=CK*(-2.0*Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/9.0);
CKTe[4][6]=CK*(Solid_b*Solid_c/Solid_a/9.0-2.0*Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/9.0);
CKTe[4][7]=CK*(-Solid_b*Solid_c/Solid_a/9.0-Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/18.0);
//第六行
CKTe[5][0]=CKTe[0][5];
CKTe[5][1]=CKTe[1][5];
CKTe[5][2]=CKTe[2][5];
CKTe[5][3]=CKTe[3][5];
CKTe[5][4]=CKTe[4][5];
CKTe[5][5]=CK*(2.0*Solid_b*Solid_c/Solid_a/9.0+2.0*Solid_a*Solid_c/Solid_b/9.0+2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[5][6]=CK*(-Solid_b*Solid_c/Solid_a/9.0-Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/18.0);
CKTe[5][7]=CK*(Solid_b*Solid_c/Solid_a/9.0-2.0*Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/9.0);
//第七行
CKTe[6][0]=CKTe[0][6];
CKTe[6][1]=CKTe[1][6];
CKTe[6][2]=CKTe[2][6];
CKTe[6][3]=CKTe[3][6];
CKTe[6][4]=CKTe[4][6];
CKTe[6][5]=CKTe[5][6];
CKTe[6][6]=CK*(2.0*Solid_b*Solid_c/Solid_a/9.0+2.0*Solid_a*Solid_c/Solid_b/9.0+2.0*Solid_a*Solid_b/Solid_c/9.0);
CKTe[6][7]=CK*(-2.0*Solid_b*Solid_c/Solid_a/9.0+Solid_a*Solid_c/Solid_b/9.0+Solid_a*Solid_b/Solid_c/9.0);
//第八行
CKTe[7][0]=CKTe[0][7];
CKTe[7][1]=CKTe[1][7];
CKTe[7][2]=CKTe[2][7];
CKTe[7][3]=CKTe[3][7];
CKTe[7][4]=CKTe[4][7];
CKTe[7][5]=CKTe[5][7];
CKTe[7][6]=CKTe[6][7];
CKTe[7][7]=CK*(2.0*Solid_b*Solid_c/Solid_a/9.0+2.0*Solid_a*Solid_c/Solid_b/9.0+2.0*Solid_a*Solid_b/Solid_c/9.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);
KTG[ibuf]=KTG[ibuf]+CKTe[loop1][loop2];
}
}
CSpecialHeat=CalcuCSpecialHeat(TempAverage,CSpecilHeat_Type);//计算单前温度的质量热容
//第一行:过镇海书,P90,式(6-23)
Cce[0][0]=CSpecialHeat*Solid_a*Solid_b*Solid_c*8.0/27.0;
Cce[0][1]=Cce[0][0]/2.0;
Cce[0][2]=Cce[0][0]/2.0;
Cce[0][3]=Cce[0][0]/4.0;
Cce[0][4]=Cce[0][0]/2.0;
Cce[0][5]=Cce[0][0]/4.0;
Cce[0][6]=Cce[0][0]/4.0;
Cce[0][7]=Cce[0][0]/8.0;
//第二行
Cce[1][0]=Cce[0][1];
Cce[1][1]=Cce[0][0];
Cce[1][2]=Cce[0][3];
Cce[1][3]=Cce[0][1];
Cce[1][4]=Cce[0][3];
Cce[1][5]=Cce[0][1];
Cce[1][6]=Cce[0][7];
Cce[1][7]=Cce[0][3];
//第三行
Cce[2][0]=Cce[0][2];
Cce[2][1]=Cce[1][2];
Cce[2][2]=Cce[0][0];
Cce[2][3]=Cce[0][1];
Cce[2][4]=Cce[0][3];
Cce[2][5]=Cce[0][7];
Cce[2][6]=Cce[0][1];
Cce[2][7]=Cce[0][3];
//第四行
Cce[3][0]=Cce[0][3];
Cce[3][1]=Cce[1][3];
Cce[3][2]=Cce[2][3];
Cce[3][3]=Cce[0][0];
Cce[3][4]=Cce[0][7];
Cce[3][5]=Cce[0][3];
Cce[3][6]=Cce[0][3];
Cce[3][7]=Cce[0][1];
//第五行
Cce[4][0]=Cce[0][4];
Cce[4][1]=Cce[1][4];
Cce[4][2]=Cce[2][4];
Cce[4][3]=Cce[3][4];
Cce[4][4]=Cce[0][0];
Cce[4][5]=Cce[0][1];
Cce[4][6]=Cce[0][1];
Cce[4][7]=Cce[0][3];
//第六行
Cce[5][0]=Cce[0][5];
Cce[5][1]=Cce[1][5];
Cce[5][2]=Cce[2][5];
Cce[5][3]=Cce[3][5];
Cce[5][4]=Cce[4][5];
Cce[5][5]=Cce[0][0];
Cce[5][6]=Cce[0][3];
Cce[5][7]=Cce[0][1];
//第七行
Cce[6][0]=Cce[0][6];
Cce[6][1]=Cce[1][6];
Cce[6][2]=Cce[2][6];
Cce[6][3]=Cce[3][6];
Cce[6][4]=Cce[4][6];
Cce[6][5]=Cce[5][6];
Cce[6][6]=Cce[0][0];
Cce[6][7]=Cce[0][1];
//第八行
Cce[7][0]=Cce[0][7];
Cce[7][1]=Cce[1][7];
Cce[7][2]=Cce[2][7];
Cce[7][3]=Cce[3][7];
Cce[7][4]=Cce[4][7];
Cce[7][5]=Cce[5][7];
Cce[7][6]=Cce[6][7];
Cce[7][7]=Cce[0][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);
CcG[ibuf]=CcG[ibuf]+Cce[loop1][loop2];
}
}
}//!!!内部单元的导热矩阵和质量热容矩阵组装完毕!!!!!
//下面计算六面体边界单元的传热矩阵和传热荷载向量
for(loop=0;loop<Border_Elemnum;loop++)
{
if(Bord_heatType==1) //1表示采用综合换热系数
{
if(*(TBorderCon+loop*(TElem_node+3)+9)==2) //常温边界条件的处理
{
CalcuBt(ChangHeatCT,iniT0);
}
if(*(TBorderCon+loop*(TElem_node+3)+9)==3) //受火边界条件的处理;过镇海书P93,(6-38b)
{
CalcuBt(ChangHeatCT,Envtemper);
}
BordEnum0=*(TBorderCon+loop*(TElem_node+3)+0); //边界单元所对应的单元号
Enode1=*(TBorderCon+loop*(TElem_node+3)+1); //边界单元所对应的节点号1
Enode2=*(TBorderCon+loop*(TElem_node+3)+2); //边界单元所对应的节点号2
Enode3=*(TBorderCon+loop*(TElem_node+3)+3); //边界单元所对应的节点号3
Enode4=*(TBorderCon+loop*(TElem_node+3)+4); //边界单元所对应的节点号4
if(*(TBorderCon+loop*(TElem_node+3)+10)==4 || *(TBorderCon+loop*(TElem_node+3)+10)==5)
{
Bord6_a=fabs(*(Node_Coor+3*Enode2+0)-*(Node_Coor+3*Enode1+0))/2.0;
Bord6_b=fabs(*(Node_Coor+3*Enode3+1)-*(Node_Coor+3*Enode1+1))/2.0;
}
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×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)
{
PBe[0]=Bord6_a*Bord6_b*ChangHeatCT*Envtemper;
PBe[1]=Bord6_a*Bord6_b*ChangHeatCT*Envtemper;
PBe[2]=Bord6_a*Bord6_b*ChangHeatCT*Envtemper;
PBe[3]=Bord6_a*Bord6_b*ChangHeatCT*Envtemper;
}
}
if(Bord_heatType==2) //2表示对流和辐射在每个时间步内综合考虑
{
BordEnum0=*(TBorderCon+loop*(TElem_node+3)+0); //边界单元所对应的单元号
Enode1=*(TBorderCon+loop*(TElem_node+3)+1); //边界单元所对应的节点号1
Enode2=*(TBorderCon+loop*(TElem_node+3)+2); //边界单元所对应的节点号2
Enode3=*(TBorderCon+loop*(TElem_node+3)+3); //边界单元所对应的节点号3
Enode4=*(TBorderCon+loop*(TElem_node+3)+4); //边界单元所对应的节点号4
if(*(TBorderCon+loop*(TElem_node+3)+10)==4 || *(TBorderCon+loop*(TElem_node+3)+10)==5)
{
Bord6_a=fabs(*(Node_Coor+3*Enode2+0)-*(Node_Coor+3*Enode1+0))/2.0;
Bord6_b=fabs(*(Node_Coor+3*Enode3+1)-*(Node_Coor+3*Enode1+1))/2.0;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -