Equation Of Motion

The equation of motion is the central differential equation to solve problems in dynamics. It is an ordinary differential equation of second order and in the general, linear, damped case it looks like the following:

(1)	$\begin{array}{l}\displaystyle m \ddot{u}(t) + c \dot{u}(t) + k u(t) = f,\end{array}$

where $\begin{array}{l}m\end{array}$ denotes the mass, $\begin{array}{l}c\end{array}$ the damping and $\begin{array}{l}k\end{array}$ the stiffness. These terms are explained more in detail in the next section. Note, that for systems with multiple degrees of freedom these quantites are matrices. $\begin{array}{l}u\end{array}$ is the unknown displacement and occurs without derivative as well as once derived (velocity) and twice derived (acceleration) with respect to time.

Via Fourier Transform it is possible, to transform the equation of motion from the time domain to the frequency domain, which looks like follows:

(2)	$\begin{array}{l}\displaystyle -\omega^2 m u(\omega) + i \omega c u(\omega) + k u(\omega) = f(\omega)\end{array}$

First, the different terms in the equation of motion are explained more in detail and afterwars the derivation of the equation of motion is presented for various problems.

Terms in the equation of motion

For a linear viscoelastic material that is subject to time-dependent variations of stress and strain, the most general relation between stress and strain is given by the following partial differential equation

(3)

$\begin{array}{l}\displaystyle \left( a_0 + a_1 \frac{\text{d}}{\text{dt}} + a_2 \frac{\text{d}^2}{\text{dt}^2} + \ldots + a_n \frac{\text{d}^n}{\text{dt}^n} + \ldots \right) f = \left( b_0 + b_1 \frac{\text{d}}{\text{dt}}+ b_2 \frac{\text{d}^2}{\text{dt}^2} + \ldots + b_n \frac{\text{d}^n}{\text{dt}^n} + \ldots \right) u . \label{eq:damping_generalized_standard_model}\end{array}$

This relation is also known as the generalized standard model, which relates forces and displacements or, in an alternative setup, stresses and strains. It holds for arbitrary orders and describes every system consisting of springs and viscous dampers.

Stiffness term

Hooke model. This model corresponds to the purely elastic Hooke's law with $\begin{array}{l}\frac{b_0}{a_0} = k\end{array}$ and all $\begin{array}{l}a_i = b_i = 0\end{array}$ for $\begin{array}{l}i \geq 1\end{array}$ .

$\begin{array}{l}\displaystyle a_0 f = b_0 u \qquad \to \quad f = k u \, .\end{array}$

In this case all time derivatives vanish.

Damping term

Newton model. This model corresponds to the viscous dashpot with $\begin{array}{l}a_0\end{array}$ , $\begin{array}{l}b_1 \neq 0\end{array}$ and all other $\begin{array}{l}a_i, b_i\end{array}$ equal to zero.

$\begin{array}{l}\displaystyle a_0 f = b_1 \frac{\text{du}}{\text{dt}} \qquad \to \quad f = c \dot{u} \, .\end{array}$

For systems with multiple degrees of freedom the damping coefficients build in the general case an arbitrary damping matrix. With respect to Modal Analysis it makes sense to make simplifying assumptions which lead to a proportional damping matrix (Rayleigh damping).

Mass term

The force acting on a mass that experiences a motion is given by the principle of d'Alembert. This force occurs due to inertia and always acts in the opposite direction than the acceleration. It is derived from Newton's 2. law:

$\begin{array}{l}\displaystyle f=m\ddot{u}\end{array}$

Equation of motion for a single degree of freedom system

The single degree of freedom (SDOF) system plays a significant role in many problems in structural dynamics. This importance stems from the fact, that many engineering systems or system components can be approximated by SDOF systems. The mechanical behaviour is described by the equation of motion, which is an ordinary differential equation.

In the following the prinicple of virtual work is presented to obtain the equation of motion for a SDOF system.

The following system is excited by a dynamic force $\begin{array}{l}f(t)\end{array}$ . Gravity forces are not considered.

The system is kinematic and has one degree of freedom. The system in the deformed state is depicted below.

The kinematic relations are:

$\begin{align} w_G &= \varphi_1 \frac{l}{2} = \varphi_2 l \quad \rightarrow \quad \frac{\varphi_1}{2} = \varphi_2 \\ w_m &= \varphi_1 l \\ w_c &= \varphi_1 \frac{l}{2} = \frac{w_m}{2} \\ w_F &= w_G = \varphi_1 \frac{l}{2} = \frac{w_m}{2} \\ w_k &= \varphi_2 \frac{l}{2} = \frac{\varphi_1}{2} \frac{l}{2} = \frac{1}{4} \varphi_1 l = \frac{w_m}{4} \, . \end{align}$

Hence the following forces act on the system:

Next, a virtual displacement is superimposed upon the deformed system.

For the system to be in equilibrium, the sum of virtual work of all forces has to vanish, i.e. $\begin{array}{l}\delta W = 0\end{array}$ . Thus,

$\begin{align} - m \ddot{w}_m \delta w_m - c \dot{w}_c \delta w_c - k w_k \delta w_k + f (t) \delta w_F &= 0 \\ m \ddot{w}_m \delta w_m + c \frac{\dot{w}_m}{2} \frac{\delta w_m}{2} + k \frac{w_m}{4} \frac{\delta w_m}{4} &= f (t) \frac{\delta w_m}{2} \, . \\ % \end{align}$

With $\begin{array}{l}\delta w_m \ne 0\end{array}$ it follows

$\begin{array}{l}\displaystyle m \ddot{w}_m + \frac{1}{4} c \dot{w}_m + \frac{1}{16} k w_m = \frac{1}{2} f(t)\end{array}$

This equation corresponds to the differential equation of the system with respect to the vertical displacement of the point mass $\begin{array}{l}m\end{array}$ . It is identical to the general equation of motion (1) with the generalized quantities $\begin{array}{l}\bar{m} = m\end{array}$ , $\begin{array}{l}\bar{c} = \frac{1}{4}c\end{array}$ , $\begin{array}{l}\bar{k} = \frac{1}{16}k\end{array}$ , $\begin{array}{l}\bar{f}=\frac{1}{2}f\end{array}$ .

Helmholtz equation

In the following, a numerical derivation of the equation of motion for wave propagation in fluids is presented using the Finite Element Method. This paragraph is taken from:

Franck, Andreas. "Finite-Elemente-Methoden, Lösungsalgorithmen und Werkzeuge für die akustische Simulationstechnik." (2009).

For linear analysis, the wave equation for an ideal fluid is given by

$\begin{array}{l}\displaystyle \Delta p = \frac{1}{c^2}\frac{\partial^2p}{\partial t^2} \qquad \text{with } c^2 = \frac{1}{\kappa \rho} \qquad \text{and } \Delta = \text{div grad},\end{array}$

where $\begin{array}{l}p\end{array}$ is the sound pressure level, $\begin{array}{l}c\end{array}$ the speed of sound, $\begin{array}{l}\rho\end{array}$ the density of the fluid and $\begin{array}{l}\kappa\end{array}$ the adiabatic compressibility of the fluid. Note, that the sound pressure level $\begin{array}{l}p=p(\omega)\end{array}$ has a complex amplitude at a circular frequency $\begin{array}{l}\omega = 2\pi f\end{array}$ . Hence, the function of $\begin{array}{l}p\end{array}$ over time and its derivatives can be written as

$\begin{array}{l}p(t)=p e^{i\omega t} \\ \frac{\partial p}{\partial t} = i \omega p e ^{i\omega t} \\ \frac{\partial^2 p}{\partial t^2} = -\omega^2 p e ^{i\omega t}.\end{array}$

With this, the wave equation can be written in the frequency domain:

(4)	$\begin{array}{l}\displaystyle \Delta p + k^2p = 0 \qquad \text{with } k = \frac{\omega}{c}.\end{array}$

The basis of the Finite Element Method is always a variational formulation of the underlying physical problem. This can be done using the method of weighted residuals. Therefore, the Helmholtz equation (4) is multiplied with a weighting function $\begin{array}{l}\bar{w}\end{array}$ and integrated over the whole domain $\begin{array}{l}\Omega\end{array}$ .

$\begin{array}{l}\displaystyle \int_\Omega \bar{w} (\Delta p + k^2p) \,d\Omega = 0.\end{array}$

