Introduction

What are Panel Methods? Panel methods are numerical schemes for solving the Prandtl-Glauert equation, which describes the behaviour of linear, inviscid, and irrotational flows around aircraft at subsonic or supersonic speeds. Essentially, it’s a simplified way to model how air flows around an aeroplane without considering the complex effects of viscosity (friction) or strong shock waves. The fundamental idea behind panel methods is to divide the aircraft surface into a collection of small, flat surfaces called panels. Each panel is then assigned a distribution of sources, doublets, or vortices, which are mathematical representations of simple flow patterns. These flow patterns are fundamental solutions to the Prandtl-Glauert equation, ensuring that the overall flow field satisfies the governing equation. The general process of a panel method can be broken down as follows:

  • Discretisation: The aircraft surface is divided into a network of panels.

  • Singularity Distributions: Each panel is assigned a distribution of sources, doublets, or vortices.

  • Influence Coefficients: Mathematical expressions are derived to calculate how the flow from each panel influences the flow at every other panel.

  • Boundary Conditions: Boundary conditions are imposed at specific points on the panels to ensure the flow adheres to the desired physical behaviour, such as no flow penetrating the aircraft surface.

  • Solving the System: The influence coefficients and boundary conditions form a system of equations that are solved to determine the strengths of the source, doublet, or vortex distributions on each panel.

  • Flow Field Calculation: Once the singularity strengths are known, the flow velocity and pressure can be calculated at any point in the flow field.

  • Force and Moment Calculation: By integrating the pressure distribution over the aircraft surface, the overall aerodynamic forces and moments can be determined.

Panel methods are highly efficient for analyzing complex aircraft configurations, especially compared to more computationally intensive methods like Navier-Stokes solvers. They find wide use in the aerospace industry for preliminary design studies, evaluating aerodynamic loads, and understanding flow behaviour around intricate geometries.

Limitations of panel methods stem from the simplifications inherent in the Prandtl-Glauert equation:

  • Neglect of viscosity: Panel methods typically don’t account for viscous effects like skin friction drag and flow separation, making them less accurate for predicting drag and stall characteristics. However, some modern panel codes like FlightStream incorporate viscous-coupling equations and integral boundary layer models to address this limitation. These methods calculate the boundary layer development and adjust the flow field accordingly.

  • Linearity assumption: The Prandtl-Glauert equation is linear, which means that panel methods struggle to accurately predict flow behaviour in scenarios with strong shock waves, such as transonic flight. However, researchers have developed hybrid methods like TRANAIR that combine panel methods with techniques like finite element methods to tackle nonlinear transonic flows. Despite their limitations, panel methods remain valuable tools for aerodynamic analysis, offering a balance between accuracy and computational efficiency.

Why FlightStream?

FlightStream is a modernized panel-method based aerodynamic solver. The original inviscid theory is enhanced with several novel features to provide a more realistic and comprehensive aerodynamic analysis tool compared to legacy panel codes. These enhancements include (but are not limited to):

  • Viscous-Coupling Equations: While traditional panel codes primarily focus on inviscid flow, FlightStream incorporates viscous-coupling equations to better model real-world aerodynamic behaviour.

  • Compressibility Models: FlightStream goes beyond incompressible flow analysis by including compressible flow models, enhancing its capability to simulate a wider range of flight conditions.

  • Advanced Boundary Conditions: Unlike earlier panel codes, FlightStream introduces advanced boundary conditions, including actuator discs for simulating propellers and jet exhausts.

  • Time-Accurate Capabilities: FlightStream boasts time-accurate simulation capabilities, enabling analysis of dynamic scenarios such as trajectory simulations and prescribed motions, a feature lacking in many legacy panel codes.

  • Relaxed Vortex-Wake Filament Method: Unlike simpler wake models, FlightStream uses a relaxed vortex-wake filament method to realistically simulate wake dynamics, including vortex stretching and disintegration. This allows for more accurate representation of wake interaction with downstream components.

  • Efficient Algorithms and Modern GUI: FlightStream leverages computationally efficient algorithms, such as the Fast Multipole Method (FMM), and a user-friendly graphical user interface (GUI) for enhanced usability.

  • Surface Vorticity: FlightStream is based on surface vorticity, which provides a more physically accurate representation of flow separation and stall characteristics.

  • Integral Boundary Layer: FlightStream includes an integral boundary layer model that accounts for the effects of viscosity on the flow. This allows for more accurate predictions of aerodynamic loads, especially at high angles

  • Unsteady Flow Solver: FlightStream has an unsteady flow solver that can be used to simulate the aerodynamics of complex, time-varying maneuvers.

The combination of these features makes FlightStream a powerful tool for modern vehicle design and analysis while ensuring ease of use, high accuracy, and computational efficiency.

Formulation

Numerical Discretization

This section talks about the order of operations in the numerical procedure.

Geometry

How is the geometry discretized? How are unstructured meshes handled?

Matrix Solver

Talk about fast multipole method (FMM) and the octree.

Subsonic Solver

The potential flow \(\phi\) equation is used to define an inviscid, irrotational, and incompressible, and flow. The potential flow Laplace equation defines the smooth behavior of the flow field:

\[ \nabla^2 \phi = 0\]
immersed boundary
Figure 1. Illustration of a body, \(S_B\), immersed in a potential flow field

For a body, \(S_B\), immersed in a potential flow field the velocity, Green’s theorem can be applied to the inner and outer surfaces to obtain the following relation:

\[\Phi_P = -\frac{1}{4\pi} \iint_{S_B + S_w} {\mu} \mathbf{n} \cdot \nabla\left(\frac{1}{r}\right) dS + \frac{1}{4\pi} \iint_{S_B} \frac{1}{r} \mathbf{n} \cdot (-\sigma) dS + \phi_\infty\]

Where the doublet and source primitives are defined as, \(\mu = \phi - \phi_i\) and \(\sigma = \frac{\partial{\phi}}{\partial{n}} - \frac{\partial{\phi_i}}{\partial{n}}\), respectively. Thus, the potential at any point in the flow field may be determined as the sum of source and doublet influences integrated across each panel. The gradient of this may also be taken to find the velocity induced at any point in the flow field.

In addition to approximating the surface by a discrete set of panels, the source and doublet strengths are assumed to follow some analytic distribution over each panel (i.e. constant, linearly-varying, etc.). Most often, panel methods incorporate a linear distribution of doublet strength and a constant distribution of source strength across each panel, or quadratic and linear distributions of source and doublet strength, respectively. This is because the source strength represents a jump in velocity across a panel, while the gradient of doublet strength represents a jump in velocity. Some panel methods instead use constant distributions for both sources and doublets; however, this is typically only done for subsonic methods. For a given body and freestream condition, the values of \(\mu\) and \(\sigma\) on each panel are not known a priori. These distributions may be written in terms of some finite set of unknown source and doublet parameters. By enforcing as many boundary conditions as there are unknown source and doublet parameters, a linear system of equations may be formed. This system of equations will be linear because the induced potentials and velocities are linear in the source and doublet strengths. Then, by solving this system of equations, the unknown source and doublet parameters may be determined. Once these are known, the flow properties at any point on the body or in the flow field are easily obtained.

Surface Blowing

A solution to the above requires that the Neumann boundary condition be satisfied on the surface:

\[ n \cdot \nabla \phi = -V_N - n \cdot V_S \; \; on \; S_B\]

Where \(V_N\) is the normal velocity on the boundary which is 0 in the absence of transpiration (later it will be shown that a non-zero transpiration velocity allows for the capture of some viscous effects on the surface). \(V_S\) is the local velocity on the boundary and may be composed of several parts due to body rotation (indeed this is how rotational flow fields are modeled in FlightStream). In the simplest steady case the above equation reduces to:

\[ n \cdot \nabla \phi = 0\]

Given a displacement thickness, \(\delta^*\), the boundary condition becomes:

\[ V_N = \frac{\partial}{\partial s} (U_e \delta^*) + V_{norm}\]

where the first term represents the rate of growth of the boundary layer, and the second term, \(V_{norm}\), is positive for outflow and negative for inflow. The boundary layer term is typically zero at the start of the calculation and is then updated during each viscous/potential iteration cycle, as described in Section Viscous Modeling.

Wake Treatment

FlightStream employs a cutting-edge relaxed vortex-wake filament method to capture the effects of the wake. At each trailing edge of the body, a wake filament is released. This has the capability to resolve around solid bodies automatically and has inbuilt vortex stretching and disintegration models. The latter enables the modeling of full volumetric effects and to study the interactions between various components while maintaining solution run times measured in minutes. Traditional free wake panel methods have issues with singularities when wakes intersect solid bodies.

The unstructured surface mesh drives the wake to also be unstructured since the trailing edges (TE) are not directionally mapped. In FlightStreams, the concept of “wake strands” (as opposed to wake sheets) emenating from mapped vertices determine the wake geometry.

Once the trailing edges are marked, the vertex pool corresponding to these edges is available from the global vertex map. These vertices become the starting point, or nodes, for the individual wake strands. Since no additional vorticity sources or sinks are added by the wake, the net vorticity of the geometry remains conserved and the wake propagates only the required vorticity to satisfy the Kutta condition at the trailing edges.

wake strands
Figure 2. The directional vorticity being shed into the wake strand at the trailing edge on a) a general node and b) a corner node.

Vorticity strength reductions downstream of the trailing edge can only occur due to viscous effects. This is modeled using the application of the Lamb–Oseen vortex model obtained from the exact solution of the Navier–Stokes equations for a laminar 2D vortex. Knowing the initial vorticity shed into the strand at the trailing-edge node, the Lamb–Oseen model provides the vorticity decay equation as:

\[ \Gamma_t = \Gamma \left( 1 - e^{-\frac{\rho_\infty r^2}{4\mu t}} \right)\]

In the above equation, the time used corresponds to the solver pseudotime and starts at zero at the trailing edge. It is marched alongside the strand (using values of local velocity and the discretization size of the relaxed wake, which is fixed to the average mesh edge length) until the termination of the wake downstream at the Trefftz plane. All strands are projected or terminated upon intersection with the Trefftz plane.

Integrated Circulation

The method of integrated circulation works by evaluating the vorticity along a series of two-dimensional cross sections encompassing the entire geometry. These cross sections can be aligned with their planar vectors directed orthogonally to the freestream and load vectors. The arbitrary distribution of vorticity on the unstructured surface mesh can therefore be converted into directional vorticity along the required load vectors. An example of such cross sections is shown in Figure 3.

The vorticity for each of the cross sections is evaluated by the application of the generalized circulation equation along the perimeter of the cross section. This is easily done once the vorticity on the surface mesh has been evaluated in the solver. Since we can evaluate the velocity induced by a vortex ring at any point in space around the body using Eq.[eq:vortex_ring_velocity]), we discretize the cross section along its rectangular perimeter (Fig.[fig:cross_section_discretization]). We then evaluate the induced velocity from all of the surface vortex rings on each of the discretized vertices of the cross sections. Once a distribution of velocity along the vertices of the cross section is known, the net, integrated circulation of that cross section is evaluated using the numerical formulation of Stokes' theorem:

\[ \Gamma_k = \int_0^L V_{induced} \cdot dl\]
integrated circulation
Figure 3. Representations of a) integrated circulation cross-sections on a sample geometry and b) equivalent Prandtl lifting line

Treated in isolation in two-dimensional space, defined by the plane of the arbitrary cross section, the aerodynamic loads from this integrated circulation could be evaluated from the Kutta–Joukowski theorem [1]. A series of these cross sections enclosing the arbitrary three-dimensional body Cross sections a) allows us to convert an unstructured distribution of bound surface vorticity into adjacent sections of two-dimensional net circulations, thereby reducing them to the Prandtl lifting-line theory Cross sections b), with all the inherent benefits in load computations as provided by this formulation.

If the cross-section planes are aligned orthogonally to the freestream and lift vectors, then their alignment is identical to a lifting-line distribution of integrated vorticity along a rectangular wing. Consequently, the integrated vorticity of each cross section is identical to the "bound" vorticity from the Prandtl theory. The downwash between each of these cross sections is then evaluated by treating each cross-section bound vorticity as shedding a horseshoe-shaped vortex aligned with its planar and freestream vectors.

Furthermore, spanwise loading distributions of arbitrary bodies can also be evaluated through this approach by plotting the loads evaluated for each cross section and plotting them in spatial locations on the abscissas. This entire formulation therefore reduces an arbitrary surface distribution of vorticity to the form of the Prandtl theory, and therefore yields the known formulations of induced lift, drag, and downwash velocities:

\[ L_{i} = \int_{-b/2}^{b/2} \rho_\infty V_\infty \Gamma_y \, dy\]
\[ D_{i} = \int_{-b/2}^{b/2} \rho_\infty w_y \Gamma_y \, dy\]
\[ w_y = \frac{1}{4\pi} \int_{-b/2}^{b/2} \left(\frac{1}{y-y_0}\right) \frac{-d\Gamma_{y_0}}{dy} \, dy_0\]

For the sake of reading convenience, each of these cross sections will hereafter be referred to as an integrated circulation loop (ICL). In numerical application, the ICLs are discretized using the average edge size of the underlying surface mesh. The spacing between ICLs is also discretized using the same reference values. It is possible to include other discretization values, and this can have a varying effect on the solution fidelity, depending on the complexity of the geometry in question.

Compressibility Correction

In the Subsonic solver, the compressibility correction is applied via the standard Prandtl-Glauert rule. This is done by scaling the freestream velocity and density by the factor \(1/\sqrt{1-M^2}\), where \(M\) is the Mach number. This correction is applied to the freestream conditions and the body, and is used to correct the flow field and forces/moments. Typically this correction is only valid up to a Mach number of about 0.6. The user can easily enable the Incompressible solver to bypass this correction. If a higher fidelity compressibility analysis is required, then the user may use on of the several available transonic/low-order supersonic/high-order supersonic/hypersonic methods.

Mesh Boundary Conditions

In FlightStream, the boundary conditions are critical for the accurate simulation of aerodynamic flows. The following sections provide an overview of the boundary conditions implemented in FlightStream.

Trailing Edges

The trailing edge boundary condition is crucial for the accurate simulation of aerodynamic flows in panel methods, such as FlightStream. Trailing edge boundaries are essential for satisfying the Kutta condition, a fundamental principle in aerodynamics dictating that the flow must leave the trailing edge of an airfoil smoothly. Failure to meet this condition leads to unrealistic pressure distributions and inaccurate lift predictions. In FlightStream, trailing edge boundaries are marked as special edges on the surface mesh, where wake vortex filaments originate. The strength of these filaments, representing the shed vorticity, is determined by the difference in flow velocity above and below the wake.

FlightStream offers two types of trailing edge boundaries: Standard-Kutta and Relaxed-Kutta. In the Standard-Kutta type, the net circulation at the trailing edge is fully shed into the wake. Conversely, in the Relaxed-Kutta type, only half of the net circulation is shed, suitable for situations where the geometry of the trailing edge is not sharp, such as a rounded fuselage. The accurate placement of trailing edges is paramount for the correctness of the simulation. FlightStream provides several ways to mark trailing edge boundaries, including automatic detection, manual selection, defining perimeters around base regions, and importing from external files.

Wake Termination Nodes

Wake termination nodes are used in conjunction with trailing edges to control the behaviour of the wake vortex filaments. These nodes mark locations where the wake filament should be clipped, preventing the formation of unrealistic tip vortices in specific situations, such as at the junction of a wing and fuselage.

Base Pressure Region

Base pressure regions are used to model the surface-pressure effects of regions of bulk separated flow. In FlightStream®, these regions are explicitly defined by the user. A uniform \(C_{p,base}\) is applied for the base region during solve time. This value can be user-defined or an empirical model can be used.

FlightStream implements a novel method to calculate the base pressure as a function of freestream Mach number. This model is based on a regression of various experimental data. The base pressure coefficient monotonically increases with decreasing Mach number in the subsonic regime and with increasing Mach number in the supersonic regime. Outside the transonic regime, there is relatively little variation in the base pressure coefficient across various configurations. However, in the transonic regime the base pressure coefficient appears to be a function of more than the freestream Mach number, Base pressure model implemented in FlightStream shows the data is not collapsing for these Mach numbers. For supersonic flows, the base pressure coefficient is bounded from below by the vacuum pressure coefficient (solid black line in the figure).

base pressure
Figure 4. Base pressure model implemented in FlightStream

To give an approximate fit to the data shown, the following base pressure model is proposed. For subsonic freestream Mach numbers, the pressure coefficient is given by 1.

