Engineering at Alberta Courses » Lagrange’s equations applied to dynamic systems

Advanced Dynamics and Vibrations: Lagrange’s equations applied to dynamic systems

Analytical Mechanics – Lagrange’s Equations

Up to the present we have formulated problems using newton’s laws in which the main disadvantage of this approach is that we must consider individual rigid body components and as a result, we must deal with interaction forces that we really have no interest in. These forces come from constraints which one part of the system puts on another.

An alternate approach referred to as analytical mechanics, considers the system as a whole and formulates the equations of motion starting from scalar quantities – the kinetic and potential energies and an expression of the virtual work associated with non-conservative forces. It relies heavily on the concept of virtual displacements.

Any set of numbers that serve to specify the configuration of the system are examples of generalized coordinates. To these we can add constraints. The number of DOF is then the number of generalized coordinates minus the number of independent equations of constraint.

e.g a point moving on the surface of a sphere of radius R with center $(x_{o},y_{0},\zeta _{0})$ must obey:

$(x-x_{0})^{2} + (y-y_{0})^{2} + (z- z _{0})^{2} = R^{2}$

There are 3 coordinates and 1 equation of constraint. Therefore it has 2 DOF. If we used two angles, say $\theta$ & $\phi$ then we would not need the constraint.

These constraints are classified as:

$\textcircled{1}$ If the constraint can be written in finite form (as above) they are called holonomic if not (then in differential form) then they are non-holonomic.

$\textcircled{2}$ If they are equalities then they are bilateral and unilateral if inequalities

$\textcircled{3}$ If there is no explicit dependence on time then they are scleronomic and if there is explicit time dependence they are rheonomic.

$x_{1}^{2} + y_{1}^{2} = a^{2}$

$(x_{2} - x_{1})^{2} + (y_{2} - y_{1})^{2} = b^{2}$

Therefore, 2 DOF holonomic, scleronomic bilateral system.

$x_{1}^{2} +y_{1}^{2} + z _{1}^{2} \leq \ell^{2}$ (unilateral)

however if we allow the point of attachment to move (in the x-y plane)

$x_{1}^{2} + (y_{1} - y_{o} \sin \omega t)^{2} = \ell^{2}$ (rheonomic)

Vertical displacement & Virtual Work

The work done during the motion is

$\overline{d W} = \underline{F} \cdot d \underline{r}$

where $\underline{F}$ & $d\underline{r}$ are vectors and $\overline{dW}$ is a scalar. For a particle $\underline{F} = m \underline{\ddot{r}}$ where $\underline{F}$ may be the sum of several forces. Therefore,

$\begin{split} \overline{dW} &= m \ddot{\underline{r}} \cdot dr \\&= m \frac{d \dot{\underline{r}}}{dt} \cdot ( \frac{d \underline{r}}{dt}dt) \\&= m \frac{d \dot{r}}{dt} \cdot (\dot{r} dt) \\&= m \dot{\underline{r}} \cdot d \dot{r} \\&= d(\frac{1}{2} m \dot{\underline{r}} \cdot \dot{\underline{r}})\end{split}$

where this is a true differential of the kinetic energy.

$T = \frac{1}{2} m \dot{\underline{r}} \cdot \dot{\underline{r}}$

If the particle moves from positive $\underline{r_{1}}$ to $\underline{r_{2}}$ then

$\int_{\underline{r}_{1}}^{\underline{r}_{2}} \underline{F} \cdot d \underline{r} = \frac{1}{2} m [\underline{\dot{r}}_{2} \cdot \underline{\dot{r}}_{2}] - \frac{1}{2} m [\underline{\dot{r}}_{1} \cdot \underline{\dot{r}}_{1}]$

if $\underline{F}$ depends only on $\underline{r}$ , $\underline{F} = \underline{F}(\underline{r})$

$dW = F \cdot d \underline{r} \equiv -dV(\underline{r})$

where $V(\underline{r})$ is a scalar function that depends only on the position and not on $\underline{\dot{r}}$ or time. $V$ is term the potential energy.

$d(T+V) = 0$

The total energy is constant.

When we have a conservative system we can use Rayleigh’s to estimate the natural frequencies.

Consider a simple problem and use visual work to solve.

Two bars of length $l$ . What is the equilibrium position

Therefore,

$P \delta x + 2mg \delta h = 0$

$x = 2 \ell \sin \theta$

$h = \frac{\ell}{2} \cos \theta$

$\delta x = 2\ell \cos \theta \delta \theta$

$\delta h = \frac{- \ell}{2} \sin \theta \delta \theta$

$P(2 \ell \cos \theta \delta \theta) = 2mg(\frac{\ell}{2} \sin \theta \delta \theta)$

$\frac{P}{mg} = \tan \theta$

Consider a rigid beam (uniform) of mass $M$ loaded by a uniform time dependent load $f_{o}(t)$ and supported by a spring & dashpot as shown

Using $VW$ determine the equation of motion.

$\sum$ Virtual work of forces applied $= \sum$ Virtual work of inertia forces

$\frac{-kl \theta}{2}\Big( \frac{\ell \delta \theta}{2} \Big) - cl \dot{\theta}(l \delta \theta) + \int_{0}^{l} f_{o}(t)dx (x \delta \theta) = \frac{Ml^{2}}{12} \ddot{\theta} \delta \theta + M \frac{l \ddot{\theta}}{2} \Big( \frac{l}{2} \delta \theta \Big)$

Therefore,

$M \ell^{2}\Big[ \frac{1}{12} + \frac{1}{4}\Big] \ddot{\theta} \delta \theta + \frac{k \ell^{2}}{4} \theta \delta \theta + c \ell^{2} \dot{\theta} \delta \theta = f_{o}(t) \frac{ \ell^{2}}{2} \delta \theta$

$\frac{M}{3} \ddot{\theta} + c \dot{\theta} + \frac{k}{4} \theta = \frac{f(t)}{2}$

We can use virtual work to calculate the energy in a multi-degree of freedom system.

For a displacement of $\delta _1$ only

$F_{1} = k_{11} \delta _{1} \hspace{1cm} F_{2} = k_{21} \delta _{1} \hspace{1cm} F_{3} = k_{31} \delta _{1}$

and similarly for $\delta _2$ , $\delta _3$ alone

$F_{1} = k_{12} \delta _{2} \hspace{1cm} F_{2} = k_{22} \delta _{2} \hspace{1cm} F_{3} = k_{32} \delta _{2}$

$F_{1} = k_{13} \delta _{3} \hspace{1cm} F_{2} = k_{23} \delta _{3} \hspace{1cm} F_{3} = k_{33} \delta _{3}$

So for the combination (assuming superposition)

$F_{i} = \sum _{j=1}^{n} k_{ij} \delta _{j}$

and

$\delta _{i} = \sum _{j=1}^{n} a_{ij} F _{j}$

The work done by force $F_{i}$ in $\delta _{i}$ is energy stored due to that force:

$V_{i} = \frac{1}{2} F_{i} \delta _{i}$

If there are n forces and the displacements are $\delta _{i}$ (i = 1,2, … n) then the total energy is:

$\begin{split} V &= \sum _{i=1}^{n} V_{i} \\&= \frac{1}{2} \sum _{i=1}^{n} F_{i} \delta _{i} \\&= \frac{1}{2} \sum _{i=1}^{n} \Bigg( \sum _{j=1}^{n} k_{ij} \delta _{j} \Bigg) \delta _{i} \end{split}$

$\begin{split} V &= \frac{1}{2} \{\delta \}^{T} [k] \{ \delta\} \\&= \frac{1}{2} \{F \}^{T} [a] \{ F\} \end{split}$

We can also show that the kinetic energy is

$T = \frac{1}{2} \{\dot{\delta}\}^{T} [m] \{\dot{\delta}\}$

$T$ & $V$ are examples of so-called quadratic forms. Quadratic forms are said to be positive definite if they are never negative and only zero when the variables are zero ( $\{\delta\} = \{0\}$ )

They are said to be positively semi-definite if they are never negative but can be zero for non-negative variables.

This is useful to know because if we have a stiffness matrix that is positive semi-definite it means there is a rigid body mode with zero natural frequency $[$ [m] is positive definite $]$

For example a submarine:

has 2 rigid body modes

Principle of Virtual Work

We will consider a system of $N$ particles moving in a 3D space and define the virtual displacements as $\delta x_{i} ^{(1)}$ for particle 1 $\delta x_{i} ^{(2)}$ for particle 2 etc.

These are not true displacements but small variations in the so-called generalized coordinates. These are imagined displacements (slight variations) which are at the same time and which satisfy the constraints given by

$g_{j} (x_{1}, y_{1}, z _{1}, x_{2}, y_{2}, z _{2} ..... t) = c \hspace{1cm} j = 1\quad \text{to} \quad m$

Then for virtual changes in these coordinates:

$g_{j} (x_{1} + \delta x_{1} , y_{1} + \delta y_{1}, z _{1} + \delta z _{1} ..... t) = c$

Expanding this in a Taylor series gives:

$g_{j}(x_{1}, y_{1}, \zeta _{1},..... t) + \sum _{i=1} ^{N} \frac{\partial g_{i}}{\partial x_{i}} \delta x_{i} + \frac{\partial g_{j}}{\partial y_{i}} \delta y_{i} + \frac{\partial g_{j}}{\partial z_{i}} \delta z_{i} = C$

Therefore,

$\sum _{i=1} ^{N} \frac{\partial g_{j}}{\partial x_{i}} \delta x_{i} + \frac{\partial g_{j}}{\partial y_{j}} \delta y_{j} + \frac{\partial g_{j}}{\partial z _{j}} \delta z _{j} = 0 \hspace{1cm} \forall _{j}$

So that there are $3N - j$ arbitrary displacements

On each of these particles the forces acting are

$\underline{R}_{i} = \underline{F}_{i} + \underline{f}_{i} \hspace{1cm} i = 1,2, ... N$

Where the resultant force $\underline{R}_{i}$ is composed of the external forces $\underline{F}_{i}$ and $\underline{f}_{i}$ the constraint forces.

When all the particles are in static equilibrium

$\underline{R}_{i} = 0 \hspace{1cm} \forall _{i}$

and

$\underline{R}_{i} \cdot \delta \underline{r}_{i} = 0$

where this represents the virtual work performed by the resultant force on the $i$ particle.

$\delta \overline{W} = \sum _{i=1} ^{N} \underline{R}_{i} \cdot \delta \underline{r}_{i}$

$\delta \overline{W} = \sum _{i=1} ^{N} \left(\underline{F}_{i} \cdot \delta \underline{r}_{i} + \underline{f}_{i} \cdot \delta \underline{r}_{i}\right) = 0$

If we rule out work by the constraint forces ( as the constraint force is perpendicular to the motion), (therefore cannot handle friction).

Consider two particles held by a rigid massless rod.

$|\underline{f} _{1}| = |\underline{f} _{2}|$

$\delta W = \underline{f}_{1} \cdot \underline{\delta} r_{1} + \underline{f}_{2} \cdot \underline{\delta} r_{2}$

But because the rod is rigid

$\delta \underline{r}_{1} \cdot \underline{e}_{r} = \delta \underline{r}_{2} \cdot \underline{e}_{r}$

$\delta W = f_{1} \underline{e}_{r} \cdot \delta \underline{r}_{1} - f_{2} \underline{e}_{r} \cdot \delta \underline{r}_{2}$

Therefore for the constraint forces

$\begin{split} \delta W &= (f_{1} - f_{2}) \underline{e}_{r} \cdot \delta r_{2} \\&= 0 \end{split}$

and the Principal of Virtual Work says

$\begin{split} \delta W &= \sum _{i=1} ^{N} \underline{F}_{i} \cdot \underline{\delta} r_{i}\\& = 0 \end{split}$

that the work performed by the applied forces through infinitesimal virtual displacements compatible with the constraints is zero.

To make the virtual work equations more useful we write it in terms of the generalized coordinates $q_k$

$\begin{split} \delta W &= \sum _{i=1} ^{N} \underline{F}_{i} \cdot \delta r_{i} \\&= \sum _{i=1} ^{N} \underline{F}_i\cdot \sum _{j=1} ^{n} \frac{\partial \underline{r}_{i}}{\partial q_{j}} \delta q_{j} \\&= \sum _{j=1} ^{n} \Bigg( \sum _{i=1} ^{N} \underline{F}_i \cdot \frac{\partial \underline{r}_{i}}{\partial q_{j}} \Bigg) \delta q_{j} \end{split}$

We call

$Q_{j} =\sum _{i=1} ^{N} \underline{F}_{i} \cdot \frac{\partial \underline{r}_{i}}{\partial q_{j}}$

the generalized forces. These are not necessarily forces as they can be moments. It is only necessary that the quantity $Q_{k} \delta q_{k}$ has the units of work.

For the situation in which we have a conservative system the work can be expressed differently.

$\begin{split} \delta W &= \sum _{i=1} ^{N} \underline{F}_{i} \cdot \delta \underline{r}_{i} \\&= - \delta V \\&= - \sum _{i=1} ^{N} \Bigg( \frac{\partial V}{\partial x_{i}} \delta _{i} + \frac{\partial V}{\partial y_{i}} \delta y_{i} + \frac{\partial V}{\partial z _{i}} \delta z _{i} \Bigg) \\&= 0 \end{split}$

as we can select the virtual displacement arbitrarily (i.e. only one non-zero at a time)

$\underline{F}_{x_{i}} = -\frac{\partial V}{\partial x_{i}} = 0, \hspace{0.5cm} \underline{F}_{y_{i}} = -\frac{\partial V}{\partial y_{i}} = 0, \hspace{0.5cm} \underline{F}_{z _{i}} = -\frac{\partial V}{\partial z _{i}} = 0$

This states that $V$ has a stationary value at the equilibrium position (can use this for stability consideration). To extend this to dynamical situations we use D’Alembert’s principle.

$\underline{F}_{i} +\underline{f}_{i} - m_{i} \underline{\ddot{r}}_{i} = 0$

Where the $m_{i} \underline{\ddot{r}}_{i}$ are the inertia forces. The virtual work for a system of particles is:

$\sum _{i=1} ^{N} (\underline{F}_{i} +\underline{f}_{i} - m_{i} \underline{\ddot{r}}_{i}) \cdot \delta \underline{r}_{i} = 0$

and since we showed the virtual work of the interned forces is zero.

$\sum _{i=1} ^{N} (\underline{F}_{i} - m_{i} \underline{\ddot{r}}_{i}) \cdot \delta \underline{r}_{i} = 0$

This is called the generalized D ‘Alembert’s principle and is really not useful as it is still in vector form. To get it in a scalar form we rewrite the second term in terms of any generalized coordinates and this leads to Lagrange’s equations of motion.

Consider a coordinate transformation where radius vector is written in terms of $q_{i}$ instead $x_{i},y_{i}, z _{i}$ .

$\underline{r}_{i} = \underline{r}_{i}(q_{1} ...... q_{n}) \hspace{1cm} i = 1,2 ... N$

$\begin{split} \underline{\dot{r}}_{i} &= \frac{\partial \underline{r}_{i}}{\partial q_{1}} \dot{q}_{1} + \frac{\partial \underline{r}_{i}}{\partial q_{2}} \dot{q}_{2} + ... \\&= \sum _{k=1} ^{n} \frac{\partial \underline{r}_{i}}{\partial q_{k}} \dot{q}_{k} \end{split}$

Note: that since $\frac{\partial \underline{r}_{i}}{\partial q_{k}}$ do not depend explicitly on time that.

$\frac{\partial \underline{\dot{r}}_{i}}{\partial \dot{q}_{k}} = \frac{\partial \underline{r}_{i}}{\partial q_{k}} --- (+)$

now consider a virtual change in $\underline{r}_{i}$

$\begin{split} \delta \underline{r}_{i} &= \frac{\partial \underline{r}_{i}}{\partial q_{1}} \delta q_{1} + ... \\&= \sum _{k=1} ^{n} \frac{\partial \underline{r}_{i}}{\partial q_{k}} \delta q_{k} \end{split}$

and apply this to the second term of D ‘Alembert’s principle

$\begin{split} \sum _{i=1} ^{N} m_{i} \underline{\ddot{r}}_{i} \cdot \delta \underline{r}_{i} &= \sum _{i=1} ^{N} m_{i} \underline{\ddot{r}}_{i} \cdot \sum _{k=1} ^{n} \frac{\partial \underline{r}_{i}}{\partial q_{k}} \delta q_{k} \\&= \sum _{k=1} ^{n} \Bigg( \sum _{i=1} ^{N} m_{i} \underline{\ddot{r}}_{i} \cdot \frac{\partial \underline{r}_{i}}{\partial q_{k}} \Bigg) \delta q_{k} \end{split}$

Consider a typical term in the bracket:

$m_i \ddot{\underline{r}_i} \cdot \frac{\partial \underline{r}_i}{\partial q_k} = \frac{\text{d}}{\text{dt}} \Big(m_i \dot{\underline{r}_i} \cdot \frac{\partial \underline{r}_i}{\partial q_k}\Big) - m_i \dot{\underline{r}_i} \cdot \frac{\text{d}}{\text{dt}} \Big( \frac{\partial \underline{r}_i}{\partial q_k}\Big)$

Now interchange the order of differentiation in the last term.

$\begin{split} m_i \ddot{\underline{r}_i} \cdot \frac{\partial \underline{r}_i}{\partial q_k} &= \frac{\text{d}}{\text{dt}} \Big( m_i \dot{\underline{r}}_i \cdot \frac{\partial \underline{\dot{r}}_i}{\partial q_k}\Big) - m_i \dot{\underline{r}_i} \cdot \frac{\partial \underline{\dot{r}}_i}{\partial \dot{q}_k} \\&= \Big[\frac{\text{d}}{\text{dt}} \Big(\frac{\partial}{\partial \dot{q_k}}\Big) - \frac{\partial}{\partial q_k} \Big]\Big(\frac{1}{2}m_i \dot{\underline{r}}_i \cdot \dot{\underline{r}}_i\Big) \end{split}$

As a result, the total second term becomes:

$\begin{split} \sum m_i \ddot{\underline{r}}_i \cdot \delta\underline{r}_i &= \sum_{k=i}^n \biggr\{ \Big[ \frac{\text{d}}{\text{dt}} \Big(\frac{\partial}{\partial \dot{q_k}}\Big) - \frac{\partial}{\partial q_k} \Big] \Big( \sum_{i=1}^N \frac{1}{2} m_i \dot{\underline{r}}_i \cdot \dot{\underline{r}}_i \Big) \biggr\} \delta q_k \\&= \sum_{k=1}^n \biggr[ \frac{\text{d}}{\text{dt}} \Big( \frac{\partial T}{\partial \dot{q_k}} \Big) - \frac{\partial T}{\partial q_k} \biggr] \delta q_k \end{split}$

As $T = \frac{1}{2} \sum_{i=1}^N m_i \dot{\underline{r}}_i \cdot \dot{\underline{r}}_i = T(q_1,q_2,...,q_n,\ \dot{q_1},\dot{q_2},...,\dot{q_n})$ is the kinetic energy of the entire system. To complete the conversion, we must rewrite the forces.

$\begin{split} \delta \tilde{W} = \sum_{i=1}^N \underline{F}_i \cdot \delta \underline{r}_i &= \sum_{i=1}^N \underline{F}_i \cdot \sum_{k=i}^n \frac{\partial \underline{r}_i}{\partial q_k} \delta q_k \\&= \sum_{k=1}^n\Big(\sum_{i=1}^N \underline{F}_i \cdot \frac{\partial \underline{r}_i}{\partial q_k} \Big) \delta q_k \end{split}$

And we call:

$\sum_{i=1}^N \underline{F}_i \cdot \frac{\partial \underline{r}_i}{\partial q_k} \equiv Q_k \quad \text{(generalized forces)}$

Therefore:

$\delta \tilde{W} = \sum_{k=1}^n Q_k \delta q_k$

Some of the forces will be conservative (derivable from a potential $V = V(q_1,q_2,...,q_n)$ and some non-conservative).

$\begin{split} \sum \tilde{W} &= \delta W_c + \delta \tilde{W_{NC}} \\&= -\delta V + \sum_{k=1}^n Q_{k_{NC}} \delta q_k \\&= -\sum_{k=1}^n \Big( \frac{\partial V}{\partial q_k} - Q_{k_{NC}} \Big) \delta q_k \end{split}$

So that the rewritten D’Alembert’s principle becomes:

$\sum_{k=1} \Big[ \frac{\text{d}}{\text{dt}} \Big(\frac{\partial T}{\partial \dot{q_k}} \Big) - \frac{\partial T}{\partial q_k} + \frac{\partial V}{\partial q_k} - Q_{k_{NC}} \Big] \delta q_k = 0$

And since the displacements are arbitrary and independent.

$\frac{\text{d}}{\text{dt}} \Big( \frac{\partial T}{\partial \dot{q_k}} \Big) - \frac{\partial}{\partial q} (T - V) = Q_{k_{NC}} \quad k = 1,2,....,n$

However, since the potential energy does not depend on the velocities, this can be written as:

$\frac{\text{d}}{\text{dt}} \Big( \frac{\partial L}{\partial \dot{q_k}} \Big) - \frac{\partial L}{\partial q_k} = Q_k \quad k = 1,2,....,n$

$L \text{(Lagrangian)} = T - V$

NOTES:

Does not have to be only for “small” oscillations.
Does not matter whether the conservative forces are included in the $Q_k$ .
Sometimes the $Q_k$ are able to be put into a “potential like” function.
Can include reaction if wanted using Lagrange multipliers.
Determine the $Q_k$ using virtual work.

Example

Find the equations of motion:

$x = r\theta$

$\begin{split} T &= \frac{1}{2} m (\dot{x})^2 + \frac{1}{2}J(\dot{\theta})^2 \\&= \frac{1}{2} \Big(m + \frac{J}{r^2}\Big)\dot{x}^2 \end{split}$

$V = \frac{1}{2}kx^2$

$\delta W\ \text{(damper)}\ = -c\dot{x} \delta x$

$Q_x = -c\dot{x}$

$\frac{\text{d}}{\text{dt}}\Big(\frac{\partial L}{\partial \dot{q}} \Big) - \frac{\partial L}{\partial q} = Q_{\text{NC}}$

Thus:

$\Big(m + \frac{J}{r^2} \Big) \ddot{x} + kx = -c\dot{x}$

NOTE: Sometimes write the damping forces as a dissipation function.

$D = \frac{1}{2}c\dot{q}^2,\quad Q = \frac{\partial D}{\partial q} = c\dot{x}$

$\frac{\text{d}}{\text{dt}} \Big(\frac{\partial L}{\partial \dot{q_k}} \Big) - \frac{\partial L}{\partial \dot{q_k}} + \frac{\partial D}{\partial \dot{q_k}} = Q_{k_{\text{NC}}}$

Example

Find the equations of motion:

$J = \frac{1}{2} M r^2$

