Estimating Wing Divergence, Control Reversal, and Flutter Onset Speeds

It is possible to make a simplified estimate of the free-stream velocity that torsional divergence, control reversal, or flutter occurs at by assuming the wing to be unswept and of high aspect ratio, that is, by assuming a simple two-degree-of-freedom dynamic model of the wing. We start from the natural frequencies of the full wing in its primary flapping and twisting modes, as calculated using FEA (or indeed, as measured in the lab, see later), and use these with a simple linearized aerodynamic model as follows: First we define a two-degree-of-freedom model where the wing motion is characterized by its vertical heave h and torsional rotation a about its elastic axis (i.e., the axis about which rotation occurs and where applied forces cause only bending) being restrained by two springs of stiffness K_{h} and K_{a}, respectively, see Figure 14.23. The wing is taken to have mass M and polar moment of inertia about the elastic axis I (these are for just the port or starboard wing and not the total for the pair). The wing is also subject to an aerodynamic lift force L_{wing} acting at its center of lift x_{ac} in front of the elastic axis, while its center of gravity lies x_{cg} behind the elastic axis. The coupled equations of motion of the system are then approximately (noting the coupling caused by the fact that

Figure 14.23 Two-degrees-of-freedom model of wing aeroelasticity.

the center of mass does not lie on the elastic axis) as follows:
and

We further assume that the two spring constants can be deduced from the natural frequencies in flapping and twisting as K_{h} = ®jj_{ap}M and K« = ®2_{wis}/- Assuming a steady-state thin airfoil behavior, which allows the effect of camber and fixed angle of attack to be superimposed later, the aerodynamic lift L_{wing} of one wing can be taken to be simply proportional to the twist as L_{wing} = «pV^{2}S/2, where p, V, and S are the density, airspeed, and the wing area of both wings, respectively, as usual. This leads to

and

If we next assume a simple harmonic solution to these equations, of frequency m and amplitude ^{0} , we get the following matrix equation:

«0

The nontrivial solutions to this equation arise when the determinant of the matrix is zero. This leads to a quadratic equation in m^{2} of the form Am^{4} + B®^{2}+ C = 0, where

The point at which the roots of this equation become imaginary is the flutter speed, while the solution for w = 0 is the divergence speed. Note that for a classical thin airfoil section we can assume C = 2л per radian. (Note: for lower aspect ratio wings we could modify this using the span efficiency.) Consequently, the flutter speed can be found from 4AC = B^{2} or

which is a quadratic equation in V^{2} of the form A'V^{4} + B'V^{2} + C = 0, where

Given the values of the natural frequencies, locations of the centers, and the inertial properties, it is then possible to solve for the flutter onset speed V_{flut}.

The divergence speed is found from ®2_{wis}/ = pV^{2}S/2, giving V_{div} = у

Notice that, if the aerodynamic center aligns with the elastic axis, x_{ac} is zero and divergence is not possible as would be expected, since the aerodynamic forces do not then cause the wing to twist; placing the main wing spar at the quarter chord point therefore helps eliminate divergence.

An essentially similar analysis can be carried out where a trailing-edge control surface is assumed to be present that alters the lift of the section in a way that is directly proportional to the control flap deflection 8; that is, we define the wing lift and moment coefficients as

and

where C_{L}/_{8} = C and C_{M}/_{8} = ^Cm. Then it may be shown that the control reversal speed

!-W~lCLZ

occurs when V = ^cC s/2 = Vev, where c is the mean aerodynamic chord used to nondi-

mensionalize C_{M}, again assuming C = .л. Note that C_{M}/_{8} will be negative. The quantities C_{L}/_{8} and C_{M}/_{8} can be estimated for the control flap in question using a range of methods or taken from past practice. Thin airfoil theory gives C_{L}/_{8} = 2(cos^{-1}(d) + V(1 - d^{2})) and

C_{M}/_{&} = 0.5(-aC_{L}/_{S} + cos^{-1}(d) - dJ(1 - d^{2})), where d is the distance of the aileron hinge behind the mid-chord point as a fraction of the semi-chord (so for a 20% chord aileron d = 0.6), and a is the distance of the aerodynamic center in front of the mid-chord point as a fraction of the semi-chord (so, typically, 0.5 for foils where the center of lift is at the quarter chord point).

If a is taken as 0.5 and d as 0.6, then the ratio C_{L}/_{&}/C_{M}/_{&} = 3.45/ - 0.64 = -5.39. If we

take this ratio and setx_{ac}to c/4, we can see that V_{rev}/V_{div} = 0.5^-C_{L}/_{S}/C_{M}/_{&} = 1.16, that is, wing divergence is likely to occur before control reversal. However, this would be applicable for an infinitely long wing with a full-length aileron; in practice C is reduced for finite-span wings, and the effectiveness of the control surface scales broadly with its extent along the wing, noting, however, that control surface effectiveness near the wing tip will be reduced by downwash effects.^{[1]}

Hoerner [10] devotes an entire chapter to the behavior of control surface flaps and shows that the flap effectiveness ratio = C_{L}/_{S} / dCt can be approximated as 4ж^/(1 - d)/2, with typical values of 40% for surfaces extending over the rear 20% of the section, that is, C_{L}/_{&} « 2.5 for d = 0.6. Experimental values of C_{L}/_{S} are shown to lie in the range 1.5-3.5/rad. Experimental values of C_{M}/_{S} are shown to be 75-80% of those predicted by thin airfoil theory and thus typically lie in the range -0.25 to -0.75 per radian, with a maximum around d = 0.48. Golland [27] shows the coefficients for an unswept wing with an aileron in the outboard half of the span as C_{L}/_{S}/C_{M}/_{&} = 2.64/ - 0.72 = -3.67. If we take this ratio, and again assume that x_{ac} is given by c/4, we can see that V_{rev}/V_{div} = 0.96, that is, control reversal is now likely to occur before divergence.

