Unkownn

My Matlab code

Mar 7th, 2021 (edited)
587
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %%Program source code: calculation of interaction diagram for the resistance of CFEHS beam-column under combined axial compression and uniaxial bending numerically
  2.  
  3. %%Program source code: calculation of interaction diagram for the resistance of CFEHS beam-column under combined axial compression and uniaxial bending numerically
  4.  
  5.  
  6.  
  7.  
  8. %%Basic properties of CFEHS
  9. %Geometrical properties of CFEHS (mm)
  10. %Member length
  11. L=2000;
  12. %Radii of major and minor axis of EHS
  13. %Major axis bending should use b as the major axis dimension
  14. as=220/2;
  15. bs=110/2;
  16. %Thickness of EHS
  17. t=6.3;
  18. %Radii of major and minor axis of concrete core
  19. ac=as-t;
  20. bc=bs-t;
  21.  
  22. %Area of concrete core and EHS
  23. Ac=pi*ac*bc;
  24. As=pi*as*bs-Ac;
  25. %Number of reinforcing bars
  26. n_r=4;
  27. %Diameter of reinforcing bars
  28. dia_r=10;
  29. %Distance of reinforcing bars to the bending axis
  30. ecc_r=bc-60;
  31. %Area and plastic section modulus of reinforcing bars
  32. Ar=n_r*pi*(dia_r/2)^2;
  33. Sr=n_r*(pi*dia_r^3/32+pi*dia_r^2/4*ecc_r^2);
  34. %Material properties of CFEHS
  35. %Strength of concrete, steel tube and reinforcing bars (Mpa)
  36. fck=30; fy=306; fr=370;
  37. %Elastic modulus of concrete, steel tube and reinforcing bars (Mpa)
  38. Ec=4700*sqrt(fck); Es=210000; Er=200000;
  39.  
  40.  
  41. %% The generation of cross-sectional resistance curve
  42. %The number of increments of shifting neutral axis
  43. N=200;
  44. %Generation of the space for the matrices
  45. 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);
  46. %Creating a shifting neutral axis
  47. %Finding the corresponding value of area and plastic section moduli
  48. %The subroutine source code of 'uni_area' and 'uni_plastic_moduli' can be find in the latter file
  49.  
  50.  
  51.  
  52. i=1;
  53. for l=0:1/N:bc-1/N
  54.  
  55. if l<ecc_r-dia_r/2
  56. Ar1(i)=Ar/2; Sr1(i)=Sr/2;
  57. else
  58. Ar1(i)=Ar/2*(1-(l-(ecc_r-dia_r/2))/dia_r);
  59. if Ar1(i)<0
  60. Ar1(i)=0;
  61.  
  62. end
  63. Sr1(i)=Sr/2*(1-(l-(ecc_r-dia_r/2))/dia_r);
  64. if Sr1(i)<0
  65. Sr1(i)=0;
  66.  
  67.  
  68.  
  69. end
  70.  
  71. end
  72.  
  73. end
  74.  
  75.  
  76. 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);
  77. Sr2(i)=-Sr1(i);
  78. 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);
  79. i=i+1;
  80.  
  81.  
  82. %Define compression in cross section is positive
  83. %Condition 1, area 1 in compression, area 2 in tension
  84. N1=Ac1*fck+As1*fy+Ar1*fr-As2*fy-Ar2*fr; M1=Sc1*fck+Ss1*fy+Sr1*fr-Ss2*fy-Sr2*fr;
  85. %Condition 2, area 2 in compression, area 1 in tension
  86. N2=Ac2*fck+As2*fy+Ar2*fr-As1*fy-Ar1*fr; M2=Sc2*fck+Ss2*fy+Sr2*fr-Ss1*fy-Sr1*fr;
  87. %By combining Condition 1 and 2, the positions of neutral axis cover the whole cross-section
  88. %According to Eurocode 4 the design moment is 0.9*M
  89. NM=[N1,0.9*M1;N2,-0.9*M2];
  90. NMu=unique(NM,'rows');
  91. %MNu is the M-N sets representing the cross-sectional resistance curve
  92. MNu=[NMu(:,2),NMu(:,1)];
  93.  
  94.  
  95.  
  96.  
  97. %%Generation of load reaction curve
  98. %Loading eccentricity
  99. e=100;
  100. %Initial imperfection
  101. w=L/1000;
  102. %Second moment of area of concrete, steel tube and reinforcement
  103. 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;
  104. %Effective Second moment of area of the composite cross-section
  105. EIeffx=0.9*(Es*Isx+0.5*Ec*Icx+Er*Irx);
  106. %Elastic critical resistance
  107. Ncrx=(pi^2*EIeffx)/L^2;
  108. %The vector of normal force
  109. Ned=0:100:max(N2);
  110. %Amplification factor for second-order effects
  111. k=1.1*((1-Ned/Ncrx).^(-1));
  112. %The vector of moment regarding the normal force
  113. Mxcol=k.*Ned*(e+w);
  114. %The curve of loading
  115. reaction=[Mxcol',Ned'];
  116. %Normalise the moment and axial force
  117. MNu_normal=[MNu(:,1)/max(MNu(:,1)) MNu(:,2)/max(MNu(:,2))]; reaction_normal=[reaction(:,1)/max(MNu(:,1)) reaction(:,2)/max(MNu(:,2))];
  118. %Find the intersection point (closest point) of the resistance andloading curve
  119. %The subroutine source code of ‘close_point’can be find in the latter file
  120. [row,col]=close_point(reaction_normal,MNu_normal);
  121. %The design resistance moment and normal force of CFEHS beam-column
  122. intersection=MNu(col,:);
  123.  
  124.  
  125. %%Figures
  126. %Interaction curves diagram figure
  127. plot(MNu(:,1),MNu(:,2),'r-'); xlabel('M');
  128. ylabel('N'); hold on plot(reaction(:,1),reaction(:,2),'b-');
  129.  
  130.  
  131.  
  132. hold on
  133. %Plot(intersection(:,1),intersection(:,2),'ko'); plot(intersection(1,1),intersection(1,2),'ko') hold off
  134. fprintf('with eCC %g and length %g, the load is %g \n',e,L,intersection(1,2));
  135. %Neutral axis
  136. li=find(intersection(1,2)==N2)/N; figure
  137. %The subroutine for drawing the ellipse is given in latter script
  138. draw_ellipse(as,bs);
  139. hold on
  140. draw_ellipse(ac,bc); hline = refline(0, li); hline.Color = 'r'; stra=num2str(li); str1={'li=' stra}; text(0,50,str1);
  141. text(40,-50,'analytical uniaxial method')
  142. hold off
  143.  
  144.  
  145. The script of function for calculation of the area of the elliptical area is given here.
  146. %Find the elliptical sector area of uniaxial bending
  147.     function [area]=uni_area(a,b,l)
  148. x=a*sqrt(1-l^2/b^2);
  149. theta=atan(l/x);
  150. area=pi*a*b/2-l*x-(a*b)*(theta-atan(((b-a)*sin(2*theta))/((b+a)+(b- a)*cos(2*theta))));
  151.  
  152.  
  153. %The script of function for calculation of the plastic moduli of the elliptical area is given here.
  154.  
  155.  
  156. %Find the elliptical sector plastic_moduli of uniaxial bending using analytical method
  157.  
  158.  
  159. function [first_moment]=uni_plastic_moduli(a,b,l)
  160.  
  161.  
  162. first_moment=(2/3)*(a/b)*(b^2-l^2)^(3/2);
  163.  
  164.  
  165. %The script of function for finding the closest points in two matrices is given here.
  166.  
  167.  
  168.  
  169.  
  170. %Find the colsest point in matrix A and B
  171. function[row,col]=close_point(A,B)
  172.  
  173.  
  174. a=size(A,1); b=size(B,1); dis=zeros(a,b);
  175.  
  176.  
  177. for i=1:a for j=1:b
  178. dis(i,j)=sqrt((A(i,1)-B(j,1))^2+(A(i,2)-B(j,2))^2);
  179.  
  180.  
  181.  
  182. end
  183.  
  184.  
  185. end
  186.  
  187.  
  188.  
  189. minMatrix = min(dis(:));
  190. [row,col] = find(dis==minMatrix);
  191.  
  192.  
  193. %The script of function for drawing ellipses is given here.
  194. %Draw ellipse with radii of a and b
  195.     function[]=draw_ellipse(a,b)
  196. t = linspace(0,2*pi); X = a*cos(t);
  197. Y = b*sin(t);
  198. plot(X,Y,'b-')
  199.  
  200. circle = rsmak('circle');fnplt(circle)
  201. circle = rsmak('circle');fnplt(circle)
  202. axis equal
  203.  
  204.  
  205.  
RAW Paste Data