【钢筋】预应力混凝土梁全过程分析（含弯矩曲率图）【含Matlab源码 4713期】.zip

clear; clf; %截面参数定义 B = 200; H = 400; Ap=354; Ep=195e3; Ec=4.43e4; As=1608; Asc=402; Npe=329200;%预应力筋有效预压力N A0=B*H; I0=B*H^3/12; e0=H/2-50; %开始进行全过程分析 sigma_p1=Npe/Ap; kesai_p1=sigma_p1/Ep; kesai_p2=(Npe/A0+Npe*e0^2/I0)/Ec; fai0=Npe*e0/(I0*Ec);%预应力作用下的初始曲率 fprintf('\n 开始进行预应力混凝土矩形梁全过程分析'); %定义计算的阶次 %以梁顶应变控制 ni=100; kesai_cu=0.006; deta_kesai=kesai_cu/ni; kesai_c=0:deta_kesai:kesai_cu; c0=0:1/ni:1; %计算X即中性轴位置 fprintf('\n 计算中性轴位置'); for i=1:ni+1 %初始化 if(i==1) cl=H; else cl=c0(i-1)-1; end fsum_old=0; c1_old=0; %中性轴判断循环 while(cl >=0&&cl<=H) %计算混凝土总合力 fc=0; for yci=0:1:H kesai_ci=(yci-cl)/cl*kesai_c(i); sigma_ci=calcu_sigma_c(kesai_ci); fc=fc+sigma_ci*1*B; end deta_p=kesai_c(i)*(H-cl-50)/cl; kesai_p=kesai_p1+kesai_p2+deta_p; sigma_p=calcu_sigma_p(kesai_p)-calcu_sigma_p(kesai_p1+kesai_p2); fp=sigma_p*Ap;%预应力钢筋水平力 kesai_s= kesai_c(i)*(H-cl-30)/cl; sigma_s=calcu_sigma_s(kesai_s); fs=sigma_s*As;%受拉钢筋水平力 kesai_sc= kesai_c(i)*(30-cl)/cl; sigma_sc=calcu_sigma_s(kesai_sc); fsc=sigma_sc*Asc;%受压钢筋水平力 %fsum=fc+fs+fsc+fp; fsum=fc+fs+fsc+fp; if(fsum==0) c0(i)=cl; fprintf('\n c0(%d)=%d',i,c0(i)); break; elseif(fsum*fsum_old<0) c0(i)=(cl+cl_old)/2; fprintf('\n c0(%d)=%d',i,c0(i)); break; else fsum_old=fsum; cl_old=cl; if(fsum>0) cl=cl+1; else cl=cl-1; end end end end fprintf('\n 计算弯矩'); %计算弯矩 for i=1:ni+1 Mc=0; for yci=0:1:H kesai_ci=(yci-c0(i))/c0(i)*kesai_c(i); sigma_ci=calcu_sigma_c(kesai_ci); Mc=Mc+sigma_ci*1*B*(yci-c0(i)); end deta_p=kesai_c(i)*(H-c0(i)-50)/c0(i); kesai_p=kesai_p1+kesai_p2+deta_p; sigma_p=calcu_sigma_p(kesai_p)-calcu_sigma_p(kesai_p1+kesai_p2); if(i==1) Mp=calcu_sigma_p(kesai_p1+kesai_p2)*Ap*(H/2-50)+Ap*sigma_p*(H-c0(i)-50); else Mp=calcu_sigma_p(kesai_p1+kesai_p2)*Ap*(H-c0(i)-50)+Ap*sigma_p*(H-c0(i)-50); end fprintf('\n Mp=%d ',Mp); kesai_s= kesai_c(i)*(H-c0(i)-30)/c0(i); sigma_s=calcu_sigma_s(kesai_s); Ms=sigma_s*As*(H-c0(i)-30); kesai_sc= kesai_c(i)*(30-c0(i))/c0(i); sigma_sc=calcu_sigma_s(kesai_sc); Msc=sigma_sc*Asc*(30-c0(i)); %Msum(i)=Mc+Mp+Ms+Msc; Msum(i)=Mc+Ms+Msc+Mp; Msum(i)=Msum(i)/1e6;%单位转换KN·m %kesai_picture=0:0.006/100:0.006; fai(i)=fai0+kesai_c(i)/c0(i); fprintf('\n c0(%d)=%d',i,c0(i)); %fprintf('\n kesai(%d)=%d',i,kesai_picture(i)); fprintf('\n r(%d)=%d',i,fai(i)); end fprintf('\n 绘制弯矩曲率曲线'); %绘制弯矩曲率曲线 %r=[0,r];Msum=[0,Msum]; a=0:fai0/2:fai0;m=0:Msum(1)/2:Msum(1); plot(kesai_c,Msum,'r','LineWidth',2); hold on; plot(a,m,'b','LineWidth',2); legend('矩形预应力混凝土梁截面弯矩-曲率曲线'); xlabel('曲率');ylabel('弯矩（kN·m）')