📄 struct2d.h
字号:
*(Zload+2)=*(RFnbanlan+2);
//下面修正总刚度系数
*(F2KZtemp+3*1+1)=*(F2KZtemp+3*1+1)-(i+1)*DelP*ComLength*ComLength/pi/pi;
*(F2KZtemp+3*2+2)=*(F2KZtemp+3*2+2)-(i+1)*DelP*ComLength*ComLength/pi/pi;
for(j=0;j<3;j++)
*(Zxewywz+j)=*(Zload+j);
//下面求解位移增量
if(agaus(F2KZtemp,Zxewywz,3)==0)
{
cout<<"按韩林海简化方法求解总刚度矩阵出现病态!\n";
WResult<<" ******按韩林海简化方法求解总刚度矩阵出现病态!******\n";
cout<<" ******按韩林海简化方法求解结束!******\n";
WResult<<" ******按韩林海简化方法求解结束!******\n";
cout<<" ******恒载升温的恒载段的第"<<(i+1)<<" 分析子步求解结束******\n";
WResult<<" ******恒载升温的恒载段的第"<<(i+1)<<" 分析子步求解结束******\n";
goto s102;
}
}
if(Aitken==1) //增量步内采用Aitken加速收敛方法
{
if(liang2==(liang2/2)*2+1) //k=1,3,....奇数
{
for(liang3=0;liang3<3;liang3++)
*(uAitken31+liang3)=*(Zxewywz0+liang3)-*(Zxewywz+liang3);
//求分子
xyytemp1=0.0;
for(liang3=0;liang3<3;liang3++)
xyytemp1=xyytemp1+*(uAitken31+liang3)*(*(Zxewywz0+liang3));
//求分母
xyytemp2=0.0;
for(liang3=0;liang3<3;liang3++)
xyytemp2=xyytemp2+*(uAitken31+liang3)*(*(uAitken31+liang3));
wAitken=xyytemp1/xyytemp2;
//增量位移加速
for(liang3=0;liang3<3;liang3++)
*(Zxewywz+j)=wAitken*(*(Zxewywz+j));
}
}
for(j=0;j<3;j++)
*(Zxewywz0+j)=*(Zxewywz+j); //截面基本未知量备份
//截面未知位移更新
e0=e0+*(Zxewywz+0);
wany=wany+*(Zxewywz+2);
wanz=wanz+*(Zxewywz+1);
//构件位移更新；
if(F2method==1) //按Lie方法
{
xewywz[0][0]=e0;
xewywz[1][0]=wanz;
xewywz[2][0]=wany;
CUy1=-(wanz*ComLength*ComLength/12.0+u0y);
CUz1=-(wany*ComLength*ComLength/12.0+u0z);
UAxial1=e0*CFirLength+(1/2.0+e0/2.0)*((CUy1+u0y)*(CUy1+u0y)+(CUz1+u0z)*(CUz1+u0z))*24.0/5.0/ComLength;
}
if(F2method==2) //按韩林海方法；
{
xewywz[0][0]=e0;
xewywz[1][0]=wanz;
xewywz[2][0]=wany;
CUy1=-(wanz*ComLength*ComLength/pi/pi+u0y);
CUz1=-(wany*ComLength*ComLength/pi/pi+u0z);
UAxial1=e0*CFirLength+(1/2.0+e0/2.0)*((CUy1+u0y)*(CUy1+u0y)+(CUz1+u0z)*(CUz1+u0z))*pi*pi/2.0/ComLength;
}
//下面形成总应变增量向量
//混凝土的总应变增量
matmpl(NCSect,F2ConElem_num,3,Zxewywz,1,ZCZstrain);
//钢筋的总应变增量
matmpl(NSSect,F2Longbar_num,3,Zxewywz,1,ZSZstrain);
//形成应力应变增量向量
for(j=0;j<F2ConElem_num;j++)
{
*(ZCStstrain+j)=*(ZCZstrain+j);
}
for(j=0;j<F2Longbar_num;j++)
{
*(ZSStstrain+j)=*(ZSZstrain+j);
}
for(j=0;j<F2ConElem_num;j++)
{
*(CElemstrain+j)=*(CElemstrain+j)+*(ZCStstrain+j);
}
for(j=0;j<F2Longbar_num;j++)
{
*(SElemstrain+j)=*(SElemstrain+j)+*(ZSStstrain+j);
}
//求混凝土和钢筋应力向量
for(j=0; j<F2ConElem_num; j++)
{
//先求单元的平均温度
xyytemp1=(*(FTnode0+*(F2CElemNode_Info+j*4))+*(FTnode0+*(F2CElemNode_Info+j*4+1))+*(FTnode0+*(F2CElemNode_Info+j*4+2))+*(FTnode0+*(F2CElemNode_Info+j*4+3)))/4.0;
*(CElemstress+j)=CStress(Csresta,*(CElemstrain+j),xyytemp1);
}
for(j=0;j<F2Longbar_num;j++)
{
xyytemp1=*(STnode0+j); //钢筋温度
*(SElemstress+j)=SStress(Ssresta,*(SElemstrain+j),xyytemp1,*(F2LbarInfo+5*j+3),*(F2LbarInfo+5*j+3)/(*(F2LbarInfo+5*j+4)));
}
//计算截面总内力
//先求混凝土对内力的贡献
matmpl(ConAc,F2ConElem_num,F2ConElem_num,CElemstress,1,JuZhenC2);
for(j=0;j<3;j++)
for(k=0;k<F2ConElem_num;k++)
*(JuZhenC1+j*F2ConElem_num+k)=*(NCSect+k*3+j);
matmpl(JuZhenC1,3,F2ConElem_num,JuZhenC2,1,FcSect);
//计算钢筋对内力的贡献
matmpl(SteAs,F2Longbar_num,F2Longbar_num,SElemstress,1,JuZhenS2);
//先将钢筋单元形函数转置
for(j=0;j<3;j++)
for(k=0;k<F2Longbar_num;k++)
*(JuZhenS1+j*F2Longbar_num+k)=*(NSSect+k*3+j);
matmpl(JuZhenS1,3,F2Longbar_num,JuZhenS2,1,JuZhens31);
for(j=0;j<3;j++)
*(FcSect+j)=*(FcSect+j)+*(JuZhens31+j);
//计算不平衡力
*(RFnbanlan+0)=(i+1)*DelP-*(FcSect+0);
*(RFnbanlan+1)=-(i+1)*DelP*CUy1-*(FcSect+1);
*(RFnbanlan+2)=-(i+1)*DelP*CUz1-*(FcSect+2);
error=fabs(RFnbanlan[0]/FcSect[0]);
if(error<fabs(RFnbanlan[1]/FcSect[1]))
error=fabs(*(RFnbanlan+1)/FcSect[1]);
if(error<fabs(RFnbanlan[2]/FcSect[2]))
error=fabs(*(RFnbanlan+2)/FcSect[2]);
cout<<" 常温下迭代误差 "<<error<<" 容许误差 "<<FNMaxErr0<<endl;
WResult<<" 常温下迭代误差 "<<error<<" 容许误差 "<<FNMaxErr0<<endl;
}while(liang2<200 && error>FNMaxErr0); //迭代终止
if(liang2<200)
{
cout<<" 恒载升温横载段的第 "<<(i+1)<<" 子步成功收敛!\n";
cout<<"\n";
WResult<<" 恒载升温横载段的第 "<<(i+1)<<" 子步成功收敛!\n";
WResult<<"\n";
}
if(liang2>=200)
{
cout<<" 恒载升温横载段的第 "<<(i+1)<<" 子步收敛失败!\n";
cout<<"\n";
WResult<<" 恒载升温横载段的第 "<<(i+1)<<" 子步收敛失败!\n";
WResult<<"\n";
cout<<" ******恒载升温的恒载段的第"<<(i+1)<<" 分析子步求解结束******\n";
WResult<<" ******恒载升温的恒载段的第"<<(i+1)<<" 分析子步求解结束******\n";
goto s102;
}
//输出截面有构件的计算结果
R2Weiyi<<(i+1)<<","<<-UAxial1*1000.0<<","<<(CUy1+u0y)*1000.0<<","<<(CUz1+u0z)*1000.0<<endl;
R2Sectwei<<(i+1)<<","<<e0<<","<<wany<<","<<wanz<<endl;
}
WResult<<"******恒载升温段分析*****"<<endl;
cout<<"*****恒载升温段分析*****"<<endl;
//****************************************************
//****************************************************
//**********************恒载升温的升温段的分析******************************
//****************************************************
for(i=0; i<Step; i++) //恒载升温的升温段的分析
{
cout<<"恒载升温的升温段的第 "<<(i+1)<<" 子步,迭代步的起始时间:"<<curtime0/60.0<<" 分钟\n";
WResult<<"恒载升温的升温段的第 "<<(i+1)<<" 子步,迭代步的起始时间:"<<curtime0/60.0<<" 分钟\n";
UAxial0=UAxial1;
//读入下一时间步的各节点的温度
for(j=0; j<TNode_num; j++)
{
Struc_Tin>>*(Ttnode+j);
//printf("节点温度%.4f !\n",*(Ttnode+j));
}
//下面计算结构网格的各节点下一时间步所对应的温度
for(j=0; j<F2Node_num; j++)
{
//温度场分析单元的a和b
xyytemp3=fabs(*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j))+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j)))*2));
xyytemp4=fabs(*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j))+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j)))*2+1));
xyytemp1=*(F2Node_Coor+j*2)-(*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j))+1)*2)+*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j)))*2))/2.0;
xyytemp2=*(F2Node_Coor+j*2+1)-(*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j))+2)*2+1)+*(Node_Coor+*(TCElemNode_Info+4*(*(F2CN_TE+j)))*2+1))/2.0;
//printf("%.5f %.5f\n",xyytemp1,xyytemp2)节点单元相对坐标及信息
*(FTnode1+j)=*(Ttnode+*(TCElemNode_Info+*(F2CN_TE+j)*4))*1.0/4.0*(1-xyytemp1/(xyytemp3/2.0))*(1-xyytemp2/(xyytemp4/2.0));
*(FTnode1+j)=*(FTnode1+j)+*(Ttnode+*(TCElemNode_Info+*(F2CN_TE+j)*4+1))*1.0/4.0*(1+xyytemp1/(xyytemp3/2.0))*(1-xyytemp2/(xyytemp4/2.0));
*(FTnode1+j)=*(FTnode1+j)+*(Ttnode+*(TCElemNode_Info+*(F2CN_TE+j)*4+2))*1.0/4.0*(1-xyytemp1/(xyytemp3/2.0))*(1+xyytemp2/(xyytemp4/2.0));
*(FTnode1+j)=*(FTnode1+j)+*(Ttnode+*(TCElemNode_Info+*(F2CN_TE+j)*4+3))*1.0/4.0*(1+xyytemp1/(xyytemp3/2.0))*(1+xyytemp2/(xyytemp4/2.0));
//printf("第%d 个节点的温度%.10f !\n.",j,*(FTnode1+j));
}
for(j=0; j<F2ConElem_num; j++)
{
xyytemp1=(*(FTnode1+*(F2CElemNode_Info+j*4))+*(FTnode1+*(F2CElemNode_Info+j*4+1))+*(FTnode1+*(F2CElemNode_Info+j*4+2))+*(FTnode1+*(F2CElemNode_Info+j*4+3)))/4.0;
*(CElemtemp1+j)=xyytemp1;
//printf("第%d 个单元的温度%.4f !\n.",j,xyytemp1);
}
//求混凝土单元节点的温度增量
for(j=0;j<F2ConElem_num;j++)
{
if(*(CElemtemp1+j)>=*(CElemtemp0+j))
*(CTempZ+j)=*(CElemtemp1+j)-*(CElemtemp0+j);
if(*(CElemtemp1+j)<*(CElemtemp0+j))
*(CTempZ+j)=0.0;
}
curtime1=*(TimeHis+i); //下一时间步的时间
//求时间增量
TimeZ=curtime1-curtime0;
for(j=0;j<F2ConElem_num;j++)
*(CTimeZ+j)=TimeZ;
for(j=0;j<F2Longbar_num;j++)
*(STimeZ+j)=TimeZ/3600.0;
//钢筋下一时间步的温度
for(j=0;j<F2Longbar_num;j++)
{
//每根钢筋所对应的点的坐标
xyytemp1=*(F2LbarInfo+5*j+0); //z坐标
xyytemp2=*(F2LbarInfo+5*j+1); //y坐标
//温度场分析单元的a和b
xyytemp3=fabs(*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j))+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j)))*2));
xyytemp4=fabs(*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j))+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j)))*2+1));
//printf("%.5f %.5f\n",xyytemp1,xyytemp2);//节点坐标及信息
//节点坐标x和y
xyytemp1=xyytemp1-(*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j))+1)*2)+*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j)))*2))/2.0;
xyytemp2=xyytemp2-(*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j))+2)*2+1)+*(Node_Coor+*(TCElemNode_Info+4*(*(F2SN_TE+j)))*2+1))/2.0;
//printf("%.5f %.5f\n",xyytemp1,xyytemp2);//节点坐标及信息
*(STnode1+j)=*(Ttnode+*(TCElemNode_Info+*(F2SN_TE+j)*4))*1.0/4.0*(1-xyytemp1/(xyytemp3/2.0))*(1-xyytemp2/(xyytemp4/2.0));
*(STnode1+j)=*(STnode1+j)+*(Ttnode+*(TCElemNode_Info+*(F2SN_TE+j)*4+1))*1.0/4.0*(1+xyytemp1/(xyytemp3/2.0))*(1-xyytemp2/(xyytemp4/2.0));
*(STnode1+j)=*(STnode1+j)+*(Ttnode+*(TCElemNode_Info+*(F2SN_TE+j)*4+2))*1.0/4.0*(1-xyytemp1/(xyytemp3/2.0))*(1+xyytemp2/(xyytemp4/2.0));
*(STnode1+j)=*(STnode1+j)+*(Ttnode+*(TCElemNode_Info+*(F2SN_TE+j)*4+3))*1.0/4.0*(1+xyytemp1/(xyytemp3/2.0))*(1+xyytemp2/(xyytemp4/2.0));
// printf("升温时间%.4f小时,第%d 个钢筋节点的温度%.4f !\n.",curtime1/3600.0,j,*(STnode1+j));
}
if(Spall_thick>=Prot_thick)
{
for(j=0;j<F2Longbar_num;j++)
{ //说明钢筋中心处的混凝土应该爆裂掉
if(*(STnode1+j)>Critical_T)
{
if(Env_heatType==1)
*(STnode1+j)=iniT0+345.0*log10(8.0*curtime0/60.0+1.0);
if(Env_heatType==2)
*(STnode1+j)=iniT0+750.0*(1.0-exp(-3.79553*sqrt(curtime0/3600.0)))+170.41*sqrt(curtime0/3600.0);
}
}
}
//求钢筋位置的温度增量
for(j=0;j<F2Longbar_num;j++)
{
if(*(STnode1+j)>=*(STnode0+j))
*(STempZ+j)=*(STnode1+j)-*(STnode0+j);
if(*(STnode1+j)<*(STnode0+j))
*(STempZ+j)=0.0;
}
//钢筋温度补偿时间的导数
for(j=0;j<F2Longbar_num;j++)
{
xyytemp3=*(SJZ_Wendubu+j);
*(SJZ_Wendubu+j)=*(SJZ_Wendubu+j)+exp(-75000.0/(*(STnode0+j)*9.0/5.0+491.0))*TimeZ/3600.0; //钢筋温度补偿时间(小时)
if(*(STempZ+j)>=0.000001)
{
*(SJZ_WendubuDT+j)=(*(SJZ_Wendubu+j)-xyytemp3)/(*(STempZ+j));
}
if(*(STempZ+j)<0.000001)
{
*(SJZ_WendubuDT+j)=0.0;
}
}
//形成混凝土的切线模量矩阵,温度切线模量矩阵
for(j=0; j<F2ConElem_num; j++)
{
//先求单元的平均温度
xyytemp1=*(CElemtemp0+j);
*(ConEstr+j*F2ConElem_num+j)=Ecst(Csresta,*(CElemstress+j),*(CElemstrain+j),xyytemp1,curtime0);
*(ConET+j*F2ConElem_num+j)=Ectempt(Csresta,*(CElemstress+j),*(CElemstrain+j),xyytemp1,curtime0);
if(*(CElemstress+j)>=0.0) //应力大于等于0
{
*(ConETime+j*F2ConElem_num+j)=Ectimet(Csresta,*(CElemstress+j),*(CElemstrain+j),xyytemp1,curtime0);
}
if(*(CElemstress+j)<0.0) //应力小于0
{
*(ConETime+j*F2ConElem_num+j)=0.0;
}
if(xyytemp1>Critical_T) //大于爆裂临界温度
{ //处于爆裂区
if(Sectwid/2.0-*(NCSect+3*j+2)<Spall_thick || *(NCSect+3*j+2)+Sectwid/2.0<Spall_thick || Sectlen/2.0-*(NCSect+3*j+1)<Spall_thick || *(NCSect+3*j+1)+Sectlen/2.0<Spall_thick)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -