marlab绘制混凝土梁弯矩-曲率曲线

clear clc format long h=1000;b=600;h0=h-40; %长度单位均为 mm fc=26.8;ft=2.39;Es=200000;Ec=32500; s=0.00001:0.00004:0.008; %自变量为受压区边缘混凝土压应变，0-0.008，步长 0.00004 n=100;q=50; %混凝土划分为n 层,工字钢q层 m=zeros(1,200); %定义数组 for i=1:200 x(1)=0;x(2)=h0;x(3)=0.5*h0; for j=3:200 %如果迭代次数超过 200 次，即认为迭代发散，终止迭代 f1=0;f2=0;f3=0;f4=0;r1=0;r2=0;r3=0;r4=0; for aa=1:6 Xbar(aa)=40+(aa-1)*184; As(1)=3.14*25*25*6/4; As(2)=3.14*25*25*2/4; As(3)=3.14*25*25*2/4; As(4)=3.14*25*25*2/4; As(5)=3.14*25*25*2/4; As(6)=3.14*25*25*6/4; strain(aa)=(x(j)-Xbar(aa))*s(i)/x(j); %钢筋应变 stress(aa)=fun1(strain(aa)); %钢筋应力，fun1为钢筋的应力应变关系 f1=f1+ stress(aa)*As(aa); %钢筋合力 r1=r1+ stress(aa)*As(aa)*(x(j)-Xbar(aa)); %钢筋合力对中性轴的力矩 end for kk=1:n hs=h/n; %混凝土分层厚度 Xc(kk)=(kk-0.5)*hs; %混凝土分层坐标 cstrain(kk)=(x(j)-Xc(kk))*s(i)/x(j); %混凝土第k层应变 cstress(kk)=fun2(cstrain(kk)); %混凝土第 k 层应力。fun2 混凝土受压时的应力应变关系 f2=f2+ cstress(kk)*hs*b; %混凝土合力 r2=r2+ cstress(kk)*hs*b*(x(j)-Xc(kk)); %混凝土合力对中性轴的力矩 end for dd=1:2 Xyy(1)=150; Xyy(2)=850; Gas(1)=400*12; Gas(2)=400*12; %工字钢翼缘面积 gstrain(dd)=(x(j)-Xyy(dd))*s(i)/x(j); %工字钢翼缘应变 gstress(dd)=fun3(gstrain(dd)); %工字钢翼缘应力。fun3 为工字钢的应力应变关系 f3=f3+ gstress(dd)*Gas(dd); %工字钢翼缘合力 r3=r3+ gstress(dd)*Gas(dd)*(x(j)-Xyy(dd)); %工字钢翼缘合力对中性轴的力矩 end for ee=1:q ghs=700/q; %工字钢腹板分层厚度 Xfb(ee)=150+ghs*(ee-0.5); %工字钢腹板分层坐标 ystrain(ee)=(x(j)-Xfb(ee))*s(i)/x(j); %工字钢腹板应变 ystress(ee)=fun3(ystrain(ee)); %工字钢腹板应力。fun3 为自定义函数，钢筋的应力应变关系 f4=f4+ ystress(ee)*ghs*12; %工字钢腹板合力 r4=r4+ ystress(ee)*ghs*12*(x(j)-Xfb(ee)); %工字钢腹板合力对中性轴的力矩 end f=f1+f2+f3+f4 if f<=1&&f>=0 v(i)=s(i)/x(j); %曲率Φ m(i)=r1+r2+r3+r4 %整个截面的力对中性轴取矩 M break elseif f<0 %二分法 x(j-1)=max(x(j-2),x(j-1)); else x(j-1)=min(x(j-2),x(j-1)); end x(j+1)=(x(j-1)+x(j))/2; end end plot(1000*v,m/1000000);grid on;xlabel('φ 单位：1/m'); ylabel('M 单位：kN*m');title('M-φ 关系图')