Introduction

Researchers continue to develop advanced computational methods to analyze challenging problems in aerospace engineering. These methods are used in areas such as aerodynamics, structural analysis, aeroelasticity, flight dynamics, and other coupled multidisciplinary problems. In particular, computational fluid dynamics, commonly referred to as “CFD,” and the finite element method, or “FEM,” have become essential tools for predicting aerodynamic behavior and structural response, respectively. This chapter does not review these techniques or their numerical algorithms in depth; it introduces the foundational principles behind them, how they are applied, and their capabilities and limitations. Nor is this chapter meant to be an all-encompassing treatment.

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 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 faster code execution, increased memory, faster data storage access, 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.

Fundamental research on modeling and algorithmic development has advanced significantly 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. Predicting flows across the entire aircraft could open new design directions, such as reconsidering supersonic transports (SSTs), as illustrated 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, 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 solution accuracy, stability, and cost. 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, researchers have used a wide range of numerical techniques to solve a hierarchy of flow equations, 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 defining and constructing appropriate computational grids to solve 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, many other options also exist 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, or 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 called the “Father of Modern Computer Architectures.” During the 1960s, CFD evolved more rapidly, in parallel with advances in computer technology. 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 with Brook Taylor^[4] in the early 1700s, with other contributions by Leonhard Euler and Carl Runge. Numerical methods also advanced significantly, and the accuracy and stability of these schemes became increasingly crucial for solving aerospace problems.

The publication of numerous technical papers and important textbooks during the 1960s and subsequent decades 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 aerospace applications. In the mid-1980s, NASA established what became the NASA Advanced Supercomputing (NAS) Division at NASA Ames, with its dedicated supercomputing facility opening in 1987. 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 mainly during the 1980s, 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 improve the accuracy and robustness of CFD solutions. During this period, CFD also began to find applications beyond aerospace, 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, making 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 flight vehicle design to weather prediction, biomedical engineering, and environmental studies. Continued software development, increased access to computational resources, and ongoing research mean 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 among fidelity, cost, and the aerodynamic quantities to be predicted.

The Navier-Stokes equations provide the most general continuum description used in aerodynamics. In conservation form, these equations may be written 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 simplifies 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, from elliptic to hyperbolic. 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 helps explain why 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. Closed-form analytical solutions generally require substantial simplifying assumptions, but the full equations can be discretized and solved numerically. For turbulent aircraft flows, however, the range of spatial and temporal scales is usually too large to resolve directly, so practical calculations require turbulence modeling or partial resolution of turbulence. A laminar-flow calculation retains molecular viscous stresses but excludes turbulent fluctuations and turbulent stresses. 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

(20)  

and in tensor notation

(21)  

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.,

(22)  

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 CFD calculations, the Navier-Stokes equations are usually written in conservation form because this form tracks the fluxes of mass, momentum, and energy through the boundaries of each computational cell. The compact conservation form was introduced earlier as Eq. 1 in the aerodynamic-equation hierarchy. In that notation, is the conserved-variable vector, while , , and are the inviscid flux vectors. The corresponding viscous fluxes are represented by , , and .

This form is especially useful for finite-volume methods because the same conservation principles can be applied directly to each cell in the computational grid. It also makes clear how simplified sets of equations are obtained. For example, neglecting the viscous fluxes yields the Euler equations, whereas retaining only the dominant wall-normal viscous derivatives yields the thin-layer Navier-Stokes approximation.

One useful reduction is the thin-layer Navier-Stokes equations, in which the viscous derivatives parallel to the surface are neglected while the dominant wall-normal viscous derivatives are retained. Symbolically, this approximation may be written as

(23)  

where denotes the local wall-normal direction and is the retained viscous flux in that direction.

The thin-layer Navier-Stokes equations are useful when viscous gradients normal to a surface dominate those in the streamwise and spanwise directions. This approximation reduces computational effort while retaining the main viscous effects needed for boundary-layer-type aerodynamic calculations.

Reynolds-Averaged Navier-Stokes (RANS) Form

The RANS equations are obtained by statistically averaging the Navier-Stokes equations. This averaging is necessary because real aircraft flows are typically turbulent, and resolving all instantaneous turbulent motion is generally too computationally expensive for engineering calculations. The key idea is to decompose each instantaneous flow quantity into a mean component and a fluctuating component.

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

