Moment curvature relationship for RCC Column section
Reference: Reinforced concrete Structures" by Park and Paulay,
especially chapters 2 and 6.
Objective: To develop Moment curvature relationship for a given RCC column
section considering
- the stress-strain relation of concrete as per confinement provided
by the transverse hoops
- the stress-strain relationship of reinforcing steel, including the
strain hardening and softening
Tools: MATLAB software
Theory: Consider any RCC section like that shown in Figure1. A
rectangular section is shown in all figures. But, the principles remain
applicable for other shapes of cross-section
To obtain the moment-curvature
relationship for the above section, 1) we need to take a specific value of
strain in extreme compression fiber of concrete, 𝜀cm,, and also a depth for neutral axis, kd
(Figure 1). 2) Using
these, the strain distribution on the whole cross-section can be obtained by
linear interpolation (assuming applicability of Bernouli's principle- Ref. page
48). 3) For
each strain level (value), the stress in steel and in concrete are obtained
using the respective stress-strain curves (Figure 2 and 3). Note that the
portion of concrete that has spalled (strain > 0.004) should be deducted
while evaluating the concrete contribution to load and moment. 4) By
multiplying this stress with the relevant area on which it is applicable, the
stress resultant on that area is obtained. 5) Similarly,
the moment of resistance about the neutral axis of the cross-section can be
obtained by multiplying the normal stress-resultant with the corresponding
lever arm.6) By summing up these stress-resultants
(normal forces and moment) for the whole cross-section, the values of
axial resistance, P, and corresponding moment resistance, M, can be
obtained, for each combination of 𝜀cm,,and kd. 7) By performing 2-dimensional
interpolation on this data, the moment-curvature relationship for any given
axial load can be obtained.
Figure 2: Stress- strain relationship of concrete in
compression under confinement by transverse reinforcement
|
Figure 3: Stress- Strain relationship of steel
(assumed same under tension as well as in compression)
Results:
When we execute the above steps (refer to the MATLAB
code given below), we get a moment -curvature plot for a given section and
materials, under a given amount of axial load. The plots for various amounts of
lateral confinement (hoop reinforcements) can be superimposed to form a graph
like the one shown in figure 4. The input data for the figure is given in the
MATLAB code. This figure can be compared with the one published in Park and
Pauley (figure 6.20 on page 235)
|
MATLAB code:
function MomentCurvature()
% Data from PArk and Pauley Pg. 232
section.b=762; % Width of section
section.h=762; % Depth of section
section.cc=38; % effective cover
section.Z=52.6;
section.ro=0.026; %check compatibility with Reinf
materials.fy=276;
materials.fc=27.6;
P=0.3*(materials.fc)*(section.b)*(section.h);
section.area=pi/4*(section.b)^2;
section.d=((section.h/20):(section.h/20):(section.h)); section.d=mean([section.d;[0 section.d(1:end-1)]]);
section.b= section.b*sqrt(1-(1-2*(section.d)/(section.b)).^2); % constant width for a rectangular section.
[M,P,phi]=P_Moment_Curvature(section, materials,P);
return
% INput File ColumnSection.m
function [M,P,phi]=P_Moment_Curvature(section, materials,P)
b=section.b;h=section.h;cc=section.cc;fc=materials.fc;fy=materials.fy;
if isfield(section,'Z');Z=section.Z;else Z=CalZ(section);end
if isfield(section,'ro');ro=section.ro;else ro=Calro(section);section.ro=ro;end
ECM=logspace(-3,-1); Flag2=ones(size(ECM));
d=(h/20:h/20:h);NA=[d(1:end-1) h*logspace(0,1)];
Asj=SmearSteel(section,d);
figure(1);col=['r'; 'g'; 'b' ;'y';'m';'k'];
for i=1:length(ECM)
ecm=ECM(i);
for k=1:length(NA)
kd=NA(k);
Sj=0;Msj=0;Cj=0;Mcj=0;dat=[];FlagStageC=1;
[fcmax,FlagC]=ConfinedRCStressStrain(ecm,fc,Z);
for j=1:length(d)
if j==1;di=h/20;else di=d(j-1)+h/20;end
e=ecm/kd*(kd-di);
%% Steel contribution to P and M
[fsj,FlagS]=SteelStressStrain(e);
Sjj=fsj*Asj(j);
Sj=Sj+Sjj;
Msj=Msj+Sjj*(h/2-di);
%% COncrete Contribution
if e>0 %ecm <0.002 && e>0
[fcj,FlagC]=ConfinedRCStressStrain(e,fc,Z);
if j==1 && e >0.004
b1=0; % Ignoring the whole width beyond core
fcmax=0;
elseif e>0.004
b1=b(j)-2*cc; %Ignoring the width beyond core
else
b1=b(j);
end
if FlagStageC<FlagC
FlagStageC=FlagC;
end
Cjj=fcj*b1*h/20;
Cj=Cj+Cjj;
Mcj=Mcj+Cjj*(h/2-di);
else
fcj=0;
end
%%
dat=[dat;[di e fcj fsj FlagStageC]];
end
%%
Cc=Cj; Pi(i,k)=Cc+Sj;
Mi=Mcj+Msj;M(i,k)=Mi;
phi(i,k)=ecm/kd;
if i==50 && k==length(NA)-60
subplot(1,2,1);patch([0;fcmax;dat(:,3);0;0],[0; 0;dat(:,1);dat(end,1);0],'y');
set(gca,'YDir','reverse');xlabel('f_c_o_n_c_r_e_t_e (N/mm^2)');
if kd<h;line([min(dat(:,3)) max(dat(:,3))],[kd kd],'Linestyle','-.','color','r');end
subplot(1,2,2);patch([0;dat(1,4);dat(:,4);0;0],[0;0; dat(:,1);dat(end,1);0],'g');
set(gca,'YDir','reverse');xlabel('f_s_t_e_e_l (N/mm^2)');
if kd<h;line([min(dat(:,4)) max(dat(:,4))],[kd kd],'Linestyle','-.','color','r');end
end
%fprintf('\n %.4f %.2f %.0f %.0f %.0f',ecm,kd,Pi/1E3,Puo/1E3, M(i)/1E3);
end
if ecm>0.003
Flag2(i)=2;
elseif ecm>0.004
Flag2(i)=3;
elseif ecm/kd*(d-kd)>0.162
Flag2(i)=4;
elseif ecm/kd*(kd-cc)>0.162
Flag2(i)=5;
else
Flag2(i)=6;
end
i
end
[X,Y]=meshgrid(ECM,NA);figure(3)
[C,hand]=contour(X,Y,Pi',[P P]);
XI=C(1,2:(1+C(2,1)));YI=C(2,2:(1+C(2,1)));
Mii=interp2(X,Y,M',XI,YI);
figure(4);plot(XI./YI,Mii);
figure(5);surfc(X,Y,M');
return
function [fc1,FlagC]=ConfinedRCStressStrain(epc,fc,Z)
ep20c=0.002+0.8/Z; % Eq.2.7 Park and Pauley
if epc<0
fc1=0;
FlagC=0;
elseif epc<=0.002
fc1=fc*(2*epc/0.002-(epc/0.002)^2); % Equ 2.8 Park and pauley
FlagC=1;
elseif epc<=ep20c
fc1=fc*(1-Z*(epc-0.002)); % Eq.2.7 Park and Pauley
FlagC=2;
else
fc1=0.2*fc;
FlagC=3;
end
return
function [fsj,FlagS]=SteelStressStrain(e)
% ref Eq 6.28, 6.29 and 6.30 Park and Pauley 1975
fy=276;Es=2E5;ey=fy/Es;esh=16*ey;esu=0.14+esh;fsu=461;
r=esu-esh;
m=(fsu/fy*(30*r+1)^2-60*r-1)/15/r^2;
if length(e)==1
if abs(e)<=ey
fsj=e*Es;
FlagS=1;
elseif abs(e)<=esh
fsj=fy*sign(e);
FlagS=2;
else
fsj=sign(e)*fy*((m*(abs(e)-esh)+2)/(60*(abs(e)-esh)+2)+(abs(e)-esh)*(60-m)/2/(30*r+1)^2);
FlagS=3;
end
else
Indx1=find(abs(e)<=ey);
Indx2=find(abs(e)<=esh);
Indx3=find(abs(e)>esh);
fsj(Indx1)=e(Indx1)*Es;
fsj(Indx2)=fy*sign(e(Indx2));
fsj(Indx3)=sign(e(Indx3))*fy.*((m*(e(Indx3)-esh)+2)./(60*(e(Indx3)-esh)+2)+(e(Indx3)-esh)*(60-m)/2/(30*r+1)^2);
FlagS=0;
end
return
function Z=CalZ(section,materials)
bb=section.bb;sh=section.sh;ros=section.ros;fc=materials.fc;
b=section.b;h=section.h;cc=section.cc;diah=section.diah;
section.ros=((2*(b-2*cc)+2*(h-2*cc))*pi/4*(diah)^2)/((b-2*cc)*(h-2*cc)*sh);
eps50h=3/4*ros*sqrt(bb/sh) %Equ.2.8 Park & Pauley
Fc=fc*145.14 % to convert N/mm2 into psi
eps50u=(3+0.002*Fc)/(Fc-1000) %Eq 2.9 Park & Pauley
Z=0.5/(eps50u+eps50h-0.002) %Eq 2.8 Park & Pauley
return
function Z=Calro(section)
Reinf=section.Reinf; As=sum(Reinf(:,2));
ros=As/(section.b)/(section.h);
return
function Asj=SmearSteel(section,d)
if isfield(section,'Reinf')
Reinf=section.Reinf;
ds=Reinf(:,1);
for k=1:length(ds)
Ind(k)=find(d-ds(k)>0,1,'first');
end
Asj=zeros(size(d));%Asj(Ind)=Reinf(:,2);
Asj(2)=Reinf(1,2);Asj(end-1)=Reinf(end,2);
Asj(3:end-2)=sum(Reinf(2:end-1,2))/(length(d)-4);
elseif isfield(section,'ro')
SteelArea=(section.ro)*(section.area);
Asj=zeros(size(d));
Asj(2)=0.40*SteelArea;Asj(end-1)=0.40*SteelArea;
Asj(3:end-2)=0.2*SteelArea/(length(d)-4);
else
error('\n longitudinal steel reinforcement data is missing');
end
return
function Puo=PureAxialLoad_Column(section,materials,ecm)
%% determination of Pure axial load capacity => kd=infy
[fsj,FlagS]=SteelStressStrain(ecm);
As=(section.ros)*(section.b)*(section.h);Sj=fsj*As;
[fcj,FlagC]=ConfinedRCStressStrain(ecm,materials.fc,section.Z);
if ecm>0.004 % Concrete in cover would have spalled
bb=section.b-2*section.cc;
dd=section.h-2*section.cc;
else
bb=section.b;dd=section.h;
end
Cj=fcj*(bb*dd-As);
Puo=Sj+Cj;
return
% Data from PArk and Pauley Pg. 232
section.b=762; % Width of section
section.h=762; % Depth of section
section.cc=38; % effective cover
section.Z=52.6;
section.ro=0.026; %check compatibility with Reinf
materials.fy=276;
materials.fc=27.6;
P=0.3*(materials.fc)*(section.b)*(section.h);
section.area=pi/4*(section.b)^2;
section.d=((section.h/20):(section.h/20):(section.h)); section.d=mean([section.d;[0 section.d(1:end-1)]]);
section.b= section.b*sqrt(1-(1-2*(section.d)/(section.b)).^2); % constant width for a rectangular section.
[M,P,phi]=P_Moment_Curvature(section, materials,P);
return
% INput File ColumnSection.m
function [M,P,phi]=P_Moment_Curvature(section, materials,P)
b=section.b;h=section.h;cc=section.cc;fc=materials.fc;fy=materials.fy;
if isfield(section,'Z');Z=section.Z;else Z=CalZ(section);end
if isfield(section,'ro');ro=section.ro;else ro=Calro(section);section.ro=ro;end
ECM=logspace(-3,-1); Flag2=ones(size(ECM));
d=(h/20:h/20:h);NA=[d(1:end-1) h*logspace(0,1)];
Asj=SmearSteel(section,d);
figure(1);col=['r'; 'g'; 'b' ;'y';'m';'k'];
for i=1:length(ECM)
ecm=ECM(i);
for k=1:length(NA)
kd=NA(k);
Sj=0;Msj=0;Cj=0;Mcj=0;dat=[];FlagStageC=1;
[fcmax,FlagC]=ConfinedRCStressStrain(ecm,fc,Z);
for j=1:length(d)
if j==1;di=h/20;else di=d(j-1)+h/20;end
e=ecm/kd*(kd-di);
%% Steel contribution to P and M
[fsj,FlagS]=SteelStressStrain(e);
Sjj=fsj*Asj(j);
Sj=Sj+Sjj;
Msj=Msj+Sjj*(h/2-di);
%% COncrete Contribution
if e>0 %ecm <0.002 && e>0
[fcj,FlagC]=ConfinedRCStressStrain(e,fc,Z);
if j==1 && e >0.004
b1=0; % Ignoring the whole width beyond core
fcmax=0;
elseif e>0.004
b1=b(j)-2*cc; %Ignoring the width beyond core
else
b1=b(j);
end
if FlagStageC<FlagC
FlagStageC=FlagC;
end
Cjj=fcj*b1*h/20;
Cj=Cj+Cjj;
Mcj=Mcj+Cjj*(h/2-di);
else
fcj=0;
end
%%
dat=[dat;[di e fcj fsj FlagStageC]];
end
%%
Cc=Cj; Pi(i,k)=Cc+Sj;
Mi=Mcj+Msj;M(i,k)=Mi;
phi(i,k)=ecm/kd;
if i==50 && k==length(NA)-60
subplot(1,2,1);patch([0;fcmax;dat(:,3);0;0],[0; 0;dat(:,1);dat(end,1);0],'y');
set(gca,'YDir','reverse');xlabel('f_c_o_n_c_r_e_t_e (N/mm^2)');
if kd<h;line([min(dat(:,3)) max(dat(:,3))],[kd kd],'Linestyle','-.','color','r');end
subplot(1,2,2);patch([0;dat(1,4);dat(:,4);0;0],[0;0; dat(:,1);dat(end,1);0],'g');
set(gca,'YDir','reverse');xlabel('f_s_t_e_e_l (N/mm^2)');
if kd<h;line([min(dat(:,4)) max(dat(:,4))],[kd kd],'Linestyle','-.','color','r');end
end
%fprintf('\n %.4f %.2f %.0f %.0f %.0f',ecm,kd,Pi/1E3,Puo/1E3, M(i)/1E3);
end
if ecm>0.003
Flag2(i)=2;
elseif ecm>0.004
Flag2(i)=3;
elseif ecm/kd*(d-kd)>0.162
Flag2(i)=4;
elseif ecm/kd*(kd-cc)>0.162
Flag2(i)=5;
else
Flag2(i)=6;
end
i
end
[X,Y]=meshgrid(ECM,NA);figure(3)
[C,hand]=contour(X,Y,Pi',[P P]);
XI=C(1,2:(1+C(2,1)));YI=C(2,2:(1+C(2,1)));
Mii=interp2(X,Y,M',XI,YI);
figure(4);plot(XI./YI,Mii);
figure(5);surfc(X,Y,M');
return
function [fc1,FlagC]=ConfinedRCStressStrain(epc,fc,Z)
ep20c=0.002+0.8/Z; % Eq.2.7 Park and Pauley
if epc<0
fc1=0;
FlagC=0;
elseif epc<=0.002
fc1=fc*(2*epc/0.002-(epc/0.002)^2); % Equ 2.8 Park and pauley
FlagC=1;
elseif epc<=ep20c
fc1=fc*(1-Z*(epc-0.002)); % Eq.2.7 Park and Pauley
FlagC=2;
else
fc1=0.2*fc;
FlagC=3;
end
return
function [fsj,FlagS]=SteelStressStrain(e)
% ref Eq 6.28, 6.29 and 6.30 Park and Pauley 1975
fy=276;Es=2E5;ey=fy/Es;esh=16*ey;esu=0.14+esh;fsu=461;
r=esu-esh;
m=(fsu/fy*(30*r+1)^2-60*r-1)/15/r^2;
if length(e)==1
if abs(e)<=ey
fsj=e*Es;
FlagS=1;
elseif abs(e)<=esh
fsj=fy*sign(e);
FlagS=2;
else
fsj=sign(e)*fy*((m*(abs(e)-esh)+2)/(60*(abs(e)-esh)+2)+(abs(e)-esh)*(60-m)/2/(30*r+1)^2);
FlagS=3;
end
else
Indx1=find(abs(e)<=ey);
Indx2=find(abs(e)<=esh);
Indx3=find(abs(e)>esh);
fsj(Indx1)=e(Indx1)*Es;
fsj(Indx2)=fy*sign(e(Indx2));
fsj(Indx3)=sign(e(Indx3))*fy.*((m*(e(Indx3)-esh)+2)./(60*(e(Indx3)-esh)+2)+(e(Indx3)-esh)*(60-m)/2/(30*r+1)^2);
FlagS=0;
end
return
function Z=CalZ(section,materials)
bb=section.bb;sh=section.sh;ros=section.ros;fc=materials.fc;
b=section.b;h=section.h;cc=section.cc;diah=section.diah;
section.ros=((2*(b-2*cc)+2*(h-2*cc))*pi/4*(diah)^2)/((b-2*cc)*(h-2*cc)*sh);
eps50h=3/4*ros*sqrt(bb/sh) %Equ.2.8 Park & Pauley
Fc=fc*145.14 % to convert N/mm2 into psi
eps50u=(3+0.002*Fc)/(Fc-1000) %Eq 2.9 Park & Pauley
Z=0.5/(eps50u+eps50h-0.002) %Eq 2.8 Park & Pauley
return
function Z=Calro(section)
Reinf=section.Reinf; As=sum(Reinf(:,2));
ros=As/(section.b)/(section.h);
return
function Asj=SmearSteel(section,d)
if isfield(section,'Reinf')
Reinf=section.Reinf;
ds=Reinf(:,1);
for k=1:length(ds)
Ind(k)=find(d-ds(k)>0,1,'first');
end
Asj=zeros(size(d));%Asj(Ind)=Reinf(:,2);
Asj(2)=Reinf(1,2);Asj(end-1)=Reinf(end,2);
Asj(3:end-2)=sum(Reinf(2:end-1,2))/(length(d)-4);
elseif isfield(section,'ro')
SteelArea=(section.ro)*(section.area);
Asj=zeros(size(d));
Asj(2)=0.40*SteelArea;Asj(end-1)=0.40*SteelArea;
Asj(3:end-2)=0.2*SteelArea/(length(d)-4);
else
error('\n longitudinal steel reinforcement data is missing');
end
return
function Puo=PureAxialLoad_Column(section,materials,ecm)
%% determination of Pure axial load capacity => kd=infy
[fsj,FlagS]=SteelStressStrain(ecm);
As=(section.ros)*(section.b)*(section.h);Sj=fsj*As;
[fcj,FlagC]=ConfinedRCStressStrain(ecm,materials.fc,section.Z);
if ecm>0.004 % Concrete in cover would have spalled
bb=section.b-2*section.cc;
dd=section.h-2*section.cc;
else
bb=section.b;dd=section.h;
end
Cj=fcj*(bb*dd-As);
Puo=Sj+Cj;
return
