Introduction

Advanced computational methods continue to be developed to analyze challenging problems in aerospace engineering. Such methods include aerodynamics, structural analysis, and coupling disciplines such as aeroelasticity and flight dynamics. In particular, computational fluid dynamics, commonly referred to as “CFD,” and the finite element method, or “FEM,” have been revolutionary computational tools for predicting aerodynamic behavior and structural response, respectively. The intention of this chapter is not to review the details of such techniques and their associated numerical algorithms in themselves, but to provide an introduction and general overview of the foundational principles and how they are applied, along with a summary of their capabilities and some of their limitations. Nor is it meant to be an all-encompassing treatment, far from it.

Traditionally, the aircraft industry has relied on faster and more cost-effective computational methods for aircraft design and flight performance prediction, with considerable success. One does not need CFD or the FEM to design and build a flight vehicle that works as intended; history speaks for itself. However, both tools can significantly help refine aircraft design to maximize performance at the lowest cost, while also improving safety and reliability. Understanding why something works (or doesn’t) is crucial to achieving commercial success and driving continued technological advancement. Wind tunnel testing has long been a critical component of flight vehicle design and will remain essential for CFD validation. Advances in computer technology, including code execution speed, memory, data storage access speed, and the availability of inexpensive, high-capacity data storage systems, continue to encourage the development of more ambitious numerical methods for use in the design process.

Significant advances in fundamental research on modeling and algorithmic development have been made over the last two decades. Aerodynamics and structural problems, once relegated to the fastest available supercomputers at national centers, can now be solved using local computational resources. For example, most universities today have powerful computing clusters, enabling research into new computational methods applicable across various technical fields. The ability to predict flows across the entire aircraft could open new design directions, such as reconsidering supersonic transports (SSTs), an example being shown in Figure 1. It has also inspired broader innovations in flight-vehicle design, enabling the investigation of concepts outside conventional wisdom with greater confidence. Thinking “out of the box” can now be quantified for the right reasons!

Continued advances in computer capabilities continue to accelerate,^[1] which have also played a significant role in motivating and accelerating the development of more integrated computational techniques that encompass several traditionally separate engineering disciplines. For example, the governing equations and solution methods used in one field are rarely the same as those used in another; therefore, some form of method coupling is required. Increasingly tight mathematical, algorithmic, and numerical coupling has enabled, for example, the simultaneous solution of the coupled aerodynamics-structural dynamics problem. If fully realized, this integration strategy will enable the prediction of an entire flight vehicle’s behavior under arbitrary flight conditions, thereby enabling more capable and efficient vehicle designs. However, significant integration challenges and advanced numerical techniques remain necessary to further develop this technology, with the ultimate goal of predicting the vehicle’s characteristics before its first flight.

Aerodynamics

Over the past three decades, significant advances have been made in understanding challenging aerodynamic problems through computational fluid dynamics (CFD). Numerical methods, such as finite-difference (FDM) and finite-volume (FVM), are applied to the governing equations (e.g., the Navier-Stokes equations) to model the flow field around an aircraft. The choice of governing equations affects the level of physics the CFD scheme can approximate, as well as the computational effort and time required to solve the problem. Not all numerical methods are created equal, and the choice of method significantly affects the accuracy, stability, and cost of the solution. The definition of “cost” comes down to how much time^[2] and effort go into obtaining a solution, and most CFD methods are considered costly by any preliminary design standard.

Developing comprehensive CFD methods is typically regarded as the pinnacle of aerodynamic modeling. So far, a wide range of numerical techniques has been used to solve a hierarchy of equations governing flow, from irrotational, inviscid potential-flow approximations to the Euler and Navier-Stokes equations. The latter set of equations can account for compressibility, thermodynamics, and turbulence. Each level in the hierarchy of governing equations increases the physical modeling capability but also requires greater computational effort and cost to solve.

A key issue in using these numerical methods is the definition and construction of appropriate computational grids for solving the governing equations, a process that requires considerable skill and effort. Other problems include the difficulty of modeling turbulent flow, for which no single method exists. Furthermore, CFD methods must be validated against measurements to achieve the high predictive confidence required for aircraft design. Commercially available tools, such as ANSYS Meshing and Pointwise, can aid in grid generation, while CFD solvers such as ANSYS Fluent and OpenFOAM can perform the actual CFD simulations. However, there are also many other options for grid generation and CFD solvers.

History

The earliest numerical solutions to viscous flow problems about bodies go back almost a century. Alexander Thom, a Carnegie Fellow at the University of Glasgow, published a paper in the Proceedings of the Royal Society in 1933 titled: “The Flow Past Circular Cylinders at Low Speeds.”^[3] In this work, the viscous flow was solved over a numerically defined grid by “arithmetical” means. Back then, there were no digital computers, so calculations were performed “by hand,” i.e., using a pencil, paper, a slide rule, and perhaps a mechanical calculator. Thom remarked in his paper that the “whole process is long,” but his computed results compared favorably with flow experiments.

The advent of digital computers in the 1950s provided the computational power to begin to tackle numerical solutions to the governing equations of fluid dynamics. Early work by pioneers such as John von Neumann and his collaborators laid the groundwork for numerical methods to solve partial differential equations. Von Neumann is often referred to as the “Father of Modern Computer Architectures.” During the 1960s, CFD began to evolve more rapidly in parallel with advances in computer technologies. Developing algorithms based on finite-difference, finite-element, and finite-volume methods enabled more accurate and efficient solutions to various forms of the governing equations of fluid dynamics. Finite-difference techniques were already well established by then, starting from Brook Taylor^[4] in the early 1700s, with other contributions by Leonhard Euler and Carl Runge. There was also much new work in numerical methods, and the accuracy and stability of such schemes were increasingly crucial for solving aerospace problems.

The publication of numerous technical papers and essential textbooks during the 1960s helped disseminate these numerical methods and apply them to fluid-flow problems. Textbooks by Brian Spalding^[5] and Suhas Patankar^[6] were used to teach CFD to students as part of courses in fluids. NASA and other aerospace laboratories worldwide have played a crucial role in advancing CFD for practical applications in aerospace engineering. In 1987, NASA established the Advanced Supercomputing Division (ASD) at NASA Ames. Additionally, the U.S. Army Aeronautical Research Laboratory played a crucial role in advancing the field of CFD. The first commercial CFD software packages emerged during this period, making the technology accessible to a broader range of engineers and researchers.

The 1980s and 1990s saw substantial improvements in numerical techniques, computational power, digital memory, and storage devices. Desktop computers were still somewhat limited in their capabilities, and most serious computational work was done on “mainframe” computers. Remotely accessible supercomputers, such as the CRAY, became available, and desktop computers also became increasingly capable, enabling more complex local CFD simulations. Advances in turbulence models, grid generation techniques, and numerical methods continued to contribute to more accurate and robust CFD solutions. During this period, CFD also started to find applications beyond the aerospace field, including automotive, chemical, and civil engineering.

CFD has continued to advance since the early 2000s, driven by the development of more sophisticated algorithms, increased computational power, and the emergence of parallel computing. High-performance computing (HPC) and graphics processing units (GPUs) have become increasingly commonplace, enabling the simulation of more complex, large-scale aerospace problems. Additionally, CFD has benefited from advances in data visualization techniques, which have made it easier for engineers to interpret and analyze large datasets and understand complex flows. Today, CFD is an essential tool in engineering analysis, with applications ranging from the design of flight vehicles to weather prediction, biomedical engineering, and environmental studies. The continuous development of software, increased accessibility of computational resources, and ongoing research mean that CFD will remain at the forefront of future scientific and engineering innovation.

Hierarchy of Aerodynamic Equations

Computational aerodynamic methods are built on a hierarchy of governing equations. Each level in the hierarchy is obtained from a more general level by making specific physical assumptions, such as neglecting viscosity, assuming irrotational flow, assuming small perturbations, or neglecting compressibility. These assumptions reduce the mathematical complexity and computational cost, but they also remove physical effects from the model. Therefore, the choice of governing equation is always a compromise between fidelity, cost, and the aerodynamic quantities that must be predicted.

The Navier-Stokes equations provide the most general continuum description used in aerodynamics. In conservation form, these equations may be written schematically as

(1)  

where is the vector of conserved variables, , , and are the inviscid flux vectors, and , , and are the viscous flux vectors. These equations express conservation of mass, momentum, and energy for a viscous fluid. In principle, they include compressibility, vorticity, viscous shear stress, heat conduction, shock waves, boundary layers, and flow separation. In practice, however, aircraft flows are usually turbulent, so the Navier-Stokes equations must be solved using some level of turbulence modeling or turbulence resolution, such as Reynolds-averaged Navier-Stokes (RANS), large-eddy simulation (LES), or, in highly restricted research problems, direct numerical simulation (DNS).

If the viscous terms are neglected, then

(2)  

and the Navier-Stokes equations reduce to the Euler equations, i.e.,

(3)  

The Euler equations describe inviscid compressible flow. They can represent pressure waves, compressibility effects, shock waves, and rotational inviscid motion. However, because viscosity has been removed, they cannot directly predict boundary-layer development, skin-friction drag, heat transfer at a wall, or viscous separation. For high-Reynolds-number flows, Euler calculations can still yield useful pressure loads over much of an aircraft. However, they must be interpreted with caution, as the missing viscous physics may still govern the actual flow.

A further reduction is possible if the inviscid flow is also assumed to be irrotational. In that case,

(4)  

and a velocity potential can be introduced such that

(5)  

The governing equation can then be written in terms of a single scalar unknown, , rather than in terms of the velocity components and thermodynamic variables separately. For compressible, inviscid, irrotational flow, this leads to the full-potential equation. In two dimensions, one common steady form can be written as

(6)  

where

(7)  

and is the local speed of sound. This equation is nonlinear because the local speed of sound depends on the local velocity magnitude through the isentropic relation. Full-potential methods were historically important because they retained nonlinear compressibility effects while being much less expensive than solving the Euler or Navier-Stokes equations. However, the irrotational assumption limits their ability to describe rotational, separated, and vortical flows unless additional modeling is introduced.

The full-potential equation can be simplified further when the flow consists of a dominant freestream plus small perturbation velocities. For example, for a freestream in the direction, the potential may be written as

(8)  

where is a perturbation potential, and the perturbation velocities are small compared to . With appropriate scaling and retention of the leading nonlinear compressibility terms, this reduction leads to the small-disturbance equations. Near transonic speeds, the resulting transonic small-disturbance (TSD) equation retains the essential nonlinear behavior associated with mixed subsonic and supersonic regions. A representative two-dimensional steady TSD form is

(9)  

The TSD equation is much simpler than the full-potential or Euler equations. However, it still captures the important transonic feature that the governing equation may change type locally between elliptic and hyperbolic behavior. This makes it useful for attached transonic airfoil and wing flows with relatively small disturbances. Its limitations are equally important: it is not appropriate for large-disturbance flows, bluff bodies, strong separation, highly vortical flows, or configurations where the geometry violates the small-disturbance assumptions.

If the perturbation equation is linearized further, then the classical linearized potential-flow equation is obtained, i.e.,

(10)  

This equation shows explicitly how compressibility changes the character of the flow. For subsonic flow, where , the equation is elliptic. For supersonic flow, where , the coefficient of changes sign, and the equation becomes hyperbolic. This change in mathematical type is one reason compressible aerodynamic flows behave so differently below and above Mach 1.

Finally, if the flow is assumed to be steady, incompressible, inviscid, and irrotational, then the velocity potential satisfies Laplace’s equation, i.e.,

(11)  

or, in Cartesian coordinates,

(12)  

This equation is linear and forms the basis of classical incompressible potential-flow methods, including source, doublet, vortex, and panel methods. Such methods are computationally efficient and remain useful for estimating pressure distributions, lift, induced effects, and attached-flow behavior. However, because the equations are inviscid and incompressible, they cannot directly predict viscous drag, boundary-layer growth, shock waves, heat transfer, or flow separation.

Therefore, the hierarchy may be summarized as:

Each step removes some physical content from the governing equations. Moving down the hierarchy reduces computational cost and often improves mathematical transparency, but it also narrows the range of flows that can be predicted reliably. The best aerodynamic model is not necessarily the most complex; it is the simplest model that captures the dominant physics of the problem being studied.

The natural starting point for this hierarchy is the Navier-Stokes equations because many lower-order aerodynamic models can be understood as simplifications of them. Even when a practical calculation uses an Euler, full-potential, boundary-layer, or panel method, the assumptions behind that method are best understood by asking what has been retained from, and what has been removed from, the full viscous-flow equations. For this reason, the discussion begins with the Navier-Stokes equations, not because they are always the most practical design tool, but because they provide the most complete continuum statement of the aerodynamic problem.

Navier-Stokes Equations

Most aerodynamic problems with flight vehicles involve viscous effects and turbulence. These effects are significant in turbulent boundary layers on wings and the airframe, as well as in a range of component-interaction phenomena, such as wing/fuselage and powerplant/wing interference. Modeling flow separation is also essential, as its onset can define the boundaries of an aircraft’s operational flight envelope. The Navier-Stokes equations are the fundamental equations that govern fluid dynamic behavior. They encompass the principles of conservation of mass, momentum, and energy, as well as the changes in these quantities within a fluid.

The Navier-Stokes equations are a set of coupled, highly nonlinear partial differential equations. This means they can only be solved by making simplifying assumptions, which in some cases may be rather sweeping. One possible simplification is to assume laminar flow, in which case viscous shear stresses are still present, but turbulent fluctuations and turbulent stresses are absent. However, flows over aircraft are typically turbulent at operational Reynolds numbers, so practical aerodynamic calculations usually require some form of turbulence modeling or turbulence resolution.

Tensor Form

The tensor form is often used to express the Navier-Stokes equations. The conservation of mass (continuity equation) is given by

(13)  

and in tensor form by

(14)  

where is the flow density and , , is the velocity vector. For an incompressible flow, where , then

(15)  

Conservation of momentum (momentum equation) gives

(16)  

and in tensor form by

(17)  

where is the pressure, is the viscous stress tensor, and is a body force per unit mass, e.g., acceleration under gravity.

Rounding out the conservation laws is the conservation of energy. Written in terms of the internal energy per unit mass, one common form of the energy equation is

(18)  

and in tensor notation

(19)  

where is the internal energy per unit mass, is temperature, is the thermal conductivity, and is the external heat-addition rate per unit mass. The term involving accounts for compression or expansion work, while the term involving accounts for viscous dissipation. Finally, the ideal gas equation (equation of state) can be used to relate the thermodynamic quantities, i.e.,

(20)  

where is the gas constant. Recall that the equation of state is one of the few truly “handy” equations in aerodynamics.

Conservation Form

For practical use, the Navier-Stokes equations are usually rewritten more compactly in the so-called conservation form, compared to the tensor differential form given previously, i.e.,

(21)  

where is the conserved variable vector and , , and are the flux vectors. The flux vectors express the rate at which mass, momentum, and energy are transported at any flow point. This equation is often expanded into a form with the viscous terms on the right-hand side to give

(22)  

where , , and are the fluxes resulting from the viscosity of the flow.

This form of the governing equations allows a mathematical reduction to a simpler form called the thin-layer Navier-Stokes equations, where the viscous derivatives with respect to the and directions are neglected to give

(23)  

for a suitably defined viscous flux, , that accounts for the retained derivatives.

The thin-layer Navier-Stokes equations are particularly useful when the viscous gradients normal to a surface dominate those in the streamwise and spanwise directions, thereby reducing computational effort and modeling complexity. Engineers often apply these simplified equations in boundary-layer-type analyses and other situations where viscous diffusion is primarily important in the wall-normal direction.

Reynolds-Averaged Navier-Stokes (RANS) Form

The RANS equations are a form of the Navier-Stokes equations used to describe turbulent flows through statistical averaging. The key idea behind the RANS equations is to decompose the instantaneous quantities (such as velocity and pressure) into their mean (time-averaged) and fluctuating components.

The Navier-Stokes equations for an incompressible fluid can be written as

(24)  

and

(25)  

where represents the external body force per unit mass.

In the Reynolds decomposition, the instantaneous velocity and pressure are decomposed into mean and fluctuating components, e.g., for the component of velocity, then

(26)  

and for the static pressure

(27)  

Here, and are the mean (time-averaged) components, while and are the fluctuating components, as illustrated in Figure 2.

When the Navier-Stokes equations are averaged over time, the nonlinear term introduces additional terms that account for the Reynolds stresses, which are the stresses created by the turbulent fluctuations. The averaged form of the continuity and momentum equations then becomes

(28)  

and

(29)  

The term represents the divergence of the Reynolds stress tensor, which arises from the turbulent fluctuations. This term is often modeled using various turbulence models, which will be discussed later, because it cannot be computed directly from the mean-flow variables.

The RANS equations provide a CFD method that balances computational cost and accuracy, making them a popular choice for many practical fluid dynamics problems. RANS methods are widely used in engineering and scientific applications to predict and analyze turbulent flows across a range of contexts, including aircraft and automobile design, ship hull design, weather forecasting, and environmental engineering, such as pollutant dispersion. Several CFD software packages are available to solve the RANS equations, including ANSYS Fluent, OpenFOAM, COMSOL Multiphysics, and STAR-CCM+.

Gridding of the Computational Domain

Numerous algorithmic advances have been made in the numerical solution of the Navier-Stokes equations. A fundamental issue in any flow problem involving a body is generating a grid structure that represents the body’s geometry for solving the governing equations. This means that the computational domain should be as small as possible for cost considerations, yet large enough to avoid nonphysical effects that can lead to excessive memory requirements.

Several types of grid systems are appropriate for analyzing flow problems; a summary is shown in Figure 3. Grids are typically classified by their geometric shape, such as O-, C-, and H-type grids. O-grids and C-grids are often used to analyze flows around wings because they conform to the curved airfoil surface, a property known as body-fitted grids. H-type grids are also used; they allow a compromise between the advantages of a Cartesian grid geometry and the ability to refine the grid near a surface to resolve the boundary layer. Nodes are the points of intersection of the grid, forming cells. The numerical solution to the flow properties is then obtained using finite-difference or finite-volume methods applied to the nodes and cells.

As shown in Figure 4, creating a CFD grid for a wing or other lifting surface involves accurately defining the geometry and refining the grid to capture flow details. It must include fine surface meshes to resolve boundary-layer flows, denser grids near the wing tip to accurately capture the formation and roll-up of wing-tip vortices, and smooth cell transitions to prevent numerical discretization issues. Adaptive Mesh Refinement (AMR) may better capture transient (time-dependent) phenomena. This approach ensures an accurate representation of aerodynamic characteristics.

Overset or chimera grid systems are locally refined grids superimposed on a primary or background grid, as shown in Figure 5, analogous to adding a grid patch where and when needed. These overlapping grids help resolve local flow phenomena without gridding the entire flow domain at high resolution, thereby saving computational effort. They can also be used to address specific difficulties in dealing with a domain best represented by separate grids that move relative to one another, such as propellers and rotors with respect to an airframe.

With chimera grids, interpolation methods are required to transfer flow-property information between grids. The accuracy of these procedures must be carefully examined to ensure the overall numerical solution remains accurate. Such grids can also be very labor-intensive to create. Unstructured grids, composed of irregularly distributed cells or elements, as illustrated in Figure 6, are increasingly being used. In two-dimensional calculations, these elements are often triangular or quadrilateral; in three-dimensional calculations, they may include tetrahedral, hexahedral, prismatic, or polyhedral cells. They have the advantage that semi-automatic mesh-generation software can create this grid type around complex surface geometries. Such grids also have the advantage of achieving a more efficient distribution of grid points, thereby accelerating the numerical solution.

Turbulence Modeling

Modeling turbulence in most flow problems often requires introducing empirical “closure” models, many of which have been developed. Indeed, resolving all turbulence scales requires extremely fine grids with precise numerical discretization, demanding large amounts of memory and the fastest computers. Calculating turbulence remains complex even for the most straightforward flows and is impractical for flows around a complete airplane.

The viscous and turbulent shear stresses can be expressed in the form

(30)  

where is the dynamic viscosity. This formulation implies that

(31)  

Here, is known as the eddy viscosity, representing an effective viscosity resulting from the effects of turbulence. Therefore, turbulence models that aim to simulate variations in eddy viscosity are called eddy-viscosity models. In these models, higher levels of turbulent stresses are simulated by augmenting the effects of molecular viscosity with a turbulent viscosity coefficient, i.e.,

(32)  

This preceding result assumes that the turbulent stress tensor is proportional to the mean strain-rate tensor, with a proportionality constant, an approach known as the Boussinesq hypothesis. It also assumes that the exchange of turbulent energy in the eddy-cascading process is analogous to the effects of molecular viscosity, as illustrated in Figure 7. Many experiments and observations of turbulent flows support this behavior. Therefore, the Boussinesq hypothesis implies that turbulent flow can be modeled as an effective total viscosity that varies with local flow conditions.

Therefore, with the addition of turbulence in the form of an eddy viscosity, the Navier-Stokes equations can be modified into the form

(33)  

Eddy viscosity models belong to a class of so-called “turbulence models” or “turbulence closure models,” which provide a mathematical basis for solving the mean flow field and the turbulent flow quantities. Numerous turbulent flow models are available, particularly for use in computational fluid dynamics (CFD), and selecting the most suitable model requires careful consideration of the specific problem. The most appropriate way to represent turbulent viscosity and select a turbulence model depends on several factors, including the specific flow application, the desired level of numerical accuracy, and the available computational resources.

Deciding which turbulence models are likely to yield good predictive performance depends on understanding the length and time scales of the specific flow problem at hand. It should be noted that turbulent flows exhibit a wide range of spatial scales, each of which plays a vital role in the flow’s overall behavior. The spatial scales of turbulence can be broadly classified into three categories: large, intermediate, and small eddies, as illustrated in Figure 8. Therefore, there is no “one size fits all” turbulence model. The different scales of turbulence interact in complex ways, and a wide range of phenomena, including vorticity, turbulence intensity, and Reynolds stresses, characterize their behavior. Modeling these scales is crucial for predicting the overall behavior of turbulent flows, and it remains a key challenge in aerodynamic modeling.

1. Large-scale eddies are the most significant structures in a turbulent flow. They transport bulk momentum and energy, and can cause substantial fluctuations in the local flow properties.
2. Intermediate-scale eddies: These eddies are smaller than the large-scale eddies, but larger than the smallest scales. They are essential in transferring energy from the large-scale eddies to the smallest scales, which takes place through the “energy cascade” process.
3. Small-scale eddies are the smallest structures, with linear scales typically on the order of millimeters or less. They dissipate energy, converting the flow’s kinetic energy into thermal energy through viscous shear.

The Boussinesq hypothesis is commonly employed to develop turbulence models for solving the Reynolds-averaged Navier-Stokes (RANS) equations and simulating turbulent flows in practical engineering applications. This hypothesis offers a computationally efficient method for incorporating turbulence effects into the Navier-Stokes equations, thereby enabling the prediction of the mean flow and associated turbulence statistics. However, in more complex turbulent flows, the Boussinesq hypothesis may need to be modified or replaced with more advanced turbulence models, given that turbulence is often anisotropic.

Various turbulence models are available, each with its own assumptions and levels of complexity. Some commonly used turbulence models for RANS simulations include:

1. The Spalart-Allmaras (SA) turbulence model is a one-equation model that directly solves for the eddy viscosity. It is known for its simplicity and computational efficiency and is widely used, especially as a benchmark.
2. k-epsilon or – models solve for two transport equations, one for the turbulent kinetic energy and another for the turbulent dissipation rate . Variations of this model include the RNG – and the Realizable – models.
3. k-omega or – models also solve two equations but use different variables: turbulent kinetic energy, , and the specific dissipation rate, . Variations include the SST – and BSL – models.
4. Reynolds Stress Models (RSMs) solve modeled transport equations for the Reynolds stresses, rather than using an eddy-viscosity assumption to relate them directly to the mean strain rate. They can provide a more detailed representation of anisotropic turbulence and are often used for flows with strong curvature, separation, or recirculation. However, RSMs are computationally more expensive than simpler eddy-viscosity models and still require closure assumptions.

In addition, two other methods for solving the Navier-Stokes equations are worth mentioning: Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS).

1. Large Eddy Simulation (LES) is not a RANS model, but it is worth mentioning. LES focuses on capturing the larger, energy-containing eddies while modeling the effects of the smaller, dissipative eddies. It is computationally expensive but yields more accurate results for complex flows.
2. DNS stands for Direct Numerical Simulation, a computational technique used in fluid dynamics to simulate turbulent flows. In DNS, all scales of turbulence are resolved numerically without the need for a turbulence closure model. This means that even the smallest eddies and flow fluctuations are explicitly simulated. However, these simulations are limited to elementary geometries and small geometric scales.

