OpenVSP's Airfoil Analysis capability and VSPAERO's strip-wise parasite drag model both use libNeuralFoil. libNeuralFoil is a C++ library implementation of Peter Sharpe's NeuralFoil. NeuralFoil is a neural network model of airfoil performance trained on a large set of XFoil runs. For details on NeuralFoil and XFoil, visit the respective sites linked above.

OpenVSP and VSPAERO both use the `xlarge` NeuralFoil model. NeuralFoil has ben trained on incompressible XFoil runs. Several compressibility corrections have been applied. These corrections are similar, but not identical to the corrections applied in Peter Sharpe's AeroSandbox.

## Airfoil Representation

NeuralFoil (and therefore libNeuralFoil) requires eighth order Class Shape Transformation (CST) curve coefficients for the upper and lower curves, as well as terms for the leading edge camber and trailing edge angle. NeuralFoil also requires the airfoils to be normalized in a particular manner and for the split between the upper and lower curves to be the most forward point of the normalized airfoil.

Internally, OpenVSP represents airfoils with piecewise cubic Bezier curves that have been fit to the definition of the airfoil (whether from NACA formulas, files, CST, etc.). Although OpenVSP provides an option for modeling airfoils with a CST parameterization, that definition does not perfectly match NeuralFoil's definition and even those airfoils are processed from the Bezier curves.

In OpenVSP's airfoil representation, the piecewise Bezier curve starts at the trailing edge (TE) of the lower surface and continues around the leading edge (LE) to the TE of the upper surface. When preparing an airfoil for NeuralFoil, the TE point is defined as the average of the first and last points of the airfoil curve. Next, OpenVSP searches for the airfoil point furthest from the TE; this point is designated as the LE. This definition of the LE may not coincide with OpenVSP or the airfoil's native definition of the LE.

The upper and lower curves are evaluated at 100 points with equal arc length spacing. This results in 199 points around the airfoil with the LE represented exactly and shared by both curves.

The airfoil is normalized (rotated and scaled) such that the LE is at (0, 0) and the TE is at (1, 0). The rotation angle and scale factor are recorded such that they can later be applied as corrections to the angle of attack and the chord as needed. Eighth order CST curves are fit to the points through a least squares fit.

When preparing to run VSPAERO, several other geometric parameters are calculated at the same time as the CST coefficients. These include the airfoil's thickness to chord ratio (t/c) and the three-dimensional leading and trailing edge lines of the lifting surfaces. All of this data is written to an `*.airfoilstack` file.

When VSPAERO is running in analysis mode, it will read this data from the file. However, when VSPAERO is running in optimization mode, it can access this information in memory via the API.

## Compressibility Corrections

In OpenVSP, when the compressible analysis mode is selected, the incompressible results from NeuralFoil are corrected for compressibility effects using the Prandtl-Glauert transformation, a buffet correction, and a transonic drag rise model. In VSPAERO, these models are applied whenever they are applicable.

### Sweep Analogy

For a swept wing, only the velocity component normal to the sweep line drives compressibility effects. The local Mach number for an infinite swept wing is:

M l o c=M∞cos⁡(Λ c/4)M_{loc} = M_\infty \cos(\Lambda_{c/4})

where Λ c/4\Lambda_{c/4} is the quarter-chord sweep angle. All subsequent compressibility corrections use M l o c M_{loc} rather than M∞M_\infty.

When using the compressible analysis mode in OpenVSP, the user specifies the freestream Mach number (M∞M_\infty) and the sweep angle (Λ c/4\Lambda_{c/4}). Conversely, VSPAERO calculates these quantities on the fly based on the three dimensional leading and trailing edges of the lifting surfaces, the flight condition, and any relative motion for the case. Consequently, the local sweep and Mach number in VSPAERO will account for differences in sweep caused by sideslip and changes due to rotating blades at any orientation.

### Prandtl-Glauert Lift and Moment Correction

The Prandtl-Glauert compressibility factor is:

β=1−M l o c 2\beta = \sqrt{1 - M_{loc}^2}

The compressible lift and moment coefficients are:

c l,c m p=c l,i n c β c m,c m p=c m,i n c β c_{l,cmp} = \frac{c_{l,inc}}{\beta} \qquad c_{m,cmp} = \frac{c_{m,inc}}{\beta}

### Compressible Pressure Coefficient

The incompressible pressure coefficient is recovered from the NeuralFoil boundary layer edge velocity ratio:

C p,i n c=1−(u e V∞)2 C_{p,inc} = 1 - \left(\frac{u_e}{V_\infty}\right)^2

The Prandtl-Glauert corrected pressure coefficient is:

C p,c m p=C p,i n c β C_{p,cmp} = \frac{C_{p,inc}}{\beta}

### Skin Friction Drag Correction

