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