Within the Aerospace Engineering community, there is an entire sub-discipline devoted to understanding the dynamics of a system and generating loads (propulsive, inertial, aerodynamic, etc). In smaller companies, engineers often need to wear multiple hats, and the lines between classical stress analysis and loads and dynamics analysis begin to blur. A few of our engineers, here at Structural Design and Analysis (SDA), worked with a local resident who had taken it upon himself to build a Fokker D-VII Biplane from scratch, and wanted to know how much weight he could save if he used an aluminum spar for the main wings instead of the original wooden spar design. Our engineers had to develop a finite element model (FEM) and conduct the basic loads and dynamics analysis to define the load cases for the vehicle.

Generating aerodynamic loads is relatively straight forward for aircraft with more conventional designs. Typically, a combination of 2D aerodynamic theory, as well as corrections for wings with finite span, are used to generate the loads in the early stages of the design phase. These loads are then applied to the structure at the quarter chord location of the wing. The Fokker is a biplane, which means that there are interference effects between the upper and lower wing which cannot be neglected and must be considered when determining the loads of the aircraft. The goal of this case study is to show that various approaches can be taken to solve this loads generation problem, and that the “best” approach for an engineer is dependent on their technical expertise, their available resources (time and/or money), and the desired accuracy of the results. Three different methods were selected to calculate the aerodynamic loading on the Fokker D-VII biplane, and they are listed in increasing order of required technical expertise and accuracy: The three different methods:

The first assumes each wing acts separately from one another. This type of solution is best suited to an aircraft enthusiast, or an engineer without much background in theoretical aerodynamics.

The second method applies corrections to account for the interaction between the upper and lower wing. This type of solution is best suited to an engineer with a BS in Aerospace Engineering. The third method implements NX Nastran’s Static Aeroelastic SOL 144 for the loads generation. This solution requires the least amount of effort on the user since the loads are calculated internally by Nx Nastran, but is best suited to an engineer with at least an MS in Aerospace Engineering. This article will compare the efficacy of the three methods, as well as cover the accuracy of their respective results.

Method 1 The first step in calculating the aerodynamic loads on the aircraft is to get the airfoil data. The Fokker D-VII uses a modified Goettigen GOE 418 airfoil for the upper wing. The airfoil data points were imported into XFOIL, a popular open-source 2D potential flow code, and lift coefficients were extracted for various angles of attack (AOA). The airfoil is scaled to a unit chord and plotted in XFOIL. it has a 15% thickness and a 6% camber. The airfoil shape is shown below. Using the XFOIL data, a plot of the lift polar is constructed. A trendline is added to the data to estimate the lift curve slope for the airfoil of Cl,α= 0.095/° and the zero lift AOA is -8 degrees. This value will be used to calculate the lift on the wings. To calculate the lift on the upper and lower wings, a simple approximation from Prandtl Lifting Line Theory is used which relates 3D lift coefficient to 2D lift curve slope, Aspect Ratio (AR) and AOA. The lower wing of the Fokker D-VII has 1 degree of incidence while the upper wing has no angle of incidence. The lift equation was then used to calculate the lifting force on the wings. A moment balance is conducted to calculate what tail force is needed to balance the aircraft.

The sum of the loads is equal to 5,394.32 lb or 4.31 g’s. Since this is a 4.0 g load case, the lift on the wings will need to be reduced. As the lift on the wings is reduced the pitching moment will change which, in turn, will change the required tail force to balance the aircraft. Excel’s goal seek was used to reduce the wing loading and balance the aircraft such that the total lift (including the tail) is equal to 4.0g’s. The final loads are shown below. These final loads are applied to the ¼ chord location of the wings. A rectangular spanwise lift distribution is applied to the upper and lower wings.

Method 2 By having two wings subject to the same flow, each wing interacts with the other. Effectively, each wing will interact with the other’s vortex system such that the upper wing will experience an increase in lift and the lower wing will experience a decrease in lift, denoted ΔC L U and ΔC L L respectively. The following method uses the simple biplane theory which is detailed in NACA Technical Report No. 458, Relative Loading on Biplane Wings [1]. It is shown that the values for ΔC L follow a linear relation to the overall vehicle C L in the following form:

Where K 1 and K 2 are constants relating to wing gap, stagger, decalage, overhang, and wing thickness. The change in lift for the lower wing is related to the change in lift of the upper wing by the ratio of wing areas. The methods of finding the values of K 1 and K 2 follow a graphical approach using the biplane ratios of gap to chord, stagger to chord, max thickness to chord, percent overhang, and decalage. The final lift coefficient for the upper and lower wings are found by adding the correction to the vehicle lift coefficient. The Fokker D-VII biplane configuration is shown in the table below. These values are used for calculation of the upper and lower wing lift coefficients. Using these values and the method described in NACA Report 458, the following values are calculated for K 1 and K 2 : The Max Weight of the Fokker D-VII is 1,259 lbs. For a 4.0g (n = 4) maneuver, the aircraft lift coefficient is calculated as follows: Plugging in values for K 1 , K 2 , and C L into the ΔC L equations give the following values: Using these new corrected values for the upper and lower wing lift coefficients, the total load can be calculated for the upper and lower wings. A moment balance is performed and the following loads are calculated for the aircraft: The loads are applied to the ¼ chord position of the aircraft wings. Two different spanwise lift distributions are applied to the model for this comparison study. The first assumes an elliptical lift distribution. The second uses Schrenk’s Approximation to estimate a more accurate spanwise lift distribution. These two distributions are shown below along with a reference to a rectangular distribution.

Method 3 The third and final method is to use NX Nastran’s static aeroelastic SOL 144 analysis to generate the loads using a vortex lattice formulation. Below is a picture of the potential flow model created in FEMAP to generate the aerodynamic loading for the Fokker. As with most potential flow methods, the NASTRAN aeroelasticity module models all aerodynamic lifting surfaces using flat plates, which is acceptable for any thin airfoil (less than 12% t/c). To account for effects such as aerodynamic twist, or in this case a constant camber, an initial downwash is imposed on the wings. As discussed in the first method, the airfoil used by the aircraft’s upper and lower wings has a zero-lift angle-of-attack of -8 degrees, and the lower wing has a built-in angle of incidence of 1 degree and so the enforced downwash allows those characteristics to be taken into account. One of the powerful functionalities about SOL 144 trim analysis is that given high-level information about any flight condition, Nastran can not only calculate the aerodynamic forces, but can also ensure that the vehicle is stable. Below is a picture showing the high-level information used to generate the 4G pull up case. For this case, the pitch and plunge degrees of freedom of the aircraft FEM were left free. To balance the equations of motion, there must be two aerodynamic degrees of freedom free to be used as variables to trim the aircraft’s stability. In this case, the variables were the aircraft angle of attack, and the elevator deflection angle. With just a few clicks, this load case can be modified to model any corner of the flight envelope by changing the dynamic pressure and load factor of the aircraft. Since Nastran calculates all the loads internally using this high-level flight condition information, it can save an incredible amount of time that might be devoted to calculating loads externally, bringing them in, and applying them to the structure accurately. The way the solution works is as follows: Given the flight case conditions, Nastran calculates the rigid Aerodynamic Influence Coefficient (AIC) matrix using a vortex lattice panel method This AIC matrix can then be transformed into an incremental stiffness matrix due to aerodynamic loads (K_AIC) Nastran then solves the following equation: Where q is the free-stream dynamic pressure. To those familiar, this form might appear similar to the classical linearized buckling system of equations. This solution requires the least amount of effort on the user since the loads are calculated internally by NX Nastran and then applied to the FEM through the points defined in the model. NX Nastran automatically calculates the trim condition of the aircraft and calculates the loading on the upper wing, lower wing, and horizontal tail. The resulting loads are shown below.