$T = \frac{1}{2} M \dot{x}^2 + \frac{1}{2}J\dot{\theta}^2 + \frac{1}{2}m\nu_A^2$

The potential energy is:

$V = mgy + \frac{1}{2}kx^2$

And the damping force is the same. Note that:

$x = r\theta,\quad y = \ell(1 - \cos \theta)$

$\dot{x} = r \dot{\theta}$

$\nu_A = ?$

$\underline{\nu}_A = \underline{\nu}_o + \underline{\nu}_{A/o}$

$\ell \dot{\theta} = |\underline{\nu}_{A/o}|$

$|\underline{\nu}_o| = \dot{x}$

$\begin{split} |\underline{\nu}_A| &= |\underline{\nu}_o|^2 + (\ell \dot{\theta})^2 - 2\nu_o\ell\dot{\theta}\cos\theta \\&= \dot{x}^2 + \frac{\ell^2}{r^2}\dot{x}^2 - 2\dot{x}^2\frac{\ell}{r} \cos \frac{x}{r} \end{split}$

Thus:

$T = \frac{\dot{x}^2}{2}\Big[\frac{3}{2}M + m\Big(1 + \frac{\ell^2}{r^2} - \frac{2\ell}{r}\cos\frac{x}{r}\Big)\Big]$

$V = mg\ell\Big(1-\cos\frac{x}{r}\Big) + \frac{1}{2}kx^2$

$\frac{\partial L}{\partial \dot{x}} = \dot{x}\Big[\frac{3M}{2} + m \Big(1+\frac{\ell^2}{r^2} - \frac{2\ell}{r}\cos\frac{x}{r}\Big)\Big]$

$\frac{\partial L}{\partial x} = \frac{\dot{x}^2}{2} \Big[ \frac{2\ell m}{r^2} \sin \frac{x}{r} \Big] -mg\frac{\ell}{r}\sin\frac{x}{r} - kx$

$\frac{\text{d}}{\text{dt}}\Big(\frac{\partial L}{\partial \dot x}\Big) = \ddot{x}\Big[\frac{3M}{2} + m \Big( 1 + \frac{\ell^2}{r^2} - \frac{2\ell}{r}\cos\frac{x}{r} \Big)\Big] + \dot{x}^2\frac{2m\ell}{r^2}\sin\frac{x}{r}$

The equation of motion comes from:

$\frac{\text{d}}{\text{dt}} \Big(\frac{\partial L}{\partial \dot{x}}\Big) - \frac{\partial L}{\partial x} = -c\dot{x}$

And becomes:

$\Big[ \frac{3M}{2} + m \Big(1+ \frac{\ell^2}{r^2} - \frac{2\ell}{r}\cos\frac{x}{r}\Big)\Big]\ddot{x} + \dot{x}^2\frac{m\ell}{r^2}\sin{\frac{x}{r}} + mg\frac{\ell}{r}\sin{\frac{x}{r}} + kx + c\dot{x} = 0$

For small deflection $\cos\frac{x}{r} \approx 1$ and $\sin \frac{x}{r} \approx \frac{x}{r}$ :

$\Big[\frac{3}{2}M + m\Big(1+\frac{\ell^2}{r^2}-\frac{2\ell}{r}\Big)\Big] \ddot{x} + \frac{m\ell x}{r^2}\Big[\frac{\dot{x}^2}{r}+g\Big] + kx + c\dot{x} = 0$

$\left(\text{Set} \left(\frac{\dot{x}^2}{r}\right) \approx 0 \right)$

2 DOF Example

The uniform rod of length $\ell$ and mass $m$ carries a slider of mass $m_o$ which is attached by a spring of stiffness $k$ . The spring is unstretched at $r = r_o$ .

$T = \frac{1}{2}m_o\big[(r\dot{\theta})^2 + (\dot{r})^2\big] + \frac{1}{2}\Big(\frac{m\ell^2}{3}\dot{\theta}^2\Big)$

$V = \frac{1}{2}k(r-r_o)^2-m_ogr\cos\theta - mg\frac{\ell}{2}\cos\theta$

$L = T - V \quad \text{(} \theta = \frac{\pi}{2} \text{ is the datum)}$

$\frac{\partial L}{\partial \dot{\theta}} = m_or^2\dot{\theta} + \frac{m\ell^2}{3}\dot{\theta}$

$\frac{\partial L}{\partial \theta} = -m_ogr\sin \theta - mg\frac{\ell}{2}\sin\theta$

$\frac{\text{d}}{\text{dt}} \Big(\frac{\partial L}{\partial \dot{\theta}}\Big) - \frac{\partial L}{\partial \theta} = 0$

Thus:

$m_o r [r\ddot{\theta} + 2\dot{r}\dot{\theta}] + \frac{m\ell^2}{3}\dot{\theta} +m_ogr\sin\theta + mg\frac{\ell}{2}\sin\theta = 0$

Also, note that:

$\frac{\partial L}{\partial \dot{r}} = m_o\dot{r}$

$\frac{\partial L}{\partial r} = -k(r-r_o) + m_og\cos\theta + m_or\dot{\theta}^2$

$\frac{\text{d}}{\text{dt}}\Big(\frac{\partial L}{\partial \dot{r}}\Big) - \frac{\partial L}{\partial r} = 0$

Thus:

$m_o\ddot{r} - m_or\dot{\theta}^2 + k(r-r_o) - m_og\cos\theta = 0$

The equilibrium positions are:

$\sin\theta = 0 \Rightarrow \theta = 0,\ \pi$

$k(r_\text{ST} - r_o) = \pm m_og$

We could now study the stability of these configurations.

Example with External Forces

$T = \frac{1}{2}m_1\dot{x}^2 + \frac{1}{2}m_2v_G^2 + \frac{1}{2}J_G\dot{\theta}^2$

$\nu_G^2 = \dot{x}^2 + \Big(\frac{\ell \dot{\theta}}{2}\Big)^2 - 2\dot{x}\Big(\frac{\ell\dot{\theta}}{2}\Big)\cos(180-\theta)$

Using $\cos(A-B) = \cos A \cos B + \sin A \sin B$ , we can say:

$\begin{split} \cos(180-\theta) &= \cos 180 \cos \theta + \sin 180 \sin \theta \\& = -\cos \theta \end{split}$

Thus:

$\nu_G^2 = \dot{x}^2 + \Big(\frac{\ell\dot{\theta}}{2}\Big)^2 + \dot{x}\ell\dot{\theta}\cos \theta$

Therefore:

$T = \frac{1}{2}m_1\dot{x}^2 + \frac{1}{2}m_2\Big[\dot{x}^2 + \Big(\frac{\ell\dot{\theta}}{2}\Big)^2 + \dot{x}\dot{\theta}\ell\cos\theta \Big] + \frac{1}{2}J_G\dot{\theta}^2$

$V = \frac{1}{2}kx^2 + m_2g\frac{\ell}{2}(1-\cos\theta)$

Where $\theta = 0$ is the datum. Now we must find the generalized forces. For a virtual displacement $\delta x$ only:

$\delta W = F(t)\delta x \Rightarrow Q_x = F(t)$

For a virtual displacement $\delta \theta$ only:

$\delta W = \ell \delta \theta \cos\theta F(t)$

Therefore:

$Q_\theta = F(t)\ell\cos\theta$

Then, note that:

$\frac{\text{d}}{\text{dt}}\Big(\frac{\partial L}{\partial \dot{x}}\Big) - \frac{\partial L}{\partial x} = Q_x \quad (1)$

$\frac{\text{d}}{\text{dt}} \Big(\frac{\partial L}{\partial \dot{\theta}}\Big) - \frac{\partial L}{\partial \theta} = Q_\theta \quad (2)$

$\frac{\partial L}{\partial \dot{x}} = m_1\dot{x} + m_2\dot{x} + m_2\frac{\ell}{2}\dot{\theta}\cos\theta$

$\frac{\partial L}{\partial x} = -kx$

$\frac{\partial L}{\partial \dot{\theta}} = m_2\frac{\ell^2}{4}\dot{\theta} + m_2\dot{x}\frac{\ell}{2}\cos\theta + J_G\dot{\theta}$

$\frac{\partial L}{\partial \theta} = -\frac{m_2}{2}\dot{x}\dot{\theta}\ell\sin\theta - m_2g\frac{\ell}{2}\sin\theta$

Therefore, (1) becomes:

$m_1\ddot{x} + m_2\ddot{x} + m_2\frac{\ell}{2}[\ddot{\theta}\cos\theta - \dot{\theta}^2\sin\theta]+kx = F(t)$

And (2) becomes:

$m_2 \frac{\ell^2}{4}\ddot{\theta} + m_2\ddot{x}\frac{\ell}{2}\cos\theta - m_2\dot{x}\frac{\ell}{2}\dot{\theta}\sin\theta + \frac{1}{12}m\ell^2\ddot{\theta} + m_2\dot{x}\frac{\ell}{2}\dot{\theta}\sin\theta + m_2g\frac{\ell}{2}\sin\theta = F(t)\ell\cos\theta$

$m_2\frac{\ell^2}{3}\ddot{\theta} + \frac{m_2}{2}\ddot{x}\ell\cos\theta + m_2g\frac{\ell}{2}\sin\theta = F(t)\ell\cos\theta$

For small angles these become:

$(m_1+m_2)\ddot{x} + m_2\frac{\ell}{2}\ddot{\theta} + kx = F(t)$

$m_2\frac{\ell^2}{3}\ddot{\theta} + m_2\frac{\ell}{2}\ddot{x} + m_2g\frac{\ell}{2}\theta = F(t)\ell$

If we have support motion instead of forcing,

then we cannot neccesairly use a potential, so now give the system a $\delta x$ and $\delta \theta$ . For a $\delta x$ :

$\delta W = -k(x-y)\delta x$

Therefore:

$Q_x = -k(x-y)$

$Q_\theta = 0$

$Q_\theta = 0$ as no work is done for $\delta \theta$ alone. Thus, $T$ is the same, and:

$V = m_2g\frac{\ell}{2}(1-\cos\theta)$

For small motions:

$(m_1+m_2)\ddot{x} +m_2\frac{\ell}{2}\ddot{\theta} + kx = ky(t)$

$m_2\frac{\ell^2}{3} \ddot{\theta} + m_2\frac{\ell}{2}\ddot{x} + m_2g\frac{\ell}{2}\theta = 0$

Thus:

$\begin{bmatrix} m_1 + m_2 & m_2\frac{\ell}{2} \\ m_2\frac{\ell}{2} & m_2\frac{\ell^2}{3} \end{bmatrix} \begin{Bmatrix} \ddot{x} \\ \ddot{\theta} \end{Bmatrix} + \begin{bmatrix} k & 0 \\ 0 & \frac{m_2g\ell}{2} \end{bmatrix} \begin{Bmatrix} x \\ \theta \end{Bmatrix} = \begin{Bmatrix} ky(t) \\ 0 \end{Bmatrix}$

Example

Consider the double pendulum:

$T = \frac{1}{2}m(\ell\dot{\theta_1})^2 + \frac{1}{2}m(\ell\dot{\theta_1} + \ell\dot{\theta_2})^2$

$V = mg\ell(1-\cos\theta_1) + mg\ell(1 - \cos\theta_1 + 1-\cos\theta_2)$

Consider only small angles. Thus:

$V \approx mg\ell\frac{\theta_1^2}{2} + mg\ell\Big(\frac{\theta_1^2}{2} + \frac{\theta_2^2}{2}\Big)$

Then:

$\frac{\partial L}{\partial \dot{\theta_1}} = m\ell^2\dot{\theta_1} + m\ell^2(\dot{\theta_1} + \dot{\theta_2})$

$\frac{\partial L}{\partial \dot{\theta_2}} = m\ell^2(\dot{\theta_1} + \dot{\theta_2})$

$\frac{\partial L}{\partial \theta_1} = -mg\ell\theta_1 - mg\ell\theta_1$

$\frac{\partial L}{\partial \theta_2} = -mg\ell\theta_2$

Therefore:

$\begin{bmatrix} 2 & 1 \\ 1 & 1 \end{bmatrix} \begin{Bmatrix} \ddot{\theta_1} \\ \ddot{\theta_2} \end{Bmatrix} + \frac{g}{\ell} \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix} \begin{Bmatrix} \theta_1 \\ \theta_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

$[m]^{-1} = \begin{bmatrix} 1 & -1 \\ -1 & 2 \end{bmatrix}$

$\begin{split} [m]^{-1}[m] &= \begin{bmatrix} 1 & -1 \\ -1 & 2 \end{bmatrix} \begin{bmatrix} 2 & 1 \\ 1 & 1 \end{bmatrix} \\&= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \end{split}$

$-p^2\begin{Bmatrix} X_1 \\ X_2 \end{Bmatrix} + \frac{g}{\ell}\begin{bmatrix} 1 & -1 \\ -1 & 2 \end{bmatrix} \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix} \begin{Bmatrix} X_1 \\ X_2 \end{Bmatrix} = 0$

$\implies -\frac{p^2\ell}{g} \begin{Bmatrix} X_1 \\ X_2 \end{Bmatrix} + \begin{bmatrix} 2 & -1 \\ -2 & 2 \end{bmatrix} \begin{Bmatrix} X_1 \\ X_2\end{Bmatrix} = 0$

Set $\lambda^2 = \frac{p^2\ell}{g}$ :

$\begin{vmatrix} 2 - \lambda^2 & -1 \\ -2 & 2 - \lambda^2 \end{vmatrix} = 0$

$(2 - \lambda ^2 ) ^2 - 2 = 0$

$4 - 4\lambda^2 + \lambda^4 - 2 = 0$

$\lambda^4 - 4\lambda^2 + 2 = 0$

$\lambda^2 = 2 \pm \sqrt{2}$

$(2-\lambda^2)X_1 - X_2 = 0$

$\frac{X_2}{X_1} = 2 - \lambda^2$

$\bigg(\frac{X_2}{X_1}\bigg)_1 = 2 - 2 + \sqrt{2} = \underline{\underline{\sqrt{2}}}$

$\bigg(\frac{X_2}{X_1}\bigg)_1 = 2 - 2 - \sqrt{2} = \underline{\underline{-\sqrt{2}}}$

Example

Now try a different set of generalized coordinates for the same problem.

If we consider small angles:

$T = \frac{1}{2} m(\ell \dot{\theta}_1)^2 + \frac{1}{2}m\big[\ell\dot{\theta}_1 + \ell(\dot{\theta}_1+ \dot{\theta}_3)\big]^2$

$V = mg\ell(1-\cos{\theta_1}) + mg\ell\big[ 1-\cos{\theta_1} + 1 - \cos{(\theta_1 + \theta_3)}\big]$

Therefore:

$V = -mg\ell\frac{\theta_1^2}{2} - mg\ell\bigg[ \frac{\theta_1^2}{2} + \frac{(\theta_2 + \theta_3)^2}{2} \bigg]$

$T = \frac{1}{2}m\ell^2(\dot{\theta}_1)^2 + \frac{1}{2}m\ell^2[2\dot{\theta}_1+\dot{\theta}_3]^2$

$V = -mg\ell\theta_1^2-mg\frac{\ell}{2}(\theta_1 +\theta_3)^2$

$L = T-V$

$\begin{split}\frac{\partial L}{\partial \dot\theta_1} &= m\ell^2\dot\theta_1 +\frac{1}{2}m\ell^22(2\dot\theta_1+\dot\theta_3)2 \\&= m\ell^2[5\dot\theta_1 + 2\dot\theta_3]\end{split}$

$\begin{split}\frac{\partial L}{\partial \dot\theta_3} &= \frac{1}{2}m\ell^2(2)[2\dot\theta_1+\dot\theta_3] \\&= m\ell^2[2\dot\theta_1+\dot\theta_3]\end{split}$

$\frac{\partial L}{\partial \theta_1 } = -2mg\ell\theta_1 - mg\ell(\theta_1+\theta_3)$

$\frac{\partial L}{\partial \theta_3} = -mg\ell(\theta_1 +\theta_3)$

$\frac{d}{dt}\bigg(\frac{\partial L}{\partial \dot\theta_1}\bigg) - \frac{\partial L}{\partial \theta_1} = 0$

$m\ell^2[5\ddot\theta_1 + 2 \ddot\theta_3] +3mg\ell\theta_1 + mg\ell\theta_3 = 0$

$m\ell^2[2\ddot\theta_1 + \ddot\theta_3] + mg\ell(\theta_1 + \theta_3) = 0$

Therefore:

$\begin{bmatrix} 5 & 2 \\ 2 & 1 \end{bmatrix}\begin{Bmatrix} \ddot\theta_1 \\ \ddot\theta_3 \end{Bmatrix} + \frac{g}{\ell}\begin{bmatrix} 3 & 1 \\ 1 & 1 \end{bmatrix}\begin{Bmatrix} \theta_1 \\ \theta_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

Therefore:

$- \frac{p^2\ell}{g}\begin{bmatrix} 5 & 2 \\ 2 & 1 \end{bmatrix}\begin{Bmatrix} X_1 \\ X_2\end{Bmatrix} + \begin{bmatrix} 3 & 1 \\ 1 & 1 \end{bmatrix} \begin{Bmatrix} X _1 \\ X_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

$[m]^{-1}[m] = \begin{bmatrix} 1 & -2 \\ -2 & 5 \end{bmatrix} \begin{bmatrix} 5 & 2 \\ 2 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$

Set $\lambda^2 = \frac{p^2\ell}{g}$ :

$-\lambda^2\begin{Bmatrix} X_1 \\ X_3 \end{Bmatrix} + \begin{bmatrix} 1 & -2 \\ -2 & 5 \end{bmatrix} \begin{bmatrix} 3 & 1 \\ 1 & 1 \end{bmatrix} \begin{Bmatrix} X_1 \\ X_3 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

$-\lambda^2\begin{Bmatrix} X_1 \\ X_3 \end{Bmatrix} + \begin{bmatrix} 1 & -1 \\ -1 & 3 \end{bmatrix} \begin{Bmatrix} X-1 \\ X_3 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

Therefore:

$(1 - \lambda^2)(3 - \lambda^2) - 1 = 0$

$\lambda^4 - 4\lambda^2 + 2 = 0 \quad \text{(same characteristic equation)}$

Therefore:

$\lambda^2 = 2 \pm \sqrt{2}$

For mode shapes:

$(1 - \lambda^2) X_1 - X_3 = 0$

$\implies \frac{X_3}{X_1} = 1 - \lambda^2$

$\begin{split} \bigg(\frac{X_3}{X_1}\bigg)_1 &= 1 - 2 + \sqrt{2} \\&= -1 + \sqrt{2} \end{split}$

$\begin{split} \bigg(\frac{X_3}{X_1} \bigg) _2 &= 1 - 2 - \sqrt{2} \\&= -1 - \sqrt{2} \end{split}$

These are different values of ratios than in the first case. However, when plotted:

Now use the newtonain approach to the same problem:

$+\circlearrowleft \sum M_0 = m\ell\ddot\theta_1 = - mg\ell\theta_1 + T( \theta_2 - \theta_1) \ell \quad \text{(A)}$

$\begin{split} + \downarrow \sum F_y &= m\ddot y \\&= -T\cos\theta_2 + mg \end{split}$

$\begin{split} + \rightarrow \sum F_x &= m\ddot x \\&= - T\sin{\theta_2} \end{split}$

$\ddot{x} \approx \ell ( \ddot\theta_1 + \ddot\theta_2) , \ \ddot y \approx 0$

Therefore:

$T \approx mg$

$m\ell(\ddot\theta_1 + \ddot\theta_2) + mg\theta_2 = 0 \quad \text{(B)}$

$\begin{bmatrix} m\ell & 0 \\ m\ell & m\ell \end{bmatrix} \begin{Bmatrix} \ddot\theta_1 \\ \ddot\theta_2 \end{Bmatrix} + \begin{bmatrix} 2mg & -mg \\ 0 & mg \end{bmatrix} \begin{Bmatrix} \theta_1 \\ \theta_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

Note both mass and stiffness matrices are non-symmetric!

Calculating the eigenvalues gives the same result and the mode shapes are identical with the other 2 formulations.

NOTE: These sets (first and third formulation) are linearly dependent.

Lagrange’s Equation for Small Oscillations

Our main interest is in the motion of MDOF systems in the neighborhood of equilibrium positions. WLOG we assume the equilibrium position is $q_1 = q_2 \dotsc = q_n = 0$ and that the displacements are sufficiently small that the linear force- displacement and force-velocity relations hold. This means that the generalized coordinates and their time derivatives appear in the differential equations to only the first power.

Therefore:

$\underline{\dot r} = \sum_{k = 1}^n \frac{\partial \underline{r}}{\partial q_k}\dot{q}_k$

And $\frac{\partial \underline r}{\partial q_k}$ are constants.

$\begin{split} T &= \frac{1}{2} \sum_{i = 1}^N m_i \underline{\dot r}_i \underline{\dot r}_i \\&= \frac{1}{2} \sum_{i = 1}^N m_i\bigg(\sum_{r = 1}^n \frac{\partial \underline{r}_i}{\partial q_r}\dot{q}_r\bigg)\bigg(\sum_{s = 1}^n \frac{\partial \underline{r}_i}{\partial q_s}\dot{q}_s\bigg) \\&= \frac{1}{2} \sum_{r = 1}^n \sum_{s = 1}^n\bigg(\sum_{i=1}^N m_i \frac{\partial \underline{r}_i}{\partial q_r}\frac{\partial \underline{r}_i}{\partial q_s}\bigg)\dot q_r\dot q_s \end{split}$

Set $m_{rs} = \sum_{i = 1}^N m_i \frac{\partial \underline{r}_i}{\partial q_r}\frac{\partial \underline{r}_i}{\partial q_s} = m_{sr}$ as the generalized masses.

Similarly we can write the potential energy as $V = V(q_1, q_2 \dotsc q_n)$ and consider a Taylor expansion of V about the equilibrium configuration.

$V(q_n) = V_0 \sum_{n = 1} ^n \frac{\partial V}{\partial q_i} + \frac{1}{2}\sum_{i= 1}^n \sum_{j = 1}^n \frac{\partial^2V}{\partial q_i \partial q_j}q_iq_j+ \dotsc$

Where the partial derivatives are all evaluated at $q_i = 0$ .

As the equilibrium point is a condition in which:

$\frac{\partial V}{\partial q_i} = 0$

$V$ is approximately:

$\begin{split} V &= \frac{1}{2}\sum_{r = 1}^n \sum_{s = 1}^n \frac{\partial^2 V}{\partial q_r \partial q_s}q_rq_s \\&=\frac{1}{2}\sum_{r = 1}^n \sum_{s = 1}^n k_{rs}q_rq_s\end{split}$

Where $k_{rs} = k_{sr}$ are the stiffness coefficients. If we now use Lagrange’s equation the equations of motion become:

$\frac{\partial L}{\partial \dot{q}_j} = \sum_{s = 1}^n m_{js}\dot{q}_s$

$\frac{\partial L}{\partial q_j} = \sum_{s = 1}^n k_{js}q_s$

