Mass spring damper system problem

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;

[Sam Chak]

Hi @OnurErdogan Erdogan

First, you should verify if the Input signal produces the expected response, as shown below. Since omega is 1 rad/s, you should count 10 peaks and 10 troughs within 63 seconds. Also, please note that your input signal is unbounded from a mathematical perspective. Thus, the output responses will be unbounded too, despite the system itself having negative real eigenvalues. You can check "eig(A_matrix)" or "isstable(sys)".

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)
t       = 0:0.01:63;

%% 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(t, u), grid on, xlabel('t'), title('Input Signal, u(t)'), 
xticks(0:7:63);
xlim([0 63])

[OnurErdogan Erdogan]

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]

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.