102_Matlab analysis of the mass-spring-damper system

Matlab analysis of the mass-spring-damper system (time responses based on the analytical solutions)

(1) msd class

To handle multiple configurations of mass-spring-damper systems, it is convenient to define a mass-spring-damper class ‘msd‘ to store properties of the system. In application, the user creates multiple msd instances, each of which corresponds with one configuration of interest.

The code below is the definition of the msd class. The class consists of the three parameters of a mass-spring-damper, namely, the mass \(m\), damper coefficient \(c\) and spring coefficient \(k\). The code also contains methods to calculate other dependent properties of the system, such as the undamped natural frequency \(\omega_{n}\), damping ratio \(\zeta\), transfer function, the root of the characteristic equation \(s\) and the time responses \(x(t)\) due to the step and impulse input forces \(f(t)\).

In the functions ‘Step’ and ‘Impulse’ in the below msd class definition, the time responses are calculated based on the analytical solutions of the differential equation of the system, listed in the left column on the tables in Summary of the analytical solutions 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

classdef msd
    % Mass Spring Damper Data Set Class
    properties
        m   % Mass [kg]
        c   % Damper coefficient [N/(m/s)]
        k   % Spring coefficient [N/m]
    end
    properties (Dependent)
        wn  % Undamped natural frequency [rad/s]
        wd  % Damped natural frequency [rad/s]
        zeta % Damping ratio
        s   % Roots of the characteristic equation
        TF  % Transfer function
    end
    methods
        % constructor
        function obj = msd(m, c ,k) 
            if nargin > 0
                obj.m = m;
                obj.c = c;
                obj.k = k;
            end
        end
        % dependent properties calculated on request
        function wn = get.wn(obj)
            wn = sqrt(obj.k/obj.m);    % undamped natural frequency [rad/s]
        end
        function zeta = get.zeta(obj)
             zeta = obj.c/(2*obj.m*obj.wn); % damping ratio
        end
        function wd = get.wd(obj)
            wd = obj.wn*sqrt(1-obj.zeta^2); % damped natural frequency [rad/s]
        end
        function s = get.s(obj)
            s = complex( [-obj.zeta*obj.wn + obj.wn*sqrt(obj.zeta^2-1)
                -obj.zeta*obj.wn - obj.wn*sqrt(obj.zeta^2-1)] );
            % root of the characteristic eq
        end
        function TransferFunction = get.TF(obj)
            TransferFunction = tf(1, [obj.m, obj.c, obj.k]); % transfer function
        end
        % functions that outputs time responses
        function [step_t, step_x] = Step(obj, t, step_value)
            % step_value - the magnitude of the input force [N]
            % for unit step input set it to 1
            % t - time vector at which the step response is calculated
            if obj.zeta < 1 % underdamped
                theta = atan(sqrt(1-obj.zeta^2)/obj.zeta);
                x = step_value/obj.k * (1  - exp(-obj.zeta*obj.wn.*t) .* sin(obj.wd.*t + theta) ./sqrt(1-obj.zeta^2));
 
            elseif obj.zeta == 1 % critically damped
                x = step_value/obj.k * (1- (1+obj.wn.*t).*exp(-obj.wn.*t));
 
            else
                alpha1 = obj.wn * (obj.zeta + sqrt(obj.zeta^2-1));
                alpha2 = obj.wn * (obj.zeta - sqrt(obj.zeta^2-1));
 
                x =  (step_value/obj.wn^2) * ...
                    (1 - (alpha1.*exp(-alpha2.*t) - alpha2.*exp(-alpha1.*t) )./(alpha1-alpha2) );
            end
            step_t = t;
            step_x = x;
        end
        function [impulse_t, impulse_x] = Impulse(obj, t, scale_value)
            % scale_value - the multiplier of the scaled impulse
            % for unit impulse set it to 1
            if obj.zeta < 1 % underdamped
                x = scale_value/ (obj.m * obj.wd) * exp(-obj.zeta.*obj.wn.*t) .* sin(obj.wd.*t);
 
            elseif obj.zeta == 1 % critically damped
                x = scale_value/obj.m .*t.*exp(-obj.wn.*t);
 
            else
                alpha1 = obj.wn * (obj.zeta + sqrt(obj.zeta^2-1));
                alpha2 = obj.wn * (obj.zeta - sqrt(obj.zeta^2-1));
 
                x =  (scale_value/(alpha1-alpha2)) * ...
                    ( exp(-alpha2.*t) - exp(-alpha1.*t) );
            end
            impulse_t = t;
            impulse_x = x;
        end
    end
end

(2) Example script file (for \(\zeta\) sweep with the step response of \(f(t) = 9 \) N)

The below Matlab script stands as an example of the method to plot the time responses and the root locus using the msd class defined above.

In the first section of the code, several configurations of mass-spring-dampers with changing damper coefficients \(c\) are set, then the time vector against which the time response is calculated are given, along with the magnitude of the input force after the step event \(f(t)\).

In the second section, ‘msd’ instances are created for each of the different configurations of the mass-spring-dampers and the step responses are calculated. The last section of the code allows the visual comparison of the mass-spring-damper systems. Figure (1) is the Root Locus and Figure (2) is the time response based on the functions in the msd class definition. Figure(3) outputs the same time response as Figure (2). The differences are Figure (3) utilises the transfer function of each msd instance in conjunction with the built-in function ‘step’ in the Matlab Control System Toolbox to confirm the result. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