Therefore

$\sum_{s = 1}^n[m_{js} \ddot{q}_s(t) + k_{js}q(t)] = Q_j$

$j = 1,2 \dotsc n$

While this is very straightforward, we can use this result to obtain approximate formulations of the equations of motion as continuous systems it replaces them with a $n$ degree of freedom system. It can also extend the technique to find approximate response to a specific forcing function. This is called the assumed modes method (AMM).

For discrete system it replaces them with a smaller number of DOF, while for continuous system it replaces them with an $n$ degree of freedom system. We can also easily extend it to a forced system.

Then the method assumes:

$y(x,t) = \sum_{i = 1}^N\phi_iq_i(t)$

Where $\phi_i(x)$ are trial functions and $q_i(t)$ are the generalized coordinates.

$T = \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^Nm_{ij}\dot{q}_i\dot{q}_j$

$V = \frac{1}{2}\sum_{i = 1}^N\sum_{j = 1}^Nk_{ij}q_iq_j$

$\phi_i(x)$ needs only satisfy the geometric boundary conditions.

For continuous systems:

$\begin{split}T &= \frac{1}{2} \int dm\bigg(\frac{\partial y}{\partial t}\bigg)^2 \\&= \frac{1}{2} \sum_{i = 1}^N\sum_{j = 1}^N\dot{q}_i\dot{q}_j\int_m \phi_i(x)\phi_j(x)dm\end{split}$

Therefore:

$m_{ij} = \int_m \phi_i(x)\phi_j(x)dm$

While for discrete systems

$m_{ij} = \sum_{p = 1}^N m_p\phi_i(x_p)\phi_j(x_p)$

The potential energy can be determined is a similar manner.

Consider a generalized beam:

$V = \frac{1}{2}\int_m EI \frac{\partial^2 y}{\partial x^2}dx$

$\frac{\partial^2 y}{\partial x^2} = \sum_{i=1}^N \phi''(x)(t)$

Therefore:

$V = \frac{1}{2}\sum_{i = 1}^N\sum_{j=1}^Nq_iq_j\int_LEI\phi''_i\phi''_jdx$

$k_{ij} = \int EI\phi''_i\phi''_jdx$

This will allow calculation of a consistent mass and stiffness matrices.

Example: Discrete System

$m_1 = m_2 = m_3 = m$

$[k] = \begin{bmatrix} k & \frac{k}{3} & 0 \\ \frac{k}{3}& \frac{11k}{9} & -k \\ 0 & -k & k\end{bmatrix}$

We want to estimate the lowest natural frequency using assumed modes method.

Try:

$\phi_1 = \begin{bmatrix} 1 \\ -1 \\ -\frac{3}{2} \end{bmatrix}$

Now calculate the consistent mass and stiffness using this assumed shape:

$m_{11} = m\begin{bmatrix} 1 & -1 & -\frac{3}{2} \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix} 1 \\ -1 \\ -\frac{3}{2} \end{bmatrix} = \frac{17}{4}m$

$k_{11} = k\begin{bmatrix} 1 & -1 & -\frac{3}{2} \end{bmatrix}\begin{bmatrix} 1 & \frac{1}{3} & 0 \\ \frac{1}{3} & \frac{11}{9} \ -1 \\ 0 & -1 & 1\end{bmatrix}\begin{bmatrix} 1 \\ -1 \\ -\frac{3}{2} \end{bmatrix} = \frac{29}{36}k$

$m_{11}\ddot{q}_1 + k_{11}q_1 = 0$

$p = \sqrt{\frac{29}{36}\frac{4}{17}\frac{k}{m}} = 0.435\sqrt{\frac{k}{m}}$

What if we choose $\phi_1 = \begin{bmatrix} 1 \\ -1 \\ -\frac{3}{2} \end{bmatrix}, \ \phi_2 = \begin{bmatrix} 1 \\ -2 \\-3 \end{bmatrix}$ .

Then $m_{11}$ and $k_{11}$ are the same.

$m_{22} = m(1)(1) + m(-2)(-2) +m (-3)(-3) = 14m$

$k_{22} = k\begin{bmatrix} 1 & -2 & -3 \end{bmatrix}\begin{bmatrix} 1 & \frac{1}{3} & 0 \\ \frac{1}{3} & \frac{11}{9} & -1 \\ 0 & -1 & 1\end{bmatrix}\begin{bmatrix} 1 \\ -2 \\ -3 \end{bmatrix} = \frac{14}{9}k$

$k_{12} = k_{21} = k\begin{bmatrix} 1 & -1 & -\frac{3}{2} \end{bmatrix}\begin{bmatrix} 1 & \frac{1}{3} & 0 \\ \frac{1}{3} & \frac{11}{9} & -1 \\ 0 & -1 & 1\end{bmatrix}\begin{bmatrix} 1 \\ -2 \\ -3 \end{bmatrix} = \frac{17}{18}k$

$m_{12} = m\begin{bmatrix} 1 & -2 & -3 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix} 1 \\ -1 \\ -\frac{3}{2} \end{bmatrix} = \frac{15}{2}m$

So we have replace the system by:

$m\begin{bmatrix} \frac{17}{4} & \frac{15}{2} \\ \frac{15}{2} & 14 \end{bmatrix} \begin{bmatrix} \ddot{q}_1 \\ \ddot{q}_2\end{bmatrix} + \begin{bmatrix} \frac{29}{36} & \frac{17}{18} \\ \frac{17}{18} & \frac{14}{9}\end{bmatrix} \begin{bmatrix} q_1 \\ q_2 \end{bmatrix}$

Now calculate the $p$ ‘s.

$\begin{bmatrix} \frac{29}{36}k - \frac{17}{4}mp^2 & \frac{17}{18} k - \frac{15}{2}mp^2 \\ \frac{17}{18}k - \frac{15}{2}mp^2 & \frac{14}{9} k - 14mp^2 \end{bmatrix} \begin{Bmatrix} Q_1 \\ Q_2 \end{Bmatrix} = 0$

$\bigg( \frac{29}{36}k - \frac{17}{4}mp^2\bigg)\bigg(\frac{14}{9}k - 14 mp^2 \bigg) - \bigg(\frac{17}{18}k - \frac{15}{2}mp^2\bigg)^2 = 0$

$\frac{29\cdot 14}{36\cdot 9}k^2 - kmp^2\bigg(\frac{17\cdot 14}{36} + \frac{14\cdot 29}{36}\bigg) + m^2p^4\bigg(\frac{17\cdot 14}{4}\bigg)$

$-\frac{17^2}{18^2}k^2 + 2\bigg(\frac{17\cdot 15}{18\cdot 2}\bigg)kmp^2 - \frac{225}{4} m^2p^4 = 0$

$p^4\bigg[ \frac{238}{4} - \frac{225}{4} \bigg] - \frac{k}{m} p^2 \bigg[ \frac{17\cdot 14}{36} + \frac{14\cdot 29}{36} - \frac{17\cdot 15}{18}\bigg] + \bigg[\frac{29\cdot 14}{36\cdot 9} - \frac{17^2}{18^2}\bigg]\bigg(\frac{k}{m}\bigg)^2 = 0$

$\frac{13}{4}p^4 - p^2\frac{k}{m}\bigg[ \frac{17\cdot 14 + 14\cdot 29 - 2\cdot 17\cdot 15}{36}\bigg] + \left[\frac{406-289}{324}\right]\frac{k^2}{m^2} = 0$

$\frac{13}{4}p^4 - \frac{134}{36}p^2\frac{k}{m} + \frac{117}{324} = 0$

$\frac{13}{4}p^2 - \frac{67}{18}p^2\frac{k}{m} + \frac{39}{108}\frac{k^2}{m^2} = 0$

$p^2 = \frac{\frac{67}{18} \pm \sqrt{(\frac{67}{18})^2 - 4\frac{13}{4}\cdot \frac{39}{108}}}{\frac{13}{2}}\frac{k}{m} = \frac{3.722 \pm \sqrt{3.722^2 - 4.694}}{6.5}\frac{k}{m}$

$p^2 = \frac{3.722 \pm 3.024}{6.5}\frac{k}{m}$

$p_1^2 \leq 0.1074 \frac{k}{m}$

$\implies p_1 \leq 0.3277\sqrt{\frac{k}{m}}$

The answer is dependent on the assumption of mode shape

EXAMPLE: Use static deflection as mode.

Add masses to the diagram

$F = 5mg$

$\Delta_1 = \frac{mg}k-5\frac{mg}k = -4\frac{mg}k$

$\Delta_2 = 15\frac{mg}k$

$\Delta_3 = 16\frac{mg}k$

Therefore choose $\phi_1 = \begin{bmatrix}4 \\ -15 \\ -16 \end{bmatrix}$

$m_{11} = m(4)(4) + m(-15)(-15) + m(-16)(-16) = 497m$

$\begin{split} k_{11} &= k\begin{bmatrix} 4 && -15 && -16 \end{bmatrix} \begin{bmatrix} 1 && \frac{1}{3} && 0 \\ \frac{1}{3} && \frac{11}{9} && -1 \\ 0 && -1 && 1 \end{bmatrix} \begin{bmatrix} 4 \\ 15 \\ 16 \end{bmatrix} \\ &= 27k \end{split}$

Therefore:

$\begin{split} p_1 &= \sqrt{\frac{27}{497}}\sqrt{\frac{k}{m}} \\ &\leq 0.23308\sqrt{\frac{k}{m}} \end{split}$

It is sometimes useful to consider these approximate solutions to finding the eigenvalues and eigenvectors. The UAMM gives upper bounds as do most techniques. We can get lower bounds as well.

Consider again the equations and motion.

$\begin{bmatrix} m \end{bmatrix} \begin{Bmatrix} \Ddot{q} \end{Bmatrix} + \begin{bmatrix} k \end{bmatrix} \begin{Bmatrix} q \end{Bmatrix} = \begin{Bmatrix} 0 \end{Bmatrix}$

To get to an eigenvalue problem we can do it in two ways.

Premultiply by $\begin{bmatrix} k \end{bmatrix}^{-1} = \begin{bmatrix} a \end{bmatrix}$
$\begin{split} \begin{Bmatrix} \begin{bmatrix} a \end{bmatrix} \begin{bmatrix} m \end{bmatrix} && - \frac{1}{p^2} \begin{bmatrix} I \end{bmatrix} \end{Bmatrix} \begin{Bmatrix} Q \end{Bmatrix} &= 0 \\ \begin{bmatrix} a \end{bmatrix} \begin{bmatrix} m \end{bmatrix} &:= \begin{bmatrix} D \end{bmatrix} -\text{dynamical matrix} \end{split}$
Premultiply by $\begin{bmatrix} m \end{bmatrix}^{-1}$
$\begin{split} \begin{Bmatrix} \begin{bmatrix} m \end{bmatrix}^{-1} \begin{bmatrix} k \end{bmatrix} && - p^2 \begin{bmatrix} I \end{bmatrix} \end{Bmatrix} \begin{Bmatrix} Q \end{Bmatrix} &= 0 \\ \begin{bmatrix} m \end{bmatrix}^{-1} \begin{bmatrix} k \end{bmatrix} &:= \begin{bmatrix} \Delta \end{bmatrix} -\text{system matrix} \end{split}$

either gives the same eigenvalues and equivalent eigenvectors.

One can use these to find the eigenvalues/vectors through matrix iteration.

One clever idea that gives a lower bound is Dunkerleys formula. Consider the Dynamic Matrix formulation

