Why the plot is not coming? i want to plot for every Phi value corresponding M_total value.
1 次查看(过去 30 天)
显示 更早的评论
%% 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)
3 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Stress and Strain 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!