Undefined function or variable 'sigma_sj'. Error in m_phi_code (line 110) Ps_j = sigma_sj*(​N/2)*(3.14​/4)*d^2*0.​001; why this error is showing?

[SL: edited to format code]

%% M-Phi Curve
clc;
clear all;
close all;

% Input data
fck = 20;           % Grade of concrete in N/mm^2
Fe = 415;           % Grade of steel in N/mm^2
P = 300;            % Axial compressive load in kN
% Phi = 0.00004;    % Initial slope in per meter
N = 6;              % Number of bars
D = 500;            % Depth of the section in mm
B = 250;            % Width of the section in mm
e_c = 50;           % Effective cover in mm
d = 16;             % Diameter of the bars in mm
Phi = 0.00001:0.00001:0.00004;
k = length(Phi);
M_total = zeros(1,length(k));
Phi = zeros(1,k);

%% Cross sectional area of the section
A = B*D;            % Rectangular section Area in mm^2

%% Calculate strain under only axial compressive force
sigma_x = (P*1000)/(A);    % Stress under axial compressive force in N/mm^2

% Quadratic equation solving
c = sigma_x/(0.45*fck);
a = 1;
b= -2;
delta = (b.^2)-(4*a*c);
if(delta < 0)      
    disp("Delta < 0 The equation does not have a real root");
elseif (delta == 0)
    disp('The equation has only one real roots');
    disp(-b./(2*a));   
else
    disp('The equation has two real roots');
    epsilon_1 = 0.002*((-b+sqrt(delta))./(2*a));
    epsilon_2 = 0.002*((-b-sqrt(delta))./(2*a));
end
epsilon = min(epsilon_1,epsilon_2);
disp(epsilon)