(24)  

and for the static pressure

(25)  

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, the nonlinear convective term introduces additional terms involving products of the fluctuating velocity components. These terms are associated with the Reynolds stresses, which represent the apparent stresses created by turbulent momentum transport. The averaged form of the incompressible continuity and momentum equations then becomes

(26)  

and

(27)  

The quantity is a tensor formed from products of fluctuating velocity components, such as , , and . Its divergence represents the effect of turbulent momentum transport on the mean flow. These Reynolds-stress terms cannot be determined from the mean-flow variables alone, so they must be modeled by a turbulence closure model.

RANS methods are widely used in practical aerodynamic calculations because they include viscous effects and turbulence modeling at a computational cost that is usually acceptable for engineering design. They are not exact solutions of turbulent flow; rather, they solve for the mean flow field and model the effects of unresolved turbulent fluctuations. This compromise makes RANS the most common CFD approach for many aircraft and component-level aerodynamic analyses.

Gridding of the Computational Domain

Numerous algorithmic advances have improved numerical solutions of the Navier-Stokes equations. A fundamental issue in any flow problem involving a body is generating a grid that represents the body’s geometry to solve the governing equations. This means the computational domain should be large enough to prevent artificial boundary effects from contaminating the solution, but not so large that it unnecessarily increases memory requirements and computational cost.

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 and can provide a convenient structured topology for some geometries. As with O- and C-type grids, the grid points can be clustered near a surface to resolve the boundary layer. Nodes are the intersection points 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 type of grid 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

(28)  

where is the dynamic viscosity. This formulation implies that

(29)  

Here, denotes the eddy viscosity, an effective viscosity arising from turbulence. Therefore, turbulence models that aim to simulate variations in eddy viscosity are called eddy-viscosity models. In these models, the effect of turbulent momentum transport is represented by augmenting the molecular viscosity with a turbulent viscosity coefficient. Under the Boussinesq eddy-viscosity approximation, the Reynolds stresses for incompressible flow are modeled as

(30)  

where is the turbulent or eddy viscosity, is the turbulent kinetic energy, and is the Kronecker delta. In incompressible-flow formulations, the isotropic term is often absorbed into a modified mean pressure.

This modeled form assumes that the turbulent stress contribution is proportional to the mean strain-rate tensor, an approach known as the Boussinesq hypothesis. It assumes that turbulent momentum transport can be represented by an effective eddy viscosity, analogous in form to 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.

With the eddy-viscosity approximation, a steady incompressible mean momentum equation may be written as

(31)  

where may include the isotropic turbulent-stress contribution. This term reduces to a constant effective viscosity multiplying only under additional assumptions, including spatially uniform effective viscosity.
Eddy viscosity models belong to a class of so-called “turbulence models” or “turbulence closure models,” which provide a mathematical basis for solving for the mean flow field and turbulent 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.

Choosing a turbulence model that yields useful predictions depends on the flow’s length and time scales, the degree of separation, the level of streamline curvature, and the expected anisotropy of the turbulence. Turbulent flows contain a wide range of spatial scales, each contributing to overall flow behavior. These scales are often described as large, intermediate, and small eddies, as illustrated in Figure 8. Because these scales interact in complex ways, no single turbulence model works equally well for all aerodynamic flows. Modeling turbulent stresses remains a central challenge in aerodynamic prediction.

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 are smaller than the large-scale eddies but larger than the smallest scales. They transfer energy from large-scale eddies to the smallest scales through the energy cascade.
3. Small-scale eddies are the smallest turbulent structures, with length scales that can be much smaller than the flow’s characteristic length. They dissipate 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 provides a computationally efficient way to incorporate turbulence effects into the Navier-Stokes equations, enabling prediction of the mean flow and associated turbulence statistics. However, in more complex turbulent flows, the Boussinesq hypothesis may need modification or replacement with more advanced turbulence models because 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 solves a transport equation for a modified eddy viscosity. It is widely used for external aerodynamic flows because it is relatively economical and often gives useful predictions for attached boundary layers.
2. The k-epsilon, or –, model solves two transport equations, one for the turbulent kinetic energy and another for the turbulent dissipation rate . Variants include the RNG – and Realizable – models.
3. k-omega, or –, models also solve two transport equations but use the specific dissipation rate, , instead of . Variants include the SST – and BSL – models.
4. Reynolds Stress Models (RSMs) solve modeled transport equations for the Reynolds stresses rather than relating them directly to the mean strain rate through an eddy viscosity. They can better represent 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 to RANS models, two higher-fidelity turbulence-simulation approaches are worth mentioning: Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS).

1. Large Eddy Simulation (LES) resolves the larger, energy-containing turbulent eddies while modeling the effects of the smaller dissipative scales. LES can give more detailed unsteady-flow information than RANS, but it is much more computationally expensive, especially near solid walls at high Reynolds numbers.
2. Direct Numerical Simulation (DNS) resolves all dynamically important turbulent scales without using a turbulence closure model. DNS provides valuable research data, but its computational cost is so high that it is generally limited to simple geometries and relatively low Reynolds numbers.

The choice of turbulence model depends on the flow physics, the expected level of separation, available computational resources, and the required accuracy. No turbulence model is reliable for all flows, so its assumptions, calibration range, and empirical content must be understood before interpreting the results. In practical CFD work, engineers often compare predictions from multiple turbulence models and assess them against experimental data, benchmark calculations, or other trusted reference solutions.

Numerical Methods

Numerical techniques for solving the Navier-Stokes or RANS equations usually begin by discretizing the governing partial differential equations in space and, for unsteady problems, in time. Common approaches include finite-difference, finite-volume, finite-element, and related methods. Finite-difference methods directly approximate derivatives using Taylor series expansions at the nodes of a computational grid.

Finite-Difference Methods

Leonhard Euler made important early contributions to finite-difference methods in the 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

(32)  

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. 32, 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

(33)  

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

(34)  

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

(35)  

where is some constant. If = (temperature), then Eq. 35 is the heat equation, which models the transport 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

(36)  

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

(37)  

Combining the terms leads to

(38)  

and after rearrangement, then

(39)  

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 specifying initial and boundary conditions (if any). After that, the stencil can be applied to obtain values over the entire grid; Figure 10 illustrates the idea. Applying the moving stencil to the first row of initial values yields the solution at 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, expressed in differential conservation form in Eq. 1, and treat the computational domain as a collection of control volumes or cells rather than as discrete nodes alone. 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 fluxes are evaluated over the boundaries of each control volume. 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. Integrating over each cell converts the governing equations into algebraic equations, and the resulting system is solved to obtain the flow variables.

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

(40)  

Applying the divergence theorem to the spatial term gives

(41)  

which in one dimension simplifies to

(42)  

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

(43)  

Substituting the approximations into the integral form gives

(44)  

which can be simplified to obtain

(45)  

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

(46)  

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 useful predictive accuracy at a reasonable cost. Numerically solving the RANS equations involves discretizing the averaged governing equations to obtain a system of algebraic equations. For incompressible flow, the primary mean-flow equations are the continuity and momentum equations. For compressible aerodynamic flows, the energy equation, equation of state, and the selected turbulence model equations must also be included.

1. The first step is to divide the computational domain into a suitable computational grid or mesh. The RANS equations are then discretized over the computational domain using numerical methods. For unsteady calculations, the governing equations must also be discretized in time. Steady calculations may instead use iterative or pseudo-time methods to obtain convergence.
2. Apply the governing conservation equations to each control volume, including continuity, momentum, and, for compressible flow, energy. 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 are assembled into algebraic systems, as outlined previously for the FVM. For steady calculations, the solution is usually advanced through nonlinear or pseudo-time iterations until convergence is achieved. For unsteady calculations, the equations are advanced through physical time steps, often with inner iterations at each step. Implicit methods solve coupled algebraic systems, whereas explicit methods update the solution using known values from the current time level.
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 until a steady state or a time-accurate solution is reached. For pressure-based formulations, 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 in Figure 11 illustrate RANS predictions of 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 is achieved by discretizing the partial differential equations. Consider, for example, a two-dimensional, steady, incompressible mean flow with an eddy-viscosity closure. In simplified form, the modeled momentum equation may be written as

(47)  

In a finite-volume formulation, the convective contribution to the -momentum equation is written as the net flux of -momentum through the cell faces. For a rectangular two-dimensional control volume,

(48)  

