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?

1 次查看(过去 30 天)
[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

回答(1 个)

Steven Lord
Steven Lord 2023-3-11
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.
  4 个评论

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Stress and Strain 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by