The choice of turbulence model depends on the specific characteristics of the simulated flow, the available computational resources, and the desired level of accuracy. However, the lack of generality of such models across all flow types limits their ability to predict turbulent flow properties, particularly the development of turbulent boundary layers. Developing better turbulence models remains a goal in modern fluid dynamics, and all such models must be used with awareness of their calibration range, assumptions, and empirical content. It is common for engineers and researchers to perform sensitivity analyses using different turbulence models to assess their impact on CFD results and to compare them with available benchmark solutions. These benchmarks may encompass both experiments and other theoretical approaches.

Numerical Methods

Numerical techniques for solving the Navier-Stokes equations can be based on finite-difference, finite-volume, finite-element, or related approaches, and numerous algorithms have been developed, each with its own relative advantages. Finite-difference methods provide a direct way to approximate the time and space derivatives in the equations using Taylor series expansions at each node of the computational grid, which is the intersection of grid lines.

Finite-Difference Methods

Leonhard Euler popularized the finite-difference method in the late 18th century. Consider some quantity, , representing a generic scalar field that varies with space and time. For example, could represent pressure, a velocity component, or temperature. One finite-difference approximation to a spatial derivative can be written as

(34)  

which is called a forward difference. The value of is the spatial discretization. The best way to understand this form of discretization is to use a stencil, as shown in Figure 9. For the forward-difference approximation in Eq. 34, the stencil uses the two neighboring values and . More general finite-difference stencils may also use values such as and , depending on the derivative and differencing scheme being approximated.

A backward difference approximation to the derivative is given by

(35)  

The backward and forward difference approximations are first-order accurate, i.e., the error is of the order . A central difference approximation to the derivative is

(36)  

which is second-order accurate, i.e., the error is of order . Higher-order differencing schemes can improve accuracy at the expense of increased computational cost and arithmetic operations.

Consider, for example, the numerical solution of a differential equation of the form

(37)  

where is some constant. If = (temperature), then Eq. 37 is called the heat equation, and models the movement of heat energy through space and time. It provides a useful test case for numerical methods because an exact solution can be obtained, albeit in series form.

Using a forward difference scheme for the time derivative gives

(38)  

where is the temporal discretization. Using a central difference scheme for the second spatial derivative gives

(39)  

Combining the terms leads to

(40)  

and after rearrangement, then

(41)  

Therefore, the new value of at the next time step is determined from the values at the previous time step .

The numerical process also requires the specification of initial and boundary conditions (if any). After that, the stencil can be applied to obtain the values of over the entire grid. The idea is shown in Figure 10. Applying the moving stencil to the first row of initial values will give the solution in step 1. These values can then be used to find more values in step 2. By moving the stencil to the left side of the grid, the boundary condition values, in this case = , can be populated. Applying the stencil further up the time grid and across physical space yields the values at step 4. The process can then be continued to cover the entire grid.

Finite Volume Methods

Finite volume methods operate on an integral form of the Navier-Stokes equations, derived from Eq. 21, and treat the computational domain as cells rather than nodes. The integral of mass, momentum, and energy within each cell is divided by the cell volume to obtain an average value of the conserved variable. The average fluxes are derived similarly by integrating over the cell boundary. This latter approach ensures the conservation laws are satisfied across each control volume or cell. Fluxes of the conserved quantities are calculated across the boundaries. The governing equations are converted into algebraic equations by integrating over each cell, and the resulting system is then solved to obtain the flow variables.

Consider again the heat equation given in Eq. 37. Divide the domain into control volumes, each with width . Integrating over a control volume gives

(42)  

Applying the divergence theorem to the spatial term gives

(43)  

which in one dimension simplifies to

(44)  

The gradient at the faces of the control volume can be obtained using central differences, i.e.,

(45)  

Substituting the approximations into the integral form gives

(46)  

which can be simplified to obtain

(47)  

To solve for , this latter equation must be integrated in time. For example, using an explicit forward difference method gives

(48)  

There are several compelling computational reasons to adopt a finite-volume approach rather than a finite-difference approach. One reason is that classical finite-difference methods are most naturally applied on structured grids, whereas finite-volume methods can be applied to both structured and unstructured grids. For simple structured-grid problems, finite-difference and finite-volume discretizations may yield the same algebraic stencil, as in the one-dimensional heat equation example above. In general, however, the two approaches differ in their formulations, i.e., finite-difference methods approximate derivatives directly, whereas finite-volume methods enforce conservation within each control volume.

Numerically Solving the RANS Equations

RANS methods are prevalent because they give good predictive accuracy at a reasonable cost. Solving the RANS equations numerically involves discretizing the governing equations into a set of algebraic equations. The RANS equations consist of the continuity and momentum equations, which are typically solved using finite-volume or finite-difference methods.

1. The first step is to divide the computational domain into a fine grid of finite-volume cells. The Navier-Stokes equations are then discretized over these cells using numerical methods. The numerical method must also be discretized into time steps to resolve the flow’s temporal evolution.
2. Apply the conservation equations (continuity and momentum) to each control volume. This approach involves discretizing the spatial and temporal derivatives in the equations. Typical schemes for solving the flow velocities and pressures include explicit and implicit time-stepping methods.
3. The resulting discretized equations can be assembled into a system of linear algebraic equations, as outlined previously for the FVM. This system is solved iteratively at each time step until convergence is achieved. Implicit methods solve a system of equations at each time step, whereas explicit methods update the solution using the current values.
4. Implement the selected turbulence model within the selected numerical framework. As mentioned earlier, this step involves solving additional equations for the turbulence quantities using a turbulence model of your choice.
5. Apply appropriate boundary conditions to model the problem’s physical boundaries. This approach may include specifying inlet, wall, and outlet conditions.
6. Iterate through the time steps until a steady state or a time-accurate solution is reached. Use a pressure-velocity coupling algorithm, such as the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE), or its modified form, SIMPLEC.
7. Monitor convergence and adjust the solution until the changes between iterations are within acceptable tolerances or other criteria.
8. Analyze and post-process numerical results to extract relevant information, including flow velocities, pressure distributions, Reynolds stresses, turbulent kinetic energy (TKE), and other turbulence quantities.

The examples shown in Figure 11 illustrate RANS predictions of the wing-tip vortex roll-up, computed with two different turbulence models. Although both solutions exhibit similar bulk-flow properties, they yield different quantitative details and turbulence levels. Which one is “correct” can be established only with reference to experiments, assuming measurements are available.

Converting the RANS equations to algebraic form arises from discretizing the partial differential equations. Consider, for example, a two-dimensional, steady, incompressible flow, i.e.,

(49)  

Using forward or upwind differencing for the convective term gives

(50)  

where , , , and represent the convective fluxes in the -direction at the eastern, western, northern, and southern faces of the control volume , respectively. The basic concept is shown in Figure 12.

Using central differencing for the diffusive term gives

(51)  

Using central differencing for the pressure gradient term gives

(52)  

Combining these terms for the component of the momentum equation gives

(53)  

Similar steps would be followed for the and components of the momentum equation and for other control volumes in the domain for three-dimensional flow problems.

This type of numerical process is typically implemented in commercial CFD software packages. These packages usually include preprocessing tools for mesh generation and turbulence model setup, as well as postprocessing tools for visualizing and analyzing the results. The choice of numerical method and solver depends on factors such as the complexity of the flow, available computational resources, and desired resolution and accuracy.

Flowfield Predictions

While CFD predictions of the flow around an entire aircraft remain challenging, as shown in Figure 13, confidence in such predictions for lift and drag can be high when the methods are carefully validated. The resolution of local flow properties may not always be sufficient to resolve every turbulence scale, but it can still be useful for design refinements, e.g., tailoring the wing’s shape in the presence of the engines. Remember also that many turbulence models used in CFD are partly postdictive, in the sense that they contain empirical assumptions and calibrations based on observed turbulent flows. They can be extremely useful engineering models, but they should not be treated as universal first-principles descriptions of turbulence. RANS solutions predict the mean flow field and model the effects of turbulence on that mean flow; they do not resolve the instantaneous turbulent eddies themselves. Nevertheless, remarkable details of the flows around complete aircraft can now be computed. Ultimately, however, all CFD solutions must be considered tentative and validated against experiments or another benchmark.

