弯矩曲率计算程序

                     #include<stdio.h>

#include<math.h>

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<300;i++)                    /*为数组m，c，p，d，mom，coc赋初值0*/

{m[i]=0.0;

c[i]=0.0;

p[i]=0.0;

d[i]=0.0;

}

for(i=0;i<100;i++)

{mom[i]=0.0;

coc[i]=0.0;

}

//*****Enter Data To Store In Input.dat*****

FILE *file1,*file2,*file3,*file4;

float fy,es,esh;                      /*定义钢筋屈服强度，钢筋弹性模量，钢筋的极限拉应变*/

float fc,fct;                         /*定义混凝土抗压强度，抗拉强度*/

float as1,as;                         /*定义抗压钢筋面积，抗拉钢筋面积*/

float l,a,h,b;                        /*定义梁跨长，作用点到左端距离，截面宽、高*/

float ao[2];                          /*抗压和抗拉钢筋中心到梁顶距离*/

int sn,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);

fclose(file1);

//*****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;              /*塑性铰长度的一半，分成sn个条带后每条带的高度，钢筋面积的中间变量，梁顶混凝土应变*/

float esy,da;                         /*定义钢筋的屈服应变，梁纵向弹性区段分段ln后每段长度*/

float dd,dsn,dl;                      /*均为挠度中间变量*/

float mm,mo,co,dp,tan;                /*弯矩中间变量，外力作用点处弯矩，对应曲率

                     ，外力变化量，屈服点处弯矩与曲率的比值*/

int j,k,n,ii,jj;

int jmax1=500,jmax2=0,jmax3=0;        /*定义jmax1，jmax2，jmax3三个整型变量并赋初值*/

esy=fy/es;                            /*计算钢筋的屈服应变*/

hn=h/sn;                              /*计算h高度的截面所分成的sn条带时，每条带的高度hn*/

eu=-0.004;                            /*定义混凝土破坏时的应变为-0.004*/

eo=-0.002;                             /*给条带应变赋初值*/ 

ee=-0.0001;                           /*给截面中间处应变赋初值*/

j=0;

ii=0;

do                                    /*do-while循环，用于将曲率与计算截面的弯矩一一对应*/

{j++; ii=0;                           /*监测程序运行的整个过程发现ii增长的不健康，无法给sf1动态赋值，因而在这里加一个ii=0*/                               

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<sn;i++)                     /*for循环，用于各条带混凝土产生的合力与合弯矩*/

{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<2;k++)                      /*计算截面钢筋的合力和弯矩*/

{ 

z=ao[k]-h/2;                     /*钢筋距上边缘的距离*/

e=ee+z*c[j];                     /*截面曲率为c[j]时钢筋的应变，其中ee为截面中心的应变，为保证截面合力为0，会做调整*/

if(abstr(e)<=esy) s=e*es;        /*弹性阶段钢筋应力的计算*/

else if(abstr(e)<=esh)  s=sign(e)*fy;         /*屈服阶段钢筋应力的计算，其中sign是根据应变的符号确定应力符号的函数*/

else s=sign(e)*(fy+0.01*es*(abstr(e)-esh));   /*钢筋硬化后钢筋应力的计算*/

if(z<0.0) aas=as1;               /*把受压钢筋的截面

                     面积赋给变量aas*/

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(dsf==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<=jmax1) jmax1=j;              /*如果在曲率增加次数小于500并且底部受拉钢筋屈服时，将曲率增加次数赋给jmax1*/

}

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<=jmax3;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<=c[jmax1])                   /*截面计算曲率小于钢筋屈服前的曲率时，计算作用的荷载及相应截面的挠度*/

