# How to simulate the Natural Frequencies of a Planetary Gearbox?

I have a dynamical model of a planetary gearbox of the form:

$\mathbf{M{\ddot q}+{\Omega_c}G{\dot q}+{[K_b+K_m-\Omega_c^2K_\Omega]q}=T(t)+F(t)}$

$K_b$: bearing stiffness symmetrical matrix $K_m$: mess stiffness matrix $K_\Omega$: diagonal centripetal stiffness matrix $T$: applied external torque $F$: static transmission error excitation $q$: the vector with the traslation and rotational coordinates for the sun,ring and carrier

For the free vibrations I have the equation:

$\mathbf{M{\ddot q}+{[K_b+K_m]q}=0}$

$\mathbf{-\omega_i^2M{\phi_i}+{[K_b+K_m]\phi_i}=0}$

Where $\phi_i$ is the vibration modes.

$\mathbf{\omega_i^2M{\phi_i}={[K_b+K_m]\phi_i}}$

I know that the $\Omega_c$ term for the gyroscopic effect is neglected at first and there's no excitation input because it's free vibration. But then I don't understand where the $\phi_i$ comes from.

I have basic knowledge of vibration theory for a single degree system, but not for multi DOF which is this case.

I need to determine the natural frequencies and vibrational modes of the system from this equation but I don't know how or where to start. I have initial numerical values for Stiffness, mass, preassure angle of the carrier, sun and planet gears.

Numerically how would I know how many modes of vibration I have and how to determine them with their natural frequencies? Like I said this is the first time I deal with a MDOF system and I've read a lot but this model is far to complicated and the math is confusing me.

• It would be ideal if you could you define each of the variables in the matrix equation. For example, is $\mathbf G$ a damping matrix (symmetric), a rotation matrix (antisymmetric), or a combination of damping and rotation (asymmetric)? I'll have a go at answering, but it certainly helps if you define the symbols. Mar 9 '17 at 12:46
• @SprocketsAreNotGears I learned mathjax and included all the other equations I have, I should have stated them from the beggining, the main problem I have now is that I don't know how to substitute all the values or represent the final form of the equation, from there I suppose I could just solve for the natural frequencies and vibration modes? Mar 9 '17 at 17:41

To obtain the modes shapes and resonant frequencies, you start from your equation of motion with no externally applied forces, which is indeed as you've stated.

$$\mathbf M \mathbf{\ddot q} + \mathbf K \mathbf q = \mathbf 0 \qquad (1)$$

For brevity, I've let $\mathbf K = \mathbf K_b + \mathbf K_m$.

Currently, $\mathbf q(t)$ is a function of time. If the system is vibrating in the $i$th mode, everywhere will be vibrating at frequency $\omega_i$ in the same phase, assuming there is no damping. Therefore, we can assume the following: $$\mathbf q(t) = \mathbf u^{(i)}\exp{(j\omega_i t)}\qquad(2)$$ , where $j$ is the imaginary number defined as $\sqrt{-1}$ and $\mathbf u^{(i)}$ is the $i$th mode shape, which is independent of time. (This assumption may seem bizarre as $\exp{(j\omega_i t)}$ is a complex number, but this does not matter in this scenario, and is used just to show that $\mathbf q(t)$ is sinusoidal in time with frequency $\omega_i$. Equivalently, you could instead use $\mathbf q(t) = \mathbf u^{(i)} (A \sin{\omega_i t} + B \cos{\omega_i t})$.)

Note that:

$$\mathbf {\ddot q}(t) = \frac{d}{dt}\left(\mathbf u^{(i)}\exp{(j\omega_i t)}\right)=-\omega_i^2\mathbf u^{(i)}\exp{(j\omega_i t)}\qquad (3)$$

Substituting equations (2) and (3) into (1) will give you the following, noting that the exponential functions cancel out:

$$\left(-\omega_i^2 \mathbf M + \mathbf K \right) \mathbf u^{(i)}=0 \qquad (4)$$

Equation 4 is a type of problem called the generalised eigenvalue problem, which involves solving equations of the form $\mathbf A \mathbf x = \lambda \mathbf B \mathbf x$ (compare with the standard eigenvalue problem, which involves solving $\mathbf A \mathbf x = \lambda \mathbf x$). In this case, the eigenvalue is $\omega_i^2$ and the eigenvector is $\mathbf u^{(i)}$. To obtain the eigenvalues and eigenvectors involves solving the generalised eigenvalue problem. There are two main approaches: hand calculation or numerical computation. For your problem, the latter will definitely be the recommended option, however, the former should help with smaller problems and lends some insight to how generalised eigenvalue problems are solved.