$\begin{Bmatrix} \begin{bmatrix} a \end{bmatrix} \begin{bmatrix} m \end{bmatrix} - \frac{1}{p^2}\begin{bmatrix} I \end{bmatrix} \end{Bmatrix} \begin{Bmatrix} Q \end{Bmatrix} = 0$

and expand the determinant of the coefficients

$\begin{vmatrix} D_{11}-\frac{1}{p^2} && D_{12} && ... && D_{1n} \\ D_{21} && (D_{22}-\frac{1}{p^2} && ... && D_{2n} \\ \begin{matrix} . \\ . \\ . \end{matrix} && && && \\ D_{n1} && && && (D_{nn}-\frac{1}{p^2} \end{vmatrix} = 0$

$\biggr( \frac{1}{p^2} \biggr)^n - \big( D_{11} + D_{22} + D_{33} + ... + D_{nn} \big) \biggr( \frac{1}{p^2} \biggr)^{n-1} + ... = 0$

Now consider that we have the roots of this polynomial $p_1^2, p_2^2, ... p_n^2$ . This equation is equivalent to (factored)

$\biggr( \frac{1}{p^2}-\frac{1}{p_1^2} \biggr) \biggr( \frac{1}{p^2}-\frac{1}{p_2^2} \biggr) \biggr( \frac{1}{p^2}-\frac{1}{p_3^2} \biggr) \biggr( \qquad \biggr) ... = 0$

and expand this to get

$\biggr( \frac{1}{p^2} \biggr)^n - \biggr( \frac{1}{p_1^2} + \frac{1}{p_2^2} + \frac{1}{p_3^2} + ... + \frac{1}{p_n^2} \biggr) \biggr( \frac{1}{p^2} \biggr)^{n-1} + ... = 0$

Now equate the coefficients

$tr D = \biggr( \frac{1}{p_1^2} + \frac{1}{p_2^2} + ... + \frac{1}{p_n^2} \biggr)$

However the largest term is the first one and often dominates

Therefore:

$\begin{split} \frac{1}{p_1^2} &< tr D \\ p_1^2 &> \frac{1}{tr D} \end{split}$

A special case occurs when the mass matrix is diagonal. Then tr D is

$tr D = a_{11}m_1 + a_{22}m_2 + ... + a_{nn}m_{nn}$

and we need calculate only $a_{ii}, i=1$ to $n$

In this case there is a physical interpretation of the expressions

$\frac{1}{p_1^2} \leq a_{11}m_1 + a_{22}m_2 + a_{33}m_3... + a_{nn}m_{nn}$

namely that $a_{11}m_1$ is the natural frequency of the system with only $m_1$ as the mass (all other masses set to zero).

Therefore:

$\frac{1}{p_1^2} \leq \frac{1}{p^2_{11}} + \frac{1}{p^2_{22}} + ... + \frac{1}{p^2_{nn}}$

This is also sometimes shown as Dunkerley’s formula.

Consider previous example

Calculate $a_{ii}$

Apply unit load to $m_1$
$\begin{split} \Delta_1 &= \frac{1}{k} + \frac{1}{k} = \frac{2}{k} = a_{11} \\ \Delta_2 &= -\frac{3}{k} = \Delta_3 \end{split}$

Apply unit load to $m_2$

$\begin{split} F(L) &= 1(3L) \\ F &= 3 \\ \hat{\Delta} &= \frac{3}{k} \\ a_{22} &= \frac{g}{k} \end{split}$

Apply unit load to $m_3$

$\begin{split} \hat{\Delta} &= \frac{g}{k} \\ a_{33} &= \frac{9}{k} + \frac{1}{k} = \frac{10}{k} \end{split}$

Therefore:

$\begin{split} \frac{1}{p_1^2} &< m\left(\frac{2}{k}\right) + m\left(\frac{9}{k}\right) + m\left(\frac{10}{k}\right) = \frac{21m}{k} \\ p_1 &> \sqrt{\frac{1}{21}}\sqrt{\frac{k}{m}} = \underline{0.2182\sqrt{\frac{k}{m}}} \end{split}$

Consider an example of a continuous system

$\begin{split} \phi_1 &= x^2 \\ \phi^{''}_1 &= 2 \\ \phi_2 &= x^3 \\ \phi^{''}_2 &= 6x \end{split}$

$\begin{split} m_{ij} &= \int_m \phi_i\phi_j dm \\ &= \int_0^L \phi_i\phi_J mdx \qquad \text{where } m= \text{ mass/length} \end{split}$

$\begin{split} m_{11} &= m\int_0^L x^4 dx = \frac{mL^5}{5} \\ m_{12} &= m\int_0^L x^5dx = \frac{mL^6}{6} = m_{21} \\ m_{22} &= \frac{mL^7}{7} \end{split}$

$\begin{split} k_{ij} &= \int EI\phi_i''\phi_j'' dx \\ k_{11} &= EI \int_0^L 4 dx\\ &= EI(4L) \\ k_{12} &= EI\int_0^L 12x dx \\ &= 6EIL^2 \\ k_{22} &= 12L^3EI \end{split}$

$mL^5 \begin{bmatrix} 1/5 && L/6 \\ L/6 && L^2/7 \end{bmatrix} \begin{Bmatrix} \Ddot{q}_1 \\ \Ddot{q}_2 \end{Bmatrix} + EIL \begin{bmatrix} 4 && 6L \\ 6L && 12L^2 \end{bmatrix} \begin{Bmatrix} q_1 \\ q_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

Therefore the eigenvalue problem is

$\begin{bmatrix} 4 - \frac{\lambda}{5} && 6L - \frac{\lambda L}{6} \\ 6L - \frac{\lambda L}{6} && 12L^2 - L^2 \frac{\lambda}{7} \end{bmatrix} \begin{Bmatrix} Q_1 \\ Q_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}$

$\begin{split} (4-\frac{\lambda}{5})L^2(12-\frac{\lambda}{7}) - L^2(6-\frac{\lambda}{6})^2 &= 0 \\ 48 - \frac{12}{5}\lambda - \frac{4}{7}\lambda + \frac{\lambda^2}{35} - (36-2\lambda + \frac{\lambda^2}{36}) &= 0 \\ 12 - \frac{(84+20-70)}{35}\lambda + \lambda^2(\frac{1}{35}-\frac{1}{36}) &= 0 \\ 12(35)(36) - (34)(36)\lambda + \lambda^2 &= 0 \end{split}$

$\lambda^2 - (34)(36)\lambda - 12(35)(36)$

$\begin{split} \lambda &= (34)(18) \pm \sqrt{(34)^2(18)^2 - 12(35)(36)} \\ &= 612 \pm \sqrt{359424} \\ &= 612 \pm 599.5 \\ &= 12.48, 1211.5 \end{split}$

Therefore

$\begin{split} p_1^2 &= \sqrt{12.48} \sqrt{\frac{EI}{mL^4}} \\ &= 3.53\sqrt{\frac{EI}{mL^4}} \\ p_2^2 &= 34.8 \sqrt{\frac{EI}{mL^4}} \end{split}$

$(4-\frac{\lambda}{5})q_1 + (6-\frac{\lambda}{6})Lq_2 = 0$

$\begin{split} \frac{Lq_2}{q_1} &= \frac{-(4-\frac{\lambda}{5})}{(6-\frac{\lambda}{6})} \\ &= \frac{-(4-\frac{12.48}{5})}{6-\frac{12.48}{6}} \\ &= -0.384 \end{split}$

The mode shape is

$\phi = x^2 - 0.384x^3$

What if we apply some dynamic loading to the structure? Then we must calculate the generalized forces for our model. Consider a virtual displacement

$\delta_y = \sum_{i=1}^n \phi_i(x)\delta q_i$

then

$\begin{split} \delta W &= \sum_{j=1}^m F_j \biggr( \sum_{i=1}^n \phi_i (x_j) \delta_{q_i} \biggr) \\ &= \sum_{i=1}^n \delta_{q_i} \biggr( \sum_{j=1}^m F_j \phi_i (x_j) \biggr) \\ &= \sum_{i=1}^n Q_i \delta q_i \end{split}$

Therefore

$Q_i = \sum_{j=1}^m F_j \phi_i (x_j)$

Consider the previous example

$\begin{split} Q_1 &= F(t) \phi_1 (L) = F(t) L^2 \\ Q_2 &= F(t) \phi_2 (L) = F(t) L^3\end{split}$

$\begin{bmatrix} m \end{bmatrix} \begin{Bmatrix} \Ddot{q} \end{Bmatrix} + \begin{bmatrix} k \end{bmatrix} \begin{Bmatrix} q \end{Bmatrix} = F(t) \begin{Bmatrix} L^2 \\ L^3 \end{Bmatrix}$

Whirling of Rotating Shafts

Rotating shafts tend to bow out at certain speeds and whirl in a complicated manner. Whirling is defined as the rotation of the bent shaft and the line of centers of the bearing. This phenomenon results from various causes as mass unbalance, hysteresis damping, gyroscopic effects, fluid friction, etc. This whirling can take place in the same or opposite direction as that of the rotation of the shaft. The whirling speed may or may not be equal to the rotation speed.

Let us assume a simple system

Note: $s$ is the center of the disc, $G$ is the center of gravity, and $e$ is an eccentricity.

$k$ is the stiffness of the system to lateral deflection.

We will assume the shaft to be rotating at an instant angular velocity $\omega$ and OS to be whirling at angular speed $\dot{\theta}$

$\underline{v}_G = \underline{v}_s + v_{G/s}$

Where $\phi = \omega t - \theta$

$\begin{split} v_r^G &= \dot{r} - e\omega\sin(\omega t - \theta) \\ v_\theta^G &= r\dot{\theta} + \omega e\cos(\omega t - \theta) \end{split}$

Therefore:

$\begin{split} \Big( v^G \big)^2 &= \Big[ \dot{r} - e\omega \sin(\omega t - \theta)\Big]^2 + \Big[ r\dot{\theta} + e\omega \cos(\omega t - \theta)\Big]^2 \\ &= (\dot{r})^2 - 2\dot{r}e\omega\sin(\omega t - \theta) + (r\dot{\theta})^2 + 2r\dot{\theta}e\omega\cos(\omega t - \theta) + e^2\omega^2 \\ (v_G)^2 &= (\dot{r})^2 + (r\dot{\theta})^2 + e^2\omega^2 - 2\dot{r}e\omega\sin(\omega t -\theta) + 2r\dot{\theta}e\omega\cos(\omega t -\theta) \end{split}$

Therefore:

$\begin{split} T &= \frac{1}{2} J_G \omega^2 + \frac{1}{2}mv^2_G \\ V &= \frac{1}{2}kr^2 \\ L &= T - V \end{split}$

Therefore:

$\begin{split} \frac{d}{dt}\Big(\frac{\partial L}{\partial\dot{r}}\Big) - \frac{\partial L}{\partial r} &= Q_r \\ \frac{d}{dt}\Big(\frac{\partial L}{\partial\dot{\theta}}\Big) - \frac{\partial L}{\partial\theta} &= Q_\theta \end{split}$

$\begin{split} \frac{\partial (V_G)^2}{\partial \dot{\theta}} &= 2r^2\dot{\theta} + 2re\omega\cos(\omega t - \theta) \\ \frac{\partial (V_G)^2}{\partial \dot{r}} &= 2\dot{r} - 2e\omega\sin(\omega t - \theta) \\ \frac{\partial (V_G)^2}{\partial \theta} &= 2\dot{r}e\omega\cos(\omega t -\theta) + 2r\theta e\omega\sin(\omega t - \theta) \\ \frac{\partial (V_G)^2}{\partial r} &= 2r(\dot{\theta})^2 + 2\dot{\theta}e\omega\cos(\omega t - \theta) \end{split}$

Therefore:

$\begin{split} \frac{md}{2dt}\big[2\dot{r}-2e\omega\sin(\omega t-\theta)\big] - \frac{m}{2}\big[2r(\dot{\theta}^2)+2\dot{\theta} e\omega\cos(\omega t -\theta) + kr &= -c\dot{r} \\ m\Ddot{r} - me\omega\cos(\omega t - \theta)(\omega - \dot{\theta}) - mr(\dot{\theta})^2 - m\dot{\theta}e\omega\cos(\omega t - \theta) + kr &= 0 \\ m\Ddot{r} - mr(\dot{\theta})^2 - me\omega^2\cos(\omega t -\theta) + kr &= -c\dot{r} \end{split}$

Therefore:

$\Ddot{r} + \frac{k}{m}r + \frac{c}{m}\dot{r}-r(\dot{\theta})^2 = e\omega^2\cos(\omega t - \theta)$

$\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{\theta}}\Big) - \frac{\partial L}{\partial \theta}$

$\begin{split} &\frac{m}{2}\frac{d}{dt}\big(2r^2\dot{\theta} + 2re\omega\cos(\omega t -\theta)\big) - \frac{m}{2}\big[2\dot{r}e\omega\cos(\omega t - \theta) + 2r\dot{\theta}e\omega\sin(\omega t -\omega)\big] \\ &= m\big[ 2r\dot{r}\dot{\theta} + r^2\Ddot{\theta} + \dot{r}e\omega\cos(\omega t -\theta) - re\omega(\omega - \dot{\theta})\sin(\omega t - \theta)\big] \\ &- m\big[\dot{r}e\omega\cos(\omega t - \theta) + r\dot{\theta}e\omega\sin(\omega t - \theta)\big] \end{split}$

Therefore:

$mr\Ddot{\theta} + 2m\dot{r}\dot{\theta} - me\omega^2\sin(\omega t - \theta) = -cr\Ddot{\theta}$

Therefore:

$r\Ddot{\theta} + \frac{c}{m}r\dot{\theta} + 2\dot{r}\dot{\theta} = e\omega^2\sin(\omega t - \theta)$

These are nonlinear equations and there are more than one solution.

Consider the simplest case of synchronous whirl i.e. try

$\begin{split} \dot{\theta} &= \omega \qquad \theta = \omega t - \phi \\ \Ddot{\theta} &= 0 \qquad \Ddot{r} = \dot{r} = 0 \end{split}$

Therefore:

$\begin{split} \frac{k}{m}r - r\omega^2 &= e\omega^2\cos\phi \\ \frac{c}{m}r\omega &= e\omega^2\sin\phi \\ r\Big(\frac{k}{m} - \omega^2\Big) &= e\omega^2\cos\phi \end{split}$

$\begin{split} \tan\phi &= \frac{\frac{c}{m}r\omega}{r\big(\frac{k}{m}-\omega^2\big)} \\ p &= \frac{k}{m} \\ \zeta &= 2mp \end{split}$

$\begin{split} r &= \frac{e\omega^2}{\big(\frac{k}{m}-\omega^2\big)}\cos\phi \\ &= \frac{e\omega^2}{\big(\frac{k}{m}-\omega^2\big)} \frac{\big(\frac{k}{m}-\omega^2\big)}{\sqrt{\big(\frac{k}{m}-\omega^2\big)^2+\big(\frac{c}{m}\omega\big)^2}} \end{split}$

$r = \frac{e\big(\frac{\omega}{p}\big)^2}{\sqrt{\big(1-\big(\frac{\omega}{p}\big)^2\big)^2 + \big( 2\zeta\frac{\omega}{p}\big)^2 }}$

This looks exactly the same as rotating inbalance excitation

Therefore the

So that synchronous whirl occurs when $\frac{\omega}{p} \approx 1$

Self Excited Vibrations

Almost all of the vibrations considered so far were either free or forced vibrations. A fundamentally different class are the so-called self-excited vibrations.

“In a self-excited vibration the alternating force that sustains the motion is created or controlled by the motion itself ; when the motion stops the alternating force disappears”

In a forced vibration the sustaining alternating force exists independently of the motion and persists even when the vibratory motion is stopped.

Engineered art

Use vibration phenomena for the display of classic ideas that ask the question “Why does that occur”
Search kinetic art
Naum Gabo Stationary Wave (1919)

Spinning Pendulum

Consider the stability of the SHM of the simple pendulum that is rotating with angular velocity $\omega$ . The total velocity of mass $m$ is:

$\nu^2 = (\ell\dot\theta)^2+(\ell\theta\omega)^2$

$T = \frac{1}{2}m\ell^2[\dot\theta^2+\theta^2\omega^2]$

$\begin{split}V &= mg\ell[1-\cos\theta ] \\&=\frac{1}{2}mg\ell\theta^2\end{split}$

$\frac{d}{dt}\bigg(\frac{\partial L}{\partial\dot\theta}\bigg)-\frac{\partial L}{\partial \theta} = 0$

$m\ell^2\ddot\theta- m\ell^2\theta\omega^2+mg\ell\theta=0$

$\ddot\theta+\bigg(\frac{g}{\ell}-\omega^2\bigg)\theta = 0$

$p^2 = \frac{g}{\ell} - \omega^2$

The frequency of oscillation of the pendulum is:

$f_p = \frac{1}{2\pi}\sqrt{\frac{q}{\ell}-\omega^2}$

Therefore the pendulum will oscillate ( $f>0$ ) until $\omega^2 = \frac{g}{\ell}$ at which point the pendulum will simply fly horizontally and never cross the $\nsubseteq$ .

To see the frequencies try a length of pendulum of $\ell = 8$ inches.

When $\omega = 0$ :

$\begin{split} f_p &=\frac{1}{2\pi}\sqrt{\frac{386}{8}} \\&=\underline{\underline{1.1}}\text{ Hz (66 RPM)}\end{split}$

When $\omega = 30$ RPM:

$\begin{split} f_p &=\frac{1}{2\pi}\sqrt{\frac{386}{8}- \bigg(\frac{2\pi \cdot30}{60}\bigg)^2} \\& = \underline{\underline{0.99}}\text{ Hz}\end{split}$

When $\omega = 40$ RPM:

$\begin{split} f_p &=\frac{1}{2\pi}\sqrt{\frac{386}{8}- \bigg(\frac{2\pi \cdot40}{60}\bigg)^2} \\& = \underline{\underline{0.88}}\text{ Hz}\end{split}$

When $\omega = 50$ RPM:

$\begin{split} f_p &=\frac{1}{2\pi}\sqrt{\frac{386}{8}- \bigg(\frac{2\pi \cdot50}{60}\bigg)^2} \\& = \underline{\underline{0.73}}\text{ Hz}\end{split}$

When $\omega = 60$ RPM:

$\begin{split} f_p &=\frac{1}{2\pi}\sqrt{\frac{386}{8}- \bigg(\frac{2\pi \cdot60}{60}\bigg)^2} \\& = \underline{\underline{0.47}}\text{ Hz}\end{split}$

Friction Testing (Kinetic)

$\begin{split} +\uparrow\sum F_y &= R_1 + R_2 - mg = 0 \qquad \textcircled{1} \\ +\rightarrow\sum F_x &= \mu(R_1-R_2) = m\Ddot{x} \qquad \textcircled{2} \\ +\circlearrowleft\sum M_G &= 0 = \mu(R_1-R_2)h + R_2(a-x) - R_1(a+x) \qquad \textcircled{3} \end{split}$

Find the natural frequency in terms of $a, h, \mu$

From $\textcircled{1}$

$R_1 + R_2 = mg \qquad \textcircled{4}$

From $\textcircled{3}$

$(R_1 - R_2)\mu h + (R_2 - R_1)a - (R_1 + R_2)x = 0$

Using $\textcircled{4}$

$(R_1 - R_2)\mu h + (R_2 - R_1)a - mgx = 0$

Therefore:

$\begin{split} (R_1 - R_2)(\mu h - a) &= mgx \\ (R_1 - R_2) &= \frac{mgx}{(\mu h - a)} \end{split}$

From $\textcircled{2}$

$\frac{\mu mgx}{(\mu h - a)} = m\Ddot{x}$

$\Ddot{x} + \frac{\mu g}{a - \mu h}x = 0$

Therefore:

$p^2 = \sqrt{\frac{\mu g}{a - \mu h}}$

Therefore, frequency of oscillation is

$f = \frac{1}{2\pi}\sqrt{\frac{\mu g}{a - \mu h}}$

“Friction keeps the world in harmonious motion”

$f = \frac{1}{2\pi}\sqrt{\frac{\mu g}{a - \mu h}}$

$a = 4" \qquad \mu = 0.5 \qquad h = 1"$

$f = \frac{1}{2\pi} \sqrt{\frac{(0.5)386}{4 - (0.5)}} = \underline{1.2 \text{ Hz}}$

$a = 4 \qquad \mu = 0.1 \qquad h = 1"$

$f = \frac{1}{2\pi} \sqrt{\frac{(0.1)386}{4 - (0.1)}} = \underline{0.5 \text{ Hz}}$

$f = \frac{1}{2\pi}\sqrt{\frac{\mu g}{\alpha - \mu h}}$

Let $\alpha = 3"$ , $h = 0.5$ and $\mu = 0.1$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.1)(386)}{3 - (0.1)(0.5)}} \\&= 0.576 \text{ Hz} \end{split}$

If $\mu = 0.2$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.2)(386)}{3 - (0.2)(0.5)}} \\&= 0.821 \text{ Hz} \end{split}$

If $\mu = 0.04$ (TEFLON/STEEL):

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.04)(386)}{3 - (0.04)(0.5)}} \\&= 0.361 \text{ Hz} \end{split}$

If $\mu = 0.3$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.3)(386)}{3 - (0.3)(0.5)}} \\&= 1.01 \text{ Hz} \end{split}$

If $\mu = 0.4$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.4)(386)}{3 - (0.4)(0.5)}} \\&= 1.18 \text{ Hz} \end{split}$

If $\mu = 0.5$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.5)(386)}{3 - (0.5)(0.5)}} \\&= 1.33 \text{ Hz} \end{split}$

If $\mu = 0.6$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.6)(386)}{3 - (0.6)(0.5)}} \\&= 1.47 \text{ Hz} \end{split}$

If $\mu = 0.7$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.7)(386)}{3 - (0.7)(0.5)}} \\&= 1.61 \text{ Hz} \end{split}$

If $\mu = 0.8$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.8)(386)}{3 - (0.8)(0.5)}} \\&= 1.73 \text{ Hz} \end{split}$

If $\mu = 0.9$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(0.9)(386)}{3 - (0.9)(0.5)}} \\&= 1.86 \text{ Hz} \end{split}$

If $\mu = 1.0$ :

$\begin{split} f &= \frac{1}{2\pi}\sqrt{\frac{(1.0)(386)}{3 - (1.0)(0.5)}} \\&= 1.98 \text{ Hz} \end{split}$

$V = \frac{1}{2}k(L \sin \theta)^2 - mgL(1 - \cos \theta)$

$T = \frac{1}{2}m(\dot{\theta L})^2$

$L = T-V$

$\frac{d}{dt}(\frac{\partial L}{\partial q_i}) - \frac{\partial L}{\partial q_i} = Q_i$

$mL^2 \ddot{\theta} + kL^2 \sin \theta \cos \theta - mgL \sin \theta + cL^2 \dot{\theta} = 0$

The equilibrium configurations are:

$\theta _o = 0, \pi$

and

$\sin \theta _o [kL^2 \cos \theta _o - mgL] = 0$