CFD predictions of flows around rotating-wing aircraft are much more challenging because of their greater complexity and the need to use simultaneously moving grids (for the rotors) and stationary grids (for the airframe). A CFD visualization of the flow of NASA’s six-passenger tilt-wing concept for urban air mobility in hover or “helicopter mode” and forward flight or “airplane mode” is shown in Figure 14. The results reveal tip vortices on the rotor blades, characterized by their vorticity (spin). In hover, each rotor generates intertwining helicoidal vortices with relatively low pitch, which convect beneath the rotor and quickly diffuse and break up to form a turbulent jet. The proximity of wake vortices to the rotor results in high induced power requirements to drive the rotor. In forward flight, the rotors produce helical wake vortices with a much greater pitch angle, resulting in lower induced effects and reduced power requirements.

These results reveal the complexity of the flow in a tilting multi-rotor configuration, in which flows from multiple rotors interact with one another, the wing, and the fuselage. Such RANS solutions can be beneficial, for example, in identifying issues with rotor/airframe aerodynamic interactions, which can produce highly variable aerodynamic forces and moments on the aircraft during transitions between hover and forward flight.

Euler Equations

Suppose that turbulence and viscous effects do not need to be resolved, assuming these can be treated as weak in a given application. The problem of gridding the computational domain at extremely fine resolution has been mitigated to some extent, thereby reducing computational costs by orders of magnitude. If the flow is assumed inviscid, the Navier-Stokes equations reduce to the simplified set known as the Euler equations.

By removing the viscous terms from the full Navier-Stokes equations, the Euler equations can be written in tensor differential form in Cartesian coordinates as

(54)  

(55)  

(56)  

Alternatively, by explicitly dropping the viscous terms from Eq. 22, the Euler equations can be written in conservation form as

(57)  

In practice, numerical solutions to the Euler equations have been found to provide excellent approximations to high Reynolds number flows around aircraft. Although the viscous terms are explicitly neglected, numerical solutions to the Euler equations can provide valuable insight into complex three-dimensional flows.

Vorticity in an inviscid flow, such as that governed by the Euler equations, should convect freely without physical viscous diffusion. However, a numerical solution of the Euler equations may still introduce artificial dissipation and diffusion of vorticity. This behavior arises from the discretization scheme, grid resolution, flux treatment, and any added numerical dissipation used to stabilize the calculation. For instance, a wing modeled with the Euler equations can produce lift only if circulation is established or imposed through the geometry, boundary conditions, and an effective Kutta condition. In a numerical Euler solution, some artificial viscosity or numerical dissipation is usually present because of the discretization scheme. Even a small amount of this numerical dissipation can help select the physically relevant lifting solution, although it is not equivalent to resolving the real viscous boundary layer.

Despite the limitations of the Euler equations, Euler-based CFD methods have helped the aircraft industry gain new insights into complex flow problems and have clarified the limitations of simpler models. Euler-based methods can provide a relatively straightforward approach for predicting many features of the flow around an aircraft, but their fundamental limitations should still be recognized.

Boundary Layer Equations

At the Reynolds numbers typical of many aircraft, the viscous terms in the Navier-Stokes equations are significant only in a thin shear layer, called the boundary layer, near the vehicle’s surface, as shown in Figure 15. Indeed, in many engineering problems, the flow details away from the surface are of limited interest, and only the effects of the boundary layer are needed, for example, because it generates most of the drag experienced on the surface. To model the flow in this layer, a simplified subset of the Navier-Stokes equations can be derived by exploiting the fact that the length scale of flow-variability is much smaller in directions normal to the surface than in directions parallel to it.

The advantage of a boundary-layer modeling approach is that it allows an inviscid “outer” solution, calculated, for instance, using the Euler equations or a potential-flow method, to be modified to account for the effects of viscosity near a surface. This approach leads to a class of numerical techniques known as zonal or viscous-inviscid interaction methods. Such methods have proven very useful in airfoil design and can also help predict the aerodynamics of an entire airplane.

The basic principle first solves the inviscid flow in the “outer” zone, then solves the viscous flow in the “inner” zone using the boundary-layer equations. The outer solution is then repeated after applying a correction to consider the effects of the boundary layer on the outer flow; the resulting solution yields, to first order, a high-Reynolds-number solution to the Navier-Stokes equations. Care must be taken when encountering separated flows, as the boundary-layer assumption becomes invalid under these conditions.

Consider this type of situation in two dimensions. If is in the direction parallel to the surface and is normal to the surface, then, as the Reynolds number of the flow becomes large, the Navier-Stokes equations for an incompressible flow reduce to

(58)  

Because the gradients of the fluid properties in the direction become much smaller than those in the direction, then

(59)  

Because the pressure across the boundary layer (i.e., normal to the surface) is constant to boundary-layer order (Eq. 59), the pressure at any point along the surface is taken to be the pressure just outside the boundary layer. Therefore, the inviscid momentum equation just outside the boundary layer can be written as

(60)  

where is the “edge” velocity just outside the boundary layer. At the edge of the boundary layer , where is the boundary layer thickness^[8]

The variation in pressure can now be eliminated, and the momentum equation within the boundary layer becomes

(61)  

Along with the boundary layer equations is the concept of a displacement thickness, , which is defined as

(62)  

This quantity allows the shape of the surface to be displaced away from the original surface by , yielding the same flow rate as in the viscous flow with the boundary layer present. This is equivalent to saying that the outer inviscid flow around the body can be obtained by solving the equations for a new shape obtained by adding the displacement thickness to the body’s actual shape. This approach exemplifies a so-called weak-interaction method. Although still an approximation, this approach yields more representative flow solutions than those obtained without including viscous effects.

The surface shear stress may also be required to estimate the viscous drag. This quantity can be obtained from the integral form of the boundary-layer equations, known as the momentum integral equation, usually attributed to von K\'{a}rm\'{a}n. In this case, the governing equation is

(63)  

where is the shear stress on the wall (surface) and is the momentum thickness, given by

(64)  

This equation is general because it makes no assumptions about the relationship between and the velocity gradient, and thus applies to both laminar and turbulent boundary-layer flows. The process is completed by choosing a suitable form for the velocity profile across the boundary layer, which allows the momentum integral equation to be solved to find the key variables describing the properties of the layer, i.e., the values of , , , and . While solutions to the momentum integral equation can be found in closed form in some cases, numerical solutions are usually obtained in most cases.

Potential Flow Equations

A further simplification of the governing equations of the flow can be achieved by assuming irrotationality, i.e., . The irrotationality condition is a good assumption for external flows away from surfaces. In this case, the velocity field can be expressed as a gradient of a potential function, , such that

(65)  

In potential flow theory, a fluid is assumed to be inviscid and irrotational. However, irrotationality does not imply the absence of viscosity. The potential equations can be expressed in many ways, but one form commonly used to model compressible flows is

(66)  

where . This form of the potential equation is usually called the full-potential equation because it retains the nonlinear terms that arise in compressible potential-flow theory. It is still based on important assumptions, especially inviscid, irrotational, and usually isentropic flow, so it is not equivalent to solving the full Euler equations. For many flow problems, these nonlinear terms are important if the characteristic variation of the flow properties with Mach number is to be represented properly. However, other simplifications are possible.

If small disturbances can be assumed, some terms are of higher order and can be neglected. In this regard, the transonic small disturbance (TSD) equations are a subset of the full potential equations and can be written in the form

(67)  

For steady flows, the TSD equation is

(68)  

or in terms of the generalized coefficients , , and , then

(69)  

where

(70)  

The effects of the airfoil shape, body motion, and far-field disturbances enter the transonic small-disturbance problem mainly through the boundary conditions rather than as arbitrary source terms in the differential equation. For example, for a thin airfoil with surface ordinate , the small-disturbance flow-tangency condition may be written in the approximate form

(71)  

or, if the perturbation potential has been nondimensionalized by , as

(72)  

Unsteady pitching, plunging, gusts, or other imposed disturbances can likewise be introduced through appropriate time-dependent boundary conditions and far-field conditions. In this sense, such effects force the solution through the boundary-value problem, even though they do not necessarily appear as explicit source terms in the governing TSD equation.