%% Finding depth of neutral axis
epsilon_max = input('Enter maximum stain value :');
for v = 1:k
    Curvature = Phi(1,v); 
    epsilon_min = epsilon_max-(D*Curvature*0.001);
    xu = (epsilon_max*D)/(epsilon_max-epsilon_min);     % Depth of neutral axis in mm
    if(xu > D)
        disp("Neutral axis outside the section");
    elseif(xu==D)
        disp("Neutral axis is at the edge of the section");
    elseif(0<=xu)&&(xu<=D)
        disp("Neutral axis inside the section");
    else
        disp("Neutral axis is above the section");
    end
    
    %% Concrete stress, strain and force calculation
    n = input('please enter number of concrete strips :');
    Pc_0 = 0;
    Mc_0 = 0;
    for i = 0:n-1
        for j = i
            p_j = (D/n)*j + (D/(2*n));
            epsilon_ci = (epsilon_max/xu)*(xu-p_j);
            sigma_ci = 0.45*fck*(2*(epsilon_ci/0.002)-(epsilon_ci/0.002)^2);
            P_ci = sigma_ci*B*(D/n)*0.001;
            if (1 <= (i+1))&&((i+1) <= n/2)
                M_ci = P_ci*(((D/2)-((D/(2*n))*(i+1)))/(i+1))*0.001; 
            elseif ((i+1) >= ((n/2)+1))
                M_ci = P_ci*((((D/2)*(i-2))-(D/(2*n))))*0.001*(-1)^i;
            end
        end
        Pc_total = Pc_0 + P_ci;
        Pc_0 = Pc_total;
        Mc_total = Mc_0 + M_ci;
        Mc_0 = Mc_total;
    end
    disp(Pc_total)
    disp(Mc_total)
    
    %% Steel stress, strain and force calculation
    s = input('please enter number of steel bar layers :');
    Ps_0 = 0;
    Ms_0 = 0;
    e_0 = e_c;
    fy = input('please enter the grade of steel :');
    if (fy == 415)
        for j = 0:s-1
            e_j = e_0;
            epsilon_sj = (epsilon_max/xu)*(xu-e_j);
            if(epsilon_sj<0.00144)
                sigma_sj = (288.7/0.00144)*epsilon_sj;
            elseif(0.00144<=epsilon_sj)&&(epsilon_sj<=0.00163)
                sigma_sj = 288.7 + (((306.7-288.7)/(0.00163-0.00144))*(epsilon_sj-0.00144));
            elseif(0.00163<=epsilon_sj)&&(epsilon_sj<=0.00192)
                sigma_sj = 306.7 + (((324.8-306.7)/(0.00192-0.00163))*(epsilon_sj-0.00163));
            elseif(0.00192<=epsilon_sj)&&(epsilon_sj<=0.00241)
                sigma_sj = 324.8 + (((342.8-324.8)/(0.00241-0.00192))*(epsilon_sj-0.00192));
            elseif(0.00241<=epsilon_sj)&&(epsilon_sj<=0.00276)
                sigma_sj = 342.8 + (((351.8-342.8)/(0.00276-0.00241))*(epsilon_sj-0.00241));
            elseif(0.00276<=epsilon_sj)&&(epsilon_sj<=0.00380)
                sigma_sj = 351.8 + (((360.9-351.8)/(0.00380-0.00276))*(epsilon_sj-0.00276));
            elseif(epsilon_sj>0.00380)
                sigma_sj = 360.9;
            end
            Ps_j = sigma_sj*(N/2)*(3.14/4)*d^2*0.001;
            Ms_j = Ps_j*((D/2)-(e_c))*0.001*(-1)^(j);
            e_0 = e_j + ((D-(2*e_c))/(s-1));
            Ps_total = Ps_0 + Ps_j;
            Ps_0 = Ps_total;
            Ms_total = Ms_0 + Ms_j;
            Ms_0 = Ms_total;
        end
    elseif (fy == 500)
        for j = 0:s-1
            e_j = e_0;
            epsilon_sj = (epsilon_max/xu)*(xu-e_j);
            if(epsilon_sj<0.00174)
                sigma_sj = (347.8/0.00174)*epsilon_sj;
            elseif(0.00174<=epsilon_sj)&&(epsilon_sj<=0.00195)
                sigma_sj = 347.8 + (((369.6-347.8)/(0.00195-0.00174))*(epsilon_sj-0.00174));
            elseif(0.00195<=epsilon_sj)&&(epsilon_sj<=0.00226)
                sigma_sj = 369.6 + (((391.3-369.6)/(0.00226-0.00195))*(epsilon_sj-0.00195));
            elseif(0.00226<=epsilon_sj)&&(epsilon_sj<=0.00277)
                sigma_sj = 391.3 + (((413.0-391.3)/(0.00277-0.00226))*(epsilon_sj-0.00226));
            elseif(0.00277<=epsilon_sj)&&(epsilon_sj<=0.00312)
                sigma_sj = 413.0 + (((423.9-413.0)/(0.00312-0.00277))*(epsilon_sj-0.00277));
            elseif(0.00312<=epsilon_sj)&&(epsilon_sj<=0.00417)
                sigma_sj = 423.9 + (((434.8-423.9)/(0.00417-0.00312))*(epsilon_sj-0.00312));
            elseif(epsilon_sj>0.00417)
                sigma_sj = 434.8;
            end
            Ps_j = sigma_sj*(N/2)*(3.14/4)*d^2*0.001;
            Ms_j = Ps_j*((D/2)-(e_c))*0.001*(-1)^(j);
            e_0 = e_j + ((D-(2*e_c))/(s-1));
            Ps_total = Ps_0 + Ps_j;
            Ps_0 = Ps_total;
            Ms_total = Ms_0 + Ms_j;
            Ms_0 = Ms_total;
        end
    end
    disp(Ps_total)
    disp(Ms_total)
    
    %% Check
    P_total = Pc_total + Ps_total;
    if((P-2)<=P_total)&&(P_total<= (P+2))
        disp('ok')
    elseif(P_total>P)
        disp('Redo with less epsilon_max value')
    elseif(P_total<P)
        disp('Redo with more epsilon_max value')
    end
    M_total(1,k) = Mc_total + Ms_total;
    disp(M_total)
    disp(P_total)
    
    %% Plot
    hold on
    plot(Phi,M_total(1,k))
end

Answer 1

[Steven Lord]

What value should sigma_sj have if none of the conditions in your large if / elseif / elseif / ... statement are satisfied by the computed value of epsilon_sj? You might want to either define sigma_sj before that statement or include an else clause (not elseif) at the end as a fallback value.

[Rajkumar]

%% M-Phi Curve

clc;

clear all;

close all;

% Input data

fck = 20; % Grade of concrete in N/mm^2

Fe = 415; % Grade of steel in N/mm^2

P = 300; % Axial compressive load in kN

% Phi = 0.00004; % Initial slope in per meter

N = 6; % Number of bars

D = 500; % Depth of the section in mm

B = 250; % Width of the section in mm

e_c = 50; % Effective cover in mm

d = 16; % Diameter of the bars in mm

Phi = (0.00001:0.00001:0.00004);

k = length(Phi);

M_total = zeros(1,k);

% Phi = zeros(1,k);

%% Cross sectional area of the section

A = B*D; % Rectangular section Area in mm^2

%% Calculate strain under only axial compressive force

sigma_x = (P*1000)/(A); % Stress under axial compressive force in N/mm^2

% Quadratic equation solving

c = sigma_x/(0.45*fck);

a = 1;

b= -2;

delta = (b.^2)-(4*a*c);

if(delta < 0)

disp("Delta < 0 The equation does not have a real root");

elseif (delta == 0)

disp('The equation has only one real roots');

disp(-b./(2*a));

else

disp('The equation has two real roots');

epsilon_1 = 0.002*((-b+sqrt(delta))./(2*a));

epsilon_2 = 0.002*((-b-sqrt(delta))./(2*a));

end

epsilon = min(epsilon_1,epsilon_2);

disp(epsilon)

%% Finding depth of neutral axis

epsilon_max = input('Enter maximum stain value :');

for v = 1:k

Curvature = Phi(1,v);

epsilon_min = epsilon_max-(D*Curvature*0.001);

xu = (epsilon_max*D)/(epsilon_max-epsilon_min); % Depth of neutral axis in mm

if(xu > D)

disp("Neutral axis outside the section");

elseif(xu==D)

disp("Neutral axis is at the edge of the section");

elseif(0<=xu)&&(xu<=D)

disp("Neutral axis inside the section");

else

disp("Neutral axis is above the section");

end

%% Concrete stress, strain and force calculation

n = input('please enter number of concrete strips :');

Pc_0 = 0;

Mc_0 = 0;

for i = 0:n-1

for j = i

p_j = (D/n)*j + (D/(2*n));

epsilon_ci = (epsilon_max/xu)*(xu-p_j);

sigma_ci = 0.45*fck*(2*(epsilon_ci/0.002)-(epsilon_ci/0.002)^2);

P_ci = sigma_ci*B*(D/n)*0.001;

if (1 <= (i+1))&&((i+1) <= n/2)

M_ci = P_ci*(((D/2)-((D/(2*n))*(i+1)))/(i+1))*0.001;

elseif ((i+1) >= ((n/2)+1))

M_ci = P_ci*((((D/2)*(i-2))-(D/(2*n))))*0.001*(-1)^i;

end

end

Pc_total = Pc_0 + P_ci;

Pc_0 = Pc_total;

Mc_total = Mc_0 + M_ci;

Mc_0 = Mc_total;

end

disp(Pc_total)

disp(Mc_total)

%% Steel stress, strain and force calculation

s = input('please enter number of steel bar layers :');

Ps_0 = 0;

Ms_0 = 0;

e_0 = e_c;

sigma_sj = 0;

epsilon_sj = 0;

fy = input('please enter the grade of steel :');

if (fy == 415)

for j = 0:s-1

e_j = e_0;

epsilon_sj = (epsilon_max/xu)*(xu-e_j);

if(epsilon_sj<0.00144)

sigma_sj = (288.7/0.00144)*epsilon_sj;

elseif(0.00144<=epsilon_sj)&&(epsilon_sj<=0.00163)

sigma_sj = 288.7 + (((306.7-288.7)/(0.00163-0.00144))*(epsilon_sj-0.00144));

elseif(0.00163<=epsilon_sj)&&(epsilon_sj<=0.00192)

sigma_sj = 306.7 + (((324.8-306.7)/(0.00192-0.00163))*(epsilon_sj-0.00163));

elseif(0.00192<=epsilon_sj)&&(epsilon_sj<=0.00241)

sigma_sj = 324.8 + (((342.8-324.8)/(0.00241-0.00192))*(epsilon_sj-0.00192));

elseif(0.00241<=epsilon_sj)&&(epsilon_sj<=0.00276)

sigma_sj = 342.8 + (((351.8-342.8)/(0.00276-0.00241))*(epsilon_sj-0.00241));

elseif(0.00276<=epsilon_sj)&&(epsilon_sj<=0.00380)

sigma_sj = 351.8 + (((360.9-351.8)/(0.00380-0.00276))*(epsilon_sj-0.00276));

elseif(epsilon_sj>0.00380)

sigma_sj = 360.9;

end

Ps_j = sigma_sj*(N/2)*(3.14/4)*d^2*0.001;

Ms_j = Ps_j*((D/2)-(e_c))*0.001*(-1)^(j);

e_0 = e_j + ((D-(2*e_c))/(s-1));

Ps_total = Ps_0 + Ps_j;

Ps_0 = Ps_total;

Ms_total = Ms_0 + Ms_j;

Ms_0 = Ms_total;

end

elseif (fy == 500)

for j = 0:s-1

e_j = e_0;

epsilon_sj = (epsilon_max/xu)*(xu-e_j);

if(epsilon_sj<0.00174)

sigma_sj = (347.8/0.00174)*epsilon_sj;

elseif(0.00174<=epsilon_sj)&&(epsilon_sj<=0.00195)

sigma_sj = 347.8 + (((369.6-347.8)/(0.00195-0.00174))*(epsilon_sj-0.00174));

elseif(0.00195<=epsilon_sj)&&(epsilon_sj<=0.00226)

sigma_sj = 369.6 + (((391.3-369.6)/(0.00226-0.00195))*(epsilon_sj-0.00195));

elseif(0.00226<=epsilon_sj)&&(epsilon_sj<=0.00277)

sigma_sj = 391.3 + (((413.0-391.3)/(0.00277-0.00226))*(epsilon_sj-0.00226));

elseif(0.00277<=epsilon_sj)&&(epsilon_sj<=0.00312)

sigma_sj = 413.0 + (((423.9-413.0)/(0.00312-0.00277))*(epsilon_sj-0.00277));

elseif(0.00312<=epsilon_sj)&&(epsilon_sj<=0.00417)

sigma_sj = 423.9 + (((434.8-423.9)/(0.00417-0.00312))*(epsilon_sj-0.00312));

elseif(epsilon_sj>0.00417)

sigma_sj = 434.8;

end

Ps_j = sigma_sj*(N/2)*(3.14/4)*d^2*0.001;

Ms_j = Ps_j*((D/2)-(e_c))*0.001*(-1)^(j);

e_0 = e_j + ((D-(2*e_c))/(s-1));

Ps_total = Ps_0 + Ps_j;

Ps_0 = Ps_total;

Ms_total = Ms_0 + Ms_j;

Ms_0 = Ms_total;

end

end

disp(Ps_total)

disp(Ms_total)

%% Check

P_total = Pc_total + Ps_total;

if((P-2)<=P_total)&&(P_total<= (P+2))

disp('ok')

elseif(P_total>P)

disp('Redo with less epsilon_max value')

elseif(P_total<P)

disp('Redo with more epsilon_max value')

end

% M_total = Mc_total + Ms_total;

M_total(1,k) = Mc_total + Ms_total;

disp(M_total)

disp(P_total)

%% Plot

hold on

plot(Phi,M_total(1,k))

end

[Rajkumar]

why the graph is not forming?

[Rajkumar]

If i run blank graph is showing.

[Walter Roberson]

k is a scalar. plot(Phi,M_total(1,k)) is asking to plot a scalar