Why the plot is not coming? i want to plot for every Phi value corresponding M_total value.
1 view (last 30 days)
Show older comments
%% M-Phi Curve
clc;
close all;
% Input data
fc = 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
%%
for Phi = 0 : 0.00004
%% Calculate strain under only axial compressive force
sigma_x = (P*1000)/(D*B); % Stress under axial compressive force in N/mm^2
% Quadratic equation solving
c = sigma_x/(0.45*fc);
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 :');
epsilon_min = epsilon_max-D*Phi*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 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*fc*(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 bar layers :');
Ps_0 = 0;
Ms_0 = 0;
e_0 = e_c;
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
elseif(0.00163<=epsilon_sj)&&(epsilon_sj<=0.00192)
sigma_sj = 306.7 + ((324.8-306.7)/(0.00192-0.00163))*epsilon_sj
elseif(0.00192<=epsilon_sj)&&(epsilon_sj<=0.00241)
sigma_sj = 324.8 + ((342.8-324.8)/(0.00241-0.00192))*epsilon_sj
elseif(0.00241<=epsilon_sj)&&(epsilon_sj<=0.00276)
sigma_sj = 342.8 + ((351.8-342.8)/(0.00276-0.00241))*epsilon_sj
elseif(0.00276<=epsilon_sj)&&(epsilon_sj<=0.00380)
sigma_sj = 351.8 + ((360.9-351.8)/(0.0038-0.00276))*epsilon_sj
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
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;
end
figure
hold on
plot(Phi,M_total)
Answers (0)
See Also
Categories
Find more on Stress and Strain in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!