Helicopter blades are flexible and bend significantly in flight. You can see it in high speed video like this. To analyze a helicopter in detail, this bending must be computed. It affects the aerodynamics and even the handling qualities of the helicopter. This article explains one of the simpler methods to compute blade bending, namely the use of normal mode shapes. While simple relative to other approaches, this is still not easy. This is a technical article and requires a calculus background.

The physics are based on classical beam bending theory. For small displacements, the bending moment is proportional to the curvature of the blade, i.e. $$\begin{equation} K(x)y''(x) = m(x) \label{eq:beamth}, \end{equation}$$

where \(m\) is the bending moment, \(K\) is the stiffness of the blade (depending on the construction), \(x\) is the distance along the blade, \(y(x)\) is the distance the blade has moved from the undeformed position and \(y''(x)\) is the second derivative of \(y\) with respect to \(x\) (the curvature). This equation, along with inertia and the external forces applied to the beam (aerodynamics), govern the “elastic motion” of a blade. We’ll start with a simple, static example below and gradually move on to the general scenario.

In this example, we consider a rotor blade that is not moving and cantilevered at the hub. Cantilevered means it's clamped down so that \(y(0)=y’(0)=0\), where \(x=0\) is the blade root / rotor hub location. An external force \(F\) holds the blade tip up at \(x=R\) (\(R\) denotes the length of the blade). Hence, the force creates a moment

$$\begin{equation} m(x) =F(R-x) \end{equation}$$

about any location \(x\) along the blade. To simplify the arithmetic, we'll assume the blade’s stiffness is constant here, i.e. \(K(x)=K\). Substituting this into \eqref{eq:beamth} gives $$\begin{equation} Ky^{\prime\prime}(x) = F(R-x), \end{equation}$$ which has the solution $$\begin{equation} y(x)=\frac{F}{K}(\frac{Rx^2}{2}-\frac{x^3}{6}), \end{equation}$$ after using the cantilevered boundary conditions. Hence, the tip of the blade (\(x=R\)) is displaced vertically by \(\frac{FR^3}{3K}.\)

In this example, we remove the force \(F\) that was applied to the tip of the blade in the prior example. There are no external forces acting on the blade, only its internal inertia and stress. As you might expect, the blade will not just settle into the undeformed position, but instead will oscillate about that position. The situation is much akin to a spring / mass problem, but the restoring force is proportional to blade curvature instead of displacement. Hence, you might expect an oscillatory solution where the blade vibrates about it's equilibrium shape. This will be the case, but with a slight added complication. A 1-dimensional spring/mass problem has just one characteristic frequency that it oscillates at. A helicopter blade will oscillate at various frequencies, each associated with a different "mode shape" of the rotor blade (like a stringed musical instrument). We will find that the differential equation \eqref{eq:beamth} has a solution

$$\begin{equation} y(x,t) = \Sigma_n y_n = \Sigma_n\hat{y}_n(x) f_n(t), \label{eq:freevsol} \end{equation} $$

where the \(\hat{y}_n\) are these mode shapes (only a function of \(x\), not \(t\)) and the \(f_n\) are the time-oscillatory terms (only a function of \(t\), not \(x\)). Below we will derive the \(\hat{y}_n \) and \(f_n\). We start with the second derivative of \eqref{eq:beamth} and then replace \(m\) with the stress function \(V(x,t)=\frac{\partial m}{\partial x}\) and use the fact that \(\frac{\partial V}{\partial x}\) is the inertial force per unit length \(-\rho \ddot{y}\), where \(\rho\) is the mass per unit length of the beam. This gives the following equation