Face values may be evaluated using upwind, central, or higher-order interpolation. The mass flux and transported momentum must be defined explicitly so that velocity is not counted twice. The basic concept is shown in Figure 12.

Using central differencing for the diffusive term gives

(49)  

Using central differencing for the pressure gradient term gives

(50)  

Combining the convective, pressure-gradient, and diffusive terms gives a discrete algebraic momentum equation for the cell value . Its exact coefficients depend on the definitions of the face mass fluxes, the interpolation scheme used for face velocities, and the treatment of the effective viscosity. A common assembled form is

(51)  

where the coefficients contain the discretized convection and diffusion contributions and contains the pressure-gradient and source terms.
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.

Commercial CFD software packages typically implement this type of numerical process. 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 flow complexity, 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 that they incorporate 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 turbulence effects; they do not resolve the instantaneous turbulent eddies themselves. Nevertheless, researchers can now compute remarkable details of flows around complete aircraft. 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 moving grids (for the rotors) and stationary grids (for the airframe) simultaneously. 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 a relatively low pitch, which convect beneath the rotor and quickly diffuse and break up into a turbulent jet. In hover, the rotors require relatively high induced power because they must generate thrust by accelerating air from essentially zero freestream velocity. 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, where 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. This mitigates the need to grid the computational domain at extremely fine resolution, reducing computational costs by orders of magnitude. If the flow is assumed inviscid, the Navier-Stokes equations reduce to 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

(52)  

(53)  

(54)  

where is the total energy per unit mass. Alternatively, by explicitly dropping the viscous fluxes from Eq. 1, the Euler equations can be written in conservation form as

(55)  

In practice, numerical solutions to the Euler equations provide excellent approximations of 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 the numerical formulation selects the appropriate circulation, typically through the boundary conditions and an effective Kutta condition at the trailing edge. 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 boundary-layer effects 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 normal to the surface than 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 account for the boundary layer’s effect 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 is no longer valid 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

(56)  

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

(57)  

Because the pressure across the boundary layer (i.e., normal to the surface) is approximately constant to boundary-layer order (Eq. 57), 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

(58)  

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 pressure gradient can now be eliminated using the outer-flow momentum equation, giving

(59)  

where . For steady flow, the term vanishes.

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

(60)  

This quantity allows the surface to be displaced from its original position by , yielding the same flow rate as in 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

(61)  

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

(62)  

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 useful approximation for external flows outside boundary layers, wakes, and other regions with significant vorticity. In this case, the velocity field can be expressed as a gradient of a potential function, , such that

(63)  

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

(64)  

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 to properly represent the characteristic variation of flow properties with Mach number. However, other simplifications are possible.

Under the small-disturbance assumptions, the full-potential equation can be reduced while retaining the leading nonlinear term needed to represent transonic mixed-type behavior. A representative two-dimensional steady transonic small-disturbance equation is

(65)  

The coefficient of depends on the perturbation velocity and can change sign locally, allowing the equation to represent adjoining subsonic and supersonic regions. Different nondimensionalizations produce equivalent forms with different numerical coefficients, so the scaling used for must always be stated.

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

(66)  

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

(67)  

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). An important breakthrough came in 1971, when Earll Murman and Julian Cole developed a finite-difference relaxation method for solving the mixed subsonic-supersonic TSD equation.^[9] Their method became an important foundation for subsequent numerical methods for transonic potential flow and helped establish CFD as a practical tool for aerodynamic analysis.

Naturally, the TSD equation is fundamental to the study of transonic flows, where speeds are near the speed of sound, and the coexistence of subsonic and supersonic regions in the flow field poses unique challenges. 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

(68)  

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.,

(69)  

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. Many classical panel methods are based on incompressible irrotational flow, for which the governing equation is Laplace’s equation (Eq. 69) for the velocity potential. The idea is to distribute surface singularities of unknown strength over the body surface. After determining the singularity strengths by imposing a flow-tangency condition on the surface, the surface velocities and, consequently, the inviscid aerodynamic loading on the body can be computed.

Several well-known panel-method computer codes have been developed and widely used in aerodynamic analysis. Examples include XFOIL, which combines a two-dimensional potential-flow panel method with integral boundary-layer and transition models for airfoil analysis. VSAERO is a three-dimensional subsonic potential-flow panel code that can be coupled with approximate viscous corrections. PANAIR is a higher-order panel code based on linearized potential-flow theory and can analyze complex configurations in subsonic and supersonic flow within the limits of that theory. 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 described by a unit-strength velocity potential . This means that the total potential induced by all the singularities in the flow is

