Free Vibration

By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Complex Number

$$ a+jb \qquad a+bj \qquad a+ib \qquad a+bi $$

$$\text{Complex number = vector in 2D}$$

1.1. Adding Complex Numbers

$$\begin{align*} z_1 = a_1 + b_1 i, \qquad \vec{z}_1 = \begin{bmatrix}a_1 \\ b_1 \end{bmatrix}\\ z_2 = a_2 + b_2 i, \qquad \vec{z}_2 = \begin{bmatrix}a_2 \\ b_2 \end{bmatrix}\\ \end{align*}$$

$$\begin{align*} z = z_1 + z_2 &= (a_1 + a_2) + (b_1 + b_2)i\\ \\ \vec{z} = \vec{z}_1 + \vec{z}_2 &= \begin{bmatrix}a_1 \\ b_1 \end{bmatrix} + \begin{bmatrix}a_2 \\ b_2 \end{bmatrix} = \begin{bmatrix}a_1 + a_2 \\ b_1 + b_2 \end{bmatrix} \end{align*}$$

1.2. Complex Exponential Function (Euler's Formular)

$$ \large e^{i\theta} = \cos\theta + i\sin\theta $$

$$\begin{align*} \vec{z} &= r\cos\theta + i\,r\sin\theta \\ \vec{z} &= re^{i\theta} = r\,(\cos\theta + i\sin\theta)\\ \\ r &\text{ : magnitude (length)}\\ \theta &\text{ : phase (angle)} \end{align*}$$

1.3. Multiplying/Dividing Complex Numbers

• Just like an exponential function

$$ \begin{align*} z_1 = r_1e^{i\theta_1} \\ z_2 = r_2e^{i\theta_2} \end{align*}$$

$$ \implies$$

$$\begin{align*} z_1 \cdot z_2 & = r_1r_2 e^{i(\theta_1 + \theta_2)}\\ \frac{z_1}{z_2} & = \frac{r_1}{r_2} e^{i(\theta_1 - \theta_2)} \end{align*}$$

• Pictorial explanation

2. Circular Motion Represented by Complex Number

2.1. Geometrical Meaning of $e^{i\theta}$

$$ e^{i\theta} = \cos\theta + i\sin\theta$$

$$\begin{align*} e^{i\theta} & \text{: point on the unit circle with angle of } \theta \\ e^{i\theta} &= e^{i(\theta + 2\pi)} = \cdots = e^{i(\theta + 2\pi k)} = \cdots\\ \\ &\text{static} \end{align*}$$

2.2. Geometrical Meaning of $e^{i\omega t}$

$$\begin{align*} \theta &= \omega t\\ e^{i\omega t} &: \text{ rotating on an unit circle with angular velocity of }\omega \\ \\ &\text{dynamic} \end{align*}$$

2.3. Sinusoidal Functions from Circular Motions (sin & cos)

$$\text{ projection of } e^{j \omega t} \text{ onto Re-axis}$$	$$\text{ projection of } e^{j \omega t} \text{ onto Im-axis}$$
$$\text{Re}\left(e^{i\omega t}\right) = \cos\omega t$$	$$\text{Im}\left(e^{i\omega t}\right) = \sin\omega t$$

• Animated explanation of

$$ \cos \omega t = \frac{e^{i \omega t} + e^{-i \omega t}}{2}$$

2.4. $i$ Multiplying

$$ \begin{align*} ie^{i\theta} & = {?}\\ \\ z_1 & = i = e^{i \frac{\pi}{2}}\\ z_2 & = e^{i\theta} \end{align*} $$

$$z_1 \cdot z_2 = e^{i(\frac{\pi}{2} + \theta)}$$

• $i$ multiplication $\implies$ $90^{o}$ rotation forward

Example

$$\begin{align*} z &= e^{i\theta} \\ z^n &= \left(e^{i\theta}\right)^n = e^{in\theta} \end{align*}$$ find the solutions of $z^n = 1$

2.5. Circular Motion

• Particle rotates on the circle with angular velocity of $\omega$

$$ p(t) = re^{i\omega t} $$

2.5.1. Velocity in Circular Motion

$$\begin{align*} \upsilon(t) &= \frac{dp(t)}{dt} = r\cdot i \omega e^{i\omega t}= i r \omega e^{i\omega t} \\ \\ \lvert \upsilon(t) \rvert &= r\omega\\ \angle \upsilon(t) &= \omega t + \frac{\pi}{2} \end{align*}$$

2.5.2. Acceleration in Circular Motion

$$\begin{align*} a(t) &= \frac{d \upsilon(t)}{dt} = r\omega i \cdot i \omega e^{i\omega t}= - r \omega^2 e ^{i\omega t} \\ \\ \lvert a(t) \rvert &= r\omega^2\\ \angle a(t) &= \omega t + \pi \end{align*}$$

Recap: Vector Representation in Circular Motion

$$\text{ How to represent } v \text{ and } a \text{ in a cross product form}\\$$$$\begin{align*} \vec v &= \vec \omega \times \vec r\\ \vec a &= \vec \omega \times \vec v\\ & = \vec \omega \times (\vec \omega \times \vec r ) \end{align*}$$

3. Simple Harmonic Motion

• At this moment, do not worry about the meaning of the simple harmonic motion

3.1. Uniform circular motion

• For $ z(t) = r e^{j\omega t} $ $$ \begin{align*}z(t) &= r e^{j\omega t}\\ \upsilon(t) &= \frac{dz(t)}{dt} = j\omega r e^{j \omega t} \\ a(t) & = \frac{d^2z(t)}{dt^2} = -r\omega^2 e^{j \omega t} \end{align*}$$

$$ -r\omega^2 e^{j \omega t} + \omega^2 \left( r e^{j\omega t} \right) = a(t) + \omega^2 z(t) = \frac{d^2z(t)}{dt^2} + \omega^2 z(t) = 0$$

• Show $ z(t) = r e^{-j\omega t} $ also satisfies $$\frac{d^2z(t)}{dt^2} + \omega^2 z(t) = 0$$

• Given the differential equation

$$\ddot z(t) + \omega^2 z(t) = 0$$

• Solution is a linear combination of

$$ z(t) = c_1 r e^{j\omega t} + c_2 r e^{-j\omega t} = A e^{j\omega t} + B e^{-j\omega t}$$

• $A,B$ are determined by initial conditions

3.2. Spring and mass system

$$ \begin{align*} m \ddot x(t) + k x(t) &= 0 \implies \\ \\ \ddot x(t) + \omega^2 x(t) &= 0 \qquad \text{where } \; \omega^2 = \frac{k}{m} \end{align*}$$

%%processing

Ball ball;
PImage spring;

void setup() {
  size(600, 100);
  background(255);
  PVector loc = new PVector(width/2 + 100, height/2);
  ball = new Ball(loc);
  spring = loadImage("./image_files/spring_image_02.png");
}

void draw() {
  background(255);  
  fill(0);
  rect(0, 0, 100, height);
  ball.run();
  image(spring, 96, height/2 - 20, ball.location.x - ball.d/2 - 92, 40);
  
}

class Ball {
  
  int d = 60;
  PVector location;
  PVector velocity;
  PVector acceleration;
  PVector center = new PVector(width/2, height/2);
  
  Ball(PVector loc) {
    location = loc.get();
    velocity = new PVector(0, 0);
    acceleration = new PVector(0, 0);
  }
  
  void run() {
    applyForce();
    update();
    display();
  }
  
  void display() {
    noStroke();
    fill(175);
    ellipse(location.x, location.y, d, d);
  }
  
  void applyForce() {
    PVector force = PVector.sub(location, center);
    force.normalize();
    force.mult(-1);
    force.mult(0.5);
    acceleration = force;
  }
  
  void update() {
    velocity.add(acceleration);
    velocity.limit(100);
    location.add(velocity);
  }  
}

3.2.1. Differential equation

• $2^{nd}$ order ODE (ordinary differential equation)

$$\ddot x(t) + \omega^2 x(t) = 0$$

• Note: $\omega$ is a system property (i.e., $\omega$ defines the system)

• the solution is given

$$ x(t) = \text{Re} \left\{ A e^{j\omega t} + B e^{-j\omega t} \right\} = \text{Re} \left\{ r e^{j(\omega t + \phi)} \right\} = r\cos(\omega t + \phi)$$

• two initial conditions determine the future motion (i.e., two unknown parameters $r$ and $\phi$)

$$ \begin{cases} x(0) = x_0 \\ \dot x(0) = \upsilon_0 \end{cases} $$

• No external force
  • system contains: $m$, $k$
  • spring force $-kx$ is an internal force

3.2.2. Determine unknown coefficients

• How to obtain $r, \phi$ from $x_0, \upsilon_0$

$$\begin{align*} x(t) &= r\cos (\omega t + \phi) & \quad x(0) &= r \cos \phi = x_0\\ \dot x(t) &= - r\omega\sin(\omega t + \phi) & \quad \dot x(0) &= - r \omega \sin (\phi) = - \upsilon_0 \end{align*}$$

3.2.3. Pictorially determine unknown coefficients from circle

$$\begin{align*} R\cos \phi = x_0 &\\ \upsilon_0 = \upsilon \sin\phi &= R \omega \sin \phi \end{align*}$$

Example 1: Initial Conditions (pictorial interpretation)

$$x(0) = x_0$$ $$\upsilon(0) = 0$$

Example 2: Initial Conditions (pictorial interpretation)

$$x(0) = x_0$$ $$\upsilon(0) = \upsilon_0$$

3.3. Pendulum

$$\begin{align*} -T + mg \cos \theta & = -ml\omega^2\\ -mg\sin\theta &= ma = ml\ddot\theta \end{align*}$$

From Nonlinear to Linear

$$ \begin{align*} -mg\sin\theta = ma = ml\ddot\theta \implies & \ddot \theta + \frac{g}{l}\sin \theta = 0\\ \implies & \ddot \theta + \frac{g}{l}\theta = 0, \qquad \omega^2 = \frac{g}{l} \end{align*}$$

• non linear system approximation possible?

• period is independent of mass (non-intuitive)

%%html
<iframe 
width="420" height="315" 
src="https://www.youtube.com/embed/YlSC2DfphaY" frameborder="0" allowfullscreen>
</iframe>

Note

• 'Equation of motion' means

  • Apply $F=ma$ to find ordinary differential equations
  • Its solution will describe how it moves

• 'Simple harmonic motion' means

$$ \ddot x + \omega_n^2x = 0$$

4. Simulation of Free Vibration

$$ z(t) = e^{j \omega t} = \cos \omega t + j \sin\omega t$$

• $\omega$: angular velocity, [rad/sec]
• $f$: frequency, [rev/sec = Hz]

$$\omega = 2\pi f$$

• One revolution per sec

$$ \begin{align*} \omega &= 2\pi \\ f &= 1 \end{align*}$$

t = 0:0.01:1;
f = 1;
w = 2*pi*f;
z = exp(1j*w*t);
plot(real(z),imag(z),'.'), axis equal, ylim([-1.1 1.1])

subplot(2,1,1), plot(t,real(z),'.')
subplot(2,1,2), plot(t,imag(z),'.')

5. Beats

$$ z(t) = e^{j \omega_1 t} + e^{j \omega_2 t} \quad \text{when }\;\omega_1 \approx \omega_2$$

%%html
<iframe 
width="560" height="315" src="https://www.youtube.com/embed/IYeV2Wq82fw" frameborder="0" allowfullscreen>
</iframe>

t = 0:0.005:3;
f1 = 4;
f2 = 5;
w1 = 2*pi*f1;
w2 = 2*pi*f2;

A = exp(1j*w1*t);
B = exp(1j*w2*t);

subplot(3,1,1), plot(t,real(A))
subplot(3,1,2), plot(t,real(B))
subplot(3,1,3), plot(t,real(A+B))

plot(real(A+B),imag(A+B)), axis equal, xlim([-3,3])

fs = 44100; % sampling frequency
T = 5;      % seconds
t = linspace(0,T,fs*T); % time variable
x1 = exp(1j*2*pi*440*t); % pure sine wave at 440 Hz
x2 = exp(1j*2*pi*441*t); % pure sine wave at 441 Hz

sound(real(x1),fs)

sound(real(x2),fs)

x = x1 + x2;
sound(real(x),fs)

Example: The Bell of King Seongdeok (Emile Bell)

• The bell is considered a masterpiece of Unified Silla art. It is unique among Korean bronze bells because of the presence of a small hollow tube near the hook. The whole structure, including its decorative elements, produces a wide range of sound frequencies; the tube absorbes high frequency waves, contributing to a distinctive beat.

%%html
<iframe width="420" height="315" 
src="https://www.youtube.com/embed/7YEyMLX3sD8" frameborder="0" allowfullscreen></iframe>

fs = 44100; % sampling frequency
T = 5;      % seconds
t = linspace(0,T,fs*T); % time variable
r = 0.3;
x1 = exp(-r*t).*exp(1j*2*pi*440*t); % pure sine wave at 440 Hz
x2 = exp(-r*t).*exp(1j*2*pi*441*t); % pure sine wave at 441 Hz

x = x1 + x2;
sound(real(x),fs)

6. Damped Free Vibrarion

• more realisic

6.1. Experiment first

%%html
<iframe width="560" height="315" 
src="https://www.youtube.com/embed/vDulP6vTa9g" frameborder="0" allowfullscreen>
</iframe>

• in a mathematical form (from the educated guess)
  • exponentially decaying while oscillating

$$z(t) = e^{-\gamma t} e^{j \omega t}$$

r = 0.3;

f = 1;
w = 2*pi*f;

t = 0:0.01:10;
z = exp(-1*r*t).*exp(1j*w*t);

plot(real(z),imag(z)), axis equal, ylim([-1.1,1.1])

plot(t,real(z),t,exp(-1*r*t),'r--',t,-exp(-1*r*t),'r--')

• Assume damping causes exponential decay while oscillating

$$ \begin{align*}z(t) &= e^{-\gamma t} e^{-j\omega t} = e^{-(\gamma +j\omega)t}\qquad \text{normalized for simplicity}\\ \upsilon(t) &= \frac{dz(t)}{dt} = -(\gamma + j\omega) e^{-(\gamma +j\omega)t} \\ a(t) & = \frac{d^2z(t)}{dt^2} = (\gamma + j\omega)^2 e^{j \omega t} = (\gamma ^2 - \omega^2 + j2\gamma \omega) e^{-(\gamma +j\omega)t} \end{align*}$$

$$ \left\{ (\gamma ^2 - \omega^2 + j2\gamma \omega) e^{-(\gamma +j\omega)t} \right\} + 2\gamma \left\{ \left( -(\gamma + j\omega) e^{-(\gamma +j\omega)t}\right) \right\} + \left(\gamma^2 + \omega^2 \right) \left\{ e^{-(\gamma +j\omega)t} \right\} \\ = a(t) + 2\gamma \upsilon(t) + \left(\gamma^2 + \omega^2 \right) z(t) = \frac{d^2z(t)}{dt^2} + 2\gamma \frac{dz(t)}{dt} + (\gamma^2 + \omega^2) z(t) = 0$$

• Show $ z(t) = e^{-\gamma t} e^{j\omega t} $ also satisfies

$$\frac{d^2z(t)}{dt^2} + 2\gamma \frac{dz(t)}{dt} + \left(\gamma^2 + \omega^2 \right) z(t) = 0$$

• Given the differential equation

$$\ddot z(t) + 2\gamma \, \dot z(t) + \left(\gamma^2 + \omega^2 \right) z(t) = 0$$

• Solution is a linear combination of

$$ z(t) = e^{-\gamma t} \left(A e^{j\omega t} + B e^{-j\omega t}\right)$$

• $A,B$ are determined by initial conditions

6.2. Mass, spring, and damper system

$$ \begin{align*} m \ddot x(t) + c \dot x(t) + k x(t) &= 0 \implies \\ \ddot x(t) + 2\zeta \omega_n \dot x(t) + \omega_n^2 x(t) &= 0 \qquad \text{where } \; \omega_n^2 = \frac{k}{m} = \gamma^2 + \omega^2 \end{align*}$$

• parameters

$$ \begin{align*} \omega_n^2 &= \frac{k}{m} = \gamma^2 + \omega^2 & \text{: natural angular velocity}\\ \omega^2 &= \omega_n^2 - \gamma^2 & \text{: actual angular velocity}\\ \gamma & = \zeta \omega_n &\text{: decaying factor} \\ \omega^2 &= \omega_n^2 \left(1 - \zeta^2 \right)\\ 0 \leq \zeta &= \sqrt{1-\left(\frac{\omega}{\omega_n}\right)^2} \leq 1 & \text{: damping ratio} \end{align*} $$

• solution

$$ z(t) = e^{-\zeta \omega_n t}\left( A e^{j\omega_n \sqrt{1-\zeta^2}} + B e^{-j\omega_n \sqrt{1-\zeta^2}} \right) $$

%plot -s 800,250
f = 1;
wn = 2*pi*f;
zeta = [0.1 0.2 0.5 0.8];
t = 0:0.01:10;

for i = 1:4
    r = zeta(i)*wn;
    w = wn*sqrt(1-zeta(i)^2);
    z = exp(-1*r*t).*exp(1j*w*t);
    subplot(1,4,i), plot(real(z),imag(z)), grid on
    axis equal, axis([-1.1,1.1 -1.1 1.1])
    title(['\zeta = ',num2str(zeta(i))],'fontsize',8)
end

for i = 1:4
    r = zeta(i)*wn;
    w = wn*sqrt(1-zeta(i)^2);
    z = exp(-1*r*t).*exp(1j*w*t);
    plot(t,real(z)), hold on    
end
hold off
legend(['\zeta = ',num2str(zeta(1))],['\zeta = ',num2str(zeta(2))],...
['\zeta = ',num2str(zeta(3))],['\zeta = ',num2str(zeta(4))])

%%processing

Ball ball;
PImage spring;
float damping_c = 0.02;

void setup() {
  size(600, 100);
  background(255);
  PVector loc = new PVector(width/2 + 100, height/2);
  ball = new Ball(loc);
  spring = loadImage("./image_files/spring_damper_image.png");
}

void draw() {
  background(255);  
  fill(0);  
  rectMode(CORNER);
  rect(0, 0, 100, height);  
  ball.run();
  image(spring, 98, height/2 - 20, ball.location.x - ball.d/2 - 96, 40);  
}

class Ball {
  
  int d = 60;
  PVector location;
  PVector velocity;
  PVector acceleration;
  PVector center = new PVector(width/2, height/2);
  
  Ball(PVector loc) {
    location = loc.get();
    velocity = new PVector(0, 0);
    acceleration = new PVector(0, 0);
  }
  
  void run() {
    springForce();
    damping();
    update();
    display();
  }
  
  void display() {
    noStroke();
    fill(150);
    rectMode(CENTER);
    rect(location.x, location.y, d, d);
  }
  
  void damping() {
    float normal = 1;
    PVector damping = velocity.get();
    damping.normalize();
    damping.mult(-1);
    damping.mult(damping_c*normal);
    acceleration.add(damping);
  }
  
  void springForce() {
    PVector force = PVector.sub(location, center);
    force.normalize();
    force.mult(-1);
    force.mult(0.3);
    acceleration.add(force);
  }
  
  void update() {
    velocity.add(acceleration);
    velocity.limit(100);
    location.add(velocity);
    acceleration.mult(0);
  }  
}

Example: Door lock

%%html
<iframe width="560" height="315" 
src="https://www.youtube.com/embed/C2P_5Q8XROk" frameborder="0" allowfullscreen>
</iframe>

Example: Torsional Pendulum

%%html
<iframe width="420" height="315" 
src="https://www.youtube.com/embed/pnhE_-fW0pE" frameborder="0" allowfullscreen>
</iframe>

Example: unicopter

%%html
<iframe width="560" height="315" 
src="https://www.youtube.com/embed/N8UX_TDMn-Y" frameborder="0" allowfullscreen>
</iframe>

%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')