📄 struct2d.h
字号:
void Struct2D();
void Struct2D()
{
ifstream F2comp("F2comp.txt");
ifstream F2zbar("F2zbar.txt");
ifstream Struc_Tin("Struc_Tin.txt");
ofstream R2Weiyi("R2Weiyi.csv");
ofstream R2Sectwei("R2Sectwei.csv");
ofstream WResult("WResult.txt",ios::ate);
ofstream OutInitData("OutInitData.txt");
if(F2comp.fail()) //检查输入文件
{
cerr<<"error opening the F2comp.txt!\n";
exit(0);
}
if(F2zbar.fail())
{
cerr<<"error opening the F2zbar.txt!\n";
exit(0);
}
// int Expand; //轴向位移是否处于膨胀段,1是,0否
// double MaxUAxial=-1000.0; //最大轴向位移
double *Ttnode; //每一时间步各节点的温度
double e0; //截面几何中心处的正应变
double wany,wanz; //截面绕y轴的曲率,绕z轴的曲率
double z0,y0; //截面形心处的z坐标z0,y坐标y0
double *NCSect; //截面混凝土单元的形函数矩阵
double *NSSect; //截面钢筋单元的形函数矩阵
double xewywz[3][1]; //截面基本未知量向量
double *ConAc; //混凝土单元的单元面积矩阵
double *SteAs; //钢筋单元的单元面积矩阵
double *ConEstr; //混凝土单元的应力切线模量列矩阵
double *ConET; //混凝土单元的温度切线模量列矩阵
double *ConETime; //混凝土单元的时间切线模量
double *SteEstr; //钢筋单元的应力切线模量列矩阵
double *SteET; //钢筋的温度切线模量矩阵
double *SteETime; //钢筋的时间切线模量矩阵
//截面的状态
double *CElemstress; //各混凝土单元的应力
double *CElemstrain; //各混凝土单元的应变
double *CElemtemp0; //各混凝土单元的温度
double *CElemtemp1; //混凝土单元下一迭代步的单元温度
double *SElemstress; //各钢筋单元的应力
double *SElemstrain; //各钢筋单元的应变
double curtime0; //当前时间步对应的时间
double curtime1; //下一时间步对应的时间
long *F2CN_TE; //结构分析的每个混凝土单元的节点所对应的温度场分析的单元号
double *FTnode0; //结构分析当前时间步结构网格各节点所对应的温度
double *FTnode1; //结构分析下一时间步结构网格各节点所对应的温度
double *STnode0; //钢筋所在位置点的当前时间步的温度
double *STnode1; //钢筋所在位置点的下一时间步的温度
long *F2SN_TE; //钢筋位置所对应的温度场分析时混凝土单元编号
double TimeZ; //每个时间步的时间增量
double *CTempZ; //混凝土单元每个时间步的温度增量
double *STempZ; //钢筋单元每个时间步的温度增量
double *CTimeZ; //混凝土单元的时间增量
double *STimeZ; //钢筋单元的时间增量
double e0y,e0z; //荷载绕z轴,y轴的偏心矩
double EfLength; //柱的计算长度系数
double u0y,u0z; //柱在绕z轴,y轴方向的初始缺陷
double Coepx; //柱的初始缺陷系数,取柱计算长度的1/1000.0
double CFirLength; //柱的受火高度
double ComLength; //柱计算长度
double PCentral; //柱的集中荷载大小,单位KN
int Pdiv; //采用恒载升温途径时,先恒载时划分的荷载步数
double DelP; //恒载时的荷载步增量
double *F2KZ; //简化分析时的总刚度矩阵
double *CEsmAC; //混凝土应力切线模量和混凝土面积的乘积
double *SEsmAs; //钢筋应力切线模量和钢筋面积的乘积
double *CETmAC; //混凝土的温度总弹塑性矩阵 温度切线模量和混凝土面积的乘积;
double *CETmAC0; //混凝土的温度总弹塑性矩阵备份;温度切线模量和混凝土面积的乘积
double *CETimAc; //混凝土的时间总弹塑性矩阵,时间切线模量和混凝土面积的乘积
double *CETimAc0; //混凝土的时间总弹塑性矩阵备份,时间切线模量和混凝土面积的乘积
double *SETmAs; //钢筋的温度总弹塑性矩阵,钢筋温度切线模量和钢筋面积的乘积
double *SETmAs0; //钢筋的温度总弹塑性矩阵备份,钢筋温度切线模量和钢筋面积的乘积
double *SETimAs; //钢筋的时间切线模量和钢筋面积的乘积
double *SETimAs0; //钢筋的时间总弹塑性矩阵备份;钢筋的时间切线模量和钢筋面积的乘积
double *NCmES; //混凝土单元形状函数和Es,cs的乘积;中间矩阵用于形成总刚度矩阵
double *NSmES; //钢筋单元形状函数和Es,cs的乘积 中间矩阵用于形成总刚度矩阵
double *F2KZtemp; //3*3的临时矩阵
double *JuZhenC1; //混凝土单元的临时矩阵F2ConElem_num*3
double *JuZhenC10; //混凝土单元的临时矩阵F2ConElem_num*3
double *JuZhenC2; //混凝土单元的临时矩阵F2ConElem_num*1
double *JuZhenS1; //钢筋单元的临时矩阵3*F2Longbar_num
double *JuZhenS10; //钢筋单元的临时矩阵3*F2Longbar_num
double *JuZhenS2; //钢筋单元的临时矩阵F2Longbar_num*1
double *JuZhen31; //3*1的临时矩阵
double *JuZhens31; //3*1的临时矩阵
double *Zxewywz; //截面基本未知量增量向量
double *Zxewywz0; //截面基本未知量增量向量备份
double *Zload; //荷载增量
double UAxial1; //轴向变形
double UAxial0; //前一时间步的轴向位移
double CUy0,CUz0; //当前增量步绕z轴,y轴的侧向位移;
double CUy1,CUz1; //当前增量步绕z轴,y轴的迭代过程中的侧向位移
double *ZCZstrain; //混凝土的总应变增量向量
double *ZSZstrain; //钢筋的总应变增量向量
double *ZCStstrain; //混凝土的应力应变增量向量
double *ZSStstrain; //钢筋的应力应变增量向量
double *FcSect; //截面总内力向量
double *RFnbanlan; //不平衡力向量
double FNMaxErr; //常温时不平衡力的最大误差
double FNMaxErr0; //升温时不平衡力的最大误差
int Aitken; //增量步内是否采用Aitken加速收敛方法;1表示采用,0表示不采用
double wAitken; //Aitken的加速收敛因子
double *uAitken31; //计算wAitken的临时矩阵3*1
int PHxiaozheng; //是否采用平衡校正的迭代算法
int EuAdams; //解除温度增量荷载和时间增量荷载与基本变量间的耦合性
double *DelFT; //温度增量荷载向量
double *DelFTime; //时间增量荷载向量
double *CJZ0; //临时变量
double *SJZ0; //临时变量
double *SJZ_Wendubu; //钢筋温度补偿时间(小时)
double *SJZ_WendubuDT; //钢筋温度补偿时间对温度的导数
double xyytemp1,xyytemp2,xyytemp3,xyytemp4;
double *liyishuox0,*liyishuox1; //临时变量;
double *liyishuok0,*liyishuok1; //临时变量;
double *liyishuoy0,*liyishuoy1; //临时变量;
int i,j,k,liang2,liang3;
double error=0.0;
char NoUse[200];//为输入设置的无用的变量
//***从文件F2comp.txt读取初始数据***
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>Csresta;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>Ssresta;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>sgaowenx;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>fc;
fc=fc*1.0e6;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>cp;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>ceutt;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>sfu;
sfu=sfu*1.0e6;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>e0y>>e0z;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>CFirLength;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>EfLength;
ComLength=CFirLength*EfLength; //柱的计算长度
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>Coepx;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>PCentral;
PCentral=1000.0*PCentral;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>Pdiv;
DelP=PCentral/Pdiv;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>FNMaxErr0;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>FNMaxErr;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>Aitken;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>PHxiaozheng;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>EuAdams;
F2comp.getline(NoUse,sizeof(NoUse));
F2comp>>F2method;
//***从文件F2zbar.txt读取初始数据***
F2zbar>>F2Longbar_num;
F2LbarInfo=(double *)calloc(F2Longbar_num*5,8);
for(i=0; i<F2Longbar_num; i++)
{
F2zbar.getline(NoUse,sizeof(NoUse));
for(j=0; j<5; j++)
F2zbar>>*(F2LbarInfo+5*i+j);
*(F2LbarInfo+5*i+3)=*(F2LbarInfo+5*i+3)*1000000.0;
*(F2LbarInfo+5*i+4)=*(F2LbarInfo+5*i+4)*1000000.0;
}
OutInitData<<" 从文件F2mesh.txt读取初始数据"<<endl;
OutInitData<<"矩形截面边长为: "<<Sectwid<<" "<<Sectlen<<endl;
OutInitData<<"保护层厚度为: "<<Prot_thick<<endl;
OutInitData<<"爆裂厚度为: "<<Spall_thick<<endl;
if(MeshID==1)OutInitData<<"与温度场网格划分一样"<<endl;
if(MeshID==0)OutInitData<<"与温度场网格划分不一样"<<endl;
OutInitData<<"结构单元网格划分数为: "<<div1<<" "<<div2<<" "<<divProt<<endl;
OutInitData<<"温度场单元网格划分数为: "<<Tdiv1<<" "<<Tdiv2<<endl;
OutInitData<<endl;
OutInitData<<" 从文件F2comp.txt读取初始数据"<<endl;
if(Csresta==1) OutInitData<<"采用过镇海的混凝土高温本构关系"<<endl;
if(Csresta==2) OutInitData<<"采用LIE的混凝土高温本构关系"<<endl;
if(Csresta==3) OutInitData<<"采用Kodur硅质骨料高强混凝土高温本构关系"<<endl;
if(Ssresta==1) OutInitData<<"采用过镇海的钢筋高温本构关系"<<endl;
if(Ssresta==2) OutInitData<<"采用LIE的钢筋高温本构关系"<<endl;
if(sgaowenx==0) OutInitData<<"不考虑钢筋高温蠕变"<<endl;
if(sgaowenx==1) OutInitData<<"考虑钢筋高温蠕变"<<endl;
OutInitData<<"混凝土的棱柱体强度为: "<<fc/1000000.0<<"MPa"<<endl;
OutInitData<<"混凝土受压峰值应力所对应的峰值应变 "<<cp<<endl;
OutInitData<<"混凝土极限拉应变 "<<ceutt<<endl;
OutInitData<<"钢筋极限抗拉强度 "<<sfu/1000000.0<<"MPa"<<endl;
OutInitData<<"荷载偏心矩 "<<e0y<<" "<<e0z<<endl;
OutInitData<<"受火长度 "<<CFirLength<<endl;
OutInitData<<"计算长度有效系数 "<<EfLength<<endl;
OutInitData<<"初始缺陷系数 "<<Coepx<<endl;
OutInitData<<"荷载大小 "<<PCentral/1000.0<<"KN"<<endl;
OutInitData<<"采用恒载升温途径时,先恒载时划分的荷载步数 "<<Pdiv<<endl;
OutInitData<<"不平衡力的最大误差 "<<FNMaxErr0<<" "<<FNMaxErr<<endl;
OutInitData<<"初始温度和爆裂临界温度为："<<iniT0<<" "<<Critical_T<<endl;
OutInitData<<"温度场分析时间步数 "<<Step<<endl;
if(Aitken==1) OutInitData<<"增量步内采用AitKen加速收敛方法"<<endl;
if(Aitken==0) OutInitData<<"增量步内不采用AitKen加速收敛方法"<<endl;
if(PHxiaozheng==1) OutInitData<<"采用平衡校正的迭代算法"<<endl;
if(PHxiaozheng==0) OutInitData<<"不采用平衡校正的迭代算法"<<endl;
if(EuAdams==1) OutInitData<<"采用Euler方法解除耦合"<<endl;
if(EuAdams==0) OutInitData<<"采用Euler和Adams方法联合解除耦合"<<endl;
if(F2method==1) OutInitData<<"采用LIE的简化方法"<<endl;
if(F2method==2) OutInitData<<"采用韩林海的简化方法"<<endl;
OutInitData<<endl;
OutInitData<<" 从文件F2zbar.txt读取初始数据"<<endl;
OutInitData<<"钢筋根数: "<<F2Longbar_num<<endl;
for(i=0; i<F2Longbar_num; i++){
for(j=0; j<5; j++)
OutInitData<<*(F2LbarInfo+5*i+j)<<" ";
OutInitData<<endl;
}
if(fabs(e0y)<=0.000001)
{
if(e0y>=0.0)
u0y=Coepx*ComLength; //柱在绕z轴方向的初始缺陷
if(e0y<0.0)
u0y=-Coepx*ComLength;
}
if(fabs(e0y)>=0.000001)
u0y=e0y;
if(fabs(e0z)<=0.000001)
{
if(e0z>=0.0)
u0z=Coepx*ComLength; //柱在绕z轴方向的初始缺陷
if(e0z<0.0)
u0z=-Coepx*ComLength;
}
if(fabs(e0z)>=0.000001)
u0z=e0z;
if(Sect_Type==1) //矩形截面几何中心
{
z0=Sectwid/2.0;
y0=Sectlen/2.0;
}
if(MeshID==0)
{
TNode_num=(Tdiv1+divProt*2+1)*(Tdiv2+divProt*2+1); //温度场分析时节点总数目
TConElem_num=(Tdiv1+divProt*2)*(Tdiv2+divProt*2); //温度场分析时单元总数目
F2Node_num=(div1+divProt*2+1)*(div2+divProt*2+1); //节点总数
F2ConElem_num=(div1+divProt*2)*(div2+divProt*2); //单元总数
Node_Coor=(double *)calloc(TNode_num*2,8);
F2Node_Coor=(double *)calloc(F2Node_num*2,8);
TCElemNode_Info=(long *)calloc(TConElem_num*4,4);
F2CElemNode_Info=(long *)calloc(F2ConElem_num*4,4);
for(i=0; i<TConElem_num; i++)
{
*(TCElemNode_Info+4*i+0)=i+int(i/(Tdiv1+divProt*2));
*(TCElemNode_Info+4*i+1)=i+int(i/(Tdiv1+divProt*2))+1;
*(TCElemNode_Info+4*i+2)=i+int(i/(Tdiv1+divProt*2))+(Tdiv1+divProt*2)+1;
*(TCElemNode_Info+4*i+3)=i+int(i/(Tdiv1+divProt*2))+(Tdiv1+divProt*2)+2;
}
for(i=0; i<Tdiv2+2*divProt+1; i++) //温度场坐标
{
for(j=0; j<Tdiv1+2*divProt+1; j++)//
{
//节点z坐标
if(j<divProt+1)
Node_Coor[2*i*(divProt*2+Tdiv1+1)+2*j]=Prot_thick*j/divProt;
if(j>divProt && j<Tdiv1+divProt+1)
Node_Coor[2*i*(divProt*2+Tdiv1+1)+2*j]=Prot_thick+(Sectwid-Prot_thick*2)*(j-divProt)/Tdiv1;
if(j>Tdiv1+divProt)
Node_Coor[2*i*(divProt*2+Tdiv1+1)+2*j]=Sectwid-Prot_thick+(j-divProt-Tdiv1)*Prot_thick/divProt;
//节点y坐标
if(i<divProt+1)
Node_Coor[2*i*(divProt*2+Tdiv1+1)+2*j+1]=Prot_thick*i/divProt;
if(i>divProt && i<Tdiv2+divProt+1)
Node_Coor[2*i*(divProt*2+Tdiv1+1)+2*j+1]=Prot_thick+(Sectlen-Prot_thick*2)*(i-divProt)/Tdiv2;
if(i>Tdiv2+divProt)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -