4
$\begingroup$

The Problem

I have the following system for which I would like to determine its dynamic response, taken from these lecture notes pg.20:

spring-mass_system

where the forcing functions are $u_1 = 1-U(t - t_0)$, $u_2=0$, and $u_3=0$, $U(\cdot)$ is the unit step function and $t_0 = 2\pi/\omega_1$ where $\omega_1$ is the first undamped natural frequency.


My Attempt

I know the time domain representation of the system to be of the typical format as shown in Eq.(1), and therefore the Fourier transform is as shown in Eq.(2):

\begin{gather} \mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{Kx} = \mathbf{u}(t) \tag{1} \\ % \Rightarrow [-\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}]\mathbf{X}(i\omega) = \mathbf{U}(i\omega) \tag{2} \end{gather}

where since we know that $\mathcal{L}\{1 - U(t - t_0)\} = 1 - \frac{\exp{(-st_0})}{s}$, then the Fourier transform of the right hand side of Eq.(2) is:

$$\mathbf{U}(i\omega) = \left[\left(1 - \frac{e^{-i\omega t_0}}{i\omega}\right),\, 0,\, 0,\right]^T$$

As a result, the dynamic response would be:

$$\mathbf{X}(i\omega) = [-\omega^2\mathbf{M} + i\omega\mathbf{C} + \mathbf{K}]^{-1}\mathbf{U}(i\omega) \tag{3}$$

Since the response by which I am comparing my answer is given in generalised modal coordinates $\mathbf{q}$, then the modal transformation of Eq.(3) follows where $\mathbf{V}$ is the eigenvector matrix thus we know that the modal mass, stiffness, damping matrices should be $\bar{\mathbf{M}} = \mathbf{I}$ (i.e. the identity matrix), $\bar{\mathbf{K}} = \mathbf{\Omega}^2$ (i.e. a diagonal marix containing the square of the natural frequencies), and $\bar{\mathbf{C}} = \mathbf{V}^T\mathbf{C}\mathbf{V}$ respectively. So Eq.(3) becomes:

$$\mathbf{q}(i\omega) = \left[-\omega^2\mathbf{I} + i\omega\bar{\mathbf{C}} + \mathbf{\Omega}^2 \right]^{-1}\left[\mathbf{V}^T\mathbf{U}(i\omega)\right] \tag{4}$$

Thefore, the dynamic response would be the magnitude of the imaginary numbers obtained from Eq.(4) above, for all positive $\omega$.

However, when I follow the above steps and solve the problem in MATLAB, the answer is far from the one provided in the excercise. Admitedly, the code is very simple as well (see below) so I struggle to see where I could have gone wrong:

clear ; clc

% Input
n = 3 ;         % Number of DOFs
m = 1 ;         % Mass [kg]
k = 1 ;         % Stiffnes [N/m]
c = 0.2 ;       % Viscous damping [Ns/m]
nv = [1 0 0]' ; % Load application vector (load applied to 1st DOF)

% Matrices
M = diag(m*ones(n, 1)) ;
K = diag(2*k*ones(n,1)) - diag(k*ones(n-1,1),-1) - diag(k*ones(n-1,1),1) ;
C = diag(c*ones(n,1)) ;

% Undamped Natural Frequencies & Modes
[V, Omsq] = eig(K, M) ;
Om = diag(sqrt(abs(Omsq))) ; % Natural Frequencies

% Frequency Aanalysis
t0 = 2*pi/Om(1) ;           % Time where unit step is applied
omega = linspace(0, 3, 100) ; % Frequency vector

for i = 1:numel(omega)
    
    s = 1i*omega(i) ;
    
    Us = V'*((1 - exp(-s*t0)/s)*nv) ;
    FRF(:,i) = (s^2*eye(n) + s*(V'*C*V) + Omsq)\Us ;
end

% Plotting
figure ; plot(omega, abs(FRF)) ; grid on ; set(gca, 'YScale', 'log') ; ylim([0 10]) ; 

enter image description here

$\endgroup$
4
  • $\begingroup$ Since this infolves frequency analysis and frequency response functions, I figured I might find an audience here, but your suggestion seems to make sense. Would you know if I can transfer this question there or do I need to repost from the beginning? $\endgroup$
    – kostas1335
    Oct 17 at 21:38
  • $\begingroup$ Actually you have at least a frequency analysis mistake, and a programming error here... $\endgroup$
    – Fat32
    Oct 17 at 21:43
  • $\begingroup$ Migrated here at the request of the author. $\endgroup$
    – Peter K.
    Oct 17 at 21:49
  • 1
    $\begingroup$ @PeterK. It's definetely fitting for control theory, and the closest site probably i robotics.se, but this place, or dsp.se could also work... $\endgroup$
    – Fat32
    Oct 17 at 21:53
2
$\begingroup$

As TimWescott indicated, this problem is more suitable for control theory for a fuller understanding than for digital signal processing, yet there's no control.se as far as I know. Hence you can give it a try here too since such damped-spring systems is not uncommon subject for audio / acoustic engineers etc.

First of all you have a programming error: Your for-loop variable is set to i , a potential error because there's already the complex unit $j$ in your body code in the line $$s = j \omega \tag{1}$$which you have implemented as :

s = 1i*omega(i)

which is wrong: you are using the same letter i , both for indexing integer into the angular frequency vector omega, and also for the imaginary unit $j$. I've replaced the for-loop variable with k instead.

Second, you have a signal processing mistake. Your nonzero input function is $$ u_1(t) = 1 - u(t-t_0) \tag{2}$$whose (bilateral - double sided) Laplace transform is

$$ - \frac{ e^{ -s t_0} }{s} \tag{3} $$

You can see this by noting that the input is actually: $$u_1(t) = u(t_0 - t) \tag{4}$$

And the third issue is about the diagonalized computation. I'have replaced the diagonalised frequency response calculation, with the undiagonalised one with the matrcices M,K, C. The result is improved,as shown on the plot below, but still not exactly what you have been given as the answer plot. I guess you may correct it further.

I put the corrected code and the corresponding MATLAB output below :

clc; clear variables ; close all

% Input
n = 3 ;         % Number of DOFs
m = 1 ;         % Mass [kg]
k = 1 ;         % Stiffnes [N/m]
c = 0.2 ;       % Viscous damping [Ns/m]
nv = [1 0 0]' ; % Load application vector (load applied to 1st DOF)

% Matrices
M = diag(m*ones(n, 1)) ;
K = diag(2*k*ones(n,1)) - diag(k*ones(n-1,1),-1) - diag(k*ones(n-1,1),1) ;
C = diag(c*ones(n,1)) ;

% Undamped Natural Frequencies & Modes
[V, Omsq] = eig(K,M); 
Om = diag(sqrt(abs(Omsq))) ; % Natural Frequencies

% Frequency Aanalysis
t0 = 2*pi/Om(1) ;           % Time where unit step is applied
omega = linspace(0, 2.5, 100) ; % Frequency vector

for k = 1:numel(omega)
    
    s = 1i*omega(k) ;
           
    Us = -(exp(-s*t0)/s)*nv; 
    FRF(:,k) = (s^2*M + s*C + K)\Us ;

    %Us = V'*((1 - exp(-s*t0)/s)*nv) ;
    %FRF(:,k) = (s^2*eye(n) + s*(V'*C*V) + Omsq)\Us ;

end

% Plotting
figure ; plot(omega, abs(FRF)) ; grid on ; set(gca, 'YScale', 'log') ; ylim([0.01 10]) ;

enter image description here

$\endgroup$
4
  • $\begingroup$ Thanks for the above. There's one thing I don't get though; why do you calculate the undamped natural frequencies using eig(M, C) instead of eig(K, M) ? $\endgroup$
    – kostas1335
    Oct 18 at 7:45
  • $\begingroup$ Sorry, say it's a typo there. Indeed. Btw, when I've written the answer I hadn't actually read the link you've posted. Indeed now I see that, the author gave an explicit solution of the problem you ask. Yet he doesn't provide a matlab code for the simulation? Furthermore, I see that his definition of input u1(t) is problematic: 1-u(t-t0) or u(t)-u(t-t0) ? He writes the 1st, but seems to mean the 2nd: left ambigious. Hence he should have provided the code so that we could really see it. May be he used different settings in his simulation. Also the Modal approach shoul yield the same result. $\endgroup$
    – Fat32
    Oct 18 at 15:18
  • $\begingroup$ "I see that his definition of input u1(t) is problematic: 1-u(t-t0) or u(t)-u(t-t0) ?". Yup. Particularly when you're dealing with the Fourier transform. If you're playing with the single-sided Laplace transform, it's safe to assume that "x(t) = 1" really means "x(t) = u(t)". But the Fourier transform is always two-sided, so a signal could really be 1 from $t = -\infty$ until $t = t_0$, and zero thereafter, not (as the author probably intended) $x(t) = 1$ for $0 \le t < t_0$ and zero everywhere else. In the latter case, $x(t) = u(t) - u(t - t_0)$. $\endgroup$
    – TimWescott
    Oct 18 at 15:33
  • $\begingroup$ @TimWescott exactly, if he already has $u(t)$ in the formula then it should be $u(t)-u(t-t_0)$, rather than the ambiguous $1 - u(t-t_0)$ if he means a causal input definition, etc... $\endgroup$
    – Fat32
    Oct 18 at 15:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.