The Transonic Small Disturbance (TSD) equation is significant in the history of computational fluid dynamics (CFD). It was one of the first flow equations that researchers such as Antony Jameson^[9] solved numerically using finite-difference methods, which helped establish CFD as a practical field of study. Naturally, the TSD equation is fundamental to the study of transonic flows, where speeds are near the speed of sound, posing unique challenges as a result of the coexistence of subsonic and supersonic regions in the flow field. Solving the TSD equation numerically provided insights into these complex flow regimes and laid the groundwork for more advanced CFD techniques and software used today.

If the nonlinear transonic terms are neglected, the governing equation reduces to a linear small-disturbance potential equation. In the special acoustic limit of small disturbances in a quiescent medium, this equation takes the familiar wave-equation form

(73)  

where is the freestream speed of sound. This simplified equation describes the propagation of small pressure or velocity-potential disturbances at the speed of sound. For aerodynamic small-disturbance theory about a finite freestream, however, the linearized equation also contains mean-flow convection effects. The full TSD equation includes additional nonlinear terms needed to represent the mixed subsonic-supersonic behavior characteristic of transonic flow.

A final simplification is to assume steady, incompressible flow (), in which case Laplace’s equation governs the velocity potential, i.e.,

(74)  

A feature of this latter equation is its linearity, which allows complex potential flows to be synthesized by combining two or more elementary potential-flow solutions. For example, if and are solutions to Laplace’s equation, then is also a solution. This property enables the construction of complex potential-flow patterns by combining simpler solutions, such as uniform flow, sources, sinks, and vortices. Laplace’s equation is the governing equation for incompressible, irrotational flow expressed in terms of a velocity potential and underlies surface-singularity, or “panel,” methods.

Laplace’s equation also applies to unsteady incompressible, irrotational potential flow. In this case, the unsteady nature of the flow is accounted for in the potential function , which now depends on time , while the spatial potential still satisfies at each instant. The instantaneous adaptation of the incompressible potential-flow field to changes in boundary conditions results from the idealization of an infinite speed of sound, whereby pressure disturbances propagate instantaneously throughout the flow field. This is an idealization of the incompressible model, not a property of real compressible fluids.

Surface Singularity Methods

Panel methods are widely used to calculate potential flows around bodies, particularly when the flow can be approximated as inviscid and irrotational. These methods offer valuable insights into aerodynamic characteristics, especially pressure distributions, lift, and induced effects. They are instrumental in preliminary design stages and in situations where viscous effects can be neglected or accounted for by separate boundary-layer or empirical corrections. Panel methods are based on incompressible irrotational flow, for which the governing equation is Laplace’s equation (Eq. 74) for the velocity potential. The idea is that surface singularities of unknown strength can be distributed over the body surface. After determining the strengths of the singularities by imposing a flow-tangency condition on the surface, the surface velocities and, consequently, the inviscid aerodynamic loading on the body can be computed.

Large computer codes that use panel methods are commercially available and widely used in the aerospace industry. Some well-known examples of codes that use panel methods for aerodynamic analysis include XFOIL, a widely used tool for predicting flows around airfoils of arbitrary shape. VSAERO is used to analyze subsonic, transonic, and supersonic aerodynamics using panel methods. PANAIR, developed by NASA, is widely used to predict the aerodynamic characteristics of complex configurations. These codes typically provide comprehensive functionalities for setting up geometries, applying boundary conditions, running simulations, and analyzing results, making them indispensable tools in aerodynamic analysis and design in the aerospace industry.

To illustrate the principles, consider the flow of a two-dimensional airfoil. The airfoil surface is replaced by a distribution of discrete panels, as shown in Figure 16, and each panel contains a singularity distribution that can be described by a velocity potential function of unit strength  . This means that the total potential induced by all the singularities in the flow is

(75)  

where the constants must be found subject to the boundary condition that the flow must be tangential to the surface at discrete collocation points on the surface, i.e.,

(76)  

where is the outward pointing normal unit vector at the th collocation point. Typically, collocation points are selected at the midpoints of the panels.

Suppose the panels contain vortex singularities that vary linearly in strength over the panel length. Then, with panels, there are flow-tangency equations containing unknown vortex strengths. The Kutta condition at the airfoil’s trailing edge provides the auxiliary equation required to solve the problem uniquely and define the airloads on the airfoil, i.e., , which is implemented by setting . Unless the Kutta condition is imposed in some form, the solution to the panel strengths becomes indeterminate. Therefore, the equations governing the problem can be expressed as

(77)  

where is the influence coefficient giving the normal velocity induced at collocation point by a unit-strength singularity on panel . An influence coefficient is effectively the normal component of the velocity induced at the collocation point on panel by the unit strength singularity distribution of on panel ; it depends on the relative distance between and orientation of the two panels. The assembled matrix has dimensions after the Kutta condition is included.

Notice that solving Eq. 77 involves finding the inverse of the matrix of influence coefficients, for which many efficient algorithms are available. In addition to calculating the influence coefficients, this inversion is the primary computational expense in the panel method. If the angle of attack, , is included in the vector on the right-hand side of Eqs. 77, the influence coefficients need only be calculated once. Once the singularity strengths are known, the airfoil loading immediately follows from the singularity strengths.

Panel methods have also been used to model the aerodynamics of non-lifting bodies. In such an application, the fuselage surface may be represented by small, flat quadrilateral panels, each with singularities consisting of sources, sinks, and/or doublets. The principles are very similar to those in the two-dimensional case, but the size of the numerical problem increases significantly. Assume that only source singularities are used for illustration. If the singularity strength on the th panel is , then the velocity induced by this panel at the control point on the th panel can be written as . At the th control point, the component normal to the surface of the velocity induced by all sources is then

(78)  

The boundary condition of flow tangency is imposed at the control point of each panel, allowing the unknown singularity strengths, , to be evaluated. For a uniform freestream , the component of the freestream velocity normal to panel is , where is the unit normal vector to panel . Therefore, at the th panel

(79)  

One equation of this type applies to each of the panels.

When applied to all panels, this latter approach yields a set of linear simultaneous equations that can be solved for the singularity strengths, , using standard linear-algebraic methods. The governing equations can be written in matrix form as

(80)  

which may be written formally as

(81)  

Therefore, the numerical task is to solve a dense linear system of algebraic equations derived from the influence coefficient matrix. In practice, the matrix is not usually inverted explicitly; instead, direct or iterative linear-system solvers are used. Once the values of are obtained, the velocity components tangential to the panels can be determined. Hence, the pressure distribution over the fuselage can be calculated using Bernoulli’s equation.

Typically, thousands of panels are required to accurately represent the complex shape of an entire aircraft. Care must be taken, for instance, to concentrate a sufficient number of panels in regions of higher surface curvature, as also illustrated in Figure 16. For a conventional dense direct solve, the numerical cost of solving the panel-method system scales approximately as , so their use is by no means inexpensive, even on a modern computer. However, panel methods are still considerably less computationally intensive than most CFD solutions based on the Navier-Stokes or Euler equations, and accelerated or iterative implementations can reduce the practical cost for large problems.

One advantage of panel methods is that they avoid the grid-based numerical diffusion and dissipation that can occur in finite-difference or finite-volume field methods. However, the downstream wake, as shown in Figure 17, must still be represented by an approximate wake model, so its accuracy depends on the wake discretization, convection assumptions, and any roll-up model used. If sufficient care is taken to distribute panels appropriately across the fuselage surface, with more panels assigned to regions of expected higher pressure gradients, then a three-dimensional surface-singularity method can be used to obtain credible predictions of the inviscid loading on an airframe, which is sufficient for preliminary design.

As an enhancement to the primary panel method, the solution of the momentum integral equation from boundary-layer theory enables displacement-thickness corrections to panel-method solutions, thereby modeling attached viscous boundary-layer effects, the principle being shown in Figure 18. Flow separation requires additional modeling, such as imposing a separation point and representing the separated region with a free shear layer.

Including boundary-layer displacement effects in lift predictions for an airfoil improves accuracy by accounting for the effective change in body shape caused by the viscous boundary layer. Classic potential-flow methods, which ignore boundary-layer growth, may overpredict lift or misrepresent the pressure distribution, especially as adverse pressure gradients become important. By incorporating displacement-thickness corrections, the influence of attached viscous effects on the outer inviscid flow is better accounted for, yielding more accurate lift predictions, as shown in Figure 19. However, the onset and extent of flow separation require additional modeling and cannot be predicted solely from a displacement-thickness correction.

