My Matlab code

%%Program source code: calculation of interaction diagram for the resistance of CFEHS beam-column under combined axial compression and uniaxial bending numerically

%%Program source code: calculation of interaction diagram for the resistance of CFEHS beam-column under combined axial compression and uniaxial bending numerically




%%Basic properties of CFEHS
%Geometrical properties of CFEHS (mm)
%Member length
L=2000;
%Radii of major and minor axis of EHS
%Major axis bending should use b as the major axis dimension
as=220/2;
bs=110/2;
%Thickness of EHS
t=6.3;
%Radii of major and minor axis of concrete core
ac=as-t;
bc=bs-t;

%Area of concrete core and EHS
Ac=pi*ac*bc;
As=pi*as*bs-Ac;
%Number of reinforcing bars
n_r=4;
%Diameter of reinforcing bars
dia_r=10;
%Distance of reinforcing bars to the bending axis
ecc_r=bc-60;
%Area and plastic section modulus of reinforcing bars
Ar=n_r*pi*(dia_r/2)^2;
Sr=n_r*(pi*dia_r^3/32+pi*dia_r^2/4*ecc_r^2);
%Material properties of CFEHS
%Strength of concrete, steel tube and reinforcing bars (Mpa)
fck=30; fy=306; fr=370;
%Elastic modulus of concrete, steel tube and reinforcing bars (Mpa)
Ec=4700*sqrt(fck); Es=210000; Er=200000;


%% The generation of cross-sectional resistance curve
%The number of increments of shifting neutral axis
N=200;
%Generation of the space for the matrices
Ac1=zeros(floor(N*bc),1); Ac2=zeros(floor(N*bc),1); As1=zeros(floor(N*bc),1); As2=zeros(floor(N*bc),1); Ar1=zeros(floor(N*bc),1); Ar2=zeros(floor(N*bc),1); Sc1=zeros(floor(N*bc),1); Sc2=zeros(floor(N*bc),1); Ss1=zeros(floor(N*bc),1); Ss2=zeros(floor(N*bc),1); Sr1=zeros(floor(N*bc),1); Sr2=zeros(floor(N*bc),1);
%Creating a shifting neutral axis
%Finding the corresponding value of area and plastic section moduli
%The subroutine source code of 'uni_area' and 'uni_plastic_moduli' can be find in the latter file



i=1;
for l=0:1/N:bc-1/N

if l<ecc_r-dia_r/2
Ar1(i)=Ar/2; Sr1(i)=Sr/2;
else
Ar1(i)=Ar/2*(1-(l-(ecc_r-dia_r/2))/dia_r);
if Ar1(i)<0
Ar1(i)=0;

end
Sr1(i)=Sr/2*(1-(l-(ecc_r-dia_r/2))/dia_r);
if Sr1(i)<0
Sr1(i)=0;



end

end

end


Ar2(i)=Ar-Ar1(i); Ac1(i)=uni_area(ac,bc,l); Ac2(i)=Ac-Ac1(i); As1(i)=uni_area(as,bs,l)-Ac1(i); As2(i)=As-As1(i);
Sr2(i)=-Sr1(i);
Sc1(i)=uni_plastic_moduli(ac,bc,l); Sc2(i)=-Sc1(i); Ss1(i)=uni_plastic_moduli(as,bs,l)-Sc1(i); Ss2(i)=-Ss1(i);
i=i+1;


%Define compression in cross section is positive
%Condition 1, area 1 in compression, area 2 in tension
N1=Ac1*fck+As1*fy+Ar1*fr-As2*fy-Ar2*fr; M1=Sc1*fck+Ss1*fy+Sr1*fr-Ss2*fy-Sr2*fr;
%Condition 2, area 2 in compression, area 1 in tension
N2=Ac2*fck+As2*fy+Ar2*fr-As1*fy-Ar1*fr; M2=Sc2*fck+Ss2*fy+Sr2*fr-Ss1*fy-Sr1*fr;
%By combining Condition 1 and 2, the positions of neutral axis cover the whole cross-section
%According to Eurocode 4 the design moment is 0.9*M
NM=[N1,0.9*M1;N2,-0.9*M2];
NMu=unique(NM,'rows');
%MNu is the M-N sets representing the cross-sectional resistance curve
MNu=[NMu(:,2),NMu(:,1)];