The weak form is obtained by partial integration

$\begin{array}{l}\displaystyle \int_\Omega (\nabla \bar{w} \nabla p - k^2\bar{w} p) \,d\Omega - \int_\Gamma \bar{w} \frac{\partial p}{\partial \textbf{n}} \,d\Gamma = 0.\end{array}$

The second term is the integral over the boundary $\begin{array}{l}\Gamma = \partial \Omega\end{array}$ , which can be split up into two parts regarding the following boundary conditions:

Neumann boundary condition: The velocity normal to the boundary $\begin{array}{l}v_n\end{array}$ is set. This leads to

$\begin{array}{l}\displaystyle \nabla p \cdot \textbf{n} = \frac{\partial p }{\partial \textbf{n}} = - i \omega \rho_F v_n = \omega^2 \rho_F u_n \quad \forall x \in \Gamma_v.\end{array}$

$\begin{array}{l}u_n\end{array}$ denotes the normal velocity of the boundary.

Addmittance boundary condtion: the sound pressure level is set in a relation to its normal derivative (which is the normal velocity). On these boundaries $\begin{array}{l}\Gamma_A\end{array}$ the following holds:

$\begin{array}{l}\displaystyle \nabla p \textbf{n} = \frac{\partial p }{\partial \textbf{n}} = - i \omega \rho_F (A_n \cdot p) \quad \forall x \in \Gamma_A.\end{array}$

The addmitance $\begin{array}{l}A_n\end{array}$ is a material property and can be obtained from measurements.

With these boundary conditions the variational formulation can be written as

(5)	$\begin{array}{l}\displaystyle \int_\Omega (\nabla \bar{w} \nabla p - k^2\bar{w} p) \,d\Omega - \int_{\Gamma_v} \bar{w}(\omega^2 \rho_F u_n) \,d \Gamma_v + \int_{\Gamma_A} \bar{w} (i \omega \rho_F A_n p) \,d \Gamma_A = 0 .\end{array}$

Discretization

In order to apply the Finite Element Method, the domain has to be approximated using elements. The solution for the physical problem is then given by the discrete solution at the nodes of those elements and interpolated via ansatz functions for the rest of the domain. Hence, every element has as many ansatz functions as nodes. These are usually simple polynomials, which are expressed in local coordinates $\begin{array}{l}\xi, \eta, \zeta\end{array}$ . The distribution of a physical quantity $\begin{array}{l}f\end{array}$ is then given by a weighted sum over all nodal values $\begin{array}{l}f_i\end{array}$ mulitplied with their respective ansatz function $\begin{array}{l}N_i\end{array}$ :

$\begin{array}{l}\displaystyle f(\xi, \eta, \zeta) = \sum N_i (\xi, \eta, \zeta) \cdot f_i.\end{array}$

Due to the interpolating characteristic of $\begin{array}{l}N_i\end{array}$ , the following must hold:

$\begin{array}{l}N_i (\xi_j, \eta_j, \zeta_j) = \begin{cases} 1 \quad i = j \\ 0 \quad i \neq j \end{cases}.\end{array}$

This means, that at each node only the respective ansatz function is $\begin{array}{l}1\end{array}$ on the rest are $\begin{array}{l}0\end{array}$ . Hence, the nodal values are interpolated at each node.

Finally, the weighting function $\begin{array}{l}\bar{w}\end{array}$ has to be chosen. One possibility is to use $\begin{array}{l}\bar{w}_i = N_i\end{array}$ which is the so-called Galerkin-approach. Hence, the sound pressure level and the weighting functions can be written in a discretized form:

$\begin{array}{l}\displaystyle p = \sum_{i=0}^n N_i p_i \qquad \bar{w} = \sum_{j=0}^n N_j \bar{w}_j,\end{array}$

where $\begin{array}{l}n\end{array}$ denotes the number of nodes per element. Inserting this into (5) leads to

(6)	$\begin{array}{l}\displaystyle \bar{w}_i \int_\Omega \nabla N_i \nabla N_j \,d\Omega p_j - k^2 \bar{w}_i \int_\Omega N_i N_j \,d\Omega p_j - \bar{w}_i \omega^2 \rho_F \int_{\Gamma_v} N_i u_n \,d\Gamma_v + \bar{w}_i i \omega \rho_F \int_{\Gamma_A} N_iN_j A_{n,j} \,d\Gamma_A p_j = 0.\end{array}$