\[\begin{equation} C_{p_{base}} = -0.2 - \frac{0.22}{50(M_\infty - 1)^2 + 1} \end{equation}\]

For supersonic freestream Mach numbers, the pressure coefficient is given by 2.

\[\begin{equation} C_{p_{base}} = -\frac{0.42}{(M_\infty - 1)^2 + 1} \end{equation}\]

These equations both predict a base pressure coefficient of -0.42 at \(M_\infty = 1\). This model is shown in Base pressure model implemented in FlightStream as the solid red line. Note that this supersonic model can predict pressure coefficients lower than the vacuum pressure for high freestream Mach numbers. Thus, in practice, the supersonic base pressure coefficient is the greater between the vacuum pressure coefficient and that predicted by 2.

While this model yields reasonable results the user is encouraged to run experiments or consult literature for \(C_{p,base}\) values for their specific application.

Turbulent Trip Edges

In many situations, boundary layer transition is deliberately induced. On production vehicles, imperfections or intentional surface features can act as trip points. In wind tunnel tests, researchers often use physical trip strips to force transition at a desired location, ensuring consistent flow behaviour for measurements. Turbulent trip edges in FlightStream are a boundary condition used to trigger artificial boundary layer transition from laminar to turbulent flow. This is important because, in real-world scenarios and wind tunnel testing, various factors can induce early transition, impacting the overall flow behaviour and aerodynamic characteristics. Once trip edges are defined, FlightStream’s viscous solver incorporates them into the boundary layer calculations. The primary impact is the shift in boundary layer characteristics, influencing parameters like skin friction, displacement thickness, and momentum thickness, ultimately affecting the predicted lift, drag, and moment values.

Actuator Discs

Actuator disc models are instrumental in representing steady-state effects of inherently unsteady systems like propellers and jet exhausts in aerodynamics. These models simulate the impact of momentum additions in flow fields due to propulsion components, offering insights into how devices like propellers and exhaust jets influence surrounding airflow. The theory behind actuator discs provides an efficient way to estimate the aerodynamic effects of propulsion systems.

Types of Actuator Discs Actuator discs generally simulate two primary propulsion effects:

  • Propeller Actuator Disc (Conway Model): This model represents the effects of a propeller operating under steady-state conditions, often used in cruise-condition simulations to study how propellers modify the overall flow field.

  • Exhaust Jet Actuator Disc (Reichardt Model): Designed to simulate steady-state jet exhaust behavior, this model allows for examining the momentum effects of jet propulsion systems as they interact with the surrounding air.

Each type of actuator disc applies different boundary conditions and momentum additions to reflect the specific characteristics of the propulsion system.

Actuator discs possess several modifiable physical parameters:

  • Dimensions and Orientation: The radius or elliptical dimensions of the actuator disc can be adjusted, allowing for alignment with the intended propulsion component’s geometry.

  • Velocity Distributions: Actuator discs enable custom-defined velocity distributions across their surfaces. For instance, propeller discs may include swirl velocity to account for rotational flow effects in the slipstream, while exhaust jet discs can specify exhaust velocity and density.

  • Custom Thrust Profiles: Propeller discs can integrate custom thrust profiles, offering a more accurate representation of the non-uniform slipstream velocities, especially important for advanced or distributed electric propulsion systems.

Jet exhaust actuator discs feature parameters tailored to capture jet dynamics accurately:

  • Exhaust Velocity and Density: Define the jet flow’s initial conditions, providing a basis for momentum addition within the flow field.

  • Jet Spreading Rate: Governs the diffusion rate of jet momentum downstream, simulating jet mixing and its influence on the flow field.

Although actuator disc models generally offer steady-state approximations, unsteady simulations can provide deeper insights in dynamic scenarios, such as varying propeller pitch or RPM. In summary, actuator disc models are crucial for simulating propulsion effects on surrounding flow fields, enabling aerodynamic evaluations of various propulsion systems in aircraft design. Understanding the principles and parameters associated with these models allows for precise and adaptable simulation setups tailored to different aerodynamic and propulsion studies.

Thrust Distribution Examples

The following examples show the thrust distributions for various propeller blade shapes. The presented equations can be used to create a thrust profile and then applied to a Custom actuator disc in FlightStream.

Generic Thrust Distribution

Let the thrust distribution over a single blade of a propeller be given by:

\[T(r) = -A(r - 1)(r - R_h)(r + 1 - 2R_h)(r + R_h) \tag{1}\]

where \(r\) is the non-dimensional blade location (varying from 0 to 1), \(R_h\) is the non-dimensional location of the hub, and \(A\) is a scaling coefficient. This is shown in Thrust distribution for a generic propeller blade.

generic thrust
Figure 5. Thrust distribution for a generic propeller blade

This can be rearranged into the form:

\[T(r) = -A(r^4 + c_3r^3 + c_2r^2 + c_1r + c_0) \tag{2}\]

where:

\[c_3 = -2R_h \\ c_2 = -R_h^2 + 2R_h - 1 \\ c_1 = 2R_h^3 \\ c_0 = R_h^2 - 2R_h^3\]

If the total thrust, \(T\), produced by the propeller is known, then \(A\) may be determined from:

\[T = N_bR_t \int_{R_h}^{1} T(r)dr \\ T = -N_bR_tA \left[ \frac{1}{5}(1 - R_h^5) + \frac{c_3}{4}(1 - R_h^4) + \frac{c_2}{3}(1 - R_h^3) + \frac{c_1}{2}(1 - R_h^2) + c_0(1 - R_h) \right]\]
\[A = -\frac{T}{N_bR_t} \left[ \frac{1}{5}(1 - R_h^5) + \frac{c_3}{4}(1 - R_h^4) + \frac{c_2}{3}(1 - R_h^3) + \frac{c_1}{2}(1 - R_h^2) + c_0(1 - R_h) \right] \tag{3}\]

where \(N_b\) is the number of blades and \(R_t\) is the tip radius (dimensional). \(T\) and \(A\) may be either dimensional or non-dimensional coefficients, so long as they are consistent with each other.

Elliptic Thrust Distribution

The thrust distribution over a single blade of a propeller may also be given by:

\[T(r) = A \sqrt{1 - r^2} \tag{4}\]

where \(A\) is a scaling constant and \(r\) is the non-dimensional blade location. We call this the elliptic distribution.

Again, \(A\) may be determined from:

\[T = N_bR_t \int_{0}^{1} T(r)dr \\ T = N_bR_t \int_{0}^{1} A\sqrt{1 - r^2}dr\]

Note that for this distribution, we assume no hub. If we let:

\[\sqrt{1 - r^2} = \cos \theta\]

then:

\[T = N_bR_tA \int_{0}^{\pi/2} \cos^2\theta d\theta \\ T = N_bR_tA \cdot \frac{1}{2} \left( \theta + \sin(\theta)\cos(\theta) \right) \Big|_{0}^{\pi/2} \\ T = N_bR_tA \cdot \frac{1}{2} \left( \frac{\pi}{2} \right)\]

Thus:

\[A = \frac{4T}{\pi N_bR_t} \tag{5}\]

Viscous Modeling

Given the inviscid solution, FlightStream® estimates the boundary layer properties such as thickness, \(\delta\), momentum thickness, \(\theta\), and shape factor \(H\) for both laminar and turbulent flows. Additionally, the viscous effects can be used to virtually perturb the boundaries and modify the inviscid bulk flow.

Laminar / Turbulent Boundary Layer

Two models for the boundary layer are implemented in FlightStream®: laminar and turbulent. A boundary layer transition model is also implemented. All of these models are two-dimensional methods implemented along the on-body surface streamlines. These integral methods are computationally efficient and can be applied to any general three-dimensional wall boundary with exceptions of regions involving three-dimensional cross flow. The laminar boundary layer analysis method used in FlightStream® is the standard two-parameter model of Thwaites [2] integral method with the momentum integral equation written as follows:

\[U \frac{{d}}{{d \eta}} {\left[ \frac{{ \theta^2}}{{ \nu}} \right]} = 0.45 - 6 \frac{{ \theta^2}}{{ \nu}} \frac{{dU}}{{d \eta}}\]

An integral boundary layer model for compressible turbulent flows has been implemented into the inviscid flow solver. The Standen model is used based on the model’s robustness and versatility for the range of subsonic Mach numbers of interest. The final turbulent boundary layer model implemented has numerical improvements and optimizations over that of the original model developed by Standen , but philosophically follows his work quite closely, with the application to subsonic, turbulent, compressible flows along on-body streamlines. The method is essentially two-dimensional (similar to the laminar model described previously), but is implemented in a quasi-three-dimensional manner inside FlightStream® by assuming that the three-dimensional cross flow terms are not dominant. For the fluid properties outside the boundary layer, the flow is assumed to be isentropic and subsonic, but compressible. The primary differential equations are developed for the turbulent compressible momentum thickness (\(\theta_c\)), and the compressible shape factor (\(H_c\)), as described here:

\[\frac{d\theta_i}{dX} = -\frac{\theta_i}{M_e} \frac{dM_e}{dX} \left( 2 + H_{tr} \right) + \frac{dx}{dX} \frac{c_f}{2} \left(\frac{T_e}{T_0}\right)^3\]
\[\frac{dH_i}{dx} = -\frac{\left(H_i - 0.7\right)^{3.715}}{4.17} \left[\frac{F}{\theta_i} \frac{dX}{dx} - \frac{H_{\delta-\delta^*}}{M_e} \frac{dM_e}{dx} - \frac{H_{\delta-\delta^*}}{\theta_i} \frac{d\theta_i}{dx}\right]\]

The auxiliary equations needed to complete the above equations are :

\[M_e = \frac{U_e}{\sqrt{\gamma g R T_e}}\]
\[\frac{dM_e}{dx} = \frac{1}{a_e \left[ 1 + \frac{U_e (\gamma - 1) T_e^2 M_e}{2 T_e a_e T_0} \right]} \frac{dU_e}{dx}\]
\[F = 0.03268 \left( H_{\delta - \delta^*} \right)^{-1.3} + 0.006\]
\[H_{\delta - \delta^*} = 1.5359 \left( H_i - 0.7 \right)^{-2.715} + 3.4\]
\[\frac{dX}{dx} = \frac{T_e}{T_R} \left(\frac{T_e}{T_0}\right)^3 \left(\frac{\mu_R}{\mu_0}\right)^{0.268}\]

Solving these equations numerically (using 4th order Runge-Kutta methods) along the on-body streamline generates all relevant turbulent boundary layer data. This data is then used to compute the boundary layer thicknesses. The skin friction equation used is:

\[\frac{c_f}{2} = \left[0.123 e^{-1.56 H_i} \left(U_e \theta_i \frac{\rho_0}{\mu_0}\right)^{-0.268}\right] \left[\left(\frac{T_e}{T_R} \left(\frac{T_R}{T_0}\right)^{0.402}\right) \left(\frac{T_0 + 198}{T_R + 198}\right)^{0.268}\right]\]

In order to account for surface roughness a modification from Olsen is applied to the surface roughness. In his work, from experiments on rough pipes, it is shown that the overlap layer in the turbulent velocity profile maintains a logarithmic shape and that the intercept at zero velocity is shifted towards higher roughness height. This shift depends on the size of the roughness. For sufficiently large roughness heights, the shift can be shown to have a logarithmic relation with roughness height. As the outer part of the velocity profile is unaffected by the roughness, Coles law of the wake is still applicable for flow over rough surfaces. By including the shift in the intercept in Coles law of the wake as a function of roughness height, it is possible to derive the following expression for the skin friction coefficient (\(C_f\)) for rough surfaces in the same manner as for the smooth case:

\[C_f = \frac{0.90 e^{-2.40H}}{\left( \log_{10} Re_{\theta} - \log_{10} Re_k + 1.11 \right)^{2.45 - 0.15H}}\]

In this equation, \(R_\theta\) is the Reynolds number based on the local momentum thickness of the boundary layer (\(\theta\)), \(Re_k\) is the Reynolds number based on the surface roughness height (\(k\)) and \(H\) is the local shape factor of the boundary layer.

Boundary Layer Transition

The presence of both laminar and turbulent separation models in FlightStream® allow for the possibility to model complex transitional separation physics along a given surface streamline. A suitable transition model is required to determine the location where a switch between the laminar and turbulent separation models can be performed. The transition model created by Dvorak et.al. is applied in FlightStream®, based on the criteria of robustness and general applicability. The general outline of this model is detailed in this section.

As the local Reynolds number based on momentum thickness gets larger, the laminar boundary layer becomes unstable and small disturbances in the boundary layer begin to amplify rather than damp out. The amplification of disturbances eventually leads to transition to a turbulent boundary layer. The local Reynolds number at which the laminar boundary layer becomes unstable can be correlated with the local pressure gradient by means of the following empirical relationships. The local pressure gradient along the laminar boundary layer on the streamline can be defined as:

\[K = \frac{\theta^2}{\nu} \frac{dU}{d\eta}\]

Here, \(\theta\) is the momentum thickness of the laminar boundary layer; \(\nu\) is the kinematic viscosity, \(U\) is the velocity at the outer edge of the boundary layer (potential flow surface velocity in FlightStream®); and \(\eta\) is the coordinate along the surface streamline. The pressure gradient parameter \(K\) is computed as a field variable along the surface streamline for each mesh face location. The work done by Dvorak et.al. has resulted in the following set of empirical equations for \(K\) for a range of local Reynolds Number based on the momentum thickness (\(R_\theta\)):

\[K = -0.4709 + 0.11066 \ln(Re_\theta) - 0.0058591 (\ln(Re_\theta))^2 \quad (0 \leq Re_\theta \leq 650)\]
\[K = 0.69412 - 0.23992 \ln(Re_\theta) + 0.0205 (\ln(Re_\theta))^2 \quad (650 < Re_\theta \leq 10000)\]

If \(K\), as computed by the original formulation using the computed boundary layer data is greater than \(K\) as computed by these empirical equations, then the laminar boundary layer is marked as unstable. In the unstable region, an average pressure gradient parameter \(\bar{K}\) can be defined as:

\[\bar{K} = \frac{\int_{\eta_{u}}^{\eta} K d\eta}{\eta - \eta_{u}}\]

\(\eta_{u}\) is the coordinate along the surface streamline where the laminar boundary layer becomes unstable. The local Reynolds number at the transition point can be correlated with \(K\) by means of the following empirical relationships using \(R_\theta\):

\[\bar{K} = -0.0925 + 0.00007(Re_\theta) \quad (0 \leq Re_\theta \leq 750)\]
\[\bar{K} = -0.12571 + 0.000114286(Re_\theta) \quad (750 \leq Re_\theta \leq 1100)\]
\[\bar{K} = 1.59381 - 0.45543 \ln(Re_\theta) + 0.032534(\ln(Re_\theta))^2 \quad (1100 \leq Re_\theta \leq 3000)\]

When \(K_u\), as computed computation is greater than \(\bar{K}\) as computed by the empirical equations, natural transition is predicted. The transition from laminar to turbulent boundary layer is assumed to take place instantaneously at the transition point. At this point, the solver switches the separation models from laminar to turbulent. The user is provided an option to toggle between completely turbulent separation models (no laminar transition) or a transitional separation model . An additional field variable for the boundary layer transition model has been implemented into the FlightStream® post-processor to visualize these results .

Performance Enhancements

With the addition of advanced viscous models to FlightStream®, additional computational enhancements have been added to the solver to maintain an \(O(nlog(n))\) scalability to the performance of the boundary layer algorithms. To this end, spatial octree algorithms that mimic the existing Fast Multipole Method implementation (for the inviscid solver) have been implemented for the viscous flow near-wall regions. This allows extremely fast near-field sorting of points closest-to, or inside, a viscous flow region. The Viscous Spatial Tree (VST) is created in the global coordinate system of the solver and increases the memory footprint of the FlightStream® solver by <5% for all cases. However, the performance benefits for viscous modules in FlightStream® are substantial with an otherwise \(O(n \times m)\) algorithm (where \(m\) is a very large number of surface mesh faces) being reduced to \(O(n \cdot log(n))\).

Further performance benefits are obtained by specifying proximity-search boxes in the local coordinate system of the mesh faces being checked for viscous flow. These local-frame proximity boxes are computed at the time of solver initialization and are therefore available without performance penalties for all subsequent analysis work. The localized proximity-search boxes are invoked as a subset within the VST. The height of the proximity search box is the displacement thickness of the local viscous boundary layer flow.

Attachment Lines

FlightStream utilizes a robust model for computing attachment lines and critical points on the inviscid surface-velocity field. This is an enabling requirement for robust coupling with the viscous boundary layer models described previously. It should be noted that the inviscid surface streamlines in FlightStream are computed by solving the Ordinary Differential Equations (ODEs) for a three-dimensional streamline, constrained to the vehicle surface. The implementation assumes a linear variation of the surface velocity field across each surface triangle until a triangle edge is reached. A triangle-to-edge-to-triangle data structure is used to move to the adjacent triangle where the integration process is continued. The streamline ODEs are integrated in reverse, starting at a point near the base of the vehicle, and ending at the stagnation point.

Attachment lines can be located using numerical methods based on the vector field topology of such inviscid flows . The topology of a vector field consists of critical points, i.e., points where the velocity is zero, and the tangent curves (instantaneous streamlines) which connect these points. Because the velocity at a critical point is zero, the velocity field in the neighborhood of the critical point is determined by \(\nabla u\). Critical points are classified, to a first order approximation, by the eigenvalues and eigenvectors of \(\nabla u\). Common classifications include a saddle, node, spiral, and center.

Using the inviscid surface velocity field generated in FlightStream, the method developed by Kenwright is used to locate the attachment lines. Kenwright developed two automatic feature extraction techniques that locate and distinguish separation and attachment lines on solid bodies in 3D numerical flow fields: (a) Phase Plane Algorithm and (b) Parallel Vector Algorithm.

These algorithms are based on eigenvalue analysis of the velocity gradient tensor and perform a local analysis of the vector field rather than a global analysis of the entire flow field. These methods are useful for analyzing large partitioned data sets, such as those computed on distributed memory architectures, because the elements can be processed in parallel. In 2D flows, this can occur only at critical points (half-saddles) located on 1D boundaries. These critical points are just places where flows along the surface in opposing directions converge and, by continuity, the wall shear stress is zero. In fully 3D flows over a 2D boundary, there are 1D loci where flows in relatively opposing directions meet , but the criterion of zero wall shear stress does not generally hold, because there is almost always flow along the 1D curve. In fact, the wall shear stress can reach zero only at isolated points along these so-called lines of separation.

The key requirement for any implementation of an attachment line method for FlightStream® requires that the technique be local, meaning that an independent test could be applied to each element on the surface of a body. A modified Phase Plane Algorithm [10] was implemented for this effort. This method models the flow over a simple triangular element and identifies flow patterns that contain streamline asymptotes. Given that the velocity vectors are defined at the vertices, a linear vector field [9] can be constructed that passes through the triangle and satisfies the prescribed vectors at the vertices:

\[\begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} + \begin{pmatrix} b_1 & c_1 \\ b_2 & c_2 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}\]

Here, \((x,y)\) is the Cartesian coordinate vector and \((\dot{x},\dot{y})\) the tangential velocity or shear stress vector. For a linear vector field, the coefficients \((a_1,a_2)\) and those in the 2x2 (Jacobian) matrix are constants. These constants can be computed analytically by substituting the coordinates and vectors from each vertex into this equation and then solving the resulting set of simultaneous equations. By differentiating this equation with respect to time and then algebraically manipulating the resulting equations, one can produce a pair of second order nonhomogeneous ordinary differential equations. If the determinant of the Jacobian matrix is nonzero, the solution has the form:

\[\begin{pmatrix} x(t) - x_{cp} \\ y(t) - y_{cp} \end{pmatrix} = \begin{pmatrix} \xi_1 & \eta_1 \\ \xi_2 & \eta_2 \end{pmatrix} \begin{pmatrix} a e^{\lambda t} \\ b e^{\mu t} \end{pmatrix}\]

Here \(\lambda\) and \(\mu\) are the eigenvalues of the Jacobian matrix and \((\varepsilon_1, \varepsilon_2)\) and \((\eta_1, \eta_2)\) are the eigenvectors. The two column eigenvectors form the eigenmatrix. The terms \(\alpha\) and \(\beta\) are arbitrary constants that define a particular curve in the phase plane. The constants \(x_{cp}\) and \(y_{cp}\) are the coordinates of the critical point:

\[x_{cp} = \frac{a_2 c_1 - a_1 c_2}{b_1 c_2 - b_2 c_1}, \quad y_{cp} = \frac{a_1 b_2 - a_2 b_1}{b_1 c_2 - b_2 c_1}\]

Here, \(x_{cp}\) and \(y_{cp}\) translate the coordinate system such that the origin coincides with the critical point. A linear vector field constructed in most triangles will have one critical point. However, that critical point will usually lie outside the boundary of the triangle. The critical points that do lie inside triangles correspond to those found by linear vector field topology methods . Using this formulation, tangent curves can be constructed that pass through the triangle. Rather than working in the physical plane of the triangle, we can simplify matters by transforming into canonical coordinates , that is, a coordinate system where the eigenvector directions are orthogonal:

\[\begin{pmatrix} X \\ Y \end{pmatrix} = \begin{pmatrix} \xi_1 & \eta_1 \\ \xi_2 & \eta_2 \end{pmatrix}^{-1} \begin{pmatrix} x(t) - x_{cp} \\ y(t) - y_{cp} \end{pmatrix} = \begin{pmatrix} a e^{\lambda t} \\ b e^{\mu t} \end{pmatrix}\]

These equations describe all necessary equations for this method. However, a condition under which a streamline asymptote will pass through the triangle is still required, as well as to determine the points of intersection. Tangent curves of the vector field can be constructed in the \((X, Y)\) plane using the above method. This is often referred to as the Poincare phase plane . By eliminating the integration variable, \(t\), one can express the trajectories of these curves in terms of an implicit scalar function. For the case where both eigenvalues are real numbers, the solution is either:

\[\psi(X, Y) = \frac{X^\mu}{Y^\lambda} \quad \text{or} \quad \psi(X, Y) = -\frac{Y^\lambda}{X^\mu}\]

The contours of \(\psi(X,Y)\) are everywhere tangent to the vector field and may be verified using the relationship \(\nabla \psi.U=0\), where \(U\) is the image of the vector field in the phase plane. By differentiating Equation-4 with respect to \(t\), the necessary transformation can be obtained, which maps the vector field, \(U\), into the phase plane. \(\psi(X,Y)\) is an exact solution to a rotational 2D linear vector field, which, in general, will not be divergence-free. A nonzero divergence means that mass is not conserved on the surface of the triangle. The fact that this system can lose mass is physically important because this accounts for fluid that leaves the surface as the flow converges on a separation line. The system will gain mass as the flow returns to the surface and diverges from an attachment line.

The phase portrait can assume one of two orientations in the phase plane depending on whether it is an attracting node \((\mu < \lambda < 0)\) or a repelling node \((0 < \mu < \lambda)\). Specifically, streamlines asymptotically diverge from the \(Y\)-axis for a repelling node, while they converge on the \(X\)-axis for an attracting node. A separation or attachment line will pass through a triangle if the triangle straddles the \(X\)-axis or \(Y\)-axis in the phase plane. The vertices of the triangle are mapped into the phase plane using the central expression in the above formulation.

Viscous-Inviscid Coupling

To perform viscous-coupled flow solutions, the laminar and turbulent boundary layer models were integrated with the inviscid solver via displacement of the inviscid boundary equal to the displacement thicknesses of the local boundary layers. A three-dimensional methodology applicable to unstructured surface meshes has been implemented into FlightStream®. The solver coupling takes the following form:

  1. Solve the potential flow field over the body and obtain the inviscid surface vorticity distribution.

  2. Compute the integral boundary layer

  3. Modify the surface flow boundary condition for the potential flow and solve for the second iteration.

  4. Repeat steps 1-3 until both the inviscid and viscous flow solutions are converged.

For boundary layers, the displacement thickness \((\delta^*)\) indicates the extent to which the surface would have to be displaced in order to be left with the same flow rate of the viscous flow, but with an inviscid velocity profile of \(u(z) = U_e\). Therefore, the displacement thickness is suited for the measure of displacement needed for the inviscid boundary. Fig. 6 illustrates the philosophical idea for this task.

The second option for inviscid boundary displacement is to simulate the displacement of a fixed surface mesh by blowing a velocity \((V_N)\) normal to the mesh face. The magnitude of the blowing velocity can be computed by using the momentum flux equations to give:

\[V_N = -\frac{\partial (U_e \delta^*)}{\partial s}\]

Here, \((\delta^*)\) is the displacement thickness values computed from the previous integral boundary layer computations; \(U_e\) is the local streamline velocity outside the boundary layer; and \(s\) is the direction along the local surface. The transpiration (blowing) velocity is then used to modify the FlightStream® inviscid Neumann boundary condition for the mesh faces using the following equation:

\[V_N = \frac{\partial (\phi + \phi_\infty)}{\partial n}\]

The iterative formulation of the Normal Blowing Boundary approach is illustrated in Fig. 8. The transpiration velocity approach has substantial numerical and practical advantages over the morphing boundary approach in that it does not necessitate morphing the physical mesh (which can cause solver stability and robustness issues, especially in regions involving geometric concavities) or updating the geometry influence matrices, which can be computationally expensive. Based on numerical tests, the transpiration velocity model was chosen as the superior model and implemented in FlightStream®.

normal blowing scheme
Figure 6. Normal Blowing Boundary technique for inviscid slip-wall displacement

Additional reading and validations of these methods in can be found at these references [].

Separation Criteria

In the boundary layer, if pressure in the direction of the flow increases (adverse pressure gradient), the flow can slow down and even reverse. If this adverse pressure is too severe, the boundary layer separates from the surface, leading to increased drag and possibly stall in the case of aircraft wings.

This section discusses the methods used to determine where the flow separates along a surface streamline in FlightStream.

Stratford Separation Criterion

After the inviscid solution is calculated, FlightStream® uses analytical models to estimate where the flow becomes separated. FlightStream® implements a modified Stratfod separation criterion to estimate boundary layer separation location. The Stratford criterion predicts a condition where this separation can be delayed or prevented, even under strong adverse pressure gradients. Stratfod derived an expression for predicting laminar boundary layer separation. [3]