%%Generation of load reaction curve
%Loading eccentricity
e=100;
%Initial imperfection
w=L/1000;
%Second moment of area of concrete, steel tube and reinforcement
Icx=pi/4*ac*bc^3; Isx=pi/4*as*bs^3-Icx; Irx=(pi*(dia_r/2)^4/4+pi*(dia_r/2)^2*ecc_r^2)*4;
%Effective Second moment of area of the composite cross-section
EIeffx=0.9*(Es*Isx+0.5*Ec*Icx+Er*Irx);
%Elastic critical resistance
Ncrx=(pi^2*EIeffx)/L^2;
%The vector of normal force
Ned=0:100:max(N2);
%Amplification factor for second-order effects
k=1.1*((1-Ned/Ncrx).^(-1));
%The vector of moment regarding the normal force
Mxcol=k.*Ned*(e+w);
%The curve of loading
reaction=[Mxcol',Ned'];
%Normalise the moment and axial force
MNu_normal=[MNu(:,1)/max(MNu(:,1)) MNu(:,2)/max(MNu(:,2))]; reaction_normal=[reaction(:,1)/max(MNu(:,1)) reaction(:,2)/max(MNu(:,2))];
%Find the intersection point (closest point) of the resistance andloading curve
%The subroutine source code of ‘close_point’can be find in the latter file
[row,col]=close_point(reaction_normal,MNu_normal);
%The design resistance moment and normal force of CFEHS beam-column
intersection=MNu(col,:);


%%Figures
%Interaction curves diagram figure
plot(MNu(:,1),MNu(:,2),'r-'); xlabel('M');
ylabel('N'); hold on plot(reaction(:,1),reaction(:,2),'b-');



hold on
%Plot(intersection(:,1),intersection(:,2),'ko'); plot(intersection(1,1),intersection(1,2),'ko') hold off
fprintf('with eCC %g and length %g, the load is %g \n',e,L,intersection(1,2));
%Neutral axis
li=find(intersection(1,2)==N2)/N; figure
%The subroutine for drawing the ellipse is given in latter script
draw_ellipse(as,bs);
hold on
draw_ellipse(ac,bc); hline = refline(0, li); hline.Color = 'r'; stra=num2str(li); str1={'li=' stra}; text(0,50,str1);
text(40,-50,'analytical uniaxial method')
hold off


The script of function for calculation of the area of the elliptical area is given here.
%Find the elliptical sector area of uniaxial bending
    function [area]=uni_area(a,b,l)
x=a*sqrt(1-l^2/b^2);
theta=atan(l/x);
area=pi*a*b/2-l*x-(a*b)*(theta-atan(((b-a)*sin(2*theta))/((b+a)+(b- a)*cos(2*theta))));


%The script of function for calculation of the plastic moduli of the elliptical area is given here.


%Find the elliptical sector plastic_moduli of uniaxial bending using analytical method


function [first_moment]=uni_plastic_moduli(a,b,l)


first_moment=(2/3)*(a/b)*(b^2-l^2)^(3/2);


%The script of function for finding the closest points in two matrices is given here.




%Find the colsest point in matrix A and B
function[row,col]=close_point(A,B)


a=size(A,1); b=size(B,1); dis=zeros(a,b);


for i=1:a for j=1:b
dis(i,j)=sqrt((A(i,1)-B(j,1))^2+(A(i,2)-B(j,2))^2);



end


end



minMatrix = min(dis(:));
[row,col] = find(dis==minMatrix);


%The script of function for drawing ellipses is given here.
%Draw ellipse with radii of a and b
    function[]=draw_ellipse(a,b)
t = linspace(0,2*pi); X = a*cos(t);
Y = b*sin(t);
plot(X,Y,'b-')

circle = rsmak('circle');fnplt(circle)
circle = rsmak('circle');fnplt(circle)
axis equal