therefore,

$\cos \theta _o = \frac{mg}{kL}$

$V = \frac{1}{2}k(L \sin \theta)^2 - mgL(1 - cos \theta)$

$\begin{split} \frac{\partial V}{\partial \theta} &= kL^2 \sin \theta \cos \theta - mgL \sin \theta \\&= kL^2[\sin \theta \cos \theta - \frac{mg}{kL} \sin \theta] \\&= kL^2 \sin \theta \big( \cos \theta - \frac{mg}{kL}\big) \end{split}$

$\frac{\partial^2V}{\partial \theta^2} = k L^2 \cos \theta \left(\cos \theta - \frac{mg}{kL}\right) - kL^2\sin \theta$

At $\theta = \theta _ o = 0$

$= kL^2 [1(1-\frac{mg}{kL})]$

$\frac{\partial ^2 V}{\partial \theta} >0$ if $\frac{mg}{kL} < 1$ , $\frac{k}{m} > \frac{g}{L}$

at $\theta = \pi$

$kL^2\left[\left(-1\right)\left(-1-\frac{mg}{kL}\right)\right] > 0$

Therefore always stable

at $\theta _o$ when $\cos \theta _ o = \frac{mg}{kL}$

$\frac{\partial ^2 V}{\partial \theta ^2} = -kL^2 \sin \theta _o$

therefore, always unstable.

This is easy to do for a SDOF (MDOF), however for general systems it is not as straightforward. For the general case, there are mathematical techniques that are used and referred to by the work of Routh, E.J. done in ‘Rigid Dynamics’ 1897. We will consider it for the case of the ‘1 1/2” and two DOF systems only.

7.2. Mathematical Criterion of Stability

(J.P. Den Hartog, Mechanical Vibrations, 4th Ed. 1956)

For single-degree-of-freedom systems, the criterion of dynamic stability can be derived by physical rather than by mathematical means. With systems of two or three degrees of freedom, a physical conception is always very helpful but usually does not give a complete interpretation of what happens. A mathematical approach is necessary, and this involves at first the setting up of the differential equations of the problem. As long as we deal with small vibrations (and thus disregard any non-linearities that may exist), the equations are all linear and of the second order. Their solution, as usual, is found by assuming the following:

$x_1 = X_1e^{st}$

$x_2 = X_2e^{st}$

$\vdots$

$x_n = X_ne^{st}$

where $s$ is a complex number the real part of which determines the damping and the imaginary part of which is the natural frequency. Substituting into the differential equations of the free vibration transforms these equations into a set of n homogeneous, linear algebraic equations in the (complex) unknowns $X_1 \ldots X_n$ . A process of algebraic elimination is then performed with the result that one equation is obtained which does not contain any of these variables. This equation, known as the “frequency equation,” is generally of the degree 2n in s. Thus, for a two-degree-of-freedom system we obtain a quartic; for a three-degree-of freedom system we obtain a sixth-degree equation, etc.

An algebraic equation of degree 2n in the variable s has 2n roots or 2n values of s. Real roots of s would lead to terms $e^{st}$ in the solution, which rarely occur in ordinary vibrating systems. The roots of s are usually complex and then they always occur in conjugate pairs:

$s_1 = p_1 + jq_1$

$s_2 = p_1 - jq_1$

$s_3 = p_2 + jq_2$

$s_4 = p_2 - jq_2$

The solution of the first differential equation is

$x_1 = C_1 e^{s_1t} + C_2 e^{s_2t} + C_3 e^{s_3t} + ...$

We know that these terms can be combined by pairs as follows:

$C_1e^{s_1t} + C_2e^{s_2t} = e^{p_1t}(A \sin q_1t + B \cos q_1t)$

so that the imaginary part of s is the frequency, and the real part of $s$ determines the rate of damping. If the real parts of all the values of s are negative, the system is dynamically stable; but if the real part of any one of the values of s is positive, the system is dynamically unstable.

Therefore the stability can be determined by an examination of the signs of the real parts of the solutions of the frequency equation. It is not necessary to solve the equation, because certain rules exist by which from an inspection of the coefficients of the equation a conclusion regarding the stability or instability can be drawn. These rules, which were given by Routh in 1877, are rather complicated for frequency equations of higher degree, but for the most practical cases (third and fourth degree) they are sufficiently simple.

Let us consider first the cubic equation

$s^3 + A_2s^2 + A_1s + A_0 = 0 \quad (\star)$

(7.4) which occurs in the case of two degrees of freedom where one mass or spring is zero (in a sense one and one-half degrees of freedom)If its roots are $s_1,s_2$ and $s_3$ , $(\star)$ can be written

$(s-s_1)\cdot (s-s_2) \cdot (s-s_3) = 0$

or, worked out,

$s^3 - (s_1 + s_2 + s_3)s^2 + (s_1s_2 + s_1s_3+s_2s_3)s - s_1s_2s_3 = 0$

A comparison shows the following:

$A_2 = -(s_1 + s_2 + s_3)$

$A_1 = s_1s_2 + s_1s_3 + s_2s_3 \qquad (\star \star)$

$A_0 = -s_1s_2s_3$

One of the three roots of a cubic equation must always be real, and the other two are either real or conjugate complex. Separating the roots $s_1, s_2, s_3$ into their real and imaginary parts, we may write

$s_1 = p_1$

$s_2 = p_2 + jq_2$

$s_3 = p_2 - jq_2$

Substituted into $(\star \star)$ this leads to the following:

$A_2 = -(p_1 + 2p_2)$

$A_1 = 2p_1p_2 + p_2^2 + q_2^2$

$A_0 = -p_1(p_2^2 + q_2^2)$

The criterion of stability is that both $p_1$ and $p_2$ be negative. It is seen in the first place that all coefficients $A_2, A_1$ , and $A_0$ must be positive, because, if any one of them were negative, the above requires that either $p_1$ or $p_2$ , or both $p_1$ and $p_2$ , must be positive. This requirement can be proved to hold for higher degree equations as well. Hence a frequency equation of any degree with one or more negative coefficients determines an unstable motion.

Granted that the coefficients, $A_0,A_1$ and $A_2$ are all positive, the third equation requires that $p_1$ be negative. No information about $p_2$ is available as yet. However, on the boundary between stability and instability, $p_2$ must pass from a positive to a negative value through zero. Make $p_2$ = 0 and we can derive the following:

$A_2 = -p_1$

$A_1 = q_2^2$

$A_0 = -p_1q_2^2$

These relations must be satisfied on the boundary of stability. By eliminating $p_1$ and $p_2$ we find

$A_0 = A_1A_2$

We do not know yet on which side of this relation stability exists. That can be found in the simplest manner by trying out one particular case. For example let $s_1 = -1$ and $s_{2,3} = -1 \pm j$ , which obviously is a stable solution. Substitution gives

$A_2 = 3, A_1 = 4, A_0 = 2$

so that

$A_o < A_1A_2$

The complete criterion for stability of the cubic is that all coefficients A are positive and that

$A_1A_2 > A_0$

Next consider the quartic

$s^4 + A_3s^3 + A_2s^2 + A_1s + A_0 = 0$

for which the procedure is similar. Since a quartic can be resolved into two quadratic factors, we may write for the roots:

$s_1 = p_1 + jq_1$

$s_2 = p_1-jq_1$

$s_3 = p_2 + jq_2$

$s_4 - p_2 - jq_2$

and substitite in, which leads to:

$A_3 = -2(p_1 + p_2)$

$A_2 = p_1^2 + p_2^2 + q_1^2 + q_2^2 + 4p_1p_2$

$A_1 = -2p_1 \cdot (p_2^2 + q_2^2) - 2p_2(p_1^2 + q_1^2)$

$A_0 = (p_1^2 + q_1^2) \cdot (p_2^2 + q_2^2)$

The requirement for stability is that both $p_1$ and $p_2$ be negative. Substitution of negative values of $p_1$ and $p_2$ makes all four A‘s positive, so that the first requirement for stability is that all coefficients A be positive. Granted that this is so, the first equation requires that at least one of the quantities $p_1$ or $p_2$ be negative. Let $p_1$ be negative. We still need another requirement to make $p_2$ also negative. On the boundary between stability and instability, $p_2$ = 0, which substituted in gives:

$A_3 = -2p_1$

$A_2 = p_1^2 + q_1^2 + q_2^2$

$A_1 = -2p_1q_2^2$

$A_0 = (p_1^2 + q_1^2)q_2^2$

being four equations in the three variables $p_1,q_1,q_2$

Elimination of these variables leads to a relations between the $A$ ‘s:

$A_1A_2A_3 = A_1^2 + A_3^2A_0$

To find out on which side of this equality stability exists, we try out a simple stable case, for example:

$s_\text{1,2} = -1 \pm j, \quad s_\text{3,4} = -2 \pm 2j$

Which, on substitution, gives:

$A_3 = 6,\ A_2 = 18,\ A_1 = 24,\ A_0 = 16$

So that $A_1A_2A_3 > A_1^2 + A_2^3A_0$ . Thus, the complete criterion for stability of the quartic (7.10) is that all coefficients $A$ are positive and that:

$A_1A_2A_3 > A_1^2 + A_2^3A_0$

Systems with three degrees of freedom generally have a sextic for their frequency equation and in degenerated cases, a quintic. In such cases, there are three real parts of the roots $s$ , and besides the requirements of positive signs for all coefficients A, there are two other requirements, each of which is rather lengthy. For further information, the reader is referred to the original work of Routh.

Example

Investigate the stability of motion of a pendulum of length $\ell$ and mass $m$ rotating with constant angular velocity $\omega$ about its vertical axis.

Assume during the steady state rotation about the vertical axis as small angle $\theta$ is given to the pendulum and $\dot{\phi}$ is the corresponding variation in the angular velocity. Use the Lagrangian formulation to determine the equations of motion.

$v_c^2 = \ell^2 \theta^2 (\omega + \dot{\phi})^2 + \ell_2 (\dot{\theta})^2$

Therefore:

$T = \frac{1}{2} m \ell^2 (\omega^2 \theta^2 + (\dot{\theta})^2) \quad \text{(neglecting } (\dot{\phi})^2 \text{ as small)}$

$V = \frac{1}{2}mg\theta^2 \quad \text{(} 1 - \cos \theta \approx \theta^2 \text{ for } \theta << 1 \text{)}$

$\Big[ \frac{d}{dt} \Big(\frac{\partial L}{\partial \dot{q}_i}\Big) - \frac{\partial L}{\partial q_i} = 0 \Big]$

Thus:

$m \ell^2 \ddot{\theta} - m \ell^2 \omega^2 \theta + mg \theta = 0$

$\ddot{\theta} + \Big(\frac{g}{\ell} - \omega^2\Big) \theta = 0$

If there was viscous damping $c\dot{\theta}$ , the equation would be:

$\ddot{\theta} + 2c\dot{\theta} + \Big(\frac{g}{\ell} - \omega^2\Big) \theta = 0$

Therefore, assuming $\theta = ce^{st}$ with $p^2 = \Big(\frac{g}{\ell} - \omega^2\Big)$ :

$s^2 + 2cs + p^2 = 0$

$s = -c \pm \sqrt{c^2 - p^2}$

So that if $p^2$ is negative, one of the roots is positive and the motion is unstable.

Example

Look at the stability of the steady rotation of a horizontal flywheel carrying two masses $m$ attached to hub by springs and free to slide along the spokes. If $J$ is the moment of inertia of the flywheel, $\omega$ is the angular velocity, and $a$ is the upstretched lengths of the springs. Let $x$ be the small change in $r$ and $\dot{\phi}$ the small change in $\omega$ .