How to update the variables of the function after solving for a symbolic expression?

function [curvature, moment] = moment_curvature
    curvature_list = zeros(1, 1000000);
    moment_list = zeros(1, 1000000);
    h = 600;
    d = 570;
    b = 300;
    As = 800;
    fcd = 30 / 1.5;
    Ec = 12680 + 460 * fcd;
    Es = 200 * 1e3;
    ec0 = 2 * fcd / Ec;
    ecu = 0.0038;
    for ec_top = 0:5e-4:ecu
        syms c y;
        %%% tension zone: contribution from the bottom steel layers
        fs = Es * (c - d) * ec_top / c;
        Fs = fs * As * 1e-3;
        %%% compression zone: contribution from the upper concrete layers
        ec_i = ec_top * y / c;
        fc_i = fcd * ((2 * ec_i / ec0) - (ec_i / ec0) ^ 2);
        Fc = int(fc_i * b, y, 0, c) * 1e-3;
        F_total = Fs + Fc == 0;
        %%% the net force must be zero, to find c:
        c = solve(F_total,c,'Real',true); %%% I need the positive real root
        %%% after finding c, reassign the variables
        fs = subs(fs);
        Fs = subs(Fs);
        ec_i = subs(ec_i);
        fc_i = subs(fc_i);
        Fc = subs(Fc);
        %%% calculating the M and K values after finding c
        moment = Fs * (d - h / 2) * 1e-3 + int((fc_i * b * (h/2 - y)), y, 0, c) * 1e-3;
        curvature = ec_top / c;
        %%%
        curvature_list = [curvature_list, curvature];
        moment_list = [moment_list, moment];
        %%% plot the moment-curvature curve of the beam section
        plot(curvature_list, moment_list);
    end
end

Hello,

In this code, I am trying to plot the Moment-Curvature Diagram of a reinforced concrete beam section. The for loop is there to find the values of the diagram for the specific strain values. For each value, there will be different stresses, and consequently forces and moments. All these values are based on c, which is the neutral axis depth of the section. This depth determines the compression (upper) and tension (lower) zones of the beam. After assigning the values with c, I want to solve the final equation of F_total, which is the total force acting on the cross-section. It must be zero if there is no external normal force so that the c value should be defined according to this.

However, I have many errors in the code. To start with, I need to figure out how to find the real positive root, which gives the actual neutral axis depth c, of the function of F_total. When I run the code, I see this warning: Solution is not unique because the system is rank-deficient. I will attach the full error screen. I will be pleased if I could get any help.

Thank you,

Answer 1

[Walter Roberson]

Eventually you get two real roots for c. Your line dividing by c then has a problem because you used / rather than ./ and your attempt to append values onto the list has problems because you are appending two values in a column to the end of a list.

[Walter Roberson]

function [curvature_list, moment_list] = moment_curvature
    h = 600;
    d = 570;
    b = 300;
    As = 800;
    fcd = 30 / 1.5;
    Ec = 12680 + 460 * fcd;
    Es = 200 * 1e3;
    ec0 = 2 * fcd / Ec;
    ecu = 0.0038;
    ecu_list = 0:5e-4:ecu;
    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);
        syms c y;
        %%% tension zone: contribution from the bottom steel layers
        fs = Es * (c - d) * ec_top / c;
        Fs = fs * As * 1e-3;
        %%% compression zone: contribution from the upper concrete layers
        ec_i = ec_top * y / c;
        fc_i = fcd * ((2 * ec_i / ec0) - (ec_i / ec0) ^ 2);
        Fc = int(fc_i * b, y, 0, c) * 1e-3;
        F_total = Fs + Fc == 0;
        %%% the net force must be zero, to find c:
        c = solve(F_total,c,'Real',true); %%% I need the positive real root
        %%% after finding c, reassign the variables
        c(c<0) = [];
        if numel(c) ~= 1
            keyboard;
        end
        fs = subs(fs);
        Fs = subs(Fs);
        ec_i = subs(ec_i);
        fc_i = subs(fc_i);
        Fc = subs(Fc);
        %%% calculating the M and K values after finding c
        moment = Fs * (d - h / 2) * 1e-3 + int((fc_i * b * (h/2 - y)), y, 0, c) * 1e-3;
        curvature = ec_top / c;
        %%%
        curvature_list(ecu_idx) = curvature;
        moment_list(ecu_idx) = moment;
    end
    
    %%% plot the moment-curvature curve of the beam section
    plot(curvature_list, moment_list);
    drawnow;
end

[Said Bolluk]

Thank you for your time. I am not familiar with what you changed in the code. However, I get the idea. Solving for the c with the solve function might not be ideal since there are two different equations of the fc_i for specific ec_i values. The approximation is called the Hognestad stress-strain model, and I will attach a photo of it. In the above code, I assume that the strain (ec_i) is always less than 0.2%. The results stand in the linear elastic range as the stress (fc_i) is defined by the elastic strain. Naturally, the numbers would be incorrect when the e_c value reaches out more than 0.002.

So, I tried a different approach -iterating for the c and solve accordingly. In this one, the alpha and beta coefficients determine the behavior of the section under different strain values. Although they are close to being precise, integrating over the length of the beam would be the best solution. As I mentioned, it is not possible for me to perform a symbolic integration while there are two different stress formulas for particular strains. Also, this one takes too much time to return a result. Would you please view this one and maybe give some suggestions? Thank you again for your effort. I appreciate that, and as you can see, I already used some of your techniques from the previous one.

function moment_curvature
    h = 600;
    b = 300;
    %%% concrete properties:
    fcu = 20; %MPa
    %%% reinforcing steel properties:
    fy = 420; %MPa
    Es = 200e3; %MPa
    %%% steel layer properties = [area, d, y] all in mm:
    %%% d:working length, y:moment arm 
    steel1 = [1980, 560, h/2 - 570];
    steel2 = [1188, 40, h/2 - 30];
    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 = [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-6 : 2*h
            %%% 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 = Fs1 + Fs2 + Fc;
            if -0.001 < F_total && 0.001 > 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;
    K_ductility = curvature_list(9) / curvature_list(2);
    disp(K_ductility);
end

[Walter Roberson]

Use piecewise() to define symbolic formula that are valid only for particular range. Symbolic integration works for piecewise() definitions.