How to speed up the iteration of the for loop and make some further improvement on the below code?
1 回表示 (過去 30 日間)
古いコメントを表示
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
0 件のコメント
採用された回答
Alan Stevens
2020 年 9 月 10 日
The following, which makes use of fzero to find the value of c that makes F_total equal to zero, runs a lot faster. I've no idea if the results are sensible!
%%% section property
h = 600; %height
b = 300; %width
%%% concrete properties:
fcu = 30; %MPa
%%% reinforcing steel properties:
fy = 420; %MPa
Es = 200e3; %MPa
materialProps = [h,b,fcu,fy,Es];
%%% layers input = [area, effective depth, moment arm] all in mm:
steel1 = [2000, 560, h/2 - 570];
steel2 = [1000, 40, h/2 - 30];
Layers = [steel1; steel2];
%%% 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];
alphabeta = [alpha; beta];
N = length(ecu_list);
curvature_list = zeros(1, N);
moment_list = zeros(1, N);
c = h; % Initial guess
for ecu_idx = 1 : N
ec_top = ecu_list(ecu_idx);
c0 = c; % Use previous converged value as guess for next ecu_idx
c = fzero(@Ffn, c0,[],ec_top, ecu_idx,materialProps, Layers, alphabeta);
%%% 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
moment = m1 + m2 + mc; %kN.m
curvature = ec_top / c; %1/mm
moment_list(ecu_idx) = moment;
curvature_list(ecu_idx) = curvature;
end
plot(curvature_list, moment_list,'o');
drawnow;
function F_total = Ffn(c, ec_top, ecu_idx,materialProps, Layers,alphabeta)
h = materialProps(1);
b = materialProps(2);
fcu = materialProps(3);
fy = materialProps(4);
Es = materialProps(5);
steel1 = Layers(1,:);
steel2 = Layers(2,:);
alpha = alphabeta(1,:);
beta = alphabeta(2,:);
%%% 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
%%% 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
%%% concrete layer
Fc = alpha(ecu_idx) * fcu * beta(ecu_idx) * c * b * 1e-3;
%%% 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
end
10 件のコメント
Alan Stevens
2020 年 9 月 13 日
I guess so. I'm unfamiliar with the physics of this, so you will need to decide yourself!
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!