clear all
%% inputs
% specification of the mass spring dampers.  m[kg], c[Ns/m], k[N/m]
mcks = [1 0 9
    1 1.2 9
    1 3.6 9
    1 4.8 9 
    1 6.0 9
    1 7.2 9];
 
t = [0:0.05:5].'; % time vector for plot of time response
 
step_value = 9; % the magnitude of the input force [N] after step event
%% create msd instances and calc time responses
num_msds = size(mcks,1); % the number of systems
 
%initialise data store arrays
msds(num_msds,1) = msd; % initialise the object array to store instances
step_ts = cell([num_msds, 1]);
step_xs = cell([num_msds, 1]);
 
for i = 1:num_msds
    msds(i) = msd(mcks(i,1), mcks(i,2), mcks(i,3));
    %create instances of mass-spring-damper systems using 'msd' class
    [step_ts{i}, step_xs{i}] = Step(msds(i),t, step_value);
    %calculate the response to step input f = step_value [N]
end
clear i
 
legend_zetas = strings(num_msds,1); % prepare string array for plot legend
zeta_values = [msds.zeta];
for i = 1:num_msds
    legend_zetas(i) = sprintf('zeta = %.2f', zeta_values(i));
end
clear zeta_values
 
%% plots
figure(1), clf, hold on
plot([msds.s], 'x', 'MarkerSize', 8)
legend(legend_zetas, 'Location', 'SouthEast')
FormatPlot()
xlim([-5 5]), ylim([-5 5])
axis equal
title('The roots s of the characteristic equation')
xlabel('Real'), ylabel('Imaginary')
 
figure(2), clf %step response using the msd class
hold on
for i = 1:num_msds
    plot(step_ts{i}, step_xs{i})
end
clear i
legend(legend_zetas, 'Location', 'SouthEast')
grid on
title('Step response')
xlabel('Time t [s]'), ylabel('Displacement x [m]')
FormatPlot()
 
figure(3), clf, hold on % plot step response using Control System Toolbox
opt = stepDataOptions('StepAmplitude',step_value); % alter the step amplitude
for i = 1:num_msds
    step(msds(i).TF, t, opt) % plot time response using built-in function
end
legend(legend_zetas, 'Location', 'SouthEast')
FormatPlot()

The script uses a small auxiliary function FormatPlot() below, stored in the same directory, to adjust the line width and font size used in the figure windows.

1
2
3
4
5

function[]= FormatPlot()
    set(findobj(gcf,'type','line'),'LineWidth', 2);
    set(findobj(gcf,'type','axes'),'FontSize',16, 'LineWidth', 0.5);
    grid on
end

Resulting plots from the example script

In the above example script, six mass-spring-dampers with varying damper coefficient were configured, while the mass and spring coefficient were fixed at \(m=1\) kg and \(k=9\) N/m. The damper coefficients range from \(c=0\) to \(c=7.2\) Ns/m, corresponding with no damping \(\zeta = 0 \) to overdamping \(\zeta = 1.2\).

The six mass-spring-dampers share the same natural frequency of \[ \omega_{n} = \sqrt{\dfrac{k}{m}} = \sqrt{\dfrac{9}{1}} = 3 \; \text{rad/s} \]

(1) Root Locus

The figure is the root locus formed by the six mass-spring-dampers. The two conjugate roots \(s\) of the characteristic equation \(ms^2+cs+k =0\) lie on the imaginary axis in case of no damping \(\zeta = 0\). With an increase in damping, the conjugate roots shift towards the middle of the figure, while keeping the constant distance \(\omega_n = 3 \) to the origin.
At critical damping \(\zeta = 1\), the two roots meet on the real axis. With further increase in damping the system become overdamped \(\zeta > 1\). One root moves to the left on the imaginary axis, while the other root moves to the right on the same axis.

(2) Time response with step input \(f(t)= 9, t>0\)

The figure is the time responses of the six mass-spring-dampers. Since the spring coefficient is \(k = 9 \;\text{Nm}\) for all the configurations, the step force input of the amplitude \(f = 9\; \text{N}\) results in a unit steady-state response of \(x =\dfrac{f}{k} = 1 \; \text{m}\).

For no damping case \(\zeta = 0\), the oscillation does not reduce in its amplitude. With an increase in damping, the peak value reduces. At critical damping \(\zeta = 1\), there is no oscillation. With overdamping \(\zeta = 1.2\), there is still no oscillation, however, it takes a longer time to settle to the steady-state than the critical damping. It is notable that at \(\zeta = 0.8\), the system largely reaches the steady-state value quicker than the critical damping except for some minor overshoot.

(3) Time response with unit impulse input \(f(t)= \delta, t=0\)

The figure is the time responses of the same six mass-spring-dampers subject to the unit impulse input. Similar to the step input response, the oscillation does not reduce for no damping \(\zeta = 0\) case. The peak value reduces with an increase in damping. At \(\zeta = 0.8\), the steady-state is reached quicker than the critical damping \(\zeta =1\), however, has the higher first peak.