Mass spring damper system problem

7 次查看（过去 30 天）
OnurErdogan Erdogan
评论： Sam Chak 2024-6-21
Ran in:
Hello, I am designing the mass spring damper system in the image. I am creating this code for the design of a ship crane. In my road equation, there is an increase of approximately 20 meters, and the spring k_c is stretched by 20 meters. I just want the spring to follow the road, but the spring follows the road like an arrow, which is very high like 2000N/mm². Stress occurs. I think there is a mistake in my input equation, but I couldn't figure it out completely. I would be very happy if you could help me
The image of the system is also attached.
Matlab code:
% System parameters
M_c = 917; % Mass of the crane boom (kg)
M_s = 87; % Mass of the shock absorber (kg)
M_WL = 4075; % Mass of the working load (kg)
k_c = 554376; % Spring constant of the crane + crane hook ropes (N/m)
k_s = 24162; % Spring constant of the shock absorber system (N/m)
k_R = 9090572; % Spring constant of the lifting beam and lower rope (N/m)
c_s = 5800; % Damping constant of the shock absorber system (Ns/m)
A = 1.5; % Amplitude of input signal (m)
omega = 1; % Angular frequency of input signal (rad/s)
V_H1 = 0.69; % Linear term coefficient for first interval (m/s)
V_H2 = 0; % Linear term coefficient for second interval (m/s)
V_H3 = -0.1; % Linear term coefficient for third interval (m/s)
d_c = 0.026; % Diameter of spring kc (m)
d_s = 0.03; % Diameter of spring ks (m)
d_R = 0.0226; % Diameter of spring kR (m)
% Combined time vector with smaller sampling period
t = 0:0.01:45; % From 0 to 45 seconds with 0.01s sampling period
% Input signal (m)
u = zeros(size(t));
u(t <= 25) = A * sin(omega * t(t <= 25)) + V_H1 * t(t <= 25);
u(t > 25 & t <= 35) = A * sin(omega * t(t > 25 & t <= 35)) + V_H1 * 25 + V_H2 * (t(t > 25 & t <= 35) - 25);
u(t > 35) = A * sin(omega * t(t > 35)) + V_H1 * 25 + V_H3 * (t(t > 35) - 35);
% State-space matrix
A_matrix = [0 1 0 0 0 0;
-k_c/M_c -c_s/M_c k_s/M_c c_s/M_c 0 0;
0 0 0 1 0 0;
k_s/M_s c_s/M_s -(k_s+k_R)/M_s -c_s/M_s k_R/M_s 0;
0 0 0 0 0 1;
0 0 k_R/M_WL 0 -k_R/M_WL 0];
B_matrix = [0;
k_c/M_c;
0;
0;
0;
0];
C_matrix = [1 0 0 0 0 0;
0 0 1 0 0 0;
0 0 0 0 1 0];
D_matrix = [0;
0;
0];
% Create state-space model
sys = ss(A_matrix, B_matrix, C_matrix, D_matrix);
% Simulate system response
[y, t, x] = lsim(sys, u, t);
Warning: The input signal is undersampled. Use a sampling period smaller than 3.9e-03.
% Calculate spring extensions (m)
delta_kc = u'- y(:,1); % Extension of spring kc
delta_ks = y(:,1) - y(:,2); % Extension of spring ks
delta_kR = y(:,2) - y(:,3); % Extension of spring kR
% Convert spring extensions to meters (m)
delta_kc_m = delta_kc; % Assuming input signal u is in meters
delta_ks_m = delta_ks; % Assuming displacements are in meters
delta_kR_m = delta_kR; % Assuming displacements are in meters
% Calculate spring stresses (N/m²)
A_c = pi * (d_c/2)^2; % Cross-sectional area of spring kc (m²)
A_s = pi * (d_s/2)^2; % Cross-sectional area of spring ks (m²)
A_R = pi * (d_R/2)^2; % Cross-sectional area of spring kR (m²)
F_kc = k_c * delta_kc_m; % Force in spring kc (N)
F_ks = k_s * delta_ks_m; % Force in spring ks (N)
F_kR = k_R * delta_kR_m; % Force in spring kR (N)
sigma_kc = F_kc / A_c; % Stress in spring kc (N/m²)
sigma_ks = F_ks / A_s; % Stress in spring ks (N/m²)
sigma_kR = F_kR / A_R; % Stress in spring kR (N/m²)
% Convert stress to N/mm²
sigma_kc_mm2 = sigma_kc * 1e-6;
sigma_ks_mm2 = sigma_ks * 1e-6;
sigma_kR_mm2 = sigma_kR * 1e-6;
% Plot input signal and outputs
figure;
subplot(5,1,1);
plot(t, u);
title('Input Signal u');
xlabel('Time (s)');
ylabel('u (m)');
subplot(5,1,2);
plot(t, y(:,1));
title('Output y1');
xlabel('Time (s)');
ylabel('y1 (m)');
subplot(5,1,3);
plot(t, y(:,2));
title('Output y2');
xlabel('Time (s)');
ylabel('y2 (m)');
subplot(5,1,4);
plot(t, y(:,3));
title('Output y3');
xlabel('Time (s)');
ylabel('y3 (m)');
% Plot spring stresses in separate figures
subplot(5,1,5);
hold on;
plot(t, sigma_kc_mm2, 'r');
plot(t, sigma_ks_mm2, 'g');
plot(t, sigma_kR_mm2, 'b');
title('Stress in Springs');
xlabel('Time (s)');
ylabel('Stress (N/mm²)');
legend('Stress \sigma_{kc}', 'Stress \sigma_{ks}', 'Stress \sigma_{kR}');
grid on;
% Animation
figure;
% Define positions of masses and springs
y1_base = 2; % Base position of mass 1 (m)
y2_base = 1; % Base position of mass 2 (m)
y3_base = 0; % Base position of mass 3 (m)
% Plot settings
axis([-0.5 1.5 -1 30]);
hold on;
% Initial positions of masses
h_m1 = rectangle('Position',[0.4 y1_base + y(1,1) 0.2 0.2], 'Curvature', [0.1 0.1], 'FaceColor', 'blue');
h_m2 = rectangle('Position',[0.4 y2_base + y(1,2) 0.2 0.2], 'Curvature', [0.1 0.1], 'FaceColor', 'green');
h_m3 = rectangle('Position',[0.4 y3_base + y(1,3) 0.2 0.2], 'Curvature', [0.1 0.1], 'FaceColor', 'red');
% Labels for masses
text(0.65, y1_base + y(1,1) + 0.1, 'M_c', 'FontSize', 12, 'Color', 'blue');
text(0.65, y2_base + y(1,2) + 0.1, 'M_s', 'FontSize', 12, 'Color', 'green');
text(0.65, y3_base + y(1,3) + 0.1, 'M_{WL}', 'FontSize', 12, 'Color', 'red');
% Initial positions of springs
h_spring1 = plot([0.5 0.5], [3 y1_base + y(1,1) + 0.2], 'k', 'LineWidth', 2);
h_spring2 = plot([0.5 0.5], [y1_base + y(1,1) y2_base + y(1,2) + 0.2], 'k', 'LineWidth', 2);
h_spring3 = plot([0.5 0.5], [y2_base + y(1,2) y3_base + y(1,3) + 0.2], 'k', 'LineWidth', 2);
% Labels for springs and dampers
text(0.3, 2.5, 'k_c', 'FontSize', 12, 'Color', 'black');
text(0.3, 1.5, 'k_s', 'FontSize', 12, 'Color', 'black');
text(0.3, 0.5, 'k_R', 'FontSize', 12, 'Color', 'black');
text(0.7, y1_base + (y2_base - y1_base) / 2 + y(1,1) / 2 + y(1,2) / 2, 'c_s', 'FontSize', 12, 'Color', 'black');
% Initial position of the damper
h_damper = plot([0.45 0.55], [y1_base + (y2_base - y1_base) / 2 + y(1,1) / 2 + y(1,2) / 2 y1_base + (y2_base - y1_base) / 2 + y(1,1) / 2 + y(1,2) / 2], 'k', 'LineWidth', 6);
% Animation loop
for i = 1:length(t)
% Update positions of masses
if isvalid(h_m1)
set(h_m1, 'Position', [0.4 y1_base + y(i,1) 0.2 0.2]);
end
if isvalid(h_m2)
set(h_m2, 'Position', [0.4 y2_base + y(i,2) 0.2 0.2]);
end
if isvalid(h_m3)
set(h_m3, 'Position', [0.4 y3_base + y(i,3) 0.2 0.2]);
end
% Update positions of springs
if isvalid(h_spring1)
set(h_spring1, 'YData', [3 y1_base + y(i,1) + 0.2]);
end
if isvalid(h_spring2)
set(h_spring2, 'YData', [y1_base + y(i,1) y2_base + y(i,2) + 0.2]);
end
if isvalid(h_spring3)
set(h_spring3, 'YData', [y2_base + y(i,2) y3_base + y(i,3) + 0.2]);
end
% Update position of the damper
if isvalid(h_damper)
set(h_damper, 'YData', [y1_base + (y2_base - y1_base) / 2 + y(i,1) / 2 + y(i,2) / 2 y1_base + (y2_base - y1_base) / 2 + y(i,1) / 2 + y(i,2) / 2]);
end
% Pause to create animation effect
pause(0.01);
end
hold off;
  3 个评论
