Nonlinear vortex lattice method for stall prediction

. The stall behavior of an empennage is a crucial and conditioning factor for its design. Thus, the preliminary design of empennages requires a fast low-order method which reliably computes the stall behavior and which must be sensitive to the design parameters (taper, sweep, dihedral, airfoil, etc.). Handbook or semi-empirical methods typically have a narrow scope and low fidelity, so a more general and unbiased method is desired. This paper presents a nonlinear vortex lattice method (VLM) for the stall prediction of generic fuselage-empennage configurations which is able to compute complete aerodynamic polars up to and beyond stall. The method is a generalized form of the van Dam algorithm, which couples the potential VLM solution with 2.5D viscous data. A novel method for computing 2.5D polars from 2D polars is presented, which extends the traditional infinite swept wing theory to finite wings, relying minimally on empirical data. The method has been compared to CFD and WTT results, showing a satisfactory degree of accuracy for the preliminary design of empennages.


Introduction
The core of the methodology is a vortex lattice method coupled with viscous aerodynamic polars. This is achieved with a generalized version of the α-method presented by van Dam [1]. The vortex lattice method implementation employed is described in section 2. The viscous coupling algorithm is presented in section 3, under the assumption that the section polars are available. Section 4 presents a methodology for computing section polars from 2D polars. Finally, some results obtained with the present method are shown in section 5. This project has been supported entirely by the Clean Sky 2 program.

Vortex lattice method
There are several possible implementations of vortex lattice methods (VLM). Katz and Plotkin [2] present a thorough explanation of the method. For this implementation, a strip method with horseshoe elements is chosen. As the chordwise information will be provided by the viscous polars, there is no benefit to increasing the number of chordwise panels.
A wing is described by its origin position and a set of parameters as a function of its local spanwise coordinate �, namely the chord ( �), the and offset ��� ( �) and ��� ( �), the geometric twist ( �), and the local airfoil. The wing is discretized into spanwise panels. The discrete geometrical parameters are taken from the center of the panel in the local � direction. For each panel, the relevant local parameters are the chord , the dihedral , the geometric quarter-chord sweep Λ � � , and the twist . The normal vector of each panel is computed and a control point is placed at the 3 /4 line, in the middle of the panel, as shown in Figure 1. Horseshoe vortices (HSV) are used as singularity elements. The HSV head is placed at the quarter-chord line, and each leg extends chordwise downstream over the wing, then in the direction of the free stream towards infinity ( Figure 1). Thus, the effect of a HSV can be computed as the sum of five vortex filament segments. The solution for a vortex filament, as well as a code subroutine, can be found in Katz and Plotkin [3].
For a flight condition given by the angles of attack and sideslip and , the free-stream velocity is � = � (cos cos + cos sin + sin ) (1) The velocity induced by the HSV of panel on the CP of panel is , . Since the intensity of the HSVs at this stage is unknown, it is useful to use � , = , /Γ . The total velocity on the CP of panel is The Neumann boundary conditions arise from the fact that the flow cannot penetrate the wing surface, so a zero normal velocity is imposed at the CPs, ⋅ = 0, yielding: and is a function of the flow angles and seen by each panel. The HSV circulations are computed solving the linear system, and the � of each panel is found with the Kutta-Joukowski theorem. [

Viscous coupling
The VLM can be coupled with 2.5D viscous polars to model the nonlinear behavior of the wing. In this section, it will be assumed that the 2.5D polars for each wing section are available (see section 4 for a method for computing 2.5D polars from the airfoil geometry). In practice, the 2.5D polars are computed for certain key sections (such as the root and the tip) and the rest are linearly interpolated. The coupling is achieved by introducing a local angle of attack correction for each panel Δ∆ � , such that the flow angles seen by each panel are modified. This alters the [ ] vector in equation 6 and, consequently, the HSV intensities [Γ] and � distribution. The flow angles and seen by a given panel are related to the local angle of attack � and the local dihedral through: = � cos (6) = � sin (7) The modified van Dam algorithm is as follows (where the subindex indicates that the computation must be done for each panel): a) Set Δ∆ = 0 b) Find the corrected flow angles * and * , and compute [ ] with the corrected angles: (9) c) Solve the VLM. The resulting lift coefficient of each panel is � � . d) Find the effective and induced angles of attack of each strip: e) Interpolate the �� in the 2.5D polar of the corresponding panel to find � �� and � �� .
f) Perform an axes transformation to find the lift and drag coefficients in the local flow axes.
� � = � �� cos � − � �� sin � (11) g) Update the value of Δ∆ . The introduction of a smoothing factor is not strictly necessary but beneficial in practice to aid numerical convergence.
h) Find the maximum error from all panels: i) If the error is below a certain tolerance , take � � as the converged value for each panel; otherwise, repeat from step b).

Aerodynamic polars: 2D to 2.5D
The main difficulty of the described method is computing the 2.5D polars. They can be can be obtained from CFD results or wind-tunnel tests, though when the results of such costand time-extensive methods are available, the present method is evidently redundant. On the other hand, 2D polars can be easily computed with a number of airfoil solvers, such as XFOIL [4], MSES [5] or 2D CFD codes. This section presents a novel method for transforming the 2D polars into the 2.5D polars given the wing geometry.
The difference between the 2.5D and 2D polars arises from the sweep, which produces a spanwise boundary layer flow (from root to tip if the sweep is backwards, as shown in Figure  2). The main consequence of this is that the BL is thinner at the root (higher effective Reynolds, delaying stall and achieving a higher � � ) and thicker at the tip (having the opposite effect). Thus, backwards sweep promotes stall from the tip, and forward sweep promotes stall from the root, as can be seen in the tuft studies from NACA TN 2445 [6]. In the general case of arbitrarily large sideslip and dihedral angles, calculating the sweep angle of a wing or wing section becomes nontrivial. The local sweep of a wing section is a factor of both the wing geometry -a combination of geometrical sweep and dihedral-and the flow angles. The relevant sweep angle considered here is that formed between the projection of � onto the panel surface and the perpendicular to the quarter-chord line. The projection of a generic vector onto a plane defined by its normal vector can be written as: where � and are the tangential and normal components of with respect to the plane. Thus, the projection of the free-stream velocity � onto the th panel with normal vector is Using the local * axis unitary vector * , the local sweep seen by panel is (Figure 3a) In the case of an infinite swept wing, the 2D airfoil polar can be transformed into the wing section 2.5D polar with the classical sweep theory, as shown in the ESDU methodology [7]. Mariens et al. [8] present a method for the transformation of � due to sweep. In practice, and especially at low Mach numbers, some simplifications can be taken, which greatly speed up the computation process without significantly compromising the results.
Given an infinite-span wing of constant sweep Λ flying at Mach number and Reynolds number , the proposed simplified methodology is the following: a) Extract the airfoil geometry from the wing section. b) Compute the airfoil polar with a 2D aerodynamic solver at Mach number and Reynolds number , obtaining � 2� ( 2� ) and � 2� ( 2� ). c) The 2.5D polar is The effect of wing sweep, as well as the correlation with the simplified methodology, is shown in Figure 3b. It can be observed that, despite the errors in the original XFOIL polar (especially in the negative stall region) the proposed transformation is adequate. In a finite swept wing, the 2.5D polars vary along the span, having a higher � � at the root than at the tip. This behavior is discussed in a paper by Hosangadi et al. [9] using the CFD results from [10]. We propose a method for modelling the 2.5D polars of finite wings based on assigning an effective (i.e. artificial) sweep to each spanwise section, then applying the sweep correction described in section 4. The effective sweep for the th panel is where the effective sweep correction Λ i is hypothesized to be a function of only the spanwise position and aspect ratio. To account for the possibility of an increase in � � as occurs at the wing root, equation 26 is replaced by A small set of wings of various aspect ratios and sweep angles, including those presented in Hosangadi et al. [9], has been studied. The spanwise distribution of the � � of these wings indicate that the hypothesized Λ function does indeed exist. Assuming that the effect of the root and the tip are independent and can be superimposed, the proposed sweep correction function is modelled as where ��� ) and � � ( � � � ) (Figure 4a) have been designed with Bézier curves to best fit the data gathered from the aforementioned set of wings. The resulting Λ ( ) function is plotted in Figure 4b for several aspect ratios. Note that for large aspect ratios, the solution for the infinite swept wing ( Λ = 1) is recovered for an area in the middle of the wing. The formulation of Λ is valid for approximately > 5; for smaller aspect ratios, the effect of the root overpowers that of the tip and the result is nonphysical.

Results
The lift curves of the set of wings from Hosangadi et al. [9] as computed by CFD and the present method are shown in Figure 5a. It may be observed that the present method has the proper sensitivity to the wing sweep. Figure 5b shows the comparison between the present method and NACA wind tunnel tests for a swept wing of = 8. The result for the van Dam -method using 2D polars has also been plotted, to demonstrate the effect of the sweep correction function Λ . Figure 5c presents the lift curve for an isolated wing with a geometry which is representative of a modern HTP. It is compared to CFD results from the Selena project [12]. The correlation between the present method and the CFD results are good, which suggests the validity of the sweep correction function Λ for aspect ratios as low as = 5. In the scope of aircraft design, the interest is generally in the analysis of wing-fuselage or tail-fuselage geometries, rather than isolated lifting surfaces. Our proposed method for accounting for the effect of the fuselage on the tailplane consists on modelling the fuselage as a cylindrical surface of vortex rings, removing all internal wing panels. The cylindrical surface must extend a length of at least ~3 ̅ in front and behind the lifting surface to avoid unphysical effects. The following results use the present method to replicate a CFD analysis of a body-tail geometry representative of a modern airliner. Figure 6a shows the VLM representation of the geometry. The spanwise loading is shown in Figure 6b. The present method correlated well to the CFD result. The 3-dimensional lift curve of the body-tail configuration presented in Figure 5d shows an acceptable correlation of the present method to wind tunnel tests and CFD computations. This indicates a strong potential of the method for analyzing complex configurations and interactions between lifting surfaces and bodies.

Conclusions and Way Forward
The present method has been implemented as a Python library with a simple user interface which enables the user to generate wing and fuselage geometries, read 2D input polars from XFOIL, MSES and Tau2D, and compute complete 3-dimensional lift curves. Typical computation times for an isolated wing are ~4 (~25 for a wing-fuselage configuration) on a single core. Given this low computational cost, the present method has great potential for use in the preliminary design process. The comparison cases shown in this report suggest that the method is able to predict the onset of stall with reasonable accuracy for a wide design space of isolated wings and wing-body configurations with proper sensitivity to the design parameters. Notably, the method is predominantly physically-based, scarcely relying on empirical data. Only the Λ function -which is indirectly a loose model for the evolution of the boundary layer thickness due to sweep-relies on empirical data. A full validation of the method is still pending. In order to achieve this, an extensive set of wind tunnel test results of a large design space of wing geometries are required. Moreover, the cases for which the 3D stall behavior cannot be derived by augmenting the 2D data must still be identified to find the limits of applicability of the present method.