Hello,
The moment_curvature function is to analyze the ductile behavior of a reinforced concrete member. There are some significant strain values, given in ecu_list, affecting the performance of the section. To calculate the curvature and moment values specific to these strains, I need to calculate the strains and stresses by finding the neutral axis depth, c. After that, I can calculate the forces. The summation of the internal reaction forces, which are provided by the steel and concrete layers, must be equal to the total sum of the external forces -zero (0) in this case.
So, my approach is to iterate over the c value and find the exact depth. As soon as the F_total closes to zero, I accepted that value of c and performed accordingly. However, these operations take too long and are a bit away from reality. I need the higher precision, but more importantly, I should speed up the running time of this code. Thus, I am open to any suggestions, and I will appreciate if I can get some help.
Thank you.
function moment_curvature %%% section property
h = 600; %height
b = 300; %width
%%% concrete properties:
fcu = 30; %MPa
%%% reinforcing steel properties:
fy = 420; %MPa Es = 200e3; %MPa %%% layers input = [area, effective depth, moment arm] all in mm:
steel1 = [2000, 560, h/2 - 570]; steel2 = [1000, 40, h/2 - 30]; %%% strain values
ecu_list = [2.5e-4, 5e-4, 1e-3, 1.5e-3, 2e-3, 2.5e-3, 3e-3, 3.5e-3, 4e-3]; %%% alpha&beta are the cooefficients specific to the given strain values
alpha = [0.178, 0.336, 0.595, 0.779, 0.889, 0.929, 0.934, 0.925, 0.910]; beta = [0.674, 0.682, 0.700, 0.722, 0.750, 0.785, 0.818, 0.849, 0.874]; N = length(ecu_list); curvature_list = zeros(1, N); moment_list = zeros(1, N); for ecu_idx = 1 : N ec_top = ecu_list(ecu_idx); for c = 0 : 1e-7 : 2*h %the neutral axis depth where there is no stress
%%% steel layer 1
es1 = (c - steel1(2)) * ec_top / c; %mm/mm
fs1 = Es * es1; %MPa if abs(fs1) > fy if fs1 < 0 fs1 = -fy; else fs1 = fy; end end Fs1 = fs1 * steel1(1) * 1e-3; %kN
m1 = Fs1 * steel1(3) * 1e-3; %kN.m
%%% steel layer 2
es2 = (c - steel2(2)) * ec_top / c; %mm/mm fs2 = Es * es2; %MPa if abs(fs2) > fy if fs2 < 0 fs2 = -fy; else fs2 = fy; end end Fs2 = fs2 * steel2(1) * 1e-3; %kN m2 = Fs2 * steel2(3) * 1e-3; %kN.m %%% concrete layer
Fc = alpha(ecu_idx) * fcu * beta(ecu_idx) * c * b * 1e-3; mc = Fc * (h/2 - beta(ecu_idx)*c/2) * 1e-3; %kN.m %%% F_total must be equal to the sum of external forces
F_total = Fs1 + Fs2 + Fc; %%% setting F_total to zero gives the neutral axis depth, c
if -0.0001 < F_total && 0.0001 > F_total break; end moment = m1 + m2 + mc; %kN.m curvature = ec_top / c; %1/mm
moment_list(ecu_idx) = moment; curvature_list(ecu_idx) = curvature; end end plot(curvature_list, moment_list); drawnow;end
Best Answer