why in Free vibration 2DOF system with damping , the displacement is increasing with time?

% This is the code I have written


function solve_open_circuit_free_vibration
% System parameters
params.m = 0.91;               % Inner mass (kg)
params.M = 1.2;                % Outer mass (kg)
params.K = 114.259;            % Spring stiffness (N/m)
params.c_m = 0.4;             % Mechanical damping coefficient (N·s/m)
params.mu = 3.55e-7;           % Magnetic dipole moment ( in T·m^3)
params.a = 0.0171;             % Solenoid Coil radius (m)
params.b0 = 0.001;             % Initial distance b/w mass and coil (m)
% Initial conditions [x, dx, y, dy]
Z0 = [0, 0, 0.01, 0];          % Initial displacements and velocities, initial displacement is provided on outer mass
% Time span for simulation
tspan = [0, 5];                % Simulation time range (5 seconds)
% Solver options
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
% Solve using ode15s
[T, Z] = ode15s(@(t, Z) governing_equations(t, Z, params), tspan, Z0, options);
% Extract results
x = Z(:, 1);                   % Inner mass displacement
dx = Z(:, 2);                  % Inner mass velocity
y = Z(:, 3);                   % Outer mass displacement
dy = Z(:, 4);                  % Outer mass velocity
% Compute voltage
z = params.b0 + y - x;         % Magnet-to-coil relative displacement
dz = dy - dx;                  % Relative velocity
V = compute_voltage(z, dz, params);  % Voltage
% Plot results
plot_results(T, x, dx, y, dy, z, V);
end
function dZdt = governing_equations(t, Z, params)
% Extract state variables
x = Z(1); dx = Z(2);
y = Z(3); dy = Z(4);
% Relative displacement and velocity
z = params.b0 + y - x;
dz = dy - dx;
% Equations of motion
ddx = (params.K * (y - x) - params.c_m * dz) / params.m;
ddy = (-params.K * (y - x) + params.c_m * dz) / params.M;
% Return derivatives
dZdt = [dx; ddx; dy; ddy];
end
function V = compute_voltage(z, dz, params)
% Compute voltage
V = (6 * pi * params.mu * params.a^2 * z .* dz) ./ (params.a^2 + z.^2).^(5/2);
end
function plot_results(T, x, dx, y, dy, z, V)
% Create a single figure with subplots
figure;
% Subplot 1: Displacement
subplot(1, 2, 1);
plot(T, x, 'r', 'LineWidth', 1.5); hold on;
plot(T, y, 'b--', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Inner Mass (x)', 'Outer Mass (y)');
title('Displacement vs Time');
grid on;
% Subplot 2: Velocity
subplot(1, 2, 2);
plot(T, dx, 'r', 'LineWidth', 1.5); hold on;
plot(T, dy, 'b--', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Velocity (m/s)');
legend('Inner Mass (dx)', 'Outer Mass (dy)');
title('Velocity vs Time');
grid on;
end

Answer 1

[Epsilon]

Hi Parthajit,

The incorrect displacement and velocity calculations seem to be due to the wrong direction of forces, defined in the differential equations.

In this system, the forces acting on the inner mass include the spring and the damping force along a single direction. Similarly, the outer mass experiences the same set of forces but in the opposite direction.

Hence, the correct equations would be:

Rectified code:

function solve_open_circuit_free_vibration
params.m = 0.91;              
params.M = 1.2;               
params.K = 114.259;           
params.c_m = 0.4;            
params.mu = 3.55e-7;          
params.a = 0.0171;            
params.b0 = 0.001;            
Z0 = [0, 0, 0.01, 0];         
tspan = [0, 5];                
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[T, Z] = ode15s(@(t, Z) governing_equations(t, Z, params), tspan, Z0, options);
x = Z(:, 1);                   
dx = Z(:, 2);                  
y = Z(:, 3);                   
dy = Z(:, 4);                  
z = params.b0 + y - x;         
dz = dy - dx;                  
V = compute_voltage(z, dz, params);  
plot_results(T, x, dx, y, dy, z, V);
end

function dZdt = governing_equations(~, Z, params)
x = Z(1); dx = Z(2);
y = Z(3); dy = Z(4);
dz = dy - dx;
%% Equations of motion rectified %%
ddx = (params.K * (y - x) + params.c_m * dz) / params.m;
ddy = (-params.K * (y - x) - params.c_m * dz) / params.M;
dZdt = [dx; ddx; dy; ddy];
end

function V = compute_voltage(z, dz, params)
V = (6 * pi * params.mu * params.a^2 * z .* dz) ./ (params.a^2 + z.^2).^(5/2);
end

function plot_results(T, x, dx, y, dy, ~, ~)
figure;
subplot(1, 2, 1);
plot(T, x, 'r', 'LineWidth', 1.5); hold on;
plot(T, y, 'b--', 'LineWidth', 1.5)
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Inner Mass (x)', 'Outer Mass (y)');
title('Displacement vs Time');
grid on;
subplot(1, 2, 2);
plot(T, dx, 'r', 'LineWidth', 1.5); hold on;
plot(T, dy, 'b--', 'LineWidth', 1.5)
xlabel('Time (s)');
ylabel('Velocity (m/s)');
legend('Inner Mass (dx)', 'Outer Mass (dy)');
title('Velocity vs Time');
grid on;
end
solve_open_circuit_free_vibration;

Documentation link for ode15s: https://www.mathworks.com/help/matlab/ref/ode15s.html