Moment curvature relationship for RCC Column section

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 1: Section and strain profile at two different stages of loading

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)

Figure 4: Moment -curvature relation for the section

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