{

i=-1;                   

do                               /*当满足条件i<=jmax1+1时计算曲率对应的弯矩*/

{

i++;

if((c[i]<=cc)&&(cc<c[i+1]))   /*判断截面计算曲率在已经求出的曲率列表中的位置*/

mm=m[i]+((cc-c[i])*(m[i+1]-m[i]))/(c[i+1]-c[i]);   /*假设弯矩曲率在相邻值间是线性变化的，求该曲率对应的弯矩*/

} 

while(i<=jmax1+1);               /*此条件表示的意思是弯矩曲率关系要是处在弹性阶段的关系*/

r=mm/a;                          /*计算支座反力*/

p[j]=(mm*l)/(a*(l-a));           /*根据合弯矩计算作用的集中荷载*/

for(n=1;n<=ln;n++)               /*采用共轭梁法计算挠度*/

{

mom[n]=r*n*da;               /*当计算长度为n*da时，计算截面的弯矩*/

jj=-1;

do                           /*当计算长度为n*da时，计算截面的曲率*/

{

jj++;  

if((m[jj]<=mom[n])&&(mom[n]<m[jj+1]))      /*判断计算截面的弯矩在弯矩列表中的位置*/

coc[n]=c[jj]+((mom[n]-m[jj])*(c[jj+1]-c[jj]))/(m[jj+1]-m[jj]);    /*计算截面的曲率*/

}

while(jj<jmax1);                              /*此条件表示的意思是弯矩曲率关系要是处在弹性阶段的关系*/

dd=coc[n]*da*(n*da-da/2);                     /*共轭梁法计算计算截面的挠度*/

d[j]=d[j]+dd;                                 /*挠度叠加后写到挠度数组中*/

}

}

else if(cc<=c[jmax2])                 /*截面计算曲率小于钢筋最大应变曲率时，计算作用的荷载及相应截面的挠度*/

{

dsn=0.0;                            /*定义dsn初始等于0*/

i=-1;

do                                  /*当满足条件i<=jmax1+1时计算曲率对应的弯矩*/

{i++;

if((c[i]<=cc)&&(cc<c[i+1]))     /*判断截面计算曲率在已经求出的曲率列表中的位置*/

mm=m[i]+((cc-c[i])*(m[i+1]-m[i]))/(c[i+1]-c[i]);/*假设弯矩曲率在相邻值间是线性变化的，求该曲率对应的弯矩*/

}

while(i<=jmax2+1);                    /*此条件表示的意思是此时已经出现塑性铰，但是力还处于上升阶段*/

r=mm/a;                             /*计算支座反力*/

p[j]=(mm*l)/(a*(l-a));              /*根据合弯矩计算作用的集中荷载*/

da=(a-lp)/ln;                       /*减去塑性铰的一半以后a长度上的每一小段的长度*/

for(n=1;n<=ln;n++)                  /*采用共轭梁法计算挠度*/

{mom[n]=r*n*da;                  /*当计算长度为n*da时，计算截面的弯矩*/

for(j

                     j=0;jj<=jmax2;jj++)        /*计算除了塑性铰区域后的挠度*/

{if((m[jj]<=mom[n])&&(mom[n]<m[jj+1]))            /*判断计算截面的弯矩在弯矩列表中的位置*/

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<=jmax2;jj++)        /*计算塑性铰区域的挠度*/

{if((m[jj]<=mo)&&(mo<m[jj+1]))   /*判断计算截面的弯矩在弯矩列表中的位置*/

co=c[jj]+((mo-m[jj])*(c[jj+1]-c[jj]))/(m[jj+1]-m[jj]); /*计算截面的曲率*/

}

dl=co*lp*(a-lp/2);                /*塑性铰区域的挠度*/

d[j]=dsn+dl;                      /*塑性铰区域和非塑性铰区域挠度的叠加*/

}

else

{ 

dsn=0.0;                         /*定义初始值为0*/

i=-1;

do                               /*计算到了下降段以后的弯矩挠度*/

{i++;

if((c[i]<=cc)&&(cc<c[i+1])) /*判断计算截面的曲率在曲率列表中的位置*/

mm=m[i]+((cc-c[i])*(m[i+1]-m[i]))/(c[i+1]-c[i]);     /*计算截面的弯矩*/

}

while(i<jmax3);                  /*此条件表示的意思是混凝土已处于下降段*/

r=mm/a;        

p[j]=(mm*l)/(a*(1-a));

dp=p[j-1]-p[j];

da=(a-lp)/ln;

for(n=1;n<ln;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<=c[jmax3]);

file3=fopen("out2.dat","w");     /*以写的方式打开out2.dat*/

file4=fopen("out3.dat","w");

for(j=0;j<jmax3;j++)

{fprintf(file3,"%15.4f\\n",d[j]);         

fprintf(file4,"%15.4f\\n",p[j]);  /*打印d[j],p[j]在这里为了方便形成excel表格用了两个文件输出，文件3输出挠度，文件4输出外力，将数据复制到excel中绘制p-δ关系曲线*/

}

}

int sign(float num)

{  

int q;

if(num<0.0)q=-1;

if(num>=0.0)q=1;                   /*定义sign*/

return q;

}

float abstr(float x)

{ 

if(x>=0)

x=x;

else x=-x;                          /*定义abstr*/

return x;

}