Results As the complexity increases with each of the methods discussed above, so does the accuracy of the results. However, not every stage of the design requires the same precision in results. Since this discussion focusses on how the loads approach impacts the design of the aircraft, let us first compare the calculated loads for all three methods. Below is a table outlining the loads for each method and the percentage difference when compared to the aeroelastic model. When comparing the lifting force of upper and lower wings of the aircraft, the aerodynamic loading from method 1 underestimates the lift on the upper wing by 10.2% and over estimates the lift on the lower wing by 22.5%. Applying Simple Biplane Theory in method 2 captures the interference effects and estimates the wing loading much more accurately with the upper wing lift being only 2.6% less and the lower wing lift being 6.3% greater, when compared to the aeroelastic model. A more detailed way to compare the resulting design impact of the three different load generation methods is to look at the internal shear and bending moment within the spars. Below is the shear force diagram for the upper wing leading edge spar. Starting with the simplest method, the 2D non-interfering lift generates the highest shear force. As one would expect the variation in the maximum internal shear force is small (at most 7% difference) over all the models since the total lift generated by the plane was set to be constant. The differences are partially due to how the lift is distributed between the upper and lower wings. The spanwise distribution also clearly has an impact on the internal shear force. Interestingly, the model that matches the Aeroelasticity model shear force the closest is the Biplane Theory model using Schrenk’s approximation to capture span wise loading effects. The differences produced by these aerodynamic models becomes even more apparent when inspecting the internal bending moment as one might expect. The internal bending moment clearly shows how the differences in aerodynamic models can propagate. The most basic model (2D Non-interfering Lift) produces the highest bending moment, 33% higher than the Aeroelastic solution (the most accurate loads generation approach). While it is better to be over conservative than under-conservative, this kind of inaccuracy could lead to substantially higher weight in the structure, limiting the performance. Again, of all the purely analytical models, the Biplane Theory using Schrenk’s approximation does the best job at matching the Aeroelastic model with only a 2.3% higher prediction in the maximum bending moment. The Biplane Theory using an elliptical distribution produces a substantially lower maximum bending moment, 32% lower than that predicted by the Aeroelastic model. Finally, it should be noted that since the Rigid Aeroelasticity and Aeroelasticity curves are so similar for both the leading-edge spar shear force and bending moment, it can be inferred that at the max dynamic pressure of the Fokker (125 mph at sea level), the flexibility of the structure is not great enough to impact the loading. A similar narrative is replicated in the upper wing trailing edge spar. The maximum shear force produced by the 2D non-interfering lift was 30% larger, the maximum shear force produced by the Biplane theory using an elliptical lift distribution predicted was 33% lower, and the maximum shear force produced by the Biplane Theory using Schrenk’s approximation was 1.3% greater, producing the closest correlation when compared with the most accurate Aeroelasticity generated aerodynamic model. Much as the leading-edge spar exhibited disparate internal bending moments depending on the aerodynamic model used, the trailing edge spar internal bending moment varied greatly. In this case, the max internal bending moment experienced by the trailing edge spar was 49% larger using the 2D non-interfering lift, 34% lower using the biplane theory with an elliptical distribution, and 13% higher when using the Biplane theory using Schrenk’s approximation. The shear and bending moment diagrams are often excellent indicators for simple structures of the internal stress state within the structure. As a stress engineer, comparing stress plots is the most meaningful way to compare how the different aerodynamics models impact the stress throughout the vehicle. Below is a picture of the upper wing leading and trailing edge spars Von Mises stress under 4G pull up when using the Aeroelasticity model. The max spar cap Von Mises stress is 11.4 ksi. In comparison, the same stress contour is presented below, only for the case using the Biplane theory using Schrenk’s approximation, which exhibited a max Von Mises stress of 11.5 ksi, a 0.8% difference from that produced using the aeroelastic aerodynamic model. In contrast, the Von Mises stress state for the upper wing under the 2D non-interfering lift can be seen below as a gross overprediction in the stress within the wing, predicting a max Von Mises stress of 15.2 ksi, 33% higher than the aeroelastic aerodynamic model.