How to speed up the iteration of the for loop and make some further improvement on the below code?

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

Answer 1

[Alan Stevens]

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

[Said Bolluk]

Well, thank you. It worked well, and the resulting graphs are alike. I appreciate your effort and time. Your improvements look impressive, but I am not sure what you changed. The mechanics of this new code is something new for me. Could you please describe it in more detail?

[Alan Stevens]

% The first part from here
    h = 600; %height
% to here
    moment_list = zeros(1, N);
% is basically just your original code,except that I've added 
    materialProps = [h,b,fcu,fy,Es];
    Layers = [steel1; steel2];
    alphabeta = [alpha; beta];
% to collect together data that will later be passed to a function.

% Then, within your loop
    for ecu_idx = 1 : N
        ...
    end
% I've inserted the fzero function.  This finds the value of a parameter (in this case 
% the parameter c) that makes the output of a function zero. The function here is called Ffn.
% It takes in an initial guess for c and works out the corresponding value of F_total.
% It also takes in several other parameters (materialProps etc) that it needs in order to 
% calculate F_total.  Fzero then adjusts c until it finds a value such that Ffn
% returns zero for F_total (well, zero, to within the built-in tolerance - which
% can be changed if desired).
% See help fzero for details of its syntax.
% All the calculations within the function are basically your original calculations, 
% except that I've inserted calculations to extract material properties etc so the
% subsequent calculations use the same variable names that you created.
% 
% fzero returns the value of c that makes F_total zero.
% Subsequent calculations are again yours.
%
% I also changed your plot command slightly to plot points rather than lines.
%
% I hope this helps.

[Said Bolluk]

It is a bit clear now. I need to work on it though? So, this function work only when the summation of forces is zero? Is there a similar way to find the c value when F_total is not zero? I have tried a lot of approaches, including the solve function, but could not succeed. I hope I don't bother you with these extra questions. Thanks again for the help.

[Alan Stevens]

Since fzero finds the value of c that makes F_total zero, in order to get it to find the value of c that makes F_zero equal some other number, say F1, then you would replace the line

F_total = Fs1 + Fs2 + Fc;

near the end of function Ffn with

F_total = Fs1 + Fs2 + Fc - F1;

[Said Bolluk]

Well, wise move. I studied a bit more about this fzero function, and understood what you did in the initial code. But, what is the [ ] when you call for the Ffn with the other parameters?

c = fzero(@Ffn, c0,[],ec_top, ecu_idx,materialProps, Layers, alphabeta);

Also, this approach is not valid for every strain value between zero and 4%. It works only for some particular strains -the ones in ecu_list- with their peculiar coefficients (alpha and beta). The most proper way is to integrate over the c for finding the actual strain value ec_i. Consequently, the after product Fc will be a function of c. However, there are two different equations for fc_i, being a function of ec_i. But, there might be two solutions since the results are too close. For example, let say that ec_i = 0.003. The result of the first equation is about 27 MPa, while that of the second is 24 MPa. So, how can this code decide for the right c value? I will upload a photo of these equations. Is there any way to make the code possible on finding the strain ec_i out of the c value, and then the real stress fc_i?

syms y;
ec_i = ec_top * (c - y) / c;
fc_i = fcu * ((2 * ec_i / ec0) - (ec_i / ec0) ^ 2); %for 0<ec_i<0.002
%%% or
fc_i = fcu * (1 - 0.15 * ((ec_i - ec0) / (ecu - ec0))); %for 0.002<ec_i<ecu
%%% then calculate the force
Fc = int(fc_i * b, y, 0, c) * 1e-3;

Thank you so much again for all the help. I hope I'm not asking too much.

[Alan Stevens]

The [] is a placeholder for an options parameter should you wish to change some default parameters, such as the absolute tolerance, for example. I'll have a look at your other comments later today.

[Alan Stevens]

If I've interpreted your material correctly (and there's no guarantee that I have!) then the following development leads to a value for your Fc integral that combines fc1 and fc2 appropriately. Then you just need to choose a value of c to get Fc directly. If you use the following you need to check it carefully for (1) the logic and (2) errors in the derivation.

[Said Bolluk]

Thank again so much for this huge effort. I have been trying to control the validity of your derivation. The final Fc function works for the strain values more than 0.002. For 0.002, the first part of the function -the ones with only the alpha and beta coefficients- gives the rigth results. However, the strain values below 0.002 do not match with this Fc function. I might be missing some numerical values. So, should I use only the first part of the equation (ones with alpha and beta) for the strains less than 0.002, and use all of the equation (including A and B) for the higher strain values?

[Alan Stevens]

I guess so. I'm unfamiliar with the physics of this, so you will need to decide yourself!

[Said Bolluk]

Thank you so much again for these great efforts. Your suggestions will answer the purpose for sure. I appreciate that.