Given that an Abaqus model has already been built to assess the static stressing and stiffness of the structure, it is a simple matter to switch to a natural frequency calculation to ascertain ®_{flap} and ®_{twist} (noting that these are specified in units of per radian in the above analyses). All that is needed is to make sure that all the density properties of the material models are correctly specified and any lumped masses added before solving for the natural modes.

To solve for the modes, when creating the analysis step, select the “Linear perturbation” procedure type in the “Create step” window and then chose a “Frequency” analysis. This gives a form to complete that specifies the desired analysis. We choose to find the first 10 eigenvalues, leaving all other settings as default values. The natural frequencies and corresponding mode shapes can be found in the post-processing viewport by examining the displaced shapes. The natural frequency can be seen at the bottom of the viewport, and the corresponding mode shape occurring at this frequency is shown by the contour plot or by animation. Clicking the arrow icon at the top right of the viewport allows the various modes that were requested to be stepped through. The position of the elastic axis can be determined by examining the mode shape in twist; the axis lies along the nodal line in the mode shape. This can be revealed by examining a truncated contour plot on an undeformed plan view of the wing, see, for example, Figure 14.24. Note that this is not wholly along the spar line because the spar itself flexes, even in the twist mode.

The overall total mass of the FEA model and its various moments of inertia are printed to the Data File and should be used to check that correct totals have been generated with appropriate radii of gyration. Note that since only one wing is included in the preceding analysis when using the FEA analysis to generate inertia data and radii of gyration, only those parts that make

Figure 14.24 Truncated Abaqus contour plot of a first twist mode revealing the nodal line and hence the elastic axis.

up the wing under investigation should be included; while when establishing mode shapes and natural frequencies, the more complete the model the better. These results can then be combined with estimates for the aerodynamic coefficients to assess possible aeroelastic effects.

One key aspect of the natural frequency analysis is the constraint set imposed on the wings in the FEA model that is designed to represent their attachment to the fuselage. If the inboard edges of the SLS nylon parts that support the wings and tailbooms are heavily constrained, then it is likely that the flap and twist frequencies will be higher than those experienced in practice. If, alternatively, the SLS nylon parts are not restricted in motion with restraints just applied to the main spar, then the resulting frequencies can be too low. Of course, the best approach would be to model the fuselage itself and then use a contact model between the wings and fuselage, but this may involve more effort than can be justified for the aeroelastic checks at this stage of design. Our practice is to check the flap and twist modes with several differing constraint models and then use the resulting frequencies as a set to investigate possible aeroe- lastic problems. For Decode-1, the frequencies we found (along with those from subsequent experimental testing) are given in Table 14.3, while Figure 14.25 illustrates the first flap and twist mode shapes. Note that the twist frequency is relatively unaffected by these changes, and it is this frequency that dominates the aeroelastic behavior. Clearly, the unconstrained model best matches the experimental frequencies in this case.

For Decode-1, the Abaqus-calculated wing weight for just the port wing but including the entire main spar is 1.224 kg centered 0.134 m behind the wing leading edge (0.038 m or 12% of the mean chord behind the spar center), with moments of inertia about the origin of 0.867, 0.003, and 0.886 kg m^{2} and thus a twist radius of gyration of 0.049 m or 15% chord.^{[2]} Examination of the twist mode suggests that the elastic axis can be taken to approximately lie along the main spar, that is, at the quarter chord position. The ailerons are 25% of the chord and

Table 14.3 Natural frequency results (Hz) using Abaqus modal analysis for the Decode-1 airframe.

Model

First

flap mode

Second flap mode

First

twist mode

Main spar clamped and inboard SLS part restrained in X, Y, and Z

9.2

53.4

66.7

Main spar clamped and inboard SLS part restrained in Y only

9.1

51.4

64.9

Main spar clamped and inboard SLS part restrained by torque peg

7.2

43.6

63.7

Main spar clamped only

6.3

40.1

63.1

Experimental (see Chapter 16)

5.8

31.7

59.3

Figure 14.25 Abaqus plots of first flap and twist modes for Decode-1 wing.

extend over 50% of the span. Plugging these numbers plus the lower set of calculated natural frequencies into the above formulae, and assuming thin wing theory, gives onset speeds of V_{div} = 68.1 m/s, V_{rev} = 82.6 m/s, and V_{flut} = 63.6 m/s, all comfortably above the maximum dive speed of 40 m/s. Reducing d^L to allow for a finite wing span would increase these speeds slightly.

It must be stressed, however, that all these calculations are greatly simplified from a full aeroelastic analysis. In particular, the assumption of steady-state aerodynamic equations is a significant approximation. It is possible to use an unsteady potential flow solution for the thin airfoil model and get better approximations, or to use a coupled aerodynamic and structural approach, using XFLR5, Fluent (in either steady or unsteady modes), Abaqus, and so on, but such methods lie beyond the scope of this book. And in any case, unless the initial estimates given in this section suggest that the various aeroelastic onset speeds are close to the maximum dive speed, they are almost certainly not warranted for the small fixed-wing UAVs considered here. We do, nonetheless, check the natural frequencies of the wings after assembly by simple experimental tests, as will be described shortly. We also routinely measure control surface characteristics in wind tunnel tests.

[1] Note also that the aspect ratio of the control surface itself plays a part when such aspect ratios reduce below 2; veryshort control surfaces should thus be avoided.

[2] It is important to make sure that a reduced Abaqus model without the tail boom or any of the empennage is usedwhen deriving the inertial properties of the wing since the simplified aeroelastic model used to derive the onset speedsdoes not include these, even though these features will probably be included when establishing the natural frequencies.