\[\bar{C}_p = 1.0 - \left(\frac{u}{u_{max}}\right)^{2.0}\]
\[\bar{X}^2\bar{C}_p\left(\frac{dC_p}{dx}\right)^2 = 0.0104\]
\[\bar{X} = x - x'\]

Here \(\bar{X}\) is the effective length of the boundary layer and \(x'\) is the origin of the layer. Typically the point of maximum velocity. For turbulent boundary layers, Cebeci-Smith modified Stratford’s original formulation to say that the boundary layer is separated where:

\[0.4(10^{-6}Re)^{0.1} = \bar{C}_p(\bar{X}\frac{dC_p}{dx})^{\frac{1}{2}}\]

To extend this 1D method to a quasi-3D flow, the pressure gradients are calculated along surface streamlines.

Valarezo Separation Criteria

The Valarezo criteria is a method of determining separation based on the difference in pressure at the leading and trailing edges of a wing. The method uses an empirically-derived "pressure difference rule" to predict maximum lift for complex multielement wing geometries . The pressure difference rule states that at a given Reynolds number and Mach number, there is a specific pressure difference between the suction peak and trailing edge of an airfoil at the maximum lift condition. This rule applies to both single and multielement airfoils. This method tends to predict separation later than Stratford. [4]

Separation Velocity Profile

The Swafford velocity profile is used to determine the velocity profile at the separation point . The full derivation can be found in the reference. The final expression implemented in FlightStream® is:

\[\frac{u^+}{u_\tau} = \frac{S}{\kappa} \tan^{-1} \left(0.09y^+\right) + \left(u_c^+ - \frac{S\pi}{0.18}\right) \tanh\left[\frac{1}{2}\left(\frac{y}{\delta}\right)^2\right]\]

Where the parameters \(a\) and \(b\) are function of \(C_f\) and \(Re_\theta\) and \(H\).

Additionally, FlightStream sets the \(C_f\) to zero in separated flow regions.

Transonic Solver

The transonic solver in FlightStream is based on a field-panel approach to solving the nonlinear full-potential equation. In a steady motion of an inviscid, isentropic, irrotational flow it can be shown that the nonlinear full potential equation is:

\[c_*^2 \nabla^2 \Phi = \nabla \Phi \cdot (\nabla \Phi \cdot \nabla \nabla \Phi) + \frac{\gamma - 1}{2} (\nabla \Phi)^2 \nabla^2 \Phi\]

And under the assumption of irrotational, steady, and homentropic flow, the total enthalpy is constant in the entire flowfield.

\[h + \frac{1}{2} \mathbf{u} \cdot \mathbf{u} = h_* = \text{constant in the entire flowfield}\]

Mach Number

If the magnitude of \(|\mathbf{u}|\) is equal to the local acoustic speed (i.e., \(|\mathbf{u}| = c\)), the flow is said to be sonic. The sonic speed, latexmath:[c_s\), is given by

\[\begin{equation} \frac{\gamma + 1}{2} c_s^2 = c_*^2 \end{equation}\]

The sonic speed is the speed at which the local fluid speed and the speed of sound are the same. If the flow speed is greater or less than the sonic speed, the flow is called supersonic or subsonic, respectively. From Eqs. 13 and 14, \(|\mathbf{u}|^2 < c_s^2 < c^2\) in the case of subsonic flow and \(|\mathbf{u}|^2 > c_s^2 > c^2\) in the case of supersonic flow. Thus, the dimensionless parameter \(M\) is defined as

\[\begin{equation} M = \frac{|\mathbf{u}|}{c} \end{equation}\]

at any point on a streamline. This quantity is called the Mach number. Based on the inequalities given above, the Mach number will be greater than, equal to, and less than unity for supersonic, sonic, and subsonic flow, respectively. It should be noted that the value for the local Mach number and, consequently, the designation of the flow, is entirely dependent on the frame of reference chosen.

Low-Order Supersonic Solver

Local Inclination Methods

The low-order methods in FlightStream for calculating the pressure coefficient on the surface of a body depend solely upon the local inclination of the surface. The inclination of the surface is given by:

\[\begin{equation} \theta = \cos^{-1}\left(\hat{n} \cdot \frac{V_{\infty}}{\lVert V_{\infty} \rVert}\right) \frac{\pi}{2} \end{equation}\]

where \(\hat{n}\) is the local surface normal vector and \(V_{\infty}\) is the freestream velocity vector. Surfaces where \(\theta > 0\) are categorized as windward surfaces, and surfaces where \(\theta \leq 0\) are categorized as leeward surfaces.

Tangent Cone Method

The Tangent Cone method approximates the pressure coefficient as the pressure coefficient on the surface of a cone in the same flow. The solution for the cone pressure coefficient is carried out by solving the Taylor-Maccoll equations for conical flow given as \cite{anderson2003modern}:

\[\begin{equation} \frac{\partial V_r}{\partial \theta} = -V_{\theta} \end{equation}\]

and

\[\begin{equation} \frac{dV_{\theta}}{d\theta} = V_{\theta}^2 - V_r - \frac{\gamma - 1}{2} \left(1 - \frac{V_r^2 - V_{\theta}^2}{2}\right) \left(2V_r + V_{\theta} \cot(\theta) \right) \left(\frac{1 - \frac{V_r^2 - V_{\theta}^2}{2}}{1 - V_r^2 + V_{\theta}^2}\right) \end{equation}\]

This is a first-order ODE, which can be solved numerically. The solution procedure is carried out inversely. First, an arbitrary shock angle is assumed. Then the Taylor-Maccoll equations are solved, marching inwards from the assumed shock angle. The deflection angle ( \(\theta\) ) is sought where the radial velocity ( \(V_r = 0\) ). This deflection angle then becomes the cone angle for the assumed shock angle. The surface pressure coefficient ( \(C_P\) ) on the cone is the known. Oblique shock relations and then isentropic relations give the surface pressure coefficient on the cone.

Solving the Taylor-Maccoll equations is computationally expensive. FlightStream does not directly solve these equations for every run. Instead, FlightStream uses a sixth-order polynomial to represent the cone surface pressure coefficient expressed as

\[\begin{equation} C_P(\theta) = C_0 + C_1 \theta + C_2 \theta^2 + C_3 \theta^3 + C_4 \theta^4 + C_5 \theta^5 + C_6 \theta^6 \end{equation}\]

\(C_0\) through \(C_6\) were calculated and tabulated across a range of freestream Mach Numbers and specific heat ratios. When FlightStream is run with the Tangent Cone method, the pressure coefficients on a panel with a given inclination are evaluated based on the polynomials at the nearest tabulated Mach numbers and specific heat ratios to the flight conditions. These values are linearly interpolated between to come up with a pressure coefficient for each panel.

Similarly to the wedge, there is a certain cone angle above which the shock is detached and the Taylor-Maccoll equations cannot be used to predict the surface pressure on the cone. This angle \(\theta_{max}\) arises during the solution process of the Taylor-Maccoll and is stored and tabulated. When a panel inclination is above the cone \(\theta_{max}\) for given flight conditions, another pressure method must be used.

For 90% of panel inclinations, the percent error for the interpolated pressure coefficient is within 2% of the exact pressure coefficient. For panel inclinations with a large relative error, the pressure coefficient is very low. As a result, these panels with high relative error have little impact on the overall aerothermal results. This level of accuracy is judged to be of an appropriate level of fidelity for the mid-fidelity nature of the FlightStream code.

Leeward Pressure Coefficients Using Prandtl-Meyer Expansion

The method for calculating the leeward pressure coefficients used in this work is a Prandtl-Meyer Expansion from the freestream Mach number. The Prandtl-Meyer function is given by

\[\begin{equation} \nu(M) = \sqrt{\frac{\gamma + 1}{\gamma - 1}} \tan^{-1} \left(\sqrt{\frac{\gamma - 1}{\gamma + 1}(M^2 - 1)} \right) - \tan^{-1} \left(\sqrt{M^2 - 1}\right) \end{equation}\]

For conducting analysis using the Prandtl-Meyer function, the inclination angle is related to the Prandtl-Meyer function by

\[\begin{equation} \theta = \nu(M_{\infty}) - \nu(M) \end{equation}\]

This expansion is applied through \(\theta_{max} < \theta < 0\). The max turning angle is given by

\[\begin{equation} \theta_{max} = \left[\frac{\pi}{2} - \nu(1) \right] + \nu(M_{\infty}) \end{equation}\]

Beyond the max turning angle, the flow is considered as separated and the pressure is determined using the base region relations shown later. Solving the Prandtl-Meyer function results in the Mach number at an expansion surface of inclination \(\theta\). The pressure coefficient on the expansion surface is given as

\[\begin{equation} C_P = \frac{2}{\gamma M^2_{\infty}} \left[\left(1 + \frac{\gamma - 1}{2} M^2 \right)^{-\frac{\gamma}{\gamma - 1}} - 1\right] \end{equation}\]

Due to the complexity of solving the Prandtl-Meyer function for each expansion surface, the function is solved numerically at the start of each run and the solution tabulated.

Low-Order Hypersonic Solver

The original Newtonian method, also referred to as the straight Newtonian method, utilizes the following equation for pressure calculations \cite{anderson1989hypersonic}:

\[\begin{equation} C_P = 2 \sin^2 \theta \end{equation}\]

where \(C_P\) is the pressure coefficient on the surface and \(\theta\) is the inclination angle. The pressure coefficient on the surface as calculated by the straight Newtonian method is independent of Mach number. Lester \cite{grantz1994calibration} suggested a modification to the Newtonian method where \(C_P\) is scaled by relating \(sin^2 \theta\) to \(C_P/C_{P,{max}}\). \(C_{P,{max}}\) is the stagnation pressure coefficient at aft of a normal shock wave given by:

\[\begin{equation} C_{P_{\text{max}}} = \frac{2}{\gamma M^2_{\infty}} \left(\frac{\gamma + 1}{4\gamma M^2_{\infty} - 2(\gamma - 1)}\right) \left[1 - \frac{\gamma + 1}{2\gamma M^2_{\infty}} \right] \end{equation}\]

This modification adds Mach number dependence to the Newtonian method. The modified Newtonian method can now be stated as:

\[\begin{equation} C_P = C_{P_{\text{max}}} \sin^2 \theta \end{equation}\]

The Newtonian methods only apply to the windward surface of the body. Leeward surfaces are handled using Prandtl-Meyer expansion theory as described in Leeward Pressure Coefficients Using Prandtl-Meyer Expansion.

Note
Known Current Limitations

This solver does not currently capture thermal effects.

High-Order Supersonic Solver

Panel methods arise from applying the method of Green’s functions to the Prandtl-Glauert equation [5], [6]. The Prandtl-Glauert equation is given by

\[ (1 - M^2_\infty) \phi_{xx} + \phi_{yy} + \phi_{zz} = 0\]

where \(M_\infty\) is the Mach number of the undisturbed freestream (assumed to be aligned with the x-axis) and \(\phi\) is the perturbation velocity potential. The total velocity potential is given by:

\[\Phi = U (\phi_\infty + \phi)\]

where \(\Phi\) is the velocity potential (i.e. u = ∇Φ is the fluid velocity), \(U\) is the magnitude of the freestream velocity, \(\phi_\infty\) is the freestream potential, and \(\phi\) is the perturbation potential. The Prandtl-Glauert equation is derived from the continuity equation, Euler’s equations of motion, and the energy equation by assuming the flow is homentropic and irrotational, and neglecting nonlinear terms [7]. This makes the Prandtl-Glauert equation a reasonably accurate model of fluid flow for small-perturbation flows outside of the transonic and hypersonic regimes [7]. Applying the method of Green’s functions, the perturbation velocity potential at any point in the flow may be found from [7]

\[\begin{equation} \kappa \phi(r) = \int_S G \sigma dS - \int_S \mu \frac{\partial G}{\partial n} dS \end{equation}\]

where \(r = (x, y, z)\) is the point of interest, \(\kappa\) is \(4\pi\) if the flow is subsonic and \(2\pi\) if the flow is supersonic, \(G\) is the Green’s function for the Prandtl-Glauert equation, \(S\) is the surface of the configuration being analyzed, \(\sigma = \frac{\partial \phi}{\partial n} - \frac{\partial \phi_i}{\partial n}\) and \(\mu = \phi - \phi_i\) are distributions of source and doublet strength across \(S\), respectively, and \(\frac{\partial}{\partial n}\) denotes the conormal derivative, described subsequently. The supersonic Green’s function is given by [6]

\[ G = \begin{cases} \frac{-1}{\sqrt{(x - \xi)^2 + (1 - M_\infty^2)\left[(y - \eta)^2 + (z - \zeta)^2\right]}} & \text{if } \xi - x < 0 \text{ and } (x - \xi)^2 + (1 - M_\infty^2)\left[(y - \eta)^2 + (z - \zeta)^2\right] > 0, \\ 0 & \text{otherwise}. \end{cases}\]

where \(\xi, \eta, \zeta\) is the integration point on \(S\). The region where \(G \neq 0\) for a fixed point \((x, y, z)\) is called the domain of dependence of point \((x, y, z)\). Alternatively, the region where \(G \neq 0\) for a fixed point \(\xi, \eta, \zeta\) is called the domain of influence of point \(\xi, \eta, \zeta\). These domains are shown schematically for a single point.

domain dependence
Figure 7. Domains of dependence and influence for a point P in supersonic flow

The subsonic Green’s function is identical, except that it is nowhere zero. Thus, in subsonic flow, the domains of dependence and influence for any point are the entire fluid region. The conormal derivative is given by

\[ \frac{\partial}{\partial \tilde{n}} = \tilde{n}^T \nabla\]

where \(n\) is the outward normal vector.

\[ \tilde{n} = [B]n\]
\[ [B] = [I] - M_\infty^2 \hat{c} \hat{c}^T\]
\[ \hat{c} = \frac{U}{U_\infty}\]

In these equations, \(n\) is the outward normal to \(S\), and \(\tilde{n}\) is called the conormal vector. Green’s function may be solved approximately by breaking \(S\) into \(N\) discrete surfaces, called panels. Doing so and allowing for the presence of a discretized wake surface, Eq. 26 becomes [7]

\[\begin{equation} \kappa \phi(\mathbf{r}) = \sum_{i=1}^N \left[ \iint_{S_i} \left( G \sigma_i - \mu_i \frac{\partial G}{\partial \mathbf{n}} \right) dS_i \right] - \sum_{j=1}^{N_w} \iint_{S_{w_j}} \mu_j \frac{\partial G}{\partial \mathbf{n}} dS_{w_j} \end{equation}\]

where \(\sigma_i\) and \(\mu_i\) are analytic distributions of source and doublet strength across each panel. For a general configuration, \(\sigma_i\) and \(\mu_i\) are not known a priori. The method used here for determining \(\sigma_i\) and \(\mu_i\) will be discussed later. In Eq. 27, the doublet influence on the velocity potential has the opposite sign from that of the source influence. Because of this, panel methods traditionally define a surface distribution of doublets as inducing flow through the surface in the direction opposite the normal vector (e.g. see [19, 46]). This is different from the standard definition of a doublet where flow is induced in the same direction as the normal vector. For consistency with the literature, the traditional panel method definition is maintained in FlightStream.

Sub- and Superinclined Panels

The nature of linearized flow theory on which FlightStream is based leads to some peculiar practical considerations. The most important of these considerations is the difference between subinclined and superinclined panels.

One basic assumption built into the Prandtl-Gluert equation is that the Mach number is everywhere constant. This then means that there is a constant characteristic angle (the Mach angle) throughout the flow, and all shocks/expansions have that same angle relative to the freestream. In addition, the cone bounding the domain of influence of any given point also has this angle.

Subinclined panels are inclined to the freestream at an angle less than the Mach angle. Superinclined panels have an inclination angle greater than the Mach angle. The implementations of these two types of panels are fundamentally different, and currently, only subinclined panels can be used in FlightStream. Future work will include implementing superinclined panels. However, their use will be limited due to the nature of linearized theory.

sub- and superinclined panels
Figure 8. Illustration of sub- and superinclined panels

Application of HOSS

Based on the described methodology, the solution of the supersonic panel method must satisfy Mach cone requirements.

Note
Known Current Limitations

The development of the HOSS is ongoing. As such, certain aspects are still being improved, and the HOSS may not produce good results in certain instances. Some known instances are:

  • Wake Modeling: The HOSS does not currently support wake modeling. It is anticipated that this should not be an issue for many users, as in supersonic flow, influences propagate only downstream, meaning most (if not all) of the surface solution is not affected by the wake or lack thereof, depending on the configuration. However, this should be kept in mind when modeling configurations with potential wake influences as solution accuracy can be degraded.

  • High-Inclination Panels: The HOSS is built on linearized theory, which assumes that any perturbations in the flow due to the presence of the configuration being modeled are small relative to the freestream. Thus, it has been seen in many instances that flows with large perturbations (i.e. due to panels steeply inclined to the freestream) often have poor localized results. For instance, on configurations with steeply inclined panels at higher Mach numbers, the pressure coefficient can erroneously drop and even become negative, which is of course non-physical. This can lead to a severe under-prediction of the total pressure drag on the configuration.

These limitations will be addressed in ongoing development. It is anticipated that future versions of FlightStream will include improvements to significantly enhance the accuracy and usability of the HOSS.

Hypersonic Solver

This section will talk about the hypersonic solver.

Postprocessing

Loads Calculation

Lorem ipsum dolor sit amet, consectetur adipiscing elit. Integer egestas

Drag Integration

FlightStream employs a dual approach for calculating drag:

  • Vorticity-Based Induced Drag Calculation: FlightStream offers a unique method for calculating induced drag that relies on user-defined trailing-edge boundary conditions. This method computes integrated circulation and downwash at the trailing edges to determine a downwash-based drag, conceptually similar to lifting line theory for arbitrary sectional load distributions. This feature, stemming from FlightStream’s vorticity-based formulation, offers a distinct perspective on induced drag estimation. Due to the closed-loop nature of the vorticity, this method is less mesh sensitive than the surface pressure integration method.

  • Surface Pressure Integration: This method leverages the Bernoulli equation to determine the pressure coefficient (Cp) at the surface of the body using the computed surface velocities. This method is applicable to all body types, including fuselages and bluff bodies. Only the pressure integration method is capable of capturing separation effects and off-body shocks. This integration lumps pressure and induced drag together.

The surface pressure integration method, explained in source, is based on applying Bernoulli’s equation to the body’s surface. The pressure coefficient (Cp) is calculated using the formula:

\[C_p = 1 - \left(\frac{V}{V_\infty}\right)^2\]

Nomenclature

English letter variables

Variable Description

\(C_f\)

Skin friction coefficient for rough surfaces

\(F\)

Auxiliary variable for turbulent boundary layer model

\(H\)

Shape factor

\(H_c\)

Compressible shape factor

\(H_{\delta-\delta}\)

Shape factor: Difference between boundary layer thickness and displacement thickness

\(M_e\)

Mach number at the edge of the boundary layer

\(Re_\theta\)

Reynolds number based on momentum thickness

\(Re_k\)

Reynolds number based on surface roughness height

\(S_B\)

Surface of the immersed body

\(T_e\)

Temperature at the edge of the boundary layer

\(T_R\)

Recovery temperature

\(T_0\)

Reference temperature

\(U\)

Velocity at the outer edge of the boundary layer

\(V_N\)

Normal blowing velocity

\(V_S\)

Local velocity on the boundary

\(X\)

Transformed coordinate in the phase plane

\(Y\)

Transformed coordinate in the phase plane

\(a_e\)

Speed of sound at the edge of the boundary layer

\(c_f\)

Skin friction coefficient

\(k\)

Surface roughness height

\(n\)

Unit normal vector

\(r\)

Distance between two points in the flow field

\(s\)

Coordinate along the local surface

\(x_{cp}\)

x-coordinate of the critical point

\(y_{cp}\)

y-coordinate of the critical point

Greek letter variables

Variable Description

\(\Phi_P\)

Potential at point \(P\)

\(\bar{K}\)

Average pressure gradient parameter

\(\lambda\)

Eigenvalue of the Jacobian matrix

\(\mu\)

Doublet strength

\(\mu\)

Eigenvalue of the Jacobian matrix

\(\mu_R\)

Dynamic viscosity at recovery temperature

\(\nabla^2\)

Laplace operator

\(\nabla \phi\)

Velocity field obtained from the potential function

\(\phi\)

Potential flow function

\(\phi_i\)

Potential flow function at a given interior point

\(\phi_\infty\)

Freestream potential

\(\sigma\)

Source strength

\(\theta\)

Momentum thickness

\(\theta_c\)

Compressible momentum thickness

\(\xi_1, \xi_2\)

Eigenvector components in the phase plane

\(\eta\)

Coordinate along the surface streamline

\(\eta_1, \eta_2\)

Eigenvector components in the phase plane

\(\nu\)

Kinematic viscosity

\(\delta\)

Boundary layer thickness

\(\delta^*\)

Displacement thickness

Appendix-A

Lorem ipsum dolor sit amet, consectetur adipiscing elit. Integer egestas vulputate consequat. Proin libero mi, auctor sed orci nec, lobortis blandit massa. Nullam ut tellus eros. Nulla posuere diam in viverra imperdiet. Duis id ex facilisis, lacinia velit in, posuere nisl. In quis magna suscipit, tempus risus a, lacinia nulla. Morbi vulputate magna nisi, sed laoreet lectus scelerisque sed. Suspendisse bibendum dignissim erat, in egestas ipsum lobortis non. Praesent purus felis, porttitor et quam vel, consectetur ullamcorper quam. Etiam ut scelerisque massa. Duis maximus erat vel orci accumsan suscipit.

Glossary

Abbreviations

Variable Description

UND

University of North Dakota

References

[1] W. M. Kutta, “Lift Forces in Fluid Flow,” Illustrierte Aeronautische Mitteilungen, p. 133, 1902.

[2] F. M. White and J. Majdalani, Viscous Fluid Flow, 4th ed. McGraw-Hill Higher Education, 2022.

[3] B. Stratford, “The Prediction of Separation of the Turbulent Boundary Layer,” Defense Technical Information Center, ADA091204, 1980. [Online]. Available: https://apps.dtic.mil/sti/tr/pdf/ADA091204.pdf.

[4] W. O. Valarezo and V. D. Chin, “A CPmax-Based Flow Separation Model for Application to Wings with High-Lift Devices,” Journal of Aircraft, vol. 31, no. 1, pp. 103–111, 1994, doi: 10.2514/3.46461.

[5] L. Morino, “Boundary Integral Equations in Aerodynamics,” Applied Mechanics Reviews, vol. 46, 1993.

[6] G. N. Ward, Linearized Theory of Steady High-Speed Flow. Great Britain: Cambridge University Press, 1955.

[7] C. D. Goates, “Development, Implementation, and Optimization of a Modern, Subsonic/Supersonic Panel Method,” PhD thesis, Utah State University, 2023.