(70)  

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.,

(71)  

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

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 only source singularities 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

(72)  

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

(73)  

One such equation applies to each of the panels.

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

(74)  

which may be written formally as

(75)  

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 tangential velocity components on 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 enough 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 as , so its use is by no means inexpensive, even on a modern computer. However, panel methods remain considerably less computationally intensive than most CFD solvers 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, solving the momentum integral equation from boundary-layer theory enables displacement-thickness corrections to panel-method solutions, thereby modeling attached viscous boundary-layer effects, as 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 airfoil lift predictions 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 lift predictions 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 densities, 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.,

(76)  

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, molecular collisions tend to drive the gas’s velocity distribution toward thermodynamic 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

(77)  

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, on average, more molecules move faster. The average speed of all molecules is

(78)  

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

(79)  

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 flows or is 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

(80)  

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 aerodynamic heating of spacecraft during atmospheric entry, which can differ from those from 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 portions of the spacecraft’s high-altitude atmospheric entry into the Martian atmosphere. 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, engineers use elements such as beams, trusses, triangles, or hexahedra. Like CFD, FEM uses meshes, but it usually classifies elements by topology, such as beam, shell, triangular, quadrilateral, tetrahedral, or hexahedral. 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 helps improve the representation of stress concentrations and 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. Engineers then solve these equations numerically by assembling a global system from individual element contributions, 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 design.

History

The origins of FEM trace back to structural analysis work in the 1940s. Richard Courant’s 1943 work is often cited as an important precursor to the Finite Element Method (FEM). 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 Edward Wilson’s leadership, is a notable early FEM program.

FEM was extended beyond structural mechanics into other fields such as heat transfer, fluid mechanics, and electromagnetics. This era also saw the development of general-purpose FEM software, such as NASTRAN (NASA Structural Analysis), which 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, as illustrated in Figure 24. If an element has nodes, the displacement within that element is approximated using shape functions and nodal displacements , i.e.,

(81)  

For each element, the stiffness matrix is given by

(82)  

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. 82 is implied.

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

(83)  

The global force-displacement relationship is given by

(84)  

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 linearly superimposing the element stiffness matrices. Solving the global system of equations efficiently, especially for large-scale problems, involves using various numerical methods.

2-D Triangular Element

A simple triangular element provides an opportunity to deepen understanding of the FEM. The method’s systematic, discretized nature makes it well-suited to computer implementation, making the steps readily programmable. 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,

(85)  

(86)  

(87)  

where the signed area is

(88)  

and the physical element area is . The shape-function expressions use in their denominators unless the nodes are explicitly required to be numbered counterclockwise, in which case .
The constants , , and are determined by the coordinates of the nodes, i.e.,

(89)  

(90)  

(91)  

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

(92)  

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

(93)  

Therefore, the matrix can be written as

(94)  

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

(95)  

where is the area of the triangular element. If unit thickness is assumed, then . This simplification is one reason constant-strain triangular elements are pedagogically useful, 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 according to the node connectivity. 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. Next, the appropriate finite element type (e.g., rectangular, triangular, tetrahedral, hexahedral) is selected, and the geometry is discretized into these 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

FEM has become a standard tool for predicting structural response to applied loads in engineering; Figure 27 shows an FEM solution for a wing and the spar connected to the fuselage. NASTRAN is a widely used FEM solver for structural analysis. Companion software typically handles preprocessing tasks such as meshing and applying boundary conditions, as well as postprocessing and visualization. Advances in computer technology have significantly expanded its capabilities, enabling 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 presented as separate computational methods, but many important aerospace problems require their integration. 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 distributions over the aircraft’s surface. These aerodynamic loads are then transferred to the structural finite element model. The FEM solver computes the resulting structural deflections, strains, and stresses. 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.

Figure 28 shows an example of a simulation-based aeroelastic assessment for a high-aspect-ratio transport-aircraft wing. 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 to perform load 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 structural and aeroelastic analyses must be integrated carefully with aerodynamic design; 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 to better understand 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 and loads, while material properties are specified separately to represent the structural behavior of the components.