📄 3.txt
字号:
//Lshape 200*700*500
#include "iostream.h"
#include "iomanip.h"
#include "fstream.h"
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
const int M=70;
const int H=20;
const int N=50;
const int P=N-H;
const int Q=M-H;
const int C=3;
const double fc=30.0;
const double d=2.0;//纵筋直径
const double E=0.002;
const double B=16;
const double s=6;//箍筋间距
const double ES=200000;
const double e0=0.002;
const double pi=3.1415926;
const double F=fc*0.79*145.138;
//***********************
//约束混凝土区的本构关系
double strc(double e,double p, double s)
{
double a,a0,m,nn,l,z;
double e50u,e50h;
e50u=(3+0.003*1.2*F)/(1.2*F-1000);
e50h=0.75*p*sqrt(B/s);
z=0.5/(e50u+e50h-0.0030);
a0=F;
m=1.2*F*((-2*e)/(1.5*E)-(e*e)/(1.5*1.5*E*E));
nn=1.2*F*(1-z*(-e-0.003));
l=0.24*F;
if(e>=-0.002&&e<0) a=-m/171.44;
else if(e>0) a=0;
else
{
if(nn>l)
a=-nn/171.44;
else
a=-l/171.44;
}
return(a);
}
//***********************
//非约束保护层混凝土区的本构关系
double strc1(double e)
{
double a,a0;
a0=0.67*fc;
if(e>=-0.002&&e<=0) a=-a0*((-2*e)/e0-(e*e)/(e0*e0));
else if(e>=-0.004&&e<-0.002) a=-a0*(1-0.8*(fabs(e)-0.002)/0.002);
else a=0;
return(a);
}
//***************************
//钢筋本构关系
double strs(double e)
{
double as,fy;
fy=400;
if(e<fy/ES&&e>-fy/ES)
as=ES*e;
else if(e<=-fy/ES&&e>=-0.016)
as=-fy;
else if(e<-0.016) as=-(fy+0.01*ES*(fabs(e)-0.016));
else if(e>=fy/ES&&e<=0.016) as=fy;
else as=fy+0.01*ES*(e-0.016);
return(as);
}
//************************
//主程序
void main()
{
int i,j,k;
int kk,kkk;
int z0=0;
kk=1;kkk=1;
double fy=400.0;
k=0;
double dans[14][2],danc[2000][2];
//***********************
//形成钢筋单元坐标数组
dans[0][0]=C+d/2;dans[0][1]=C+d/2;
dans[1][0]=H-C-d/2;dans[1][1]=C+d/2;
dans[2][0]=C+d/2;dans[2][1]=H-C-d/2;
dans[3][0]=H-C-d/2;dans[3][1]=H-C-d/2;
dans[4][0]=M-C-d/2;dans[4][1]=C+d/2;
dans[5][0]=M-C-d/2;dans[5][1]=H-C-d/2;
dans[6][0]=C+d/2;dans[6][1]=N-C-d/2;
dans[7][0]=H-C-d/2;dans[7][1]=N-C-d/2;
dans[8][0]=H+Q/3-C;dans[8][1]=C+d/2;
dans[9][0]=H+Q/3-C;dans[9][1]=H-C-d/2;
dans[10][0]=H+2*Q/3-C;dans[10][1]=C+d/2;
dans[11][0]=H+2*Q/3-C;dans[11][1]=H-C-d/2;
dans[12][0]=C+d/2;dans[12][1]=H+P/2-C;
dans[13][0]=H-C-d/2;dans[13][1]=H+P/2-C;
for(i=0;i<14;i++)
cout<<dans[i][0]<<" "<<dans[i][1]<<" ";
//形成混凝土单元节点数组
//**********************1
for(i=0;i<M;i++)
{for(j=0;j<3;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
//**********************2
for(i=0;i<3;i++)
{for(j=3;j<N-3;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
//**********************3
for(i=0;i<H;i++)
{for(j=N-3;j<N;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
//**********************4
for(i=M-3;i<M;i++)
{for(j=3;j<H;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
//**********************5
for(i=H-3;i<M-3;i++)
{for(j=H-3;j<H;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
//**********************6
for(i=H-3;i<H;i++)
{for(j=H;j<N-3;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
z0=k;
//形成约束区内混凝土单元的坐标数组
for(i=3;i<M-3;i++)
{for(j=3;j<H-3;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
for(i=3;i<H-3;i++)
{for(j=H-3;j<N-3;j++)
{danc[k][0]=i+0.5;
danc[k][1]=j+0.5;
k++;}
}
//************************
//混凝土单元数组形成完毕
for(i=0;i<2000;i++)
cout<<danc[i][0]<<" "<<danc[i][1]<<" ";
//输入原始数据
double n,m,p,alf,as1,as2,x0,y0,d0;
double l1,l2;
double as3;
double rou;
double zhou=0.40;
n=zhou*fc*0.67*200000;
x0=27.5;y0=17.5;alf=pi*1.75;
d0=1.0;//箍筋直径
as1=100*d*d*pi/4.0;
as2=100*d0*d0*pi/4.0;
as3=12*12*pi/4.0;
l1=10*((H-2*C)*2+(M-2*C)*2);
l2=10*((H-2*C)*2+(N-2*C)*2);
p=as2*(l1+l2)/((M+N-H-2*C)*(H-2*C)*s*1000.0);
rou=(8*as1+4*as3)/180000;
double esb,sb;
sb=s/d;
esb=44200*pow(sb,-0.412)/1000000.0;
cout<<esb<<" ";
double r,ri,r1,dr;
double fai,dfai,sita,dsita,fai1,sita1;
double nc,ns1,nsc,np;
double mcx,mcy;
double msx,msy,mscx,mscy;
double msx1,msy1;
double mx,my;
double eic,eis[14];
//***************************
double eicm,eism;
int q,u,a,w,aa;
char quf;
int end=0;
double miu;
double shi;
int tt=0;
double mfai[10000][6],mmax=0;
double fac[10000],fag[10000];
double facg;
//***************************************
//strs1数组储存钢筋屈服时各钢筋单元的应变
//strs2数组储存第一批混凝土单元屈服时各钢筋单元的应变
//strs3数组储存极限状态个钢筋单元的应变
double strs1[14],strs2[14],strs3[14];
double *ang1,*ang2;
double *dist1,*dist2;
double anga,angb,dista,distb;
double angle,dist;
int number;
number=10000;
//申请动态内存
ang1=new double[number];
ang2=new double[number];
dist1=new double[number];
dist2=new double[number];
fai=10e-7;
dfai=2e-7;
sita1=alf;
dsita=0.0005;
r1=-5;
dr=0.005;
int num1=0;
do{
fai=fai+dfai;
q=1;k=0;u=0;
do{
sita=sita1+u*q*dsita;
a=0;w=1;aa=0;
do{
eicm=0;eism=0;
r=r1+aa*w*dr;
np=0;nc=0;
ns1=0;nsc=0;
mcx=0;mcy=0;
msx=0;msy=0;mscx=0;mscy=0;
msx1=0;msy1=0;
mx=0;my=0;
for(i=0;i<z0;i++)
{ri=r-danc[i][0]*cos(sita)-danc[i][1]*sin(sita);
eic=ri*fai*10.0;
if(eic<eicm)
eicm=eic;
{
nc=nc+strc1(eic)*100;
mcx=mcx+strc1(eic)*100.0*(danc[i][1]-y0)*10.0;
mcy=mcy+strc1(eic)*100.0*(danc[i][0]-x0)*10.0;
}
}
for(i=z0;i<2000;i++)
{ri=r-danc[i][0]*cos(sita)-danc[i][1]*sin(sita);
eic=ri*fai*10.0;
nc=nc+strc(eic,p,s)*100;
mcx=mcx+strc(eic,p,s)*100.0*(danc[i][1]-y0)*10.0;
mcy=mcy+strc(eic,p,s)*100.0*(danc[i][0]-x0)*10.0;
}
for(i=0;i<8;i++)
{
ri=r-dans[i][0]*cos(sita)-dans[i][1]*sin(sita);
eis[i]=fai*ri*10.0;
if(eis[i]>eism)
eism=eis[i];
ns1=ns1+strs(eis[i])*as1;
nsc=nsc+strc(eis[i],p,s)*as1;
msx1=msx1+strs(eis[i])*as1*(dans[i][1]-y0)*10.0;
msy1=msy1+strs(eis[i])*as1*(dans[i][0]-x0)*10.0;
mscx=mscx+strc(eis[i],p,s)*as1*(dans[i][1]-y0)*10.0;
mscy=mscy+strc(eis[i],p,s)*as1*(dans[i][0]-x0)*10.0;
}
for(i=8;i<14;i++)
{
ri=r-dans[i][0]*cos(sita)-dans[i][1]*sin(sita);
eis[i]=fai*ri*10.0;
if(eis[i]>eism)
eism=eis[i];
ns1=ns1+strs(eis[i])*as3;
nsc=nsc+strc(eis[i],p,s)*as3;
msx1=msx1+strs(eis[i])*as3*(dans[i][1]-y0)*10.0;
msy1=msy1+strs(eis[i])*as3*(dans[i][0]-x0)*10.0;
mscx=mscx+strc(eis[i],p,s)*as3*(dans[i][1]-y0)*10.0;
mscy=mscy+strc(eis[i],p,s)*as3*(dans[i][0]-x0)*10.0;
}
np=nc+ns1-nsc;
msx=msx1-mscx;
msy=msy1-mscy;
a++;
aa=(a+1)/2;
w=-w;
}while((fabs((np+n)/n))>=0.02);
r1=r;
k++;
u=(k+1)/2;
q=-q;
mx=mcx+msx;
my=mcy+msy;
shi=atan2(mx,my);
}while((shi-alf)>=0.005);
sita1=sita;
m=sqrt(mx*mx+my*my)/1000000.0;
fai1=1000*fai;
//**************************
//求混凝土达到0.0033时的曲率
for(i=0;i<2000;i++)
{
ri=r-danc[i][0]*cos(sita1)-danc[i][1]*sin(sita1);
eic=ri*fai*10.0;
if(eic<-0.0033)
{
fac[kkk]=fai;
ang2[kkk]=sita1;
dist2[kkk]=r;
cout<<eic<<" "<<danc[i][0]<<" "<<danc[i][1]<<" ";
kkk++;
break;
}
}
//***********************
//求钢筋达到屈服时的曲率
for(i=0;i<14;i++)
{
ri=r-dans[i][0]*cos(sita1)-dans[i][1]*sin(sita1);
eis[i]=ri*fai*10.0;
if(eis[i]>fy/ES)
{fag[kk]=fai;
ang1[kk]=sita;
dist1[kk]=r;
cout<<eis[i]<<" "<<dans[i][0]<<" "<<dans[i][1]<<" ";
kk++;
break;
}
}
//**********************
//计算纵筋压曲时的end值
for(i=0;i<14;i++)
{
ri=r-dans[i][0]*cos(sita1)-dans[i][1]*sin(sita1);
eis[i]=ri*fai*10.0;
if(eis[i]<=-esb)
{end=1;
cout<<end<<" "<<eis[i]<<" "<<dans[i][0]<<" "<<dans[i][1]<<" ";
break;}
}
tt++;
mfai[tt][0]=m;
mfai[tt][1]=fai1;
mfai[tt][2]=eicm;
mfai[tt][3]=eism;
mfai[tt][4]=r;
mfai[tt][5]=sita1;
mmax=5.0;
for(i=1;i<=tt;i++)
{
if(mfai[i][0]>mmax)
mmax=mfai[i][0];
}
} while((mfai[tt][0]>=0.7*mmax)&&(end!=1));
//*****************************
//跳出循环～～～～～～～～～～·
angle=sita1;
dist=r;
anga=ang1[1];
dista=dist1[1];
angb=ang2[1];
distb=dist2[1];
//释放堆内存空间
delete[]ang1;
delete[]ang2;
delete[]dist1;
delete[]dist2;
int num2,num3,num4;
num2=num3=num4=0;
int conc[2000],stel[14],stelp[14];
//*******************************
//计算混凝土压碎的单元，1表示压碎
for(i=0;i<2000;i++)
{ri=r-danc[i][0]*cos(sita1)-danc[i][1]*sin(sita1);
eic=ri*fai*10.0;
if(eic<-0.0033)
{
num2++;
conc[i]=1;
}
else
conc[i]=0;
}
//*****************************************
//计算达到极限曲率时钢筋屈服的单元，1表示屈服
for(i=0;i<14;i++)
{
ri=r-dans[i][0]*cos(sita1)-dans[i][1]*sin(sita1);
eis[i]=ri*fai*10.0;
//求解极限状态时各钢筋单元的应变
strs3[i]=eis[i];
if(eis[i]>fy/ES)
{
num3++;
stel[i]=1;
}
else
stel[i]=0;
//**********************************************
//计算达到极限曲率时纵筋压曲的钢筋单元，1表示压曲
if(eis[i]<-esb)
{
num4++;
stelp[i]=1;
}
else
stelp[i]=0;
}
//***************************
// 求解屈服曲率
if(fac[1]>=fag[1])
{facg=fag[1]*1000;
quf='g';
}
else
{facg=fac[1]*1000;
quf='h';
}
miu=fai1/facg;
//***********************************
//求解钢筋开始屈服时各钢筋单元的应变
for(i=0;i<14;i++)
{ri=dista-dans[i][0]*cos(anga)-dans[i][1]*sin(anga);
strs1[i]=fag[1]*ri*10.0;
}
//***********************************
//求解混凝土压应变达到0.0033时各钢筋单元的应变
for(i=0;i<14;i++)
{ri=distb-dans[i][0]*cos(angb)-dans[i][1]*sin(angb);
strs2[i]=fac[1]*ri*10.0;
}
ofstream out;
out.open("d10-60-315.dat");
out<<setw(10)<<"轴压比n＝" <<zhou<<endl;
out<<"esb="<<esb<<endl;
out<<"配箍率p＝" <<p<<endl;
out<<"配筋率ρ="<<rou<<endl;
out<<"荷载角alf="<<alf<<endl;
out<<"end="<<end<<endl;
out<<num2<<" "<<num3<<" "<<num4<<endl;
out<<"屈服曲率："<<facg<<" "<<"极限曲率："<<fai1<<"\n";
out<<"最大承载力："<<mmax<<"截面延性："<<miu<<"\n";
out<<"最先屈服单元是："<<quf<<endl;
out<<setw(15)<<'M'<<setw(15)<<"fai"<<endl;
//****************************************
for(i=1;i<=tt;i++)
out<<setw(15)<<mfai[i][0]<<setw(15)<<mfai[i][1]<<endl;
out<<setw(15)<<"eicm"<<setw(15)<<"eism"<<endl;
for(i=1;i<=tt;i++)
out<<setw(15)<<mfai[i][2]<<setw(15)<<mfai[i][3]<<endl;
out<<setw(15)<<'R'<<setw(15)<<"sita"<<endl;
for(i=1;i<=tt;i++)
out<<setw(15)<<mfai[i][4]<<setw(15)<<mfai[i][5]<<endl;
for(i=0;i<2000;i++)
{ out<<conc[i]<<" ";
if(i%10==0)
out<<endl;}
out<<endl;
for(i=0;i<14;i++)
out<<stel[i]<<" "<<stelp[i]<<endl;
//******************************
out<<"混凝土应变达到0.0033时的曲率为:"<<endl;
out<<fac[1]*1000<<endl;
out<<"R="<<distb<<endl;
out<<"中和轴法线角度："<<angb<<endl;
//**********************************
out<<"钢筋屈服时截面的曲率为:"<<endl;
out<<fag[1]*1000<<endl;
out<<"R="<<dista<<endl;
out<<"中和轴法线角度："<<anga<<endl;
//******************************
out<<"钢筋屈服时各钢筋单元的应变为:" <<endl;
for(i=0;i<14;i++)
out<<strs1[i]<<endl;
out<<"第一批混凝土单元屈服时各钢筋单元的应变为:" <<endl;
for(i=0;i<14;i++)
out<<strs2[i]<<endl;
out<<"极限破坏时各钢筋单元的应变为:" <<endl;
for(i=0;i<14;i++)
out<<strs3[i]<<endl;
out.close();
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -