% MP.m
% elastic-plastic moment-curvature relationship for a rectangular cross section
% elastic-plastic force-deflection relationship for a fixed-fixed beam
% http://www.duke.edu/~hpgavin/ce130/MP.m

b = 2.0;		% width of beam			cm	
d = 0.1524;		% depth of beam			cm
Sy = 23e3;		% yield stress			N/cm^2
E  = 20e6;		% elastic modulus 		N/cm^2
L = 15;			% length of beam		cm
I = b*d^3/12;		% bending moment of inertia	cm^4
S = b*d^2/6;		% elastic section modulus	cm^3
Z = b*d^2/4;		% plastic section modulus	cm^3

% --------- moment - curvature for a rectangular cross section

My = S*Sy;		% yield moment			N.cm
Mp = Z*Sy;		% plastic moment		N.cm

Y_curve = My/(E*I);					% yield curvature
curve   = [ Y_curve/20 : Y_curve/10 : 5*Y_curve ];	% curvatures


M = My * ( Z/S - 0.5*Y_curve^2 ./ curve.^2 );	% elastic-plastic moments
index = find(M < My);				% points where moments are elastic
M(index) = E*I*curve(index);			% elastic moments

figure(1)
Title1 =sprintf('\\phi_y=%6.4f;    M_y=%4.0f N.cm;    M_p=%4.0f N.cm', Y_curve,My,Mp);
subplot(2,1,1)
 plot( curve, M  )				% plot Moment-Curvature
 xlabel('curvature, 1/cm')
 ylabel('moment, N.cm')
 title(Title1);
 grid

% --------- force - deflection for a fixed-fixed beam

Fy = 2*My/L;					% yield force   (elastic limit)
Fp = 2*Mp/L;					% plastic force (ultimate)
Dy = Fy/(12*E*I/L^3);				% yield displacement	

delta_x = L/500;
x       = [ 0 : delta_x : L ];			% x-axis
pts = length(x);				% points in the x-axis

m = 1 - 2*x/L;					% moment diagram	

Mc = [ 0:0.02:0.98 0.999 0.9994 0.9997 0.9999:1e-5:0.99999 ]*Mp;

F = 0;					% Force
D = 0;					% Deflection

figure(2)
Title2 = sprintf('\\Delta_y=%5.3f cm;   D_y=%5.3f cm;    F_y=%3.1f N;    F_p=%3.1f N', Dy, Mp*L^2/(6*E*I), Fy,Fp);

for i = 1:length(Mc)

 M = Mc(i);

 curve = ( Y_curve/sqrt(2) ) ./ sqrt ( Z/S - M*abs(m)/My ) .* sign(m);
 index = find(abs(M*m) < My);		% points where moments are elastic
 curve(index) = M*m(index)/(E*I);	% elastic moments
 slope = cumtrapz(curve)*delta_x;
 displ = cumtrapz(slope)*delta_x;

figure(2)
 subplot(3,1,1)
   plot (x, curve );
   xlabel('x, cm');
   ylabel('curvature, 1/cm');
   axis([0 L -0.05 0.05])		% plot axis range
   grid
 
 subplot(3,1,2)
   plot (x, slope);
   xlabel('x, cm');
   ylabel('slope');
   axis([0 L 0 0.2])			% plot axis range
   grid 
 
 subplot(3,1,3)
   plot (x, displ);
   xlabel('x, cm');
   ylabel('deflection, cm');
   axis([0 L 0 3.5])			% plot axis range
   grid 

 drawnow
 
 pause(0.1)

 F(i) = 2*M/L;
 D(i) = max(displ);


figure(1)
subplot(2,1,2)
 plot(D,F)
 xlabel('deflection, cm')
 ylabel('force, N')
 axis([0 3 0 130])			% plot axis range
 title(Title2)
 grid 

end

% ---------------------------------------------------- MP