显示 1更早的评论
OnurErdogan Erdogan
OnurErdogan Erdogan 2024-6-21
移动：Sam Chak 2024-6-21
Hello, thank you very much for your quick and valuable feedback.
I made the improvements you mentioned.
Could you have an opinion on the subject below?
The lower end of the k_c spring (the position where the u input enters) remains stable. And then, as the spring k_c stretches up to 20 meters, very high tensions occur. I expected that this end of the arc would follow the road data.
You can see lates version of my data
% System parameters
M_c = 917; % Mass of the crane boom (kg)
M_s = 87; % Mass of the shock absorber (kg)
M_WL = 4075; % Mass of the working load (kg)
k_c = 554376; % Spring constant of the crane + crane hook ropes (N/m)
k_s = 24162; % Spring constant of the shock absorber system (N/m)
k_R = 9090572; % Spring constant of the lifting beam and lower rope (N/m)
c_s = 5800; % Damping constant of the shock absorber system (Ns/m)
A = 1.5; % Amplitude of input signal (m)
omega = 1; % Angular frequency of input signal (rad/s)
V_H1 = 0.69; % Linear term coefficient for first interval (m/s)
V_H2 = 0; % Linear term coefficient for second interval (m/s)
V_H3 = -0.1; % Linear term coefficient for third interval (m/s)
d_c = 0.026; % Diameter of spring kc (m)
d_s = 0.03; % Diameter of spring ks (m)
d_R = 0.0226; % Diameter of spring kR (m)
% Combined time vector with smaller sampling period
t = 0:0.003:63; % From 0 to 63 seconds with 0.003s sampling period
% Input signal (m)
u = zeros(size(t));
u(t <= 25) = A * sin(omega * t(t <= 25)) + V_H1 * t(t <= 25);
u(t > 25 & t <= 35) = A * sin(omega * t(t > 25 & t <= 35)) + V_H1 * 25 + V_H2 * (t(t > 25 & t <= 35) - 25);
u(t > 35) = A * sin(omega * t(t > 35)) + V_H1 * 25 + V_H3 * (t(t > 35) - 35);
% Plot input signal
figure;
plot(t, u);
title('Input Signal, u(t)');
xlabel('Time (s)');
ylabel('u (m)');
grid on;
xticks(0:7:63);
xlim([0 63]);
% State-space matrix
A_matrix = [0 1 0 0 0 0;
-k_c/M_c -c_s/M_c k_s/M_c c_s/M_c 0 0;
0 0 0 1 0 0;
k_s/M_s c_s/M_s -(k_s+k_R)/M_s -c_s/M_s k_R/M_s 0;
0 0 0 0 0 1;
0 0 k_R/M_WL 0 -k_R/M_WL 0];
B_matrix = [0;
k_c/M_c;
0;
0;
0;
0];
C_matrix = [1 0 0 0 0 0;
0 0 1 0 0 0;
0 0 0 0 1 0];
D_matrix = [0;
0;
0];
% Create state-space model
sys = ss(A_matrix, B_matrix, C_matrix, D_matrix);
% Simulate system response
[y, t, x] = lsim(sys, u, t);
% Calculate spring extensions (m)
delta_kc = u' - y(:,1); % Extension of spring kc
delta_ks = y(:,1) - y(:,2); % Extension of spring ks
delta_kR = y(:,2) - y(:,3); % Extension of spring kR
% Convert spring extensions to meters (m)
delta_kc_m = delta_kc; % Assuming input signal u is in meters
delta_ks_m = delta_ks; % Assuming displacements are in meters
delta_kR_m = delta_kR; % Assuming displacements are in meters
% Calculate spring stresses (N/m²)
A_c = pi * (d_c/2)^2; % Cross-sectional area of spring kc (m²)
A_s = pi * (d_s/2)^2; % Cross-sectional area of spring ks (m²)
A_R = pi * (d_R/2)^2; % Cross-sectional area of spring kR (m²)
F_kc = k_c * delta_kc_m; % Force in spring kc (N)
F_ks = k_s * delta_ks_m; % Force in spring ks (N)
F_kR = k_R * delta_kR_m; % Force in spring kR (N)
sigma_kc = F_kc / A_c; % Stress in spring kc (N/m²)
sigma_ks = F_ks / A_s; % Stress in spring ks (N/m²)
sigma_kR = F_kR / A_R; % Stress in spring kR (N/m²)
% Convert stress to N/mm²
sigma_kc_mm2 = sigma_kc * 1e-6;
sigma_ks_mm2 = sigma_ks * 1e-6;
sigma_kR_mm2 = sigma_kR * 1e-6;
% Plot input signal and outputs
figure;
subplot(5,1,1);
plot(t, u);
title('Input Signal u');
xlabel('Time (s)');
ylabel('u (m)');
subplot(5,1,2);
plot(t, y(:,1));
title('Output y1');
xlabel('Time (s)');
ylabel('y1 (m)');
subplot(5,1,3);
plot(t, y(:,2));
title('Output y2');
xlabel('Time (s)');
ylabel('y2 (m)');
subplot(5,1,4);
plot(t, y(:,3));
title('Output y3');
xlabel('Time (s)');
ylabel('y3 (m)');
% Plot spring stresses in separate figures
subplot(5,1,5);
hold on;
plot(t, sigma_kc_mm2, 'r');
plot(t, sigma_ks_mm2, 'g');
plot(t, sigma_kR_mm2, 'b');
title('Stress in Springs');
xlabel('Time (s)');
ylabel('Stress (N/mm²)');
legend('Stress \sigma_{kc}', 'Stress \sigma_{ks}', 'Stress \sigma_{kR}');
grid on;
% Check stability of the system
is_stable = isstable(sys);
disp(['System is stable: ', num2str(is_stable)]);
eigenvalues = eig(A_matrix);
disp('Eigenvalues of A_matrix:');
disp(eigenvalues);
% Animation
figure;
% Define positions of masses and springs
y1_base = 2; % Base position of mass 1 (m)
y2_base = 1; % Base position of mass 2 (m)
y3_base = 0; % Base position of mass 3 (m)
% Plot settings
axis([-0.5 1.5 -1 30]);
hold on;
% Initial positions of masses
h_m1 = rectangle('Position',[0.4 y1_base + y(1,1) 0.2 0.2], 'Curvature', [0.1 0.1], 'FaceColor', 'blue');
h_m2 = rectangle('Position',[0.4 y2_base + y(1,2) 0.2 0.2], 'Curvature', [0.1 0.1], 'FaceColor', 'green');
h_m3 = rectangle('Position',[0.4 y3_base + y(1,3) 0.2 0.2], 'Curvature', [0.1 0.1], 'FaceColor', 'red');
% Labels for masses
text(0.65, y1_base + y(1,1) + 0.1, 'M_c', 'FontSize', 12, 'Color', 'blue');
text(0.65, y2_base + y(1,2) + 0.1, 'M_s', 'FontSize', 12, 'Color', 'green');
text(0.65, y3_base + y(1,3) + 0.1, 'M_{WL}', 'FontSize', 12, 'Color', 'red');
% Initial positions of springs
h_spring1 = plot([0.5 0.5], [3 y1_base + y(1,1) + 0.2], 'r', 'LineWidth', 2);
h_spring2 = plot([0.5 0.5], [y1_base + y(1,1) y2_base + y(1,2) + 0.2], 'g', 'LineWidth', 2);
h_spring3 = plot([0.5 0.5], [y2_base + y(1,2) y3_base + y(1,3) + 0.2], 'b', 'LineWidth', 2);
% Labels for springs and dampers
text(0.3, 2.5, 'k_c', 'FontSize', 12, 'Color', 'red');
text(0.3, 1.5, 'k_s', 'FontSize', 12, 'Color', 'green');
text(0.3, 0.5, 'k_R', 'FontSize', 12, 'Color', 'blue');
text(0.7, y1_base + (y2_base - y1_base) / 2 + y(1,1) / 2 + y(1,2) / 2, 'c_s', 'FontSize', 12, 'Color', 'black');
% Initial position of the damper
h_damper = plot([0.45 0.55], [y1_base + (y2_base - y1_base) / 2 + y(1,1) / 2 + y(1,2) / 2 y1_base + (y2_base - y1_base) / 2 + y(1,1) / 2 + y(1,2) / 2], 'k', 'LineWidth', 6);
% Animation loop
for i = 1:length(t)
% Update positions of masses
if isvalid(h_m1)
set(h_m1, 'Position', [0.4 y1_base + y(i,1) 0.2 0.2]);
end
if isvalid(h_m2)
set(h_m2, 'Position', [0.4 y2_base + y(i,2) 0.2 0.2]);
end
if isvalid(h_m3)
set(h_m3, 'Position', [0.4 y3_base + y(i,3) 0.2 0.2]);
end
% Update positions of springs
if isvalid(h_spring1)
set(h_spring1, 'YData', [3 y1_base + y(i,1) + 0.2]);
end
if isvalid(h_spring2)
set(h_spring2, 'YData', [y1_base + y(i,1) y2_base + y(i,2) + 0.2]);
end
if isvalid(h_spring3)
set(h_spring3, 'YData', [y2_base + y(i,2) y3_base + y(i,3) + 0.2]);
end
% Update position of the damper
if isvalid(h_damper)
set(h_damper, 'YData', [y1_base + (y2_base - y1_base) / 2 + y(i,1) / 2 + y(i,2) / 2 y1_base + (y2_base - y1_base) / 2 + y(i,1) / 2 + y(i,2) / 2]);
end
% Pause to create animation effect
pause(0.01);
end
hold off;
Sam Chak
Sam Chak 2024-6-21
Hi @OnurErdogan Erdogan
I didn't suggest any improvements because I'm not an expert in springs. I simply reminded you to verify if the input signal was plotted correctly as intended. However, I don't understand why you strongly assert that the stress of the spring (kc) should follow the road profile. Perhaps what you plotted previously is correct.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Aerospace Applications 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by