FlightStream Theory Manual

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

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

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

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:

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

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.

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.

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:

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.

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:

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:

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.

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.

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 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 regions are used to model the surface-pressure effects of regions of bulk separated flow. In FlightStream^®, these regions are defined by the user. And a uniform \(C_{p,base}\) is specified for the base region. 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.

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.

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 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:

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:

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

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.

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:

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:

The auxiliary equations needed to complete the above equations are :

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:

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:

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.

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:

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\)):

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:

\(\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\):

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 .

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.

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:

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:

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:

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:

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:

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.

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:

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:

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:

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^®.

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

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.

FlightStream employs a dual approach for calculating drag:

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:

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.

[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] P. M. Sinclair, “A three-dimensional field-integral method for the calculation of transonic flow on complex configurations — theory and preliminary results,” The Aeronautical Journal, vol. 92, no. 916, pp. 235–241, 1998, doi: 10.1017/S0001924000016183.