Moreover, the boundary layer affects the pressure distribution over the airfoil by altering the effective shape seen by the outer inviscid flow, with the displacement thickness providing one measure of this effect. These changes influence the flow field and, consequently, the lift generated. Including boundary-layer displacement effects improves the predicted lift by better reflecting the interaction between the viscous boundary layer and the outer pressure distribution. If a transition model or a prescribed transition location is also included, laminar-to-turbulent transition effects can be represented, further refining lift and drag predictions.

Rarefied Gas Flows

Detailed information on low-density or rarefied-flow environments is crucial for spacecraft, spanning design through mission operations. Models of such environments must be developed for various applications, including assessing orbital decay and aerodynamic heating during re-entry into Earth’s atmosphere or during atmospheric entry on other planets. A rarefied gas, characteristic of the Earth’s outer atmosphere near space, consists of molecules that are relatively far apart and collide infrequently compared with ordinary continuum-flow conditions. The continuum hypothesis can become inaccurate in these conditions, necessitating alternative flow descriptions based on statistical mechanics. Rarefied-gas behavior typically occurs at very low pressures or high temperatures, where the mean free path of the molecules is no longer negligible compared with the body’s characteristic dimensions, as illustrated in Figure 20. The value of the Knudsen number, , helps determine the degree of rarefaction, i.e.,

(82)  

where is the mean free path, and is a characteristic length for the problem. If the value of is sufficiently large, continuum assumptions begin to break down. In practice, slip-flow effects may appear for , transition-flow effects occur over intermediate values of , and free-molecular behavior is approached for .

Over time, if there are no external forces, the distribution of these rarefied gas molecules will settle into a predictable pattern, i.e., a form of equilibrium. In equilibrium, the velocities of gas molecules are described by the Maxwell-Boltzmann distribution, an idealized yet realistic probability distribution, as illustrated in Figure 21. This model describes the probability of finding an individual molecule at a given temperature with a particular speed. The speed most molecules have is given by

(83)  

where is the Boltzmann constant, is the absolute temperature, and is the mass of the gas molecule.

The Boltzmann constant relates the average kinetic energy of particles in a gas to the gas’s temperature; its value is approximately Joules per Kelvin (J/K).  At lower temperatures, the molecules have less energy. Therefore, molecular speeds are lower, and the distribution is narrower. As the molecular temperature increases, the probability distribution broadens and shifts toward higher molecular speeds because the molecules have higher average kinetic energy. So more of the molecules, on average, are moving faster. The average speed of all molecules is

(84)  

and the corresponding root mean square (rms) speed, a measure of the typical speed of the molecules, is

(85)  

In thermal equilibrium, the velocities of gas molecules follow the Maxwell-Boltzmann distribution, characterized by the most probable, average, and root-mean-square speeds. The distribution changes when the gas is not in equilibrium, for example, when it is flowing or subjected to external forces. A simplified relaxation-time form of the Boltzmann equation, often called the BGK model, describes how this distribution evolves and can be written as

(86)  

where is the molecular velocity distribution function, is the molecular velocity, is the equilibrium distribution, and is the relaxation time that characterizes the rate at which molecular collisions drive the distribution back toward equilibrium. In a rarefied gas, collisions occur less frequently because the molecules are farther apart. Collisions between molecules drive the gas toward equilibrium, and models help account for these effects.

Rarefied-gas models are often used to predict the aerodynamics and kinetic heating generated by spacecraft during atmospheric entry, which differ from those predicted by continuum models. An example of their application is shown in Figure 22.^[10] These models are essential for understanding the behavior of gases in low-density environments where traditional continuum-based descriptions of fluid dynamics are inapplicable.

In general, rarefied-gas models provide insights into the complex interactions between gas molecules and the spacecraft. These are critical for accurately predicting aerodynamic forces, heat transfer rates, and the resulting thermal protection requirements. For instance, during the Mars Science Laboratory mission, rarefied-gas models were used to simulate the effects of the Martian atmosphere on the spacecraft as it descended and landed. Such models have enabled engineers to design effective heat shields and aerodynamic surfaces, ensuring a safe landing.

Aerospace Structures

The Finite Element Method (FEM) is a widely used numerical technique in engineering for analyzing structures and predicting their response to various physical forces, including aerodynamic forces. It involves discretizing complex structures into smaller, discrete elements connected at nodes. Depending on the problem geometry and complexity, elements such as beams, trusses, triangles, or hexahedra are employed. Similar to CFD, FEM uses meshes, but the elements are usually classified by their topology, such as beam, shell, triangular, quadrilateral, tetrahedral, or hexahedral elements. Mesh refinement is particularly useful in areas requiring higher resolution, including around structural cutouts such as windows and access panels, as shown in Figure 23. This approach ensures an accurate representation of stress concentrations and stress deformation patterns.

In practice, the FEM formulates the governing equations of the structure based on the principles of energy minimization or virtual work, thereby deriving the element stiffness matrix and the system of equations to be solved. These equations are then solved numerically by assembling a global system from the contributions of individual elements, while accounting for boundary conditions and applied loads. The versatility of the FEM enables engineers to optimize designs, assess structural integrity, and simulate real-world conditions with fewer physical prototypes. This makes it indispensable in aerospace engineering because it can handle complex geometries, diverse materials, and realistic loading conditions. FEM has become a major tool for improving the strength-to-weight ratio of aerospace structures during the design process.

History

The origins of FEM can be traced to work in structural analysis in the 1940s. Richard Courant’s 1943 work is often cited as one of the earliest instances of the Finite Element Method (FEM) being employed. Courant used piecewise-linear approximations over triangular subregions to solve torsion problems, thereby laying the groundwork for the method. The technique later gained traction through the efforts of engineers and mathematicians. The 1956 paper by Turner, Clough, Martin, and Topp on the stiffness and deflection analysis of complex structures was a foundational contribution to what became the finite element method. Ray Clough, a structural engineer, is generally credited with introducing the term “finite element method” in his 1960 paper on plane-stress analysis. Along with John Argyris and other researchers, the method was extended for more complex structural analysis problems.

During this decade, the theoretical foundation of FEM was solidified. Clough’s work, along with that of others like Olgierd Zienkiewicz, played a crucial role. Zienkiewicz’s book, “The Finite Element Method in Structural and Continuum Mechanics,” first published in 1967, has become a cornerstone textbook. Early FEM software started to appear. The Structural Analysis Program (SAP) series, developed at the University of California under the leadership of Edward Wilson, is a notable early FEM program.

FEM was extended beyond structural mechanics into other fields such as heat transfer, fluid mechanics, and electromagnetics. The development of general-purpose FEM software, such as NASTRAN (NASA Structural Analysis), also began in this era and has since become one of the most widely used FEM tools in aerospace engineering. The method continued to evolve, improving theoretical understanding and computational techniques. Commercial FEM software packages, such as ANSYS and ABAQUS, emerged and gained popularity across various engineering disciplines.

However, there are many other types of FEM, including research codes developed at universities. The overarching objective of the FEM is to incorporate efficient algorithms for large-scale problems, which inevitably requires methods optimized for various types of analyses. Today, FEM methods use parallel processing and advanced numerical techniques to enhance computational efficiency. Many can handle linear and nonlinear problems, as well as dynamic analysis, thermal analysis, and other complex tasks.

General Approach

The structural domain is divided into elements with nodes, as illustrated in Figure 24. The displacement within each element is approximated using shape functions and nodal displacements , i.e.,

(87)  

For each element, the stiffness matrix is given by

(88)  

where is the strain-displacement matrix, is the material property matrix (e.g., containing Young’s modulus, , and Poisson’s ratio, ), and is the element’s spatial domain. The use of matrix multiplication in Eq. 88 is implied.

The global stiffness matrix is then assembled from the element stiffness matrices, i.e.,

(89)  

The global force-displacement relationship is given by

(90)  

where is the global displacement vector and is the global force vector. Boundary conditions are then applied to this system of equations. For an assembly of such elements, the global stiffness matrix is constructed by linear superposition of the element stiffness matrices. Solving the global system of equations efficiently, especially for large-scale problems, involves the use of various numerical methods.

2-D Triangular Element

A simple triangular element provides an opportunity to deepen understanding of the FEM. It will become apparent that the method’s systematic, discretized nature lends itself to computer implementation, allowing the steps to be readily programmed. The shape functions are linear for a triangular element with three nodes. Let the coordinates of the nodes be , , and , as shown in Figure 25.

The shape functions, in this case,

