#include
#include
#include
float abstr(float); /*定义名为abstr 的函数，指定功能在代码最后*/ int sign(float); /*定义名为sign 的函数，指定功能在代码最后*/ void main(void)
{
float m[300],c[300],p[300],d[300]; /*定义一维单精度浮点型变量数组，数组长度为300，m:弯矩;c:曲率;p:外力;d:挠度*/
float mom[100],coc[100]; /*定义一维单精度浮点型变量数组，数组长度为300，mom:纵梁方向距梁端n*da处弯矩;coc:对应曲率*/
int i; /*定义为整型变量*/
for(i=0;i
c[i]=0.0;
p[i]=0.0;
d[i]=0.0;
}
for(i=0;i
{mom[i]=0.0;
coc[i]=0.0;
}
//*****Enter Data To Store In Input.dat*****
FILE *file1,*file2,*file3;
float fy,es,esh; /*定义钢筋屈服强度，钢筋弹性模量，钢筋的极限拉应变*/
float fc,fct; /*定义混凝土抗压强度，抗拉强度*/
float as1,as; /*定义抗压钢筋面积，抗拉钢筋面积*/
float l,a,b,h; /*定义梁跨长，作用点到左端距离，截面宽、高*/ float ao[2]; /*钢筋中心到梁顶距离*/
intsn,ln,st; /*定义截面划分条带数，a 长度上的分段数，钢筋型号*/
file1=fopen("input.dat","r"); /*从input.dat 中读取相应数据*/
fscanf(file1,"%f%f%f",&fy,&es,&esh);
fscanf(file1,"%f%f",&fc,&fct);
fscanf(file1,"%f%f",&as1,&as);
fscanf(file1,"%f%f%f%f",&l,&a,&h,&b);
fscanf(file1,"%f%f",&ao[0],&ao[1]);
fscanf(file1,"%d%d%d",&sn,&ln,&st);
//*****End of Inputing Data*****
float dc=0.0000002,de=0.00005,ee,em; /*定义曲率增量，应变增量，截面中间处应变，应变
的中间变量*/
float sf1=0.0,sf2=0.0,dsf; /*定义截面合力sf1，sf2，截面合力修正值dsf*/
float ffc=0.0,fs=0.0; /*定义混凝土合力，钢筋合力*/
float mi=0.0,mic,mis; /*定义某一条带的弯矩，混凝土截面内的合弯矩，截面钢筋合弯矩*/
float z,e,s,r; /*定义条带到截面中间的距离，应变，应力，反力*/ float eo,eu,cc=0.0; /*定义混凝土达到极限强度时的应变，破坏时的应变，cc 为计算挠度时曲率变量，只对cc 赋初值0.0*/
float lp,hh,hn,aas,etop; /*塑性铰长度的一半，hh 未定义，分成sn 个条带后每条带的高度，钢筋面积的中间变量，梁顶混凝土应变*/
float esy,da; /*定义钢筋的屈服应变，梁纵向弹性区段分段ln 后每段长度*/
float dd,dsn,dl; /*均为挠度中间变量*/
float mm,mo,co,dp,tan; /*弯矩中间变量，外力作用点处弯矩，对应曲率，
外力变化量，屈服点处弯矩与曲率的比值*/
intj,k,n,ii,jj;
int jmax1=500,jmax2=0,jmax3=0; /*定义jmax1，jmax2，jmax3三个整型变量并赋初值*/
esy=fy/es; /*计算钢筋的屈服应变*/
hn=h/sn; /*计算h 高度的截面所分成的条带数*/ eu=-0.004; /*定义混凝土破坏时的应变为-0.004*/
e=-0.002; /*给条带应变赋初值*/
ee=-0.0001; /*给截面中间处应变赋初值*/
j=0;
ii=0;
do /*do-while循环，用于将曲率与算得的截面和弯矩一一对应*/
{j++;
c[j]=c[j-1]+dc;
do /*do-while循环，计算某一曲率下的混凝土与钢
筋的合弯矩，并保证钢筋混凝土截面内合力为0*/
{ii++;
mic=0.0;
mis=0.0; /*给混凝土和钢筋的弯矩赋初值为0*/
//*****Calculate The Force Of Concrete*****
ffc=0.0; /*给混凝土合力赋初值为0*/
for(i=0;i
{z=h/2-i*hn-hn/2; /*求每一条带到截面中间的距离*/
e=ee+z*c[j]; /*以截面中间处应变为基准计算每一条带中心处应变，原代码左端为"ee" 有误*/
if(e>0.00015) s=0.0; /*if语句，用来判断每一条带中心处应变所在范围，从而选择对应公式计算该处应力*/
else if(e>0.0001) s=fct;
else if(e>0.0) s=(2*fct*e)/(e+0.0001);
else if(e>eo) s=-0.85*fc*(2*e/eo-(e/eo)*(e/eo));
else if(e>eu) s=-0.85*fc*(1-100*(eo-e));
else s=0.0;
ffc=ffc+s*b*(h/sn); /*将每一条带中心的混凝土受力叠加求合力*/
mic=mic+s*b*z*(h/sn); /*将每一条带中心的混凝土弯矩叠加求合弯矩*/ }
//*****Calculate The Force Of Steel*****
fs=0.0; /*钢筋合力初值为0.0*/
for(k=0;k
{
z=ao[k]-h/2; /*钢筋距上边缘的距离*/
e=ee+z*c[j]; /*截面曲率为c[j]时钢筋的应变，其中ee 为截面中心的应变，为
保证截面合力为0，会做调整*/
if(abstr(e)
else if(abstr(e)
变的符号确定应力符号的函数*/
else s=sign(e)*(fy+0.01*es*(abstr(e)-esh)); /*钢筋硬化后钢筋应力的计算*/ if(z
else aas=as; /*把受拉钢筋的截面面积赋给变量aas*/
if(abstr(e)>esy) dc=0.0000003; /*钢筋屈服后的曲率增量*/
fs=fs+s*aas; /*截面钢筋合力*/
mis=mis+s*aas*z; /*截面钢筋合弯矩*/
}
if(ii==1) /*从此处到计算mm 之前，是在计算判断赋给截面中心应变的增量是多少时，截面的合力接近0*/
{
sf1=ffc+fs; /*截面合力*/
ee=ee+de; /*第一次计算截面内力后给ee 一个增量，然后再次计算得到截面内力，
用两次的计算结果计算ee 的合理增量*/
em=de; /*把初定的变量de 赋给em ，参与后边的计算*/
}
else
{
sf2=ffc+fs; /*用ee+de计算的截面合力*/
dsf=sf2-sf1; /*两次计算的截面合力差*/
em=-sf2*em/dsf; /*计算合理的ee 增量，也就是使截面合力趋近0时的ee 增量*/ if(em==0.0) em=de; /*如果计算的增量为0, 则第一次赋给ee 的增量是合理的*/ ee=ee+em; /*计算合理的截面中心处的应变*/
sf1=sf2; /*将最后一次计算的合力sf2赋给上一次计算的合力sf1*/ }
mm=mic+mis; /*截面合弯矩*/
}while(abstr(sf2)>=0.1); /*当截面合力的误差大于0.1时重新再次计算截面的合力*/ if(e>esy)
{
if(j
}
etop=ee-(h/2-hn/2)*c[j]; /*顶部条带中心的最大应变*/
m[j]=mm; /*把合弯矩赋给数组m[j]*/
if(m[j]>=mi) /*如果弯矩一直增大，则把弯矩赋给m[i],用于下一次比较，同
时将j 赋给jmax2*/
{
mi=m[j];
jmax2=j;
}
}
while(etop>=1.2*eu); /*顶部应变的代数值不小于1.2*eu时，增加曲率，进行下一次计算*/
jmax3=j; /*曲率增加的总次数j 赋给jmax3*/
//*****Put Out The Answer of Moment And Curvature*****//
file2=fopen("out1.dat","w");
for(j=0;j
{
fprintf(file2,"%18.8g,%18.8g\n",c[j],m[j]);
}
fclose(file2);
if(st==1) lp=0.35*h; /*钢筋类型为1时，塑性铰一半的长度*/
if(st==2) lp=0.5*h; /*钢筋类型为2时，塑性铰一半的长度*/
tan=m[jmax1]/c[jmax1]; /*弯矩曲率图斜率*/
da=a/ln; /*长度为a 的梁分成ln 段后每段的长度*/
dc=0.00000022; /*计算荷载-挠度关系时的曲率增量*/
j=0;
do
{
cc=cc+dc;
j=j+1;
if(cc
{
i=-1;
do /*当满足条件i
{
i++;
if((c[i]c[i+1])) /*判断截面计算曲率在已经求出的曲率列表中的位置*/ mm=m[i]+((cc-c[i])*(m[i+1]-m[i]))/(c[i+1]-c[i]); /*假设弯矩曲率在相邻值间是线性变化的，
求该曲率对应的弯矩*/
}
while(i
段的关系*/
r=mm/a; /*计算支座反力*/
p[j]=(mm*l)/(a*(l-a)); /*根据合弯矩计算作用的集中荷载*/
for(n=1;n
{
mom[n]=r*n*da; /*当计算长度为n*da时，计算截面的弯矩*/ jj=-1;
do /*当计算长度为n*da时，计算截面的曲率*/ {
jj++;
if((m[jj]
}
while(jj
d[j]=d[j]+dd; /*挠度叠加后写到挠度数组中*/
}
}
else if(cc
{
dsn=0.0; /*定义dsn 初始等于0*/
i=-1;
do /*当满足条件i
if((c[i]
曲率对应的弯矩*/
}
while(i
于上升阶段*/
r=mm/a; /*计算支座反力*/
p[j]=(mm*l)/(a*(l-a)); /*根据合弯矩计算作用的集中荷载*/
da=(a-lp)/ln; /*减去塑性铰的一半以后a 长度上的每一小段的长度*/ for(n=1;n
{mom[n]=r*n*da; /*当计算长度为n*da时，计算截面的弯矩*/
for(jj=0;jj
{if((m[jj]
coc[n]=c[jj]+((mom[n]-m[jj])*(c[jj+1]-c[jj]))/(m[jj+1]-m[jj]); /*计算截面的曲率*/ }
dd=coc[n]*da*(n*da-da/2); /*共轭梁法计算计算截面的挠度*/
dsn=dsn+dd; /*挠度叠加后得到最后值即为除塑性铰区域的挠度*/ }
mo=r*a; /*截面弯矩*/
for(jj=jmax1;jj
{if((m[jj]
}
dl=co*lp*(a-lp/2); /*塑性铰区域的挠度，原代码为aa 有误*/
d[j]=dsn+dl; /*塑性铰区域和非塑性铰区域挠度的叠加*/
}
else
{
dsn=0.0; /*定义初始值为0*/
i=-1;
do /*计算到了下降段以后的弯矩挠度*/
{i++;
if((c[i]
while(i
r=mm/a;
p[j]=(mm*l)/(a*(1-a));
dp=p[j-1]-p[j];
da=(a-lp)/ln;
for(n=1;n
{mom[n]=dp*n*da; /*次阶段跟上阶段意思一样*/ coc[n]=mom[n]/tan;
dd=coc[n]*da*(n*da-da/2);
dsn=dsn+dd;
}
dl=dc*lp*(a-lp/2);
d[j]=d[j-1]-dsn+dl; /*计算此阶段的挠度总值*/
}
}while(cc
file3=fopen("out2.dat","w");/*以写的方式打开out2.dat*/
for(j=0;j
{fprintf(file3,"%15.4f,%15.4f\n",d[j],p[j]); /*打印d[j],p[j]*/
}
}
int sign(float num)
{
int q;
if(num
if(num>=0.0)q=1; /*定义sign*/ return q;
}
floatabstr(float x)
{
if(x>=0)
x=x;
else x=-x;
return x;
} /*定义abstr*/