According to the fundamental lemma of variational calculus, equation (6) must be fulfilled for an arbitrary choice of $\begin{array}{l}\bar{w}_i\end{array}$ . Hence, $\begin{array}{l}\bar{w}_i\end{array}$ can be replaced with $\begin{array}{l}1\end{array}$ and one obtains

(7)	$\begin{array}{l}\displaystyle \int_\Omega \nabla N_i \nabla N_j \,d\Omega p_j - \frac{\omega^2}{c_F^2} \int_\Omega N_i N_j \,d\Omega p_j + i \omega \rho_F \int_{\Gamma_A} N_iN_j A_{n,j} \,d\Gamma_A p_j = \omega^2 \rho_F \int_{\Gamma_v} N_i u_n \,d\Gamma_v\end{array}$

$\begin{array}{l}\displaystyle \Leftrightarrow (\textbf{K}_F - \omega^2 \textbf{M}_F + i\omega \textbf{A}_F) \textbf{p} = \textbf{f}_F.\end{array}$

This equation is identical to the equation of motion in the frequency domain (2). The matrices $\begin{array}{l}\textbf{K}_F\end{array}$ , $\begin{array}{l}\textbf{M}_F\end{array}$ and $\begin{array}{l}\textbf{A}_F\end{array}$ are computed in the following way:

$\begin{align} \text{Compressibility matrix: } \qquad K_{F,ij} &= \int_\Omega \nabla N_i \nabla N_j \,d\Omega \\ \text{Mass matrix: } \qquad M_{F,ij} &= \frac{1}{c_F^2} \int_\Omega N_i N_j \,d\Omega \\ \text{Admittance matrix: } \qquad A_{F,ij} &= \rho_F \int_{\Gamma_A} A_n N_i N_j \,d\Gamma_A \\ \text{force vector: } \qquad f_{F,i} &= \omega^2 \rho_F \int_{\Gamma_v} N_i u_n \,d\Gamma_v \\ &= \omega^2 \rho_F \int_{\Gamma_v} N_i N_j \,d\Gamma_v u{n,j} \end{align}$

Euler-Bernoulli beam

We derive the equation of motion of the Euler-Bernoulli beam via the equilibrium of forces on a differential element.

We assume the material properties to be independent of time. The dynamic equilibrium for the vibration of a beam with the bending stiffness $\begin{array}{l}EI\end{array}$ and the mass per unit length $\begin{array}{l}\mu\end{array}$ , and under the assumption of Euler-Bernoulli theory (i.e., neglecting rotational inertia and shear deformation), reads

$\begin{array}{l}\displaystyle -Q' (x,t) + c (x) \dot{w} (x,t) + \mu (x) \ddot{w} (x,t) - p (x,t) = 0 \, .\end{array}$

We assume the material and cross-sectional properties to be time-independent. The shear force is $\begin{array}{l}Q (x,t) = \left( - EI (x) w'' (x,t) \right)'\end{array}$ . Thus, we obtain

$\begin{array}{l}\displaystyle \left(EI (x) w'' (x,t) \right)' + c (x) \dot{w} (x,t) + \mu (x) \ddot{w} (x,t) - p (x,t) = 0 \, .\end{array}$

Neglecting the external loading as well as damping forces, and assuming a constant bending stiffness $\begin{array}{l}EI\end{array}$ and mass density $\begin{array}{l}\mu\end{array}$ , the partial differential equation for the free vibration is obtained.

(8)	$\begin{array}{l}\displaystyle EI w'''' (x,t) + \mu \ddot{w} (x,t) = 0\end{array}$

Solution via a separation approach

The partial differential equation (8) of the free vibration can be solved by a separation approach with space function $\begin{array}{l}W ( x )\end{array}$ and time function $\begin{array}{l}f(t)\end{array}$

(9)	$\begin{array}{l}\displaystyle w (x,t) = W (x) f(t) \, .\end{array}$

We insert (9) into (8) and obtain

(10)	$\begin{array}{l}\displaystyle EI W'''' (x) f(t) + \mu W (x) \ddot{f} (t) = 0 \, .\end{array}$

Dividing (10) by $\begin{array}{l}\mu\end{array}$ and $\begin{array}{l}W(x) f(t)\end{array}$ , it follows

(11)	$\begin{array}{l}\displaystyle \frac{\ddot{f}(t)}{f(t)} = - \frac{EI}{\mu} \frac{W'''' (x)}{W (x)} \, .\end{array}$

(11) must hold at every instance in time and space, i.e., both sides have to be equal to a constant, also called separation constant, say $\begin{array}{l}\lambda = - \omega^2\end{array}$ . This at first not obvious choice will become apparent when discussing the solution characteristics.

(12)	$\begin{array}{l}\displaystyle \frac{\ddot{f}(t)}{f(t)} = - \frac{EI}{\mu} \frac{W'''' (x)}{W (x)} = - \omega^2 \, .\end{array}$

Thus, we can write the following for both terms in (12) separately to obtain the two decoupled ordinary homogeneous differential equations

(13)	$\begin{gather} \ddot{f}(t) + \omega^2 f(t) = 0 \\ W'''' (x) - k^4 W (x) = 0 \end{gather}$

with

$\begin{array}{l}\displaystyle k^4 = \dfrac{\mu}{EI} \omega^2 \, .\end{array}$

We find the following solutions under the assumption $\begin{array}{l}\omega^2 >0\end{array}$ .

(14)	$\begin{align} W (x) & = \widehat{K}_{x1} e^{kx} + \widehat{K}_{x2} e^{-kx} + \widehat{K}_{x3} e^{i k x} + \widehat{K}_{x4} e^{- i k x} \label{eq:Biegewelle1} \\ % f (t) & = \widehat{K}_{t1} e^{i \omega t} + \widehat{K}_{t2} e^{-i \omega t} \, . \end{align}$

The solution in time now corresponds to oscillating waves as the exponent of the exponential function is complex. Now, reconsider our choice for $\begin{array}{l}\lambda\end{array}$ in (12). If the constant $\begin{array}{l}\lambda\end{array}$ was assumed to be $\begin{array}{l}\omega^2\end{array}$ for $\begin{array}{l}\omega^2>0\end{array}$ we would have obtained solutions that correspond to decaying and growing exponential functions, which are not admissible in our problems.

Eqs. (14) can alternatively be written as

(15)	$\begin{align} W (x) & = C_{x1} \cosh (k x) + C_{x2} \sinh (k x) + C_{x3} \cos (k x) - C_{x4} \sin (k x) \\ f (t) & = C_{t1} \cos (\omega t) - C_{t2} \sin (\omega t) \, . \end{align}$

Thus, the solution contains six unknown coefficients in Eqs. (14), or Eqs. (15), respectively. These can be found by inserting the spatial boundary conditions and the temporal initial conditions. Possible spatial boundary conditions are depicted in the following figure.

Through the relations $\begin{array}{l}w(x,t) = W (x) f(t)\end{array}$ , $\begin{array}{l}w' (x,t) = W' (x) f(t)\end{array}$ , $\begin{array}{l}M(x,t) = - EI W'' (x) f(t)\end{array}$ , and $\begin{array}{l}Q(x,t) = - EI W''' (x) f(t)\end{array}$ all boundary conditions can be represented by derivatives of the displacement. At each boundary of a beam two boundary conditions can be found, i.e., at each boundary two equations for the unknown coefficients are given. This results in a $\begin{array}{l}4\times 4\end{array}$ system of equations for the unknown coefficients $\begin{array}{l}C_{xi}\end{array}$ . A solution is found by excluding the trivial solution and setting the determinant of the equation system to zero. As the solution for the time function $\begin{array}{l}f(t)\end{array}$ is decoupled, the unknowns $\begin{array}{l}C_{ti}\end{array}$ are found through the initial displacement and velocities for $\begin{array}{l}t=0\end{array}$ .

References

Snowdon, John C. Vibration and shock in damped mechanical systems. J. Wiley, 1968.

Nashif, Ahid D., David IG Jones, and John P. Henderson. Vibration damping. John Wiley & Sons, 1985.

Kausel, Eduardo. Advanced structural dynamics. Cambridge University Press, 2017.