(91)  

(92)  

(93)  

where is the area of the triangle given by

(94)  

The constants , , and are determined by the coordinates of the nodes, i.e.,

(95)  

(96)  

(97)  

The strain-displacement matrix relates the nodal displacements to the strains within the element. For this triangular element, is given by

(98)  

Because the shape functions are linear, their derivatives are constant over the element, i.e.,

(99)  

Therefore, the matrix can be written as

(100)  

For a constant-strain triangular element, the displacement field varies linearly over the element. Therefore, the strain components are constant within the element, and the matrix is constant. If the material properties are uniform, then is also constant. For a plane-stress triangular element of thickness , the integral then simplifies to

(101)  

where is the area of the triangular element. If unit thickness is assumed, then . This simplification is one reason why constant-strain triangular elements are useful pedagogically, although more sophisticated elements are often needed for high-accuracy structural analysis.

Finally, as previously described, the global stiffness matrix is assembled by summing the contributions of all element stiffness matrices based on the connectivity of the nodes. Boundary conditions are also applied to ensure that the constraints, such as fixed supports or prescribed displacements, are correctly imposed.

Total Structure Mesh Generation

Creating a mesh for complex problems, such as a complete wing or airplane structure, involves several steps. First, the structure’s geometry is carefully defined. The appropriate type of finite element (e.g., rectangular, triangular, tetrahedral, hexahedral) is then selected, and the geometry is discretized into these finite elements. The mesh is refined, particularly in areas of stress concentration, such as the wing-fuselage junction, attachment fittings, cutouts, and other regions with abrupt geometric changes, as illustrated in Figure 26. The mesh quality is verified to ensure accurate results, and boundary conditions and loads are applied accordingly. Finally, the mesh is exported for use in finite element method (FEM) simulations.

Various software tools and libraries provide functionality for creating high-quality meshes for accurate FEM simulations. The meshing process for FEM differs from that for CFD, although both require careful consideration of grid resolution to accurately predict structural response and aerodynamic behavior.

Predictions from the FEM

The FEM has become a standard tool for predicting the structural response to applied loads in engineering; an example of a FEM solution is shown in Figure 27 for a wing and the spar connected to the fuselage. NASTRAN is a widely used FEM software package that includes preprocessing (meshing and boundary conditions), solution (numerical methods), and postprocessing (visualization of results). Advances in computer technology have significantly enhanced its capabilities, allowing for the simulation of increasingly complex systems.

The FEM has also been integrated with other computational techniques, including optimization algorithms and multiphysics simulations. Ongoing research focuses on improving the accuracy and efficiency of FEM, developing adaptive methods, and integrating FEM with machine learning and artificial intelligence for predictive modeling and data-driven simulations.

Coupling CFD & FEM

CFD and FEM are often introduced as separate computational methods, but many important aerospace problems require their use together. The reason is straightforward, i.e., the aerodynamic loads on a flight vehicle depend on its shape, but the shape of a flexible structure can change under those loads. This interaction between the flow and the structure is called fluid-structure interaction (FSI). In aerospace engineering, it is closely related to aeroelasticity, including phenomena such as wing bending, control-surface deformation, buffet, flutter, and structural response to gusts.

In a coupled CFD-FEM calculation, the CFD method predicts the pressure and shear-stress distribution over the aircraft’s surface. These aerodynamic loads are then transferred to the structural finite element model. The FEM solver computes the resulting deflections, strains, and stresses in the structure. If the structural deflections are large enough to change the aerodynamic shape, then the deformed geometry must be passed back to the CFD solver, and the flow solution is recomputed. This exchange may be repeated until the aerodynamic loads and structural deflections are mutually consistent.

The basic coupling process can be summarized as:

1. The CFD model computes the aerodynamic pressure and shear stress distributions over the vehicle surface.
2. These surface loads are interpolated or mapped onto the FEM structural mesh.
3. The FEM model computes the structural deflections, stresses, strains, and internal load paths.
4. The structural deflections are transferred back to the aerodynamic surface mesh.
5. The CFD mesh is updated or deformed to match the new structural shape.
6. The process is repeated until the aerodynamic and structural solutions are consistent.

This type of coupling is especially important for lightweight aerospace structures, where even modest deflections can change the aerodynamic loading. A high-aspect-ratio wing, for example, may bend and twist under lift. The twist changes the local angle of attack, altering the lift distribution and, in turn, the structural deformation. In such cases, the aerodynamic and structural problems cannot always be treated independently.

Coupled CFD-FEM calculations may be performed in different ways. In a one-way coupled calculation, the aerodynamic loads are computed first and then applied to the structural model. The structural deformation is calculated but not fed back into the CFD solution. This approach is useful when the structural deflections are small enough not to significantly alter the flow. In a two-way coupled calculation, aerodynamic loads and structural deflections are iterated between the CFD and FEM solvers until a consistent solution is obtained. Two-way coupling is more accurate for flexible structures but also more computationally expensive.

For steady problems, the goal is usually to find a static aeroelastic equilibrium, where the aerodynamic loads and structural deformation are in balance. For unsteady problems, the coupling must also account for time. The CFD solver predicts time-varying aerodynamic loads, and the FEM solver predicts the corresponding dynamic structural response. This type of analysis is needed for problems such as gust response, control-surface oscillations, buffet, and flutter. Flutter is especially important because it involves a potentially unstable energy exchange between the airflow and the elastic structure.

A simple way to think about the coupled static problem is

The process continues until the loads and deflections no longer change appreciably between iterations. In practice, this coupling requires careful interpolation between the aerodynamic and structural meshes, conservation of total force and moment during load transfer, and robust numerical methods to prevent the coupled solution from becoming unstable.

An example of a simulation-based aeroelastic assessment for a high-aspect-ratio transport-aircraft wing is shown in Figure 28. Increasing the aspect ratio can reduce induced drag and fuel consumption, but the resulting wing flexibility can lower the flutter speed. In the example, a parametric finite element model was used with loads analysis, structural optimization, and flutter checks to show how changes in wing stiffness, mass distribution, and engine-pylon dynamics can affect aeroelastic stability. This case illustrates why CFD, FEM, and aeroelastic analyses must be integrated carefully; aerodynamic shaping alone cannot determine whether a high-aspect-ratio wing is viable.^[11]

A coupled CFD-FEM analysis is powerful, but it does not guarantee an automatic predictive solution. The result depends on the fidelity of both models, the quality of the CFD and FEM meshes, the turbulence model, the structural idealization, the boundary conditions, and the accuracy of the load-transfer procedure. Therefore, validation against wind-tunnel tests, ground-vibration tests, flight tests, or other benchmark data remains essential. Like CFD and FEM individually, coupled simulations are best viewed as engineering models that must be checked against physical evidence.

The increasing use of coupled CFD-FEM methods reflects the broader trend toward multidisciplinary analysis and design optimization in aerospace engineering. Modern flight vehicles must simultaneously satisfy aerodynamic, structural, propulsion, control, manufacturing, and operational constraints. Coupled computational methods allow engineers to study these interactions earlier in the design process and to identify problems that may not be visible when each discipline is analyzed separately.

Summary & Closure

This chapter summarizes advanced computational methods for analyzing flight vehicles, focusing on the critical role of computational fluid dynamics (CFD), primarily using the Navier-Stokes equations, in understanding their aerodynamics. Despite advances, significant numerical and practical issues remain, particularly in efficiently coupling CFD with the dynamic and elastic structural motion of airframe components. Developing adaptable turbulence closure models is also crucial for a deeper understanding of practical flight problems.

The extensive computational resources and algorithmic challenges posed by modern CFD methods, especially Navier-Stokes and Euler-based approaches, still limit their routine use for fully resolved, complete-aircraft predictions. Modeling the flow around a complete aircraft is daunting, and validating these models requires considerable effort. This goal challenges experimentalists to provide high-quality measurements of both surface and off-surface flow fields. Proper validation is crucial for assigning the confidence levels required for CFD use in design.

The Finite Element Method (FEM) is employed to analyze the structural integrity of mechanical components, using structured or unstructured grids to simulate stress distributions, deformations, and potential failure. An acceptable mesh resolution is required around stress-concentration areas to accurately predict structural responses under various loads. Boundary conditions in FEM involve applying constraints, loads, and material properties to simulate realistic structural behavior, optimize designs, and ensure adequate strength, durability, and safety margins.