HAND CALCULATION

Note how equation (4) also has the form $\mathbf A \mathbf u^{(i)} = \mathbf 0$. Obviously, $\mathbf u^{(i)}=\mathbf 0$ would satisfy such an equation, but this is a trivial result that corresponds to no actual vibration. Instead, we are interested in values of $\mathbf u^{(i)}\ne\mathbf 0$ that solve the equation. In order for such a non-zero answer to exist, $\mathbf A$ must have no inverse, i.e. $\mathbf A$ is a singular matrix. (If an inverse did exist, then you could premultiply both sides of the equation with the inverse matrix, getting $\mathbf u^{(i)} = \mathbf A^{-1} \mathbf 0 = \mathbf 0$.) In order for a matrix to be singular, the determinant of the matrix must be zero. i.e.

$$\left| \mathbf A \right| = \left| -\omega_i^2 \mathbf M + \mathbf K \right| = 0$$

This will produce a polynomial in $\omega_i^2$, the roots of which will be each of the resonant frequencies squared. Take care to note that some modes will have the same resonant frequencies, in the form of repeated roots.

For example, say:

$$\mathbf M = \left[\begin{matrix} m & 0\\ 0 & m\\ \end{matrix}\right] \qquad \mathbf K = \left[\begin{matrix} 2k & -k\\ -k & 2k\\ \end{matrix}\right]$$

The determinant of interest would be:

$$\left|\begin{matrix} 2k-\omega_i^2 m & -k\\ -k & 2k-\omega_i^2 m\\ \end{matrix}\right| =0$$

By expanding out the determinant expression, we obtain the polynomial:

$$(2k-\omega_i^2 m)^2-k^2=0$$

The two roots of which are:

$$\omega_i^2 = \frac{k}{m},\frac{3k}{m}$$

So the two resonant frequencies are:

$$\omega_1 = \sqrt{\frac{k}{m}} \qquad \omega_2 = \sqrt{\frac{3k}{m}}$$

Note that you should obtain a number of (not necessarily distinct) eigenvalues equal to the number of rows or columns of your $\mathbf M$ or $\mathbf K$ matrices.

To obtain the eigenvector corresponding to a particular eigenvalue, substitute back into equation (4), to obtain a series of simultaneous equation, the solving of which should determine the eigenvector.

In our example, finding the eigenvector corresponding to the eigenvalue $\omega_1^2$, substitute $\omega_i=\omega_1$ into equation (4)...

$$\left[\begin{matrix} 2k-k & -k\\ -k & 2k-k\\ \end{matrix}\right] \left[\begin{matrix} u^{(1)}_1\\ u^{(1)}_2\\ \end{matrix}\right] = \left[\begin{matrix} k u^{(1)}_1-k u^{(1)}_2\\ -k u^{(1)}_1+k u^{(1)}_2\\ \end{matrix}\right] = \left[\begin{matrix} 0\\ 0\\ \end{matrix}\right]$$

From this, we can see that $u^{(1)}_1 = u^{(1)}_2$, and so the eigenvector is therefore:

$$\mathbf u^{(1)} = \left[\begin{matrix} 1\\ 1\\ \end{matrix}\right]$$

NUMERICAL COMPUTATION

Your problem does not seem practical to do by hand. Thankfully, there are a variety of different eigen-solver algorithms out there that are dedicated to solving this very problem.

A couple of these include:

• Power Iteration method
• Lanczos method
• Krylov-Schur method
• Jacobi-Davidson method

To explain how each of these algorithms work would take more time than I currently have, however I would definitely recommend tracking down these algorithms, reading how they work and using them to solver the generalised eigenvalue problem.

Hope this helps.

• Excellent explanation. However, for some of us who primarily do fluids and touch vibration only when someone asks "is there an ME here?", a tidbit on dampening and critical dampening would make this a very nice reference point to return to.
– Mark
Mar 9 '17 at 20:08
• @SprocketsAreNotGears This explanation is great! and is basically what I needed to get started. One question, In theory, would the natural frequencies I calculate from this analytical/numerical approach, be the same I get with a Experimental Modal Analysis? or are there other factors I should consider? Mar 13 '17 at 16:29