Compressibility reduces turbulent skin friction relative to the incompressible value at the same Reynolds number. A simplified correction factor based on the stagnation-to-static temperature ratio for isentropic flow is:

T 0 T=1+γ−1 2 M 2\frac{T_0}{T} = 1 + \frac{\gamma - 1}{2}M^2

Applied for air (γ=1.4\gamma=1.4), where the square root captures the effect of temperature in the skin friction coefficient equation.

F c=T 0 T=1+0.2 M l o c 2 F_c = \sqrt{\frac{T_0}{T}} = \sqrt{1 + 0.2\, M_{loc}^2}

The compressible drag coefficient is:

c d,c m p=c d,i n c F c c_{d,cmp} = \frac{c_{d,inc}}{F_c}

## Critical Mach Number

The critical Mach number M c r i t M_{crit} is the freestream Mach number at which the local flow first reaches sonic conditions at the point of minimum pressure on the airfoil. It is calculated from the minimum incompressible pressure coefficient C p,m i n C_{p,min} using a curve fit developed by Peter Sharpe for use with NeuralFoil and AeroSandbox:

M c r i t=(1.0083619−C p,m i n+(−0.51058894 C p,m i n)0.6553655)−0.5536965 M_{crit} = \left(1.0083619 - C_{p,min} + \left(-0.51058894\, C_{p,min}\right)^{0.6553655}\right)^{-0.5536965}

The minimum incompressible Cp is found by scanning the upper and lower surface edge velocity distributions returned by NeuralFoil.

## Drag Divergence Mach Number

The drag divergence Mach number M d d M_{dd} is defined as the Mach number at which the drag rise rate reaches d c d/d M=0.1 dc_d/dM = 0.1. Using Locke's quartic drag rise model described below, this condition yields:

M d d=M c r i t+(0.1 80)1/3 M_{dd} = M_{crit} + \left(\frac{0.1}{80}\right)^{1/3}

The increment (0.1/80)1/3≈0.108\left(0.1/80\right)^{1/3} \approx 0.108 is obtained by differentiating c d,w a v e=20(M−M c r i t)4 c_{d,wave} = 20\,(M - M_{crit})^4 and setting the result to 0.1.

## Wave Drag Model

Wave drag c d,w a v e c_{d,wave} is modeled in three regions:

Subcritical (M l o c≤M c r i t M_{loc} \leq M_{crit}):

c d,w a v e=0 c_{d,wave} = 0

Transonic rise (M c r i t<M l o c≤0.95 M_{crit} < M_{loc} \leq 0.95):

c d,w a v e=20(M l o c−M c r i t)4 c_{d,wave} = 20\,(M_{loc} - M_{crit})^4

This quartic form matches the classical drag rise behavior near the critical Mach number.

High subsonic / transonic / supersonic (M l o c>0.95 M_{loc} > 0.95):

A cubic Bezier curve blends from the quartic drag rise at M=0.95 M = 0.95 to the thin-airfoil-theory supersonic value. The supersonic anchor point is taken at M=2.0 M = 2.0, where the linearized supersonic drag is:

c d,w a v e=4(t/c)2 M 2−1|M=2=4(t/c)2 3 c_{d,wave} = \frac{4\,(t/c)^2}{\sqrt{M^2 - 1}}\bigg|_{M=2} = \frac{4\,(t/c)^2}{\sqrt{3}}

The Bezier blend ensures continuity of value and slope at M=0.95 M = 0.95, with zero slope imposed at the supersonic anchor to provide a smooth transition.

The total drag coefficient is the sum of the viscous and wave components:

c d,t o t=c d,c m p+c d,w a v e c_{d,tot} = c_{d,cmp} + c_{d,wave}

## Transonic Buffet Correction

Above the drag divergence Mach number, shock-induced separation leads to buffet onset and a loss of lift effectiveness. A correction is applied to c l c_l and c m c_m for M l o c>M d d+0.06 M_{loc} > M_{dd} + 0.06:

BFact1=M l o c−M d d−0.06 1.0−M d d−0.06\text{BFact1} = \frac{M_{loc} - M_{dd} - 0.06}{1.0 - M_{dd} - 0.06}

BFact2=(1−BFact1)+0.3 BFact1\text{BFact2} = (1 - \text{BFact1}) + 0.3\, \text{BFact1}

c l,b u f=c l,c m p⋅BFact2 c m,b u f=c m,c m p⋅BFact2 c_{l,buf} = c_{l,cmp} \cdot \text{BFact2} \qquad c_{m,buf} = c_{m,cmp} \cdot \text{BFact2}

This linearly reduces effective lift and moment from their full compressible values at M d d+0.06 M_{dd} + 0.06 to 30% of those values as M l o c M_{loc} approaches 1. The threshold of 0.06 above M d d M_{dd} provides a small margin before buffet onset degrades performance.

Drag is not modified by the buffet factor; the viscous and wave drag components are reported independently.