我正在尝试开发这个功能。 [phi, M] = moment_curvature_relationship(l, b, h, c_bot, c_top, As, Asc, del_dot, del_phi, n_fibres, epsilon_s_max, N_axis, tol, dy0, rho, f_c, fy, Es, Esp);
function [sigma_st, st_ten_yield] = stress_strain_steel(epsilon_s, fy, Es, Esp, es_dot)
% epsilon_s - strain in steel
% fy - yield strength of steel
% Es - modulus of elasticity of steel
% Esp - post-yield modulus of steel
% Initialize stress array
sigma_st = zeros(size(epsilon_s));
% Calculate Dynamic Increase Factor of yield strength
DIF_st = (es_dot / (10^-4)) .^ (0.074 - 0.04 * fy / 414);
% Define yield strain
epsilon_y = DIF_st * fy / Es;
% Check if the steel strain exceeds yield strain
if abs(epsilon_s) >= abs(epsilon_y)
st_ten_yield = true; % Steel has yielded
else
st_ten_yield = false; % Steel has not yet yielded
end
% Calculate stress-strain behavior of steel bar in compression (-ve)
% and tension (+ve)
% Elastic compression
steel_comp_elastic = (epsilon_s <= 0) & (epsilon_s >= - epsilon_y);
sigma_st(steel_comp_elastic) = Es .* epsilon_s(steel_comp_elastic);
% Compression steel hardening
steel_comp_hard = (epsilon_s < - epsilon_y);
sigma_st(steel_comp_hard) = - (DIF_st .* fy + Esp .* (abs(epsilon_s(steel_comp_hard)) - (DIF_st .* fy / Es)));
% Elastic tension
steel_tens_elastic = (epsilon_s >= 0) & (epsilon_s <= epsilon_y);
sigma_st(steel_tens_elastic) = Es .* epsilon_s(steel_tens_elastic);
% Tension steel hardening
steel_tens_hard = (epsilon_s > epsilon_y);
sigma_st(steel_tens_hard) = DIF_st .* fy + Esp .* (abs(epsilon_s(steel_tens_hard)) - (DIF_st .* fy / Es));
end
function [phi, M] = moment_curvature_relationship(l, b, h, c_bot, c_top, As, Asc, del_dot, del_phi, n_fibers, epsilon_s_max, N_axial, dy0, tol, rho, f_c, fy, Es, Esp)
% Curvature rate conversion from mid-span loading, del_dot
phi_dot = del_dot * 12 / l^2;
% Preallocate arrays for output
phi = []; % Curvature
M = []; % Moment
% Initialize variables
c = 0; % Initial curvature condition
termination = false; % Analysis termination flag
% Create condition for loop termination
while ~termination
%Assume neutral axis depth
y0 = h/2;
% Assume large difference initially for the loop
diff = inf;
% Initilize Moment in fiber elements for neutral axis
bending_moment = 0;
% Limit max number of iterations
max_iter = 10; % Limit to number of iterations
iter_count = 0; % Initialize iteration count
% Convergence based on the equilibrium of normal forces
% Iterating neutral axis position
while abs(diff) > tol && iter_count < max_iter
% Initialize concrete fiber element variables
y = zeros(n_fibers, 1); % Array for position of fiber
epsilon = zeros(n_fibers, 1); % Array for strain in fiber element
e_dot = zeros(n_fibers, 1); % Array for strain rate in fiber element
sigma = zeros(n_fibers, 1); % Array for stress in fiber element
N_c = zeros(1, n_fibers); % Array to store force in each fiber element
M_c = zeros(1, n_fibers); % Array tp store moment in each fiber element
% Iterate over each fiber element for stress and force calculation
for i = 1:n_fibers
y(i) = (i - 0.5) * h / n_fibers; % fiber position
% Strain in concrete fiber element
epsilon(i) = c .* (y0 - y(i));
% Strain rate in concrete fiber element
e_dot(i) = abs(y0 - y(i)) * phi_dot;
% Stress in fiber
sigma(i) = stress_strain_concrete(rho, f_c, epsilon(i), e_dot(i));
% Force contribution from concrete fiber
A_c = b * h / n_fibers; %Area of each fiber
N_c(i) = sigma(i) * A_c;
% Moment contribution from concrete fiber
M_c(i) = N_c(i) * (y0 - y(i));
end
% Calculate strain of steel fiber
ys = c_bot; % Positions of steel rebars (tension, compression)
% If Asc (compression steel area) is greater than 0, include compression steel
if Asc > 0
ys = [c_bot, c_top]; % Positions of tensile and compression steel
end
% Initialize array for storing variables
epsilon_s = zeros(length(ys), 1);
es_dot = zeros(length(ys), 1);
sigma_st = zeros(1, length(ys));
N_s = zeros(1, length(ys)); % Total comtribution of steel fiber to axial force
M_s = zeros(1, length(ys)); % Total moment contribution of steel fiber
st_ten_yield_total = false; % Flag to track tensile steel yielding
for i = 1:length(ys)
epsilon_s(i) = c * (y0 - ys(i)); % Strain in steel rebars
es_dot(i) = abs(ys(i) - y0) * phi_dot; % Strain rate in steel rebar
[sigma_st(i), st_ten_yield(i)] = stress_strain_steel(epsilon_s(i), fy, Es, Esp, es_dot(i));
% If this is the tensile steel (i == 1) and it yields, print phi
if i == 1 && st_ten_yield(i) && ~st_ten_yield_total
fprintf('Curvature phi at which epsilon_s(1) reaches yield strain: %.6f\n', c);
st_ten_yield_total = true; % Ensure this is only printed once
end
% Loop for calculating Area of steel fiber
if i == 1
A_s = As; % Area of tensile steel fiber
elseif i == 2
A_s = Asc; % Area of compression steel fiber
end
N_s(i) = sigma_st(i) .* A_s; % Force contribution from steel fiber
M_s(i) = N_s(i) .* (y0 - ys(i)); % Moment contribution from steel fiber
end
% Total force contibution
N_total = sum(N_c) + sum(N_s); % Total axial force
% Calculate force residual
diff = N_total - N_axial;
% Adjust neutral axis position (y0)
% Numerical derivative for Newton-Raphson
y0_pert = y0 + dy0; % Perturbed neutral axis
% Recompute N_total for perturbed y0
% (Repeating calculation for N_total using y0_pert in place of y0)
% Reinitialize accumulated variables
% Re-Initialize concrete fiber element variables for perturbed
% state
y_pert = zeros(n_fibers, 1); % Array for position of fiber
epsilon_pert = zeros(n_fibers, 1); % Array for strain in fiber element
sigma_pert = zeros(n_fibers, 1); % Array for stress in fiber element
e_dot_pert = zeros(n_fibers, 1); % Array to store strain rate in fiber element
N_c_pert = zeros(1, n_fibers); % Array to store force in each fiber element
% Iterate over each fiber element for stress and force calculation
for i = 1:n_fibers
y_pert(i) = (i - 0.5) * h / n_fibers; % fiber position
% Strain in concrete fiber element
epsilon_pert(i) = c .* (y0_pert - y(i));
% Strain rate in concrete fiber element
e_dot_pert(i) = abs(y0_pert - y(i)) * phi_dot;
% Stress in fiber
sigma_pert(i) = stress_strain_concrete(rho, f_c, epsilon_pert(i), e_dot_pert(i));
% Force contribution from concrete fiber
A_c = b * h / n_fibers; %Area of each fiber
N_c_pert(i) = sigma_pert(i) * A_c;
end
% Re-Initialize array for storing variables for perturbed state
epsilon_s_pert = zeros(length(ys), 1);
es_dot_pert = zeros(length(ys), 1);
sigma_st_pert = zeros(length(ys), 1);
N_s_pert = zeros(1, length(ys)); % Reset total steel axial force for perturbed calculation
for i = 1:length(ys)
epsilon_s_pert(i) = c * (y0_pert - ys(i)); % Strain in steel rebars
es_dot_pert(i) = abs(ys(i) - y0_pert) * phi_dot; % Strain rate in steel rebar
sigma_st_pert(i) = stress_strain_steel(epsilon_s_pert(i), fy, Es, Esp, es_dot_pert(i));
% Loop for calculating Area of steel fiber
if i == 1
A_s = As; % Area of tensile steel fiber
elseif i == 2
A_s = Asc; % Area of compression steel fiber
end
N_s_pert(i) = sigma_st_pert(i) .* A_s; % Force contribution from steel fiber
end
% Total force contibution
N_total_pert = sum(N_c_pert) + sum(N_s_pert); % Total axial force
% Compute derivative of N with respect to y0
dN_dy0 = (N_total_pert - N_total) / dy0;
% Newton-Raphson update of y0
y0 = y0 - (N_total - N_axial) / dN_dy0;
% Exit if y0 goes beyond bounds
if y0 <= c_bot || y0 >= h
termination = true;
fprintf('The value of y0: %.2f mm goes beyond bounds. \n', y0);
break;
% Exit if maximum iterations reached without convergence
elseif iter_count == max_iter
fprintf('Newton-Raphson did not converge after %d iterations. \n, ', max_iter);
break;
end
% Increment iteration count
iter_count = iter_count + 1;
end
% Calculate bending moment contribution from concrete and steel
% and any Axial force component (considered here at section mid
bending_moment = bending_moment + sum(M_c) + sum(M_s) + N_axial * h / 2;
bending_moment = bending_moment / 1000; % Converted to kNmm
% Store the outputs for curvature and moment
phi(end + 1) = c;
M(end + 1) = bending_moment;
% Check termination condition for moment-curvature (Here: Fracture strain of tensile steel)
if max(abs(epsilon_s)) > epsilon_s_max
termination = true;
end
% Increment in curvature
c = c + del_phi;
end
end
功能运行正常，只是这个问题我无法解决。
问题在于，在此函数中，st_ten_yield 的打印行被打印多次。有什么解决办法吗？
出现此问题的原因是，在循环内检查并更新了
st_ten_yield_totalst_ten_yield_total