弯矩-曲率曲线

clear 
clc 
tic;
h=500;b=250;h0=h-40;                          %长度单位均为mm 
fc=20.1;ft=-2.01;                             %混凝土受压受拉强度
cs0=0.00164;                                   %峰值压应变，按混规取值
cs1=-0.000095;                                 %峰值拉应变，按混规取值
Es=200000;Ec=30000;                           %钢筋、混凝土弹性模量
t=500;                                        %混凝土划分为t层
M=zeros（1501）;                               %定义数组
phi=0:0.0000001:0.00005;                      %给定曲率Φ
for k=1:501
s（1）=0;s（2）=0.004;s（3）=0.002;                 %受压区边缘混凝土压应变
for l=3:400
xc=s（l）/phi（k）;
N1=0;N2=0;M1=0;M2=0;
for  i=1:t 
        hs=h/t;                               %混凝土分层厚度                          
        yi（i）=（i-0.5）*hs;                     %混凝土分层坐标 
        cstrain（i）=（xc-yi（i））*s（l）/xc;        %混凝土第k层应变
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%混凝土应力应变全曲线，按混凝土结构设计规范建议的本构关系取值
        if cstrain（i）<0
          at=1.25;                            %C30受拉段对应下降段参数值，按混规取值
          pt=ft/（Ec*cs1）;
          v=cstrain（i）/cs1;
          if v<=1
          dt=1-pt*（1.2-0.2*v^5）;
          else
          dt=1-pt/（at*（v-1）^1.7+v）;
          end
          cstress（i）=（1-dt）*Ec*cstrain（i）;
        else cstrain（i）>=0
          ac=1.36;                            %C30受压段对应下降段参数值，按混规取值
          pc=fc/（Ec*cs0）;         
          m=Ec*cs0/（Ec*cs0-fc）;
          vv=cstrain（i）/cs0;
          if vv<=1
          dc=1-pc*m/（m-1+vv^m）;
          else
          d