$$\begin{equation} Ky_n'''' = m'' = V' = -\rho \ddot{y_n}. \end{equation}$$

We consider a single term \(y_n\) in the sum \eqref{eq:freevsol} since the sum of such solutions is a solution. Notice that \( y_n''''=\hat{y}_n'''' f_n \) and \(\ddot{y}_n = \hat{y}_n \ddot{f}_n\) and therefore $$\begin{equation} K\hat{y}_n'''' f_n = -\rho \hat{y}_n \ddot{f}_n \end{equation} $$ or $$\begin{equation} \frac{K\hat{y}_n''''}{\rho\hat{y}_n} = \frac{-\ddot{f}_n }{f_n} \label{eq:freevsid} \end{equation} $$

Since the left side of this equation is only a function of \(x\) and the right side is only a function of \(t\), both sides must be constant. This provides two new differential equations by equating the numerator of each side to a multiple of the denominator. The right side shows that the second derivative of \(f_n\) is a multiple of \(f_n\) which is a well known oscillatory system with (real) solutions $$ \begin{equation} A_n\cos \omega_n t + B_n \sin \omega_n t \end{equation}$$

where the numerators are \(\omega_n^2\) times the denominators. The left side takes more work. Let's first collect the constants into \(\beta_n^4 = \frac{\rho \omega_n^2}{K}\) so that the equation is

$$\begin{equation}\hat{y}''''=\beta_n^4 \hat{y}. \end{equation}$$

This differential equation has (real) solutions

$$\begin{equation} \hat{y}_n = C_ne^{\beta_n x} + D_ne^{-\beta_n x} + E_n\cos\beta_n x + F_n\sin\beta_nx \label{eq:freeva}. \end{equation} $$

We can drop the \(C_n\) since we can scale \(\hat{y}_n\) as we wish (ultimately it will be scaled by \(f_n\) in the full equation, where the scale matters). We can apply the boundary conditions (\(y(0)=y'(0)=y''(R)=0\)) to each \(\hat{y}_n\) individually, even though the real solution is a sum of these terms. This is because we could start with the beam in any particular "mode shape" \(\lambda \hat{y}_n\) and apply the same logic / equations as the beam would always be shaped by a multiple \(f_n (t)\hat{y}_n\) of \(\hat{y}_n\) by definition. With a bit of arithmetic, the first two equations (at the blade root) imply that \(1+D_n+E_n=0\) and \(1-D_n+F_n=0\), respectively. Using this in the last equation \(\hat{y}''(R)=0\) gives

$$\begin{equation} D_n = \frac{e^{\beta_nR} + \cos\beta_nR+\sin\beta_nR}{\sin\beta_nR-e^{-\beta_nR}-\cos\beta_nR}, \end{equation}$$

and therefore

$$\begin{equation} \hat{y}_n=e^{\beta_nx}+D_ne^{-\beta_nx}-(1+D_n)\cos\beta_nx+(D_n-1)\sin\beta_nx . \label{eq:freevb} \end{equation}$$

What's left is to solve for the \(\beta_n\) and therefore the frequency \(\omega_n\) associated with \(\hat{y}_n\). This can be done using the last boundary condition \(y'''(R)=0\). We can apply this to \(\hat{y}_n\) as done above to get the equation \(\hat{y}_n'''(R)=0\). Taking the third derivative of Equation \eqref{eq:freevb} and dividing by \(R^3\) gives

$$\begin{equation} e^{\beta_nR}-D_ne^{-\beta_nR} - (1+D_n)\sin\beta_nR-(D_n-1)\cos\beta_nR = 0. \end{equation}$$

This equation has an infinite number of solutions (associated with an infinite number of mode shapes) and can be solved numerically. The smallest three positive solutions are \(\beta \approx 1.875/R, 4.694/R, 7.855/R\). From the definition of \(\beta_n\) above we have \(\omega_n^2 =\frac{\beta_n^4 K}{\rho} \), so the first frequency is \(\omega_1 \approx 3.516\sqrt{ \frac{K}{R^4\rho} }\). The mode shapes (from Equation \eqref{eq:freevb}) associates with these frequencies are plotted below.

If forces push a blade into the first mode shape \(\hat{y}_1\) shown above, the blade will be shaped by a multiple \(f(t) \hat{y}_1\) as long as no further forces are applied. The blade will oscillate about the equilibrium position with the frequency \(\omega_1 \approx 3.516\sqrt{ \frac{K}{R^4\rho} }\). The same is true for any mode shape \(\hat{y}_n\), but with a different frequency \(\omega_n\). In general, a blade will be in some linear combination of all the mode shapes, but we can typically approximate the shape well with just the first 10 or so modes, say \(\lambda_1\hat{y}_1+\lambda_2\hat{y}_2+...+\lambda_{10}\hat{y}_{10}\).

We now have all the pieces of \eqref{eq:freevsol} except the constants \(A_n\) and \(B_n\) in \(f_n\). These constants may be found from the initial conditions of a given problem. For example, if the blade starts out with no velocity \(\dot{y}=0\) then all the \(B_n\) must be 0.

We estimated static blade deflection and then found the frequencies (with mode shapes) of a freely vibrating blade. In order to be useful for real-world problems, we’ll need to compute the motion of nonuniform blades (\(K(x)\), \(\rho (x)\) not constant) with arbitrary forces acting on them. Such forces include aerodynamic forces \(F(x)\) acting along the blade, inertial forces due to hub motion, and centrifugal forces due to rotor rotation. The math gets tedious and is solved numerically on a computer. We won’t get into the tedious math here, but will outline how this is done in many software models of rotorcraft. If interested, the standard reference for this work with details we’ve omitted is NACA Report 1346 by Houbolt and Brooks (this does maintain several assumptions including no plastic deformation, small angles of elastic twist and neglects some inertial forces).

Like Houbolt and Brooks, we’ll consider models that use the mode shapes of free vibration \(\hat{y}_n\), which we computed above. Although the blade won’t be vibrating freely here, this approach works well if we replace the \(f_n\) with appropriate time varying quantities called “mode participation factors,” denoted \(\delta_n\). The \(\delta_n\) are a function of the forces applied to the blades. Only a finite number \(N\) of modes are used (the ones with the \(N\) smallest frequencies \(\omega\)), so that the blade shape at time \(t\) becomes

$$\begin{equation} y(x,t) = \Sigma_{n=1}^N \delta_n(t)\hat{y}_n(x). \end{equation}$$ Increasing \(N\) may increase fidelity, but reduces efficiency. Depending on the context, from 2 to 20 modes are typically used in the software.

The blade will also be discretized radially into a number of segments as shown in the diagram below. Typically 5 to 25 segments will be used. The user must input blade properties like \(K\) and \(\rho\) at each segment. The software will estimate the forces acting on said segments and use those forces to compute the mode participation factors \(\delta_n\).

For simplicity, our examples above included only one degree of freedom (DOF) along the blade - \(y\) had a single value for each pair \(x,t\). Models for real-world application typically include three DOF: out-of-plane, in-plane and twist motion. The motions are all coupled, so that a mode shape consists of the three values at each radial section (but still a single mode participation factor that scales all 3 values by the same amount).

The question now is how to compute the \(\delta_n\). We won't derive this, but as you might expect the \(\delta_n\) are “accelerated” by applied forces. They are driven by a second order differential equation

$$\begin{equation} \frac{d^2\delta_n}{d t^2}+ 2\zeta_n\omega_n \frac{d\delta_n}{dt}+\omega_n^2\delta_n=\frac{F_n}{I_n}, \end{equation}$$

where \(I_n\) is a “generalized inertia” associated with the mode shape \(\hat{y}_n\), \(F_n\) is the “virtual work” done on \(\hat{y}_n\) and \(\zeta_n\) is a sometimes added artificial damping term. These values may be computed under various assumptions, e.g. assuming certain hub motions or blade pitch acceleration are negligible. We won’t go into details here, but roughly speaking the virtual work is the inner product of the forces and the mode shape \(\int_0^R F(x)\hat{y}_n(x)dx\). You may think of this as the force applied to mode \(n\).

Stay tuned, we plan on adding an example here and even posting source code to solve such problems.

Examples of comprehensive rotorcraft analyses that use some or all of these concepts include

- C81 / COPTER from Bell Flight,
- C-60 / Tech-02 from Boeing Vertol
- Y-200 / SIMVIB from Sikorsky Aircraft,
- 2GCHAS / RCAS from the US Army,
- CHARM from Continuum Dynamics,
- Flightlab from ART,
- Dymore from the Georgia Institute of Technology and
- CAMRAD from Johnson Aeronautics.

A survey of these and other models was done by Wayne Johnson: A History of Rotorcraft Comprehensive Analyses.