Introduction to Advanced Computational Methods

J. Gordon Leishman

54 Introduction to Advanced Computational Methods

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, usually called “CFD,” and the finite element method, or “FEM,” have been revolutionary computational tools for predicting aerodynamics and structural loads, respectively. The intention of this chapter is not to review the details of such techniques and their associated numerical algorithms per se, but to give an introduction and general overview of the foundational principles and how they are used, along with a summary of their capabilities and some of their limitations. Nor is it meant to be an all-encompassing treatment.

Traditionally, the aircraft industry has relied on quicker and more economical computational methods for designing aircraft and predicting their flight performance, and successfully, too. One does not need CFD or the FEM to design and make a flight vehicle work in the way it was intended; history speaks for itself. However, both tools can help significantly refine the aircraft design to extract its maximum performance at the lowest cost, with the byproduct of improving safety and reliability. Knowing why something works (or doesn’t) is vital to commercial success and continued technological advancement. Wind tunnel testing has always been a critical part of the flight vehicle design process and will continue to be essential for measurements that can be used for CFD validation. Advances in computer technology, such as code execution speed, memory and data storage access speed, and the availability of inexpensive mass data storage systems, continue to encourage the development of more ambitious numerical methods that can be used during design.

Much new fundamental research work in modeling and algorithmic development has occurred 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 computer clusters, allowing for research into new computational methods for use in a broad range of technical fields. As shown in the image below, the ability to predict flows about complete aircraft can open up opportunities in new design directions, such as the reconsideration of supersonic transports (SSTs). It has also inspired broader innovations in flight vehicle designs, allowing concepts outside the “domain of conventional wisdom” to be investigated confidently. Thinking “out of the box” can now be quantified for the right reasons!

An advanced CFD solution around a proposed SST aircraft can be used to explore configurational changes at the design stage. (NASA public domain image.)

Continued advances in computer capabilities, which continue to accelerate,^[1] has also played a large part in motivating and accelerating the development of more integrated computational techniques, encompassing 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, so some form of coupling of methods is required. Increasingly tighter mathematical, algorithmic, and numerical coupling have allowed, for example, between aerodynamics and structural dynamics to be solved simultaneously. If fully realized, this integration strategy will eventually predict the behavior of an entire flight vehicle in an arbitrary flight condition and so allow more capable and efficient vehicles to be designed. However, significant integration problems and more ambitious numerical techniques are still needed to push forward this technology before they can fully predict the vehicle’s characteristics before their first flights.

Learning Objectives

Understand the foundational elements of modern computational aerodynamic methods starting from the governing Navier-Stokes equations.
Appreciate the hierarchy of such aerodynamic methods in terms of their relative predictive advantages and disadvantages.
Know about surface singularity methods and their capabilities and limitations for aerodynamic prediction.
Understand why boundary layers significantly affect aerodynamic behavior and why their modeling is essential.
Appreciate the fundamentals of the structural finite element method and its uses.
Learn about the field of aeroelasticity and why it is essential in the design of flight vehicles.

Aerodynamics

Over the past three decades, advances have been made toward understanding challenging problems in aerodynamics by using CFD. Numerical methods, such as finite-difference methods (FDM) or finite-volume methods (FVM), applied to the governing flow equations (e.g., Navier-Stokes), are used to model the flow field about an aircraft. The choice of which governing equations to use affects the level of physics that can be approximated by the CFD scheme and the computational effort and time taken to solve it. Not all numerical methods are created equally, and the choice of method affects the solution’s accuracy, stability, and cost. The definition of “cost” comes down to how much time^[2] and effort goes into obtaining a solution, and most CFD methods are considered costly by any design standard.

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

A key issue in using these numerical methods is the definition and creation of appropriate computational grids over which to solve the governing equations, requiring considerable skill and effort. Other problems include the difficulties in modeling flow turbulence, for which no single method exists. Furthermore, CFD methods must be validated against measurements to realize the high predictive confidence levels needed for aircraft design work. Commercially available tools such as ANSYS Meshing and Pointwise can aid in grid generation, while CFD solvers like ANSYS Fluent and OpenFOAM can run 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 the calculations were performed “by hand,” i.e., using pencil, paper, 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 called 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 of the 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 important for solving aerospace problems.

The publication of many technical papers and essential textbooks during the 1960s helped to disseminate these numerical methods and apply them to fluid flow problems. Textbooks by Brian Spalding ^[5] and Suhas V. Patankar ^[6] were used to teach CFD to students as part of courses in fluids. NASA and other aerospace laboratories worldwide played a crucial role in advancing CFD for practical applications in aerospace engineering. In 1987, NASA founded the Advanced Supercomputing (NAS) Division at NASA Ames. In addition, the U.S. Army Aeronautical Research Laboratory played an essential role in advancing the CFD field. The first commercial CFD software packages also emerged during this period, making the technology accessible to more 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, allowing for more complex CFD simulations at a local level. 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 from the early 2000s until today, with the development of more sophisticated algorithms, increased computational power, and the advent of parallel computing. High-performance computing (HPC) and graphics processing units (GPUs) have become commonplace, and have enabled the simulation of more complex and larger-scale aerospace problems. Additionally, CFD has benefited from advancements in data visualization techniques, making it easier for engineers to interpret and analyze large quantities of data and understand complex flows. Today, CFD is an essential tool in engineering analysis, with applications ranging from designing 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.

Navier-Stokes Equations

Most aerodynamic problems with flight vehicles involve viscous effects and turbulence. These effects are significant in the turbulent boundary layers found on the wings and airframe, as well as with a host of component interaction phenomena, such as wing/fuselage interference, powerplant/wing interference, etc. Modeling flow separation is also essential, as the onset can set boundaries for an aircraft’s operational flight envelope. The Navier-Stokes equations are the fundamental equations that govern fluid dynamic behavior. They encompass the principles of conserving mass, momentum, and energy and their interchange 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. The most common practical simplification is to avoid treating turbulence and the associated viscous shearing effects by assuming a laminar flow. While this approach simplifies the Navier-Stokes equations, it should be noted that the flows will be turbulent even at the lowest Reynolds numbers found on aircraft.

Tensor Form

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

(1) $\begin{equation*} \frac{\partial \varrho}{\partial t} + \nabla \bigcdot (\varrho \, \mathbf{u}) = 0 \end{equation*}$

and in tensor form by

(2) $\begin{equation*} \frac{\partial \varrho}{\partial t} + \frac{\partial (\varrho \, u_i)}{\partial x_i} = 0 \end{equation*}$

where $\varrho$ is the flow density and $\mathbf{u} = (u_i)$ , $i = 1, 2, 3$ , is the velocity vector. For an incompressible flow, where $\varrho$ = 0, then

(3) $\begin{equation*} \nabla \bigcdot \mathbf{u} = 0 \quad \mbox{and} \quad \frac{\partial u_i}{\partial x_i} = 0 \end{equation*}$

Conservation of momentum (momentum equation) gives

(4) $\begin{equation*} \varrho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \bigcdot \nabla) \mathbf{u} \right) = -\nabla p + \nabla \bigcdot \boldsymbol{\tau} + \varrho \, \mathbf{f} \end{equation*}$

and in tensor form by

(5) $\begin{equation*} \varrho \left( \frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} \right) = - \frac{\partial p}{\partial x_i} + \frac{\partial \tau_{ij}}{\partial x_j} + \varrho \, f_i \end{equation*}$

where $p$ is the pressure, $\tau_{ij}$ is the stress tensor, and ${f_i}$ is a body force per unit mass, e.g., acceleration under gravity.

Modeling Viscous Effects

The Navier-Stokes equations are general and apply to all flows, including flows with turbulence. The flow about an aircraft, at least at lower Mach numbers, can be readily approximated as an isotropic Newtonian fluid, for which the relationship between stress and strain in the Navier-Stokes equations is given by

(6) $\begin{equation*} \tau_{ij} = \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) + \lambda (\nabla \bigcdot \mathbf{u}) \delta_{ij} \end{equation*}$

where $\delta_{ij}$ is the Kronecker delta function, i.e., $\delta_{ij} = 1$ for $i = j$ , and $\delta_{ij} = 0$ otherwise, and $\mu$ is the coefficient of dynamic viscosity. The parameter $\lambda$ is called the second coefficient of viscosity or bulk viscosity, usually given as $\lambda = -(2/3) \mu$ , which is referred to as Stokes’s hypothesis.^[7] The bulk viscosity, which can be measured, as can $\mu$ , is essential in describing the compression and energy dissipation in fluids, both gasses and liquids, such as in sound wave propagation, shock waves, and other phenomena.

Notice from either Eq. 1 2, for an incompressible flow where $\nabla \bigcdot \mathbf{v} = 0$ , then the term involving $\lambda$ vanishes, so the stress tensor depends only on the dynamic viscosity, $\mu$ . For a compressible flow, however, the term $\nabla \bigcdot \mathbf{v}$ represents the rate at which the fluid’s density changes under compression or expansion. Therefore, the term $\lambda ( \nabla \bigcdot \mathbf{v})$ in the stress tensor accounts for the viscous stresses arising from these volume changes.

Rounding out the conservation laws is the conservation of energy (energy equation), which can be expressed as

(7) $\begin{equation*} \varrho \left( \frac{\partial e}{\partial t} + \mathbf{u} \bigcdot \nabla e \right) = \nabla \bigcdot (k \, \nabla T) + \tau_{ij} \, \frac{\partial u_i}{\partial x_j} + \varrho \, q \end{equation*}$

and in tensor notation

(8) $\begin{equation*} \varrho \left( \frac{\partial e}{\partial t} + u_i \frac{\partial e}{\partial x_i} \right) = \frac{\partial}{\partial x_i} \left( k \, \frac{\partial T}{\partial x_i} \right) + \tau_{ij} \, \frac{\partial u_i}{\partial x_j} + \varrho \, q \end{equation*}$

where $e$ is the internal energy per unit mass, $T$ is temperature, $k$ is the bulk thermal conductivity, and $q$ is the external heat per unit mass. Finally, the ideal gas equation (equation of state) can be used to relate the thermodynamic quantities, i.e.,

(9) $\begin{equation*} p = \varrho \, R \, T \end{equation*}$

where $R$ 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.,

(10) $\begin{equation*} \frac{\partial {\bf Q}}{\partial t} + \frac{\partial {\bf E}}{\partial x} + \frac{\partial {\bf F}}{\partial y} + \frac{\partial {\bf G}}{\partial z} = 0 \end{equation*}$

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

(11) $\begin{equation*} \frac{\partial {\bf Q}}{\partial t} + \frac{\partial {\bf E}}{\partial x} + \frac{\partial {\bf F}}{\partial y} + \frac{\partial {\bf G}}{\partial z} = \frac{\partial {\bf E_v}}{\partial x} + \frac{\partial {\bf F_v}}{\partial y} + \frac{\partial {\bf G_v}}{\partial z} \end{equation*}$

where ${\bf E_v}$ , ${\bf F_v}$ , and ${\bf G_v}$ 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 $x$ and $y$ directions are neglected to give

(12) $\begin{equation*} \frac{\partial {\bf Q}}{\partial t} + \frac{\partial {\bf E}}{\partial x} + \frac{\partial {\bf F}}{\partial y} + \frac{\partial {\bf G}}{\partial z} = \frac{\partial {\bf S}}{\partial z} \end{equation*}$

for a suitably defined viscous flux, ${\bf S}$ , that accounts for the neglected $z$ derivatives.

The thin-layer Navier-Stokes equations are particularly useful in scenarios where the flow can be approximated as two-dimensional or quasi-two-dimensional, simplifying computational efforts and modeling complexities. Engineers often apply these simplified equations in boundary layer analyses, fluid film lubrication, and other situations where the flow dynamics are primarily in the $x$ and $y$ directions.

Reynolds-Averaged Navier-Stokes (RANS) Form

The Reynolds-Averaged Navier-Stokes (RANS) equations are a form of the Navier-Stokes equations to describe the flow of incompressible fluids. These equations are particularly useful in the study of turbulent flows. 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

(13) $\begin{equation*} \nabla \bigcdot \vec{V} = 0 \quad \hbox{(Mass conservation)} \end{equation*}$

and

(14) $\begin{equation*} \frac{\partial \vec{V}}{\partial t} + \vec{V} \bigcdot \nabla \vec{V} = -\frac{1}{\varrho} \nabla p + \nu \, \nabla^2 \vec{V} + \varrho \vec{f}_{b} \quad \hbox{(Momentum conservation)} \end{equation*}$

where $f_b$ represents the external body forces per unit mass.

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

(15) $\begin{equation*} u(t) = \overline{u} + u' \end{equation*}$

and for the static pressure

(16) $\begin{equation*} p(t) = \overline{p} + p' \end{equation*}$

Here, $\overline{u}$ and $\overline{p}$ are the mean (time-averaged) components, while $u'$ and $p'$ are the fluctuating components, as illustrated in the figure below.

The classic Reynolds decomposition splits the turbulent flow properties into a statistical mean (average) and a statistically fluctuating part.

When the Navier-Stokes equations are averaged over time, the nonlinear term $\vec{V} \bigcdot \nabla \vec{V}$ 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 become

(17) $\begin{equation*} \nabla \bigcdot \vec{V} = 0 \quad \hbox{(Mass conservation)} \end{equation*}$

and

(18) $\begin{equation*} \frac{\partial \overline{\vec{V}}}{\partial t} + \overline{\vec{V}} \bigcdot \nabla \overline{\vec{V}} = -\frac{1}{\varrho} \nabla \overline{p} + \nu \, \nabla^2 \overline{\vec{V}} - \nabla \bigcdot \overline{\vec{V}' \vec{V}'} + \varrho \overline{\vec{f}_{b}} \quad \hbox{(Momentum conservation)} \end{equation*}$

The term $-\nabla \bigcdot \overline{\vec{V}' \vec{V}'}$ 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 directly calculated from the mean flow variables.

The RANS equations provide a CFD method that balances computational cost and accuracy well, 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 in aircraft and automobile design, ship hull design, water flow in pipes, weather forecasting, and environmental engineering, including pollutant dispersion. Several CFD software packages are available to solve the RANS equations, including ANSYS Fluent, OpenFOAM, COMSOL Multiphysics, and STAR-CCM+.

Griding of the Computational Domain

There have been many algorithmic developments in the numerical solution of the Navier-Stokes equations. A fundamental issue in any flow problem involving a body is generating the grid structure representing the body’s geometry on which to solve the governing equations. This means that the size of the computational domain should be as small as possible for economy, yet large enough to ensure that any nonphysical issues within the computational domain are avoided, making computer memory requirements very large.

Several types of grid systems are appropriate for analyzing flow problems; a summary is shown in the figure below. Grids are usually categorized according to their geometric shape, e.g., O-grids, C-grids, and H-type grids. The O-grids and C-grids are often used to analyze the flows about wings as they can conform to the curved airfoil shape of the surface, which are called body-fitted grids, as shown in the figure below. 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 intersection points on the grid, and the nodes make up a cell. The numerical solution to the flow properties is then based on finite-difference or finite-volume solutions applied over the nodes and cells.

Types of elementary grids used as a basis for solving problems in computational fluid dynamics.

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

Creating CFD grids for a wing and wing tip requires considerable care to compute the flow correctly.

Overset or chimera grid systems are locally refined grids superimposed upon a primary or background grid, as shown in the figure below, which is like adding a grid patch where and when needed. These overlapping grids help resolve local flow phenomena without griding the entire flow domain at high resolution, thereby saving computational effort. They can also be used to overcome specific difficulties in dealing with a domain that is most easily meshed using separate grids that move relative to each other, such as propellers and rotors relative to an airframe.

Chimera grids, which are overlapping “patch” grids, allow for local grids without requiring grid refinement of the entire CFD domain. (NASA image.)

With chimera grids, interpolation methods are required to transfer flow property information between the grids. The accuracy of these procedures needs to be carefully examined to ensure the overall accuracy of the numerical solution is maintained. Such grids can also be very labor-intensive to create. Unstructured grids, composed of irregularly distributed tetrahedral elements, as shown in the figure below, are increasingly being used. 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 reaching a more optimal number of grid points, thereby speeding up the numerical solution and reducing cost.

An unstructured “irregular” grid is composed of random tetrahedral elements.

Turbulence Modeling

Modeling turbulence in most flow problems requires introducing empirical “closure” models, of which many types 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 complete aircraft.

The viscous stresses can be expressed in the form

(19) $\begin{equation*} \nu \frac{\partial^2 \overline{u}}{\partial y^2} - \frac{\partial}{\partial y}(\overline{u'v'}) = \frac{\partial}{\partial y} \left( \nu \frac{\partial \overline{u}}{\partial y} - (\overline{u'v'}) \right) = \frac{\partial \left( \tau_l + \tau_t \right)}{\partial y} \end{equation*}$

where $\nu = \dfrac{\mu}{\varrho}$ . This latter formulation implies that

(20) $\begin{equation*} \tau_t = -\overline{u'v'} = \mu_t \left( \frac{\partial \overline{u}}{\partial y} \right) \end{equation*}$

Here, $\mu_t$ 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 in the flow are called eddy viscosity models. In these models, higher levels of turbulent stresses are simulated by augmenting the effects of molecular viscosity alone with a turbulent viscosity coefficient, i.e.,

(21) $\begin{equation*} \tau_e = \left( \mu + \mu_t \right) \nabla \overline{\vec{V}} \end{equation*}$

This preceding result assumes that the turbulent stress tensor is proportional to the mean strain rate tensor multiplied by a constant, which is called the Boussinesq hypothesis. It also assumes that the exchange of turbulent energy in the cascading process of eddies is analogous to the effects of molecular viscosity, as illustrated in the figure below. Many experiments andons of turbulent flows support this other observati behavior. Therefore, the Boussinesq hypothesis means that the turbulent flow can be modeled as an effective total viscosity that varies with the local flow conditions.

The concept of an “eddy” viscosity is used to augment the laminar molecular viscosity from the effects of mixing caused by turbulence.

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

(22) $\begin{equation*} \varrho \left( \overline{\vec{V}} \bigcdot \nabla \, \overline{\vec{V}}\right) = -\nabla \overline{p} + \left( \mu + \mu_t \right) \nabla^2 \,\overline{\vec{V}} \end{equation*}$

Eddy viscosity models are in a class of so-called “turbulence models” or “turbulence closure models,” which provide a mathematical basis for solving the mean flow field, as well as the turbulent flow quantities. Many different turbulent flow models are available, especially for use in CFD, and choosing the most suitable turbulence model requires careful consideration regarding the problem. The best way to represent the turbulent viscosity and make the decision on the most suitable turbulence model depends on several factors, including the specific flow application, the desired level of numerical accuracy, and the available computational resources.

Deciding on the types of turbulence models that can be expected to give good predictive capability depends on understanding the length and time scales of the specific flow problem at hand. It should be appreciated that turbulent flows exhibit a wide range of spatial scales, each of which plays a vital role in the overall behavior of the flow. The spatial scales of turbulence can be broadly classified into three categories: large, intermediate, and small eddies, as illustrated in the figure below. Therefore, there is no “one size fits all” turbulence model. The different scales of turbulence interact with each other in a complex manner, and a wide range of phenomena, such as vorticity, turbulence intensity, and Reynolds stresses, all characterize their behavior. Modeling these scales is crucial for predicting the overall behavior of turbulent flows, which remains a critical challenge in aerodynamic modeling.

Eddies of different scales and velocities will affect the local flow properties and cause positive and negative fluctuations.

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.
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.
Small-scale eddies are the smallest structures, with linear scales typically millimeters or less. They are responsible for energy dissipation and convert the kinetic energy in the flow into thermal energy through viscous shear.

The Boussinesq hypothesis is commonly used to develop turbulence models for solving the RANS equations and simulating turbulent flows in practical engineering applications. This hypothesis provides a computationally efficient way to incorporate the effects of turbulence into the Navier-Stokes equations, thereby allowing for the prediction of the mean flow and its turbulence statistics. However, in more complex turbulent flows, the Boussinesq hypothesis may need to be modified or replaced by more advanced turbulence models, especially considering that turbulence is often anisotropic.

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

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.
k-epsilon or $k$ – $\epsilon$ models solve for two transport equations, one for the turbulent kinetic energy $k$ and another for the turbulent dissipation rate $\epsilon$ . Variations of this model include the RNG $k$ – $\epsilon$ and Realizable $k$ – $\epsilon$ .
k-omega or $k$ – $\omega$ models also solve for two equations but use different variables, turbulent kinetic energy, $k$ , and specific dissipation rate, $\omega$ . Variations include the SST $k$ – $\omega$ and BSL $k$ – $\omega$ models.
Reynolds Stress Models (RSMs) resolve the Reynolds stresses directly, providing more detailed information about the turbulence field. They are designed for use in flows with anisotropic turbulence, flow separation, and recirculation. However, RSMs are computationally more expensive than other models.

In addition, two other methods of solving the Navier-Stokes equations are worth mentioning, namely LES and DNS.

Large Eddy Simulation (LES) is not a RANS model but is worth mentioning. LES focuses on capturing the larger, energy-containing eddies, while modeling the effects of the smaller, dissipative eddies. It is a computationally expensive method, but provides more accurate results for complex flows.
DNS stands for Direct Numerical Simulation, a computational technique used in fluid dynamics to simulate turbulent flows. In DNS, all the scales of turbulence are resolved numerically without the specification of a turbulence closure model. This means that even the smallest eddies and fluctuations in the flow 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 with such models to all types of flows limits the ability to predict turbulent flow properties, especially the development of turbulent boundary layers. Developing better turbulence models is still a goal in modern fluid dynamics, although these models must always be considered postdictive (i.e., after the fact) rather than predictive. It is common for engineers and researchers to perform a sensitivity analysis using different turbulence models, which can be used to assess the impact on the CFD results and compare them to benchmark solutions if available. These benchmarks may include both experiments and other theories.

Numerical Methods

The numerical techniques used to solve the Navier-Stokes equations may be based on either a finite-difference or a finite-volume approach, and numerous algorithms have been developed, each with relative advantages. Finite-difference methods are the most common. They approximate the time and space derivatives in the equations using Taylor series approximations at each node on the computational grid; a node can be considered the intersection between grid lines.

Finite-Difference Methods

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

(23) $\begin{equation*} \frac{\partial \phi}{\partial x} \approx \frac{\phi_{i+1}^{j} - \phi_{i}^{j}}{\Delta x} \end{equation*}$

which is called a forward difference method. The value of $\Delta x$ is the spatial discretization. The best way to understand this form of discretization is to use a stencil, as shown in the figure below. The nodes shown represent the four values in Eq. 23, i.e., $\phi_{i}^{j}$ , $\phi_{i-1}^{j}$ , $\phi_{i+1}^{j}$ and $\phi_{i}^{j+1}$ .

Example of a computational stencil, which can be applied over a grid to determine the scalar field values.

A backward difference approximation to the derivative is given by

(24) $\begin{equation*} \frac{\partial \phi}{\partial x} \approx \frac{\phi_{i }^{j} - \phi_{i-1}^{j}}{\Delta x} \end{equation*}$

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

(25) $\begin{equation*} \frac{\partial \phi}{\partial x} \approx \frac{\phi_{i+1}^{j} - \phi_{i-1}^{j}}{2 \, \Delta x} \end{equation*}$

which is second-order accurate, i.e., the error is of order $(\Delta x)^2$ . Higher-order differencing schemes can improve accuracy at the expense of more arithmetic operations and higher computer costs.

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

(26) $\begin{equation*} \frac{\partial \phi}{\partial t} = \alpha \left( \frac{\partial^2 \phi}{\partial t^2} \right) \end{equation*}$

where $\alpha$ is some constant. If $\phi$ = $T$ (temperature), then Eq. 26 is called the heat equation, and models the movement of heat energy through space and time. It forms a good test case for numerical methods because an exact solution can also be obtained, although it is still expressed as a series approximation.

Using a forward difference scheme for the time derivative gives

(27) $\begin{equation*} \frac{\partial \phi}{\partial t} \approx \frac{\phi_{i}^{j+1} - \phi_{i}^{j}}{\Delta t} \end{equation*}$

where $\Delta t$ is the temporal discretization. Using a central difference scheme for the spatial derivative gives

(28) $\begin{equation*} \frac{\partial^2 \phi}{\partial x^2} \approx \frac{\phi_{i+1}^{j} - 2\phi_i^j + \phi_{i-1}^j}{2 \Delta x } \end{equation*}$

Combining the terms leads to

(29) $\begin{equation*} \frac{\phi_i^{j+1} - \phi_i^j}{\Delta t} = \alpha \left( \frac{\phi_{i+1}^j - 2\phi_i^j + \phi_{i-1}^j}{2 \Delta x } \right) \end{equation*}$

and after rearrangement, then

(30) $\begin{equation*} \phi_i^{j+1} = \phi_i^j + \frac{\alpha \, \Delta t}{(\Delta x)^2} \bigg(\phi_{i+1}^j - 2\phi_i^j + \phi_{i-1}^j \bigg) \end{equation*}$

Therefore, the new value of $\phi$ at the next time step $j+1$ is determined from the values at the previous time step $j$ .

The numerical process also requires specifying initial and boundary conditions (if any). After that, the stencil can be applied to obtain the values of $\phi$ over the entire grid. The idea is shown in the figure below. 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. Moving the stencil to the left side of the grid, the values along the boundary condition, in this case $x$ = $a$ , can be populated. Applying the stencil further up the grid in time and across the physical space allows the values at step 4 to be obtained. The process can then be continued to cover the entire grid.

While a relatively simple example shows the basic ideas of how CFD methods work by applying a computational stencil and marching forward in time.

Finite Volume Methods

Finite volume methods operate on an integral form of the Navier-Stokes equations, derived from Eq. 10, and treat the computational domain as cells rather than nodes. The integral of the mass, momentum, and energy within each cell is divided by the cell volume to get an average conserved variable. The average fluxes are derived similarly by integrating over the boundary of the cell. 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, then solving the resulting system of algebraic equations to obtain the values of the flow variables.

Consider again the heat equation given in Eq. 26. Divide the domain into $N$ control volumes, each with width $\Delta x$ . Integrating over a control volume gives

(31) $\begin{equation*} \int_{x_{i - 1/2}}^{x_{i + 1/2}} \frac{\partial \phi}{\partial t} \, dx = \alpha \int_{x_{i - 1/2}}^{x_{i + 1/2}} \frac{\partial^2 \phi}{\partial x^2} \, dx \end{equation*}$

Applying the divergence theorem to the spatial term gives

(32) $\begin{equation*} \frac{\partial}{\partial t} \int_{V_i} \phi \, dV = \alpha \int_{\partial V_i} \nabla T \bigcdot \mathbf{n} \ dA \end{equation*}$

which in one dimension simplifies to

(33) $\begin{equation*} \frac{\partial}{\partial t} (\phi_i \Delta x) = \alpha \left( \frac{\partial \phi}{\partial x} \bigg|_{x_{i + 1/2}} - \frac{\partial \phi}{\partial x} \bigg|_{x_{i - 1/2}} \right) \end{equation*}$

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

(34) $\begin{equation*} \frac{\partial \phi}{\partial x} \bigg|_{x_{i + 1/2}} \approx \frac{\phi_{i+1} - \phi_i}{\Delta x} \quad \mbox{and} \quad \frac{\partial \phi}{\partial x} \bigg|_{x_{i - 1/2}} \approx \frac{\phi_i - \phi_{i-1}}{\Delta x} \end{equation*}$

Substituting the approximations into the integral form gives

(35) $\begin{equation*} \frac{\partial \phi_i}{\partial t} \Delta x = \alpha \left( \frac{\phi_{i+1} - \phi_i}{\Delta x} - \frac{\phi_i - \phi_{i-1}}{\Delta x} \right) \end{equation*}$

which can be simplified to obtain

(36) $\begin{equation*} \frac{\partial \phi_i}{\partial t} = \frac{\alpha}{ (\Delta x)^2} \bigg( \phi_{i+1} - 2\phi_i + \phi_{i-1} \bigg) \end{equation*}$

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

(37) $\begin{equation*} \phi_i^{n+1} = \phi_i^n + \Delta t \frac{\alpha}{( \Delta x)^2} \bigg( \phi_{i+1}^n - 2\phi_i^n + \phi_{i-1}^n \bigg) \end{equation*}$

There are many good computational reasons to adopt a finite-volume rather than a finite-difference approach. One reason is that finite-difference methods can only operate on structured grids, whereas finite-volume methods can be used on both structured and unstructured grids. However, both approaches are ultimately equivalent in the mathematical sense.

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.

The first step is to divide the computational domain into a fine grid of finite volumes or 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.
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.
The resulting discretized equations can be developed into a system of linear algebraic equations, such as the approach outlined previously with the FVM. This system is solved iteratively at each time step until convergence is achieved. Implicit methods involve solving a system of equations at each time step, while explicit methods update the solution based on the current values.
Implement the selected turbulence model within the selected numerical framework. As previously explained, this step involves solving additional equations for the turbulence quantities using a favorite turbulence model.
Apply appropriate boundary conditions to represent the physical boundaries of the problem. This approach may include specifying inlet, wall, and outlet conditions.
Iterate through the time steps until a steady state or a time-accurate solution is reached. Use a pressure-velocity coupling algorithm, such as SIMPLE (Semi-Implicit Method for Pressure Linked Equations); the SIMPLEC is a modified form.
Monitor convergence and adjust the solution until the changes between iterations are within acceptable tolerances or other criteria.
Analyze and post-process the numerical results to extract relevant information, such as flow velocities, pressure distributions, Reynolds stresses, TKE, and other turbulence quantities.

The examples below show RANS predictions of the wing tip vortex roll-up made using two different turbulence models. Even though both solutions give similar bulk flow properties, they obtain different detailed quantitative properties and turbulence levels. Which one is “correct” can only be established with reference to experiments, assuming that measurements are available.

RANS predictions of the rollup and formation of the wing tip vortex using two different turbulence models. The differences may have to be reconciled with results from experiments.

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

(38) $\begin{equation*} \overline{\vec{V}} \bigcdot \nabla \overline{\vec{V}} = -\left(\frac{1}{\varrho}\right) \nabla \overline{p} + \nu \nabla^2 \overline{\vec{V}} \end{equation*}$

Using forward or upwind differencing for the convective term gives

(39) $\begin{equation*} \left(\overline{u} \frac{\partial \overline{u}}{\partial x} + \overline{v} \frac{\partial \overline{u}}{\partial y}\right)_P \approx \frac{1}{\Delta x}(F_E - F_W) + \frac{1}{\Delta y}(F_N - F_S) \end{equation*}$

where $F_E$ , $F_W$ , $F_N$ , and $F_S$ represent the convective fluxes in the $x$ -direction at the eastern, western, northern, and southern faces of the control volume $P$ , respectively. The basic concept is shown in the figure below.

The concept of a finite-volume differencing scheme.

Using central differencing for the diffusive term gives

(40) $\begin{equation*} \nu \nabla^2 \overline{u} \approx \nu \left(\frac{\overline{u}_E - 2\overline{u}_P + \overline{u}_W}{\Delta x^2} + \frac{\overline{u}_N - 2\overline{u}_P + \overline{u}_S}{\Delta y^2}\right) \end{equation*}$

Using central differencing for the pressure gradient term gives

(41) $\begin{equation*} -\frac{1}{\varrho} \left( \frac{\partial \overline{p}}{\partial x} \right) \approx -\frac{1}{\varrho}\left(\frac{\overline{p}_E - \overline{p}_W}{2 \Delta x}\right) \end{equation*}$

Combining these terms for the $x$ component of the momentum equation gives

(42) $\begin{equation*} \overline{u}_P \left(\frac{1}{\Delta x}(F_E - F_W) + \frac{1}{\Delta y}(F_N - F_S)\right) = -\frac{1}{\varrho} \left(\frac{\overline{p}_E - \overline{p}_W}{2\Delta x}\right) + \nu \left(\frac{\overline{u}_E - 2\overline{u}_P + \overline{u}_W}{(\Delta x)^2} + \frac{\overline{u}_N - 2\overline{u}_P + \overline{u}_S}{(\Delta y)^2}\right) \end{equation*}$

Similar steps would be followed for the $z$ component of the momentum equation and 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 also include preprocessing tools for mesh generation and turbulence model setup and postprocessing tools for visualization and analysis of the results. The choice of numerical method and solver depends on factors such as the complexity of the flow, the available computational resources, and the desired resolution and accuracy outcomes.

Flowfield Predictions

While the CFD prediction of the flow about an entire aircraft may still be beyond the state-of-the-art, as shown in the figure below, the confidence in such predictions in terms of lift and drag is still 90% or better. The resolution of local flow properties may not be as good, but certainly of the fidelity that can be used for design refinements, e.g., tailoring the wing’s shape in the presence of the engines. Remember also that the turbulence models used for CFD are postdictive, i.e., after-the-fact models, so they do not lead to any extra understanding of the mechanisms of turbulence. Often, solving the flow field numerically using RANS leads to solutions that generally look like turbulence, but are just an artifact of how mathematics is used in its simulation. Nevertheless, remarkable details of the flows around complete aircraft can now be computed. In the end, however, all CFD solutions must be considered tentative, and must be validated against experiments or another benchmark, one way or another.

CFD solution for the flow about a complete aircraft. The confidence in prediction may be as high as 90%. (DLR image.)

CFD predictions of the flows about rotating wing aircraft are much more challenging because of their higher 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 the images below. The results show the tip vortices from the rotor blades, which are identified in terms of their vorticity (spin). In hover, each rotor lays down intertwining helicoidal vortices with relatively low pitch, which convect below the rotor and quickly diffuse and break up to form a turbulent jet. The proximity of the wake vortices to the rotor leads to high induced power requirements to drive the rotors. In forward flight, the rotors produce helical wake vortices with much greater pitch, resulting in lower induced effects and lower power for flight.

RANS calculation of the flow of NASA’s six-passenger tilt-wing concept for urban air mobility in hover. In this case, confidence in prediction may only be about 70%. (NASA public domain images/Patricia Ventura Diaz.)

These results reveal the flow’s complexity for a tilting multi-rotor configuration, where the flows from many rotors interact with each other, the wing, and the fuselage. Such RANS solutions can be beneficial, for example, in determining the issues with rotor/airframe aerodynamic interactions, which can produce highly variable aerodynamic forces and moments on the aircraft when transitioning from hover to forward flight and vice versa.

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 to extremely fine resolution is now prevented to some extent, allowing computational costs to be reduced by orders of magnitude. If the flow is assumed to be inviscid, then the Navier-Stokes equations reduce to a simplified set called 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

(43) $\begin{equation*} \frac{\partial \varrho}{\partial t} + \frac{\partial (\varrho \, u_j)}{\partial x_j} = 0 \quad \hbox{(Mass conservation)} \end{equation*}$

(44) $\begin{equation*} \frac{\partial (\varrho \, u_j)}{\partial t} + \frac{\partial(\varrho \, u_i \, u_j)}{\partial x_j} = -\frac{\partial p}{\partial x_i} \quad \hbox{(Momentum conservation)} \end{equation*}$

(45) $\begin{equation*} \frac{\partial(\varrho \, h)}{\partial t} + \frac{\partial (\varrho \, h \, u_j)}{\partial x_j} = \frac{\partial p}{\partial t} + u_j \left( \frac{\partial p}{\partial x_j} \right) \quad \hbox{(Energy conservation)} \end{equation*}$

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

(46) $\begin{equation*} \frac{\partial {\bf Q}}{\partial t} + \frac{\partial {\bf E}}{\partial x} + \frac{\partial {\bf F}}{\partial y} + \frac{\partial {\bf G}}{\partial z} = 0 \end{equation*}$

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

Vorticity in an inviscid flow, such as that governed by the Euler equations, should convect freely without any dissipation or diffusion. However, discretizing the Euler equations still leads to the creation of nonphysical numerical dissipation and diffusion of vorticity. This behavior is an artifact of discretizing the spatial derivatives using Taylor’s series approximations, which always results in some residual error that mimics the effects of viscosity. For instance, a wing modeled using the Euler equations should not generate lift because this process occurs through viscosity. In practice, however, the numerical solution of the Euler equations introduces an artificial source of viscosity, even a tiny amount of which will allow the essentially viscous flow problem over a wing to be modeled using an inviscid set of equations.

Despite the limitations of using the Euler equations, Euler-based CFD methods have helped the aircraft industry gain new insight into many complex flow problems, and have allowed the limitations of simpler models to be better understood. Although Euler-based methods will eventually be adopted as the most straightforward set of equations that can yield predictions of the flow around an aircraft, the fundamental limitations of the approach 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 the figure below. 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, e.g., because it creates the majority of drag experienced on the surface. To model the flow in this layer, a simplified subset of the Navier-Stokes equations can be derived, which exploits the fact that the length scale of the variation of the flow variables is very much smaller in a direction normal to the surface than it is parallel to the surface.

The boundary layer is a region of reduced velocity near a surface profoundly affected by viscosity.

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 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 proved very useful in airfoil design but can also help predict the aerodynamics of entire aircraft.

The basic principle first solves the “outer” zone as a problem and then solves an inviscid flo the “inner” zone in terms of 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, though, where separated flows are encountered because, under these conditions, the assumption of the boundary layer approach becomes invalidated.

Consider this type of situation in two dimensions. If $x$ is in the direction parallel to the surface and $y$ 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

(47) $\begin{eqnarray*} &&\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \quad \hbox{(continuity)}\\[8pt] &&u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = -\frac{1}{\varrho} \left( \frac{\partial p}{\partial x} \right) + \nu \frac{ \partial^2 u}{\partial y^2} \quad \hbox{($x$ component of momentum)} \end{eqnarray*}$

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

(48) $\begin{equation*} -\frac{\partial p}{\partial y } = 0 \end{equation*}$

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

(49) $\begin{equation*} \frac{\partial U_e}{\partial t} + U_e \frac{\partial U_e}{\partial x} = -\frac{1}{\varrho} \left( \frac{\partial p}{\partial x} \right) \end{equation*}$

where $U_e(x, t)$ is the “edge” velocity just outside the boundary layer. At the edge of the boundary layer $y=\delta$ , where $\delta$ is the boundary layer thickness^[8]

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

(50) $\begin{equation*} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = U_e \frac{\partial U_e}{\partial x} + \frac{\mu}{\varrho} \left( \frac{\partial^2 u}{\partial y^2}\right) \end{equation*}$

Along with the boundary layer equations is the concept of a displacement thickness, $\delta^*$ , which is defined as

(51) $\begin{equation*} \delta^{*} = \int_0^\delta \left( 1 - \frac{u}{U_e} \right) dy \end{equation*}$

This quantity allows the shape of the surface to be displaced away from the original surface by $\delta^{*}$ to yield the same flow rate as in the viscous flow with the boundary layer present. This is the same as saying that the outer inviscid flow about the body can be obtained by solving the equations about a new shape obtained by adding the displacement thickness to the body’s actual shape. This approach is an example of a so-called weak interaction method. Although still an approximation, this approach yields more representative flow solutions than if the viscous effects were not included.

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

(52) $\begin{equation*} \varrho \, \frac{d}{dx} (U_e^2 \, \theta ) + \varrho \, \delta^* U_e \left( \frac{d U_e}{d x} \right) = \tau_w \end{equation*}$

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

(53) $\begin{equation*} \theta = \int_0^\delta \left( 1 - \frac{u}{U_e} \right) \left( \frac{u}{U_e} \right) dy \end{equation*}$

This equation is very general because there is no assumption about the relationship between $\tau_w$ and the velocity gradient, so it will apply to both laminar or 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 $\delta$ , $\delta^*$ , $\theta$ , and $\tau_w$ . While in some cases, solutions to the momentum integral equation can be found in closed form, in most cases, numerical solutions are usually obtained.

Potential Flow Equations

A further simplification to the governing equations of the flow can be achieved by assuming that the flow is irrotational, i.e., the condition of irrotationality being $\nabla \times \mathbf{V} = 0$ . 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, $\phi$ , such that

(54) $\begin{equation*} \mathbf{V} \equiv \vec{V} = \nabla \phi \end{equation*}$

In potential flow theory, a fluid is considered 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

(55) $\begin{equation*} \nabla^2 \phi - \frac{1}{a^2} \left[ \frac{\partial^2 \phi}{\partial t^2} + \frac{\partial}{\partial t}(V^2) + \vec{V} \bigcdot \nabla \left( \frac{V^2}{2} \right) \right] = 0 \end{equation*}$

where $V = |\mathbf{V}| = |\vec{V}|$ . This form of the potential equation is usually called the full-potential equation because no terms have been neglected in its derivation from the Euler equation. For many flow problems, all the terms will be important if the characteristic nonlinear variation of the flow properties with Mach number is to be represented properly. However, other assumptions 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

(56) $\begin{equation*} \left( 1 - M_\infty^2 \right) \frac{\partial^2 \phi}{\partial x^2} + 2 M_\infty^2 \frac{\partial \phi}{\partial x} \frac{\partial^2 \phi}{\partial x \partial y} + \left( 1 + M_\infty^2 \frac{\partial \phi}{\partial x} \right) \frac{\partial^2 \phi}{\partial y^2} - \frac{1}{a_\infty^2} \frac{\partial^2 \phi}{\partial t^2} = 0 \end{equation*}$

For steady flows, then the TSD equation is

(57) $\begin{equation*} \left(1 - M_\infty^2 \right) \frac{\partial^2 \phi}{\partial x^2} + 2 M_\infty^2 \frac{\partial \phi}{\partial x} \left( \frac{\partial^2 \phi}{\partial x \, \partial y} \right) + \left(1 + M_\infty^2 \frac{\partial \phi}{\partial x}\right) \frac{\partial^2 \phi}{\partial y^2} = 0 \end{equation*}$

or in terms of the generalized coefficients $A$ , $B$ , and $C$ , then

(58) $\begin{equation*} A \frac{\partial^2 \phi}{\partial x^2} + B \frac{\partial^2 \phi}{\partial x \, \partial y} + C \frac{\partial^2 \phi}{\partial y^2} = 0 \end{equation*}$

where

(59) $\begin{equation*} A = 1 - M_\infty^2 \quad \mbox{,} \quad B = 2M_\infty^2 \left( \frac{\partial \phi}{\partial x} \right) \quad \mbox{and} \quad C = 1 + M_\infty^2 \left( \frac{\partial \phi}{\partial x} \right) \end{equation*}$

Further simplifying by collecting and isolating the cross-terms and low-order derivatives, then

(60) $\begin{equation*} A \frac{\partial^2 \phi}{\partial t^2} + B \frac{\partial^2 \phi}{\partial x \, \partial t} = F_x + \frac{\partial^2 \phi}{\partial z^2} + C \frac{\partial^2 \phi}{\partial y^2} \ + \text{~cross terms and low-order derivatives} \end{equation*}$

where $F_x$ represents a forcing term that refers to contributions that account for external influences or boundary conditions that modify the behavior of the flow.

The forcing terms can arise from the following factors:

The shape of the airfoil or body causes a disturbance to the freestream flow, generating a source of forcing in the equation. For a slender airfoil in transonic flow, the body surface imposes a boundary condition that relates the velocity potential to the body geometry. For example, the boundary condition on the surface of an airfoil can be written as
$\frac{\partial \phi}{\partial y} = \frac{\partial h}{\partial x}$

where $h(x)$ represents the airfoil’s thickness as a function of $x$ .
If the body is in unsteady motion (e.g., pitching, plunging, or otherwise oscillating), these motions contribute additional time-dependent forcing terms to the TSD equation. For instance, for a body oscillating with frequency $\omega$ , the unsteady potential can include terms like
$\phi(x, y, t) = \phi_s(x, y) + \phi_t(x, y, t)$

where $\phi_t$ accounts for the unsteady forcing because of the body’s motion.
In some cases, boundary conditions at infinity or disturbances propagating from other regions in the flow can introduce forcing terms.
Variations in the freestream flow (e.g., gusts, compressibility changes) can also act as a forcing term

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^[9] solved numerically using finite-difference methods, which helped establish CFD as a practical field of study. Naturally, the TSD equation is particularly important in the study of transonic flows where the flow speeds are around the speed of sound, causing unique challenges because of the presence of both subsonic and supersonic flow 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 all the higher-order terms are neglected, the TSD equation reduces to a wave equation of the form

(61) $\begin{equation*} \frac{1}{a_{\infty}^2} \left( \frac{\partial^2 \phi}{\partial t^2} \right) = \nabla^2 \phi \end{equation*}$

where $a_{\infty}$ is the free-stream speed of sound. In this simplified form, the equation describes the propagation of small disturbances in a compressible flow. This reduction to a wave equation is useful for understanding the basic characteristics of wave propagation in a fluid, and forms the basis for more complex analyses in transonic flow regimes. The full TSD equation includes additional nonlinear terms to account for the effects of transonic flow, but this simplified form captures the essence of wave propagation.

A final simplification is to assume steady, incompressible flow ( $a_{\infty} \rightarrow \infty$ ), in which case Laplace’s equation governs the velocity potential, i.e.,

(62) $\begin{equation*} \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = \frac{\partial^2 \phi}{\partial z^2} = \nabla^2 \phi = 0 \end{equation*}$

A feature of this latter equation is that it is linear, which allows complex flows to be synthesized by combining two or more elementary flows. For example, if $\phi_1$ and $\phi_2$ are solutions to the Laplace equation, then $\phi_1 + \phi_2$ is also a solution. This property allows complex flow patterns to be constructed by combining simpler solutions, such as uniform flow, sources, sinks, and vortices. The Laplace equation is the governing equation for incompressible flows and underlies what are known as surface singularity or “panel” methods.

The Laplace equation also applies to unsteady incompressible flows because of the incompressibility assumption implying an infinite speed of sound. The unsteady nature of the flow is accounted for in the potential function $\phi(x, y, z, t)$ , which now depends on time $t$ . The instantaneous adaptation of the flow field to changes in boundary conditions is a consequence of the idealization of the infinite speed of sound in incompressible flow, meaning pressure disturbances propagate instantly to all parts of the flow field.

Surface Singularity Methods

Panel methods are widely used for calculating potential flows around bodies, especially in cases where the flow remains incompressible and inviscid. These methods provide valuable insights into aerodynamic characteristics, such as lift and drag. They are particularly useful in preliminary design stages and in situations where the flow can be approximated as incompressible. Panel methods are based on the assumption of small disturbances and incompressible irrotational flow, for which the governing equation is Laplace’s equation (Eq. 62) 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 hence the aerodynamic loading on the body, can be calculated.

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, which is a popular tool for flow predictions about 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 $N$ discrete panels, and each panel $i$ contains a singularity distribution that can be described by a velocity potential function $\phi_i$ . This means that the total potential induced by all the singularities in the flow is

(63) $\begin{equation*} \phi = \sum_{i=1}^N \gamma_i \, \phi_i \end{equation*}$

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

(64) $\begin{equation*} \big( \nabla \phi(x_j) + V_\infty \big) \bigcdot \vec{n_j} = 0 \end{equation*}$

where $\vec{n_j}$ is the outward pointing normal unit vector at the $j$ th collocation point. Usually, the collocation points are selected as the midpoints of the panels.

A panel method solution based on a linear vorticity distribution over each panel.

Suppose the panels contain vortex singularities that vary linearly in strength over the panel length^[10], then with $N$ panels, there are $N$ equations containing $N+1$ unknowns. 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., $\gamma (TE) =0$ , which is implemented by setting $\gamma_{1} + \gamma_{N+1} = 0$ . 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

(65) $\begin{equation*} \left\{ \begin{array}{c} \gamma_1 \\[6pt] \gamma_2 \\[6pt] \bigcdot \\[6pt] \bigcdot \\[6pt] \gamma_N \\[6pt] \gamma_{N+1} \end{array} \right\} = \left[ \begin{array}{ccccc} A_{1,1} & A_{1,2} & \bigcdot & \bigcdot & A_{1,N+1} \\[6pt] A_{2,1} & A_{2,2} & \bigcdot & \bigcdot & A_{2,N+1} \\[6pt] \bigcdot & \bigcdot & \bigcdot & \bigcdot & \bigcdot \\[6pt] \bigcdot & \bigcdot & \bigcdot & \bigcdot & \bigcdot \\[6pt] A_{N,1} & A_{N,2} & \bigcdot & \bigcdot & A_{N,N+1} \\[6pt] 1 & 0 & \bigcdot & 0 & 1 \end{array} \right]^{-1} \left\{ \begin{array}{c} -V_{\infty}\sin \alpha \bigcdot n_1 \\[6pt] -V_{\infty}\sin \alpha \bigcdot n_2 \\[6pt] \bigcdot \\[6pt] \bigcdot \\[6pt] -V_{\infty}\sin \alpha \bigcdot n_N \\[6pt] 0 \end{array} \right\} \end{equation*}$

where ${\bf A} = A_{i,j} = \nabla \phi_i(x_j) \bigcdot \vec{n_j}$ is a $(N+1)\times (N+1)$ matrix of influence coefficients. An influence coefficient $A_{i,j}$ is effectively the normal component of the velocity induced at the collocation point on panel $j$ by the unit strength singularity distribution of $\phi_i$ on panel $i$ ; it depends on the relative distance between and orientation of the two panels.

Notice that solving Eq. 65 involves finding the inverse of the matrix of influence coefficients, for which many efficient algorithms are available. Besides calculating the influence coefficients themselves, this inversion is the primary computational expense involved in the panel method. If the angle of attack, $\alpha$ , is included in the vector on the right-hand side of Eqs. 65, the influence coefficients need only be calculated once. Once the singularity strengths $\gamma_i$ are known, the loading on the airfoil follows immediately 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 might be represented by $N$ small, flat, quadrilateral panels onto which singularities in the form of sources, sinks, and/or doublets are placed. The principles are very similar to the two-dimensional case, but the size of the numerical problem becomes much larger. Assume that only source singularities are used for illustration. If the singularity strength on the $i$ th panel is $\sigma_i$ , then the velocity induced by this panel at the control point on the $j$ th panel can be written as $A_{i,j} \; \sigma_i$ . At the $j$ th control point, the component normal to the surface of the velocity induced by all $N$ sources is then

(66) $\begin{equation*} V_j = \sum_{i=1}^{N} A_{i,j} \; \sigma_i {\quad} \mbox{for $i \ne j$} \end{equation*}$

The boundary condition of flow tangency is imposed at the control point of each panel to allow the unknown singularity strengths, $\sigma_i$ , to be evaluated. So, for instance, if $\vec{V}_{{\infty}_j}$ is the velocity at panel $j$ as the component of $\vec{V}_{{\infty}_j}$ normal to panel $j$ will be $\vec{V}_{{\infty}_j} \bigcdot \vec{n}_j$ , where $\vec{n}_j$ is the unit normal vector to panel $j$ . Therefore, at the $j$ th panel

(67) $\begin{equation*} \sum_{i=1}^{N} A_{i,j} \; \sigma_i + \vec{V}_{{\infty}_j} \bigcdot \vec{n}_j = 0. \end{equation*}$

One equation such as this applies to each of the $N$ panels.

When applied to all panels, this latter approach leads to a set of linear simultaneous equations that can be used to solve for the singularity strengths, $\sigma_i$ , using standard numerical methods of linear algebra. The governing equations can be written in matrix form as

(68) $\begin{equation*} [ A_{ij} ] \{ \sigma_{i} \} = {-}\{\vec{V}_{\infty} \bigcdot n_{j} \}, \end{equation*}$

which has the solution

(69) $\begin{equation*} \left\{ \sigma_{i} \right\} = \left[ A_{i,j} \right]^{-1} \left\{ -\vec{V}_{\infty} \bigcdot n_{j} \right\} \end{equation*}$

Therefore, the solution involves finding the inverse of the influence coefficient matrix. Once the values of $\sigma_{i}$ are obtained, the velocity components tangential to the panels can be obtained, and hence, the pressure distribution over the fuselage follows using Bernoulli’s equation.

Usually, thousands of panels are required to resolve adequately the complex shape of an entire aircraft, and care must be taken, for instance, to concentrate a sufficient number of panels in regions of higher surface curvature, as shown in the figure below. Because the numerical cost of panel methods is of order $N^3$ , their use is by no means inexpensive, even on a modern computer. However, they are considerably less computer-intensive than most CFD solutions based on Navier-Stokes or Euler methods. One of the advantages of panel methods is that they do not suffer the numerical diffusion and dissipation associated with finite-difference methods. This means that the downstream wake, shown in the image below, is self-preserving despite its relatively coarse approximated form. Suppose care is taken to make a suitable distribution of panels on the surface of the fuselage, with more panels being assigned in regions of expected higher pressure gradients. In that case, a three-dimensional surface singularity method can be used to obtain entirely credible predictions of the loading on an airframe, certainly to the fidelity required for preliminary design.

A panel method solution for an entire aircraft, including the downstream wake developments.

As an enhancement to the primary panel method, the solution of the momentum integral equation from boundary layer theory allows displacement corrections to be made to solutions obtained from the panel method, thereby allowing viscous effects (but not flow separation) to be modeled. However, the effects of flow separation can be included, albeit approximately, by extending the surface singularities off the airframe as free shear layers.

Displacement thickness corrections can enhance the primary panel method and improve its predictive ability.

Including boundary layer displacements in lift predictions for an airfoil enhances accuracy by accounting for viscous effects and flow separation. Classic potential flow methods, which ignore these factors, often overpredict lift. By incorporating boundary layer displacements, the viscous drag and its impact on the overall aerodynamic performance are considered, resulting in more precise lift predictions. Additionally, the boundary layer can influence the onset and extent of flow separation, significantly reducing lift at higher angles of attack.

Corrections for boundary layer displacement effects allow aerodynamic predictions to be extended to a higher angle of attack.

Moreover, the boundary layer affects the pressure distribution over the airfoil and alters its effective shape resulting from displacement thickness. These changes influence the flow field and, consequently, the lift generated. Including boundary layer displacements modifies the predicted lift to better match experiments by reflecting the actual pressure distribution and flow characteristics around the airfoil. This approach also considers the transition effects between laminar and turbulent flows, further refining the lift calculations for more representative predictions of aerodynamic performance.

Rarefied Gas Flows

Detailed information about low-density or rarefied flow environments is crucial for spacecraft, spanning from design to mission operations. Models of such environments must be developed for various applications, including assessing orbital decay and aerodynamic heating during re-entry to the Earth’s atmosphere or 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. In these conditions, the continuum hypothesis becomes invalid, 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 large compared to the dimensions of the interacting body, as illustrated in the schematic below. The value of the Knudsen number, $Kn$ , helps determine if a gas is rarefied or not, i.e.,

(70) $\begin{equation*} Kn = \frac{\lambda}{L} \end{equation*}$

where $\lambda$ is the mean free path, and $L$ is a characteristic length for the problem. If the value of $Kn$ exceeds 1.0, the gas can be considered rarefied.

Key concepts used to understand behavior of molecules in a rarefied gas.

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 the gas molecules follow are described by the Maxwell-Boltzmann distribution, which is an idealized but realistic probability assumption, as shown in the figure below. This model describes how likely it is to find an individual molecule with a particular speed at a given temperature. The speed most molecules have is given by

(71) $\begin{equation*} V_p = \sqrt{\frac{2 \, k_B \, T}{m}} \end{equation*}$

where $k_B$ is the Boltzmann constant, $T$ is the absolute temperature, and $m$ is the mass of the gas molecule.

The Maxwell-Boltzmann distribution gives the probability density of finding molecules at different speeds at a given temperature.

The Boltzmann constant relates the average kinetic energy of particles in a gas to the temperature of the gas; its value is approximately $1.380649 \times 10^{-23}$ Joules per Kelvin (J/K). At lower temperatures, the molecules have less energy. Therefore, the molecular speeds are lower and the distribution is over a smaller range. As the temperature of the molecules increases, the probability distribution flattens out because the molecules have greater energy at higher temperatures. So more of the molecules, on average, are moving faster. The average speed of all molecules is

(72) $\begin{equation*} \langle V \rangle = \sqrt{\frac{8 \, k_B \, T}{\pi \, m}} \end{equation*}$

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

(73) $\begin{equation*} V_{\text{rms}} = \sqrt{\frac{3 \, k_B \, T}{m}} \end{equation*}$

In thermal equilibrium, the velocities of gas molecules follow the Maxwell-Boltzmann distribution, characterized by most probable, average, and root mean square speeds. The distribution changes when the gas is not in equilibrium, such as when it is flowing or experiencing external forces. The Boltzmann equation describes how this distribution evolves and can be written as

(74) $\begin{equation*} \frac{\partial f}{\partial t} + \vec{V} \bigcdot \nabla f = -\frac{f - f_{\text{eq}}}{\tau} \end{equation*}$

where $f$ is the distribution function of molecules, $\vec{V}$ is their velocity, $f_{\text{eq}}$ is the equilibrium distribution, and $\tau$ is the relaxation time, i.e., the average time between molecular collisions. In a rarefied gas, this time is longer because the molecules are farther apart. Collisions between molecules drive the gas toward an equilibrium state, 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 is different from what would be predicted with a continuum model. An example of their application is shown below.^[11] These models are essential for understanding the behavior of gases in low-density environments where traditional continuum-based descriptions of fluid dynamics are inapplicable.

An atmospheric entry simulation of a spacecraft using a rarefied gas model. (NASA image.)

In general, rarefied gas models provide insights into the complex interactions between the gas molecules and the spacecraft, which 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 employed to simulate the Martian atmosphere’s effect on the spacecraft as it descended and landed on the surface. Such models helped engineers design effective heat shields and aerodynamic surfaces to ensure a safe landing.

Aerospace Structures

The Finite Element Method (FEM) is a numerical technique used extensively in engineering for analyzing structures and predicting their response to various physical forces such as aerodynamics. It involves discretizing complex structures into smaller, discrete elements connected at nodes. Depending on the geometry and complexity of the problem, elements like beams, trusses, triangles, or hexahedrons are employed. Similar to the process of CFD, the FEM utilizes grids, such as O-grids and H-grids, which are particularly useful in areas where grid refinement is crucial, like around structural cutouts such as windows or access panels, as shown in the example below. This approach ensures an accurate representation of stress concentrations and stress deformation patterns.

A mesh (grid) for use in a FEM, which will involve grid refinement in regions of higher structural loads and stress.

In practice, the FEM formulates the governing equations of the structure based on principles of energy minimization or virtual work, which is used to derive the element stiffness matrix and the system of equations to be solved. These equations are then solved numerically by assembling a global system of equations from contributions of individual elements, 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 without physical prototyping, making it indispensable in aerospace engineering for its accuracy, efficiency, and ability to handle complex geometries and different materials. The success in maximizing the strength-to-weight ratio of aerospace structures is mainly because of the use of the FEM in the process of design.

History

The origins of FEM can be traced back to the work done in the field of structural analysis during the 1940s. The work of Richard Courant in 1943 is often cited as one of the first instances where methods of FEM were employed. Courant used piecewise linear approximations over triangular subregions to solve torsion problems, which laid the groundwork for the method. The method later gained traction through the efforts of engineers and mathematicians. Ray Clough, a structural engineer, is credited with coining the term “finite element method” in a paper published in 1956. However, his significant contributions began in the late 1950s. 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 those 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, became 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 one of the notable early FEM programs.

FEM was extended beyond structural mechanics to other fields like 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 it became 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 like ANSYS, ABAQUS, and others 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 include efficient algorithms to handle large-scale problems, which inevitably means the use of different methods that are 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, dynamic analysis, thermal analysis, and more.

General Approach

The structural domain is divided into $N_e$ elements with $n_n$ nodes, the idea being shown in the figure below. The displacement within each element is approximated using shape functions $N_i(x, y, z)$ and nodal displacements $\mathbf{u}_i$ , i.e.,

(75) $\begin{equation*} \mathbf{u}(x, y, z) = \sum_{i=1}^{n_n} N_i(x, y, z) \, \mathbf{u}_i \end{equation*}$

For each element, the stiffness matrix $\mathbf{K}^e$ is given by

(76) $\begin{equation*} \mathbf{K}^e = \int_{\Omega_e} \mathbf{B}^T \, \mathbf{D} \, \mathbf{B} \, d\Omega_e \end{equation*}$

where $\mathbf{B}$ is the strain-displacement matrix, $\mathbf{D}$ is the material property matrix (e.g., containing Young’s modulus, $E$ , and Poisson’s ratio, $\nu$ ), and $\Omega_e$ is the element’s spatial domain. The use of matrix multiplication in Eq. 76 is implied.

The concept of the FEM is to discretize the structure into basic elements and combine material properties, forces, and displacements to represent the entire structure.

The global stiffness matrix $\mathbf{K}$ is then assembled from the element stiffness matrices, i.e.,

(77) $\begin{equation*} \mathbf{K} = \sum_{e=1}^{N_e} \mathbf{K}^e \end{equation*}$

The global force-displacement relationship is given by

(78) $\begin{equation*} \mathbf{K} \, \mathbf{U} = \mathbf{F} \end{equation*}$

where $\mathbf{U}$ is the global displacement vector and $\mathbf{F}$ 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 $\mathbf{K}$ is constructed by linear superposition of the element stiffness matrices. Solving the global system of equations $\mathbf{K} \, \mathbf{U} = \mathbf{F}$ efficiently, especially for large-scale problems, involves the use of various numerical methods.

2-D Triangular Element

A simple triangular element offers an opportunity to understand the FEM better. It will become apparent that the systematic discretized nature of the method lends itself naturally to the computer so that the steps can be readily programmed. The shape functions $N_i(x, y)$ are linear for a triangular element with three nodes. Let the coordinates of the nodes be $(x_1, y_1)$ , $(x_2, y_2)$ , and $(x_3, y_3)$ , as shown in the figure below.

The shape functions, in this case are given by straight lines, i.e.,

(79) $\begin{equation*} N_1(x, y) = \frac{1}{2A} (a_1 + b_1 x + c_1 y) \end{equation*}$

(80) $\begin{equation*} N_2(x, y) = \frac{1}{2A} (a_2 + b_2 x + c_2 y) \end{equation*}$

(81) $\begin{equation*} N_3(x, y) = \frac{1}{2A} (a_3 + b_3 x + c_3 y) \end{equation*}$

where $A$ is the area of the triangle given by

(82) $\begin{equation*} A = \frac{1}{2} \bigg| x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) \bigg| \end{equation*}$

The constants $a_i$ , $b_i$ , and $c_i$ are determined by the coordinates of the nodes, i.e.,

(83) $\begin{equation*} a_1 = x_2 y_3 - x_3 y_2, \quad b_1 = y_2 - y_3, \quad c_1 = x_3 - x_2 \end{equation*}$

(84) $\begin{equation*} a_2 = x_3 y_1 - x_1 y_3, \quad b_2 = y_3 - y_1, \quad c_2 = x_1 - x_3 \end{equation*}$

(85) $\begin{equation*} a_3 = x_1 y_2 - x_2 y_1, \quad b_3 = y_1 - y_2, \quad c_3 = x_2 - x_1 \end{equation*}$

The strain-displacement matrix $\mathbf{B}$ relates the nodal displacements to the strains within the element. For this triangular element, $\mathbf{B}$ is given by

(86) $\begin{eqnarray*} \mathbf{B} = \begin{bmatrix} \dfrac{\partial N_1}{\partial x} & 0 & \dfrac{\partial N_2}{\partial x} & 0 & \dfrac{\partial N_3}{\partial x} & 0 \\[8pt] 0 & \dfrac{\partial N_1}{\partial y} & 0 & \dfrac{\partial N_2}{\partial y} & 0 & \dfrac{\partial N_3}{\partial y} \\[12pt] \dfrac{\partial N_1}{\partial y} & \dfrac{\partial N_1}{\partial x} & \dfrac{\partial N_2}{\partial y} & \dfrac{\partial N_2}{\partial x} & \dfrac{\partial N_3}{\partial y} & \dfrac{\partial N_3}{\partial x} \end{bmatrix} \end{eqnarray*}$

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

(87) $\begin{equation*} \frac{\partial N_i}{\partial x} = \frac{b_i}{2A}, \quad \frac{\partial N_i}{\partial y} = \frac{c_i}{2A} \end{equation*}$

Therefore, the matrix $\mathbf{B}$ can be written as

(88) $\begin{equation*} \mathbf{B} = \frac{1}{2A} \begin{bmatrix} b_1 & 0 & b_2 & 0 & b_3 & 0 \\[6pt] 0 & c_1 & 0 & c_2 & 0 & c_3 \\[6pt] c_1 & b_1 & c_2 & b_2 & c_3 & b_3 \end{bmatrix} \end{equation*}$

For plane stress conditions (although plane strain is an option), the material property matrix $\mathbf{D}$ is

(89) $\begin{equation*} \mathbf{D} = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\[4pt] \nu & 1 & 0 \\[4pt] 0 & 0 & \dfrac{1 - \nu}{2} \end{bmatrix} \end{equation*}$

where $E$ is Young’s modulus and $\nu$ is Poisson’s ratio. Poisson’s ratio^[12] is a material property that describes the ratio of the transverse strain to the axial strain when a material is stretched.

The element stiffness matrix $\mathbf{K}^e$ is given by

(90) $\begin{equation*} \mathbf{K}^e = \int_{\Omega_e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, d\Omega_e \end{equation*}$

For a linear triangular element, this integral simplifies because $\mathbf{B}$ and $\mathbf{D}$ are constant over the element, i.e.,

(91) $\begin{equation*} \mathbf{K}^e = \mathbf{B}^T \mathbf{D} \mathbf{B} \ A \end{equation*}$

Substituting $\mathbf{B}$ and $\mathbf{D}$ gives

(92) $\begin{equation*} \mathbf{K}^e = \left( \frac{1}{2A} \begin{bmatrix} b_1 & 0 & b_2 & 0 & b_3 & 0 \\[6pt] 0 & c_1 & 0 & c_2 & 0 & c_3 \\[6pt] c_1 & b_1 & c_2 & b_2 & c_3 & b_3 \end{bmatrix} \right)^T \left( \dfrac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\[6pt] \nu & 1 & 0 \\[6pt] 0 & 0 & \dfrac{1 - \nu}{2} \end{bmatrix} \right) \left( \frac{1}{2A} \begin{bmatrix} b_1 & 0 & b_2 & 0 & b_3 & 0 \\[6pt] 0 & c_1 & 0 & c_2 & 0 & c_3 \\[6pt] c_1 & b_1 & c_2 & b_2 & c_3 & b_3 \end{bmatrix} \right) \ A \end{equation*}$

Multiplying out these matrices and integrating them over the area gives

(93) $\begin{equation*} \mathbf{K}^e = \frac{E}{4A(1 - \nu^2)} \begin{bmatrix} b_1^2 + \nu c_1^2 & (1-\nu)b_1 c_1 & b_1 b_2 + \nu c_1 c_2 & (1-\nu)b_2 c_1 & b_1 b_3 + \nu c_1 c_3 & (1-\nu)b_3 c_1 \\[8pt] (1-\nu)b_1 c_1 & \nu b_1^2 + c_1^2 & (1-\nu)b_1 c_2 & b_1 b_2 + \nu c_1 c_2 & (1-\nu)b_1 c_3 & b_1 b_3 + \nu c_1 c_3 \\[8pt] b_1 b_2 + \nu c_1 c_2 & (1-\nu)b_1 c_2 & b_2^2 + \nu c_2^2 & (1-\nu)b_2 c_2 & b_2 b_3 + \nu c_2 c_3 & (1-\nu)b_3 c_2 \\[8pt] (1-\nu)b_2 c_1 & b_1 b_2 + \nu c_1 c_2 & (1-\nu)b_2 c_2 & \nu b_2^2 + c_2^2 & (1-\nu)b_2 c_3 & b_2 b_3 + \nu c_2 c_3 \\[8pt] b_1 b_3 + \nu c_1 c_3 & (1-\nu)b_1 c_3 & b_2 b_3 + \nu c_2 c_3 & (1-\nu)b_2 c_3 & b_3^2 + \nu c_3^2 & (1-\nu)b_3 c_3 \\[8pt] (1-\nu)b_3 c_1 & b_1 b_3 + \nu c_1 c_3 & (1-\nu)b_3 c_2 & b_2 b_3 + \nu c_2 c_3 & (1-\nu)b_3 c_3 & \nu b_3^2 + c_3^2 \end{bmatrix} \end{equation*}$

Finally, as previously described, the global stiffness matrix $\mathbf{K}$ is then assembled by summing up the contributions of all of the element stiffness matrices $\mathbf{K}^e$ based on the connectivity of the nodes. Boundary conditions are also applied to ensure that the constraints (e.g., 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 geometry of the structure 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, especially in stress concentration areas, such as between the wing and the fuselage or around engine intakes, as shown in the figure below. The mesh quality is checked to ensure accurate results, and boundary conditions and loads are applied. Finally, the mesh is exported for use in FEM simulations.

Diagram of a structural FEM model of an IAR 99 aircraft. — A structural FEM model of an aircraft. Like cockpit windows, regions with cutouts need particular attention to avoid local stress concentrations.

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

Predictions from the FEM

The FEM has become a standard tool for structural load predictions in engineering; an example of a solution from a FEM is shown below for a wing and the spar connected to the fuselage. NASTRAN is a widely used FEM software, which includes pre-processing (meshing and boundary conditions), solving (numerical methods), and post-processing (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, such as optimization algorithms and multi-physics 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.

Aeroelasticity

Aeroelasticity is a specialized field in aerospace engineering that deals with the interactions between aerodynamic forces, structural dynamics, and elastic deformations of aircraft or other structures. It examines how these interactions influence the flexibility, stability, and safety of aircraft and other flight vehicles subjected to aerodynamic loads. While flutter (self-excited oscillations caused by interactions between aerodynamic forces and structural vibrations) can be predicted analytically in some cases, tightly integrating CFD and the FEM is essential for a comprehensive aeroelastic analysis of an entire flight vehicle. The complication is that there are many possible modes of aeroelastic behavior on an aircraft, not just limited to the wings.

History

The origins of aeroelasticity can be traced to aircraft developments in the early 1900s. The Wright brothers and other early aviation pioneers encountered aeroelastic phenomena such as wing flutter, a dynamic instability caused by the interaction of aerodynamic forces and structural dynamics. However, they did not understand these effects then. In 1902, Samuel Langley tried to fly an aircraft, but his monoplane’s wood and fabric wings were not stiff enough and broke off, i.e., a form of aeroelastic divergence. The Wright brothers experimented with different structures and knew how to make a more rigid structure in the form of a biplane, free of these aeroelastic effects. The catastrophic failure of the horizontal tail on a Handley Page bomber in 1916 from flutter suggested much was still to be learned about aircraft design.

Flutter became a critical concern as aircraft designs evolved, particularly with the transition to monoplanes. Monoplanes, with their single-wing configuration and lower structural stiffness, highlighted the importance of understanding and mitigating flutter during aircraft design. “The Blue Max” is a 1966 film set in WWI that tells the story of a German fighter pilot, Bruno Stachel, who aims to earn the prestigious Pour le Mérite, also known as the “Blue Max,” for outstanding acts of bravery in aerial combat. While wing flutter is not a central theme in the movie, the film depicts the increased dangers of structural failure that pilots faced when transitioning from flying biplanes to early monoplanes. The rapid advancement in aircraft design and increased flight speeds after this period led to more frequent and severe aeroelastic problems. Engineers began systematically studying these issues, formulating aeroelastic theories and methods for predicting flutter.

The development of high-speed, high-altitude aircraft during WWII and the subsequent Cold War era necessitated a deeper understanding of structural aeroelasticity and flutter. The introduction of jet engines and swept wings presented increased challenges in predicting aeroelastic effects at higher airspeeds and Mach numbers. This period saw significant theoretical and experimental advancements, including establishing formal mathematical criteria to predict the onset of flutter and divergence, including unsteady aerodynamic effects. The work of Theodore Theodorsen, among others, laid the groundwork^[13] for modern aeroelastic analysis. The space race and the development of supersonic and hypersonic vehicles further pushed the boundaries of aeroelastic research. More comprehensive computational methods emerged, allowing for more sophisticated and accurate analyses. Wind tunnel testing and the use of aeroelastically “Froude” scaled models remained crucial for validation.

With the advent of powerful computers and advanced numerical methods, such as the FEM, aeroelastic modeling has become highly sophisticated. Modern software tools enable detailed simulation of aeroelastic phenomena, including flutter, divergence, and control reversal. These tools are used extensively in the design of modern aircraft and spacecraft. Aeroelasticity has become a mature field that integrates advanced computational techniques with experimental methods to ensure the safety and performance of all aerospace flight vehicles.

Method of Analysis

The field of aeroelastics is not an easy one to master, but understanding the basics, and especially the underlying physics in a qualitative context, is the place to start. In a fairly general way, the structural dynamics of an aeroelastic system can be described by

(94) $\begin{equation*} M \,\overbigddot{q} + C \, \overbigdot{q} + K \, q = F_a \end{equation*}$

where $M$ is the mass matrix, $C$ is the damping matrix, $K$ is the stiffness matrix, $q$ is the displacement vector, and $F_a$ is the aerodynamic force vector. Aerodynamic forces can be computed using CFD, although panel methods or unsteady airfoil theory might also be used while appreciating their limitations. The aerodynamic force $F_a$ on the structure using CFD or panel methods is usually extracted from the pressure distribution, $p(x, y, z, t)$ , that acts over the surface, i.e.,

(95) $\begin{equation*} F_a (t) = \int_S p(t) \, d\vec{S} \end{equation*}$

where $\vec{S}$ is the unit normal vector area. Other methods, however, might be used to approximate aerodynamics, including linear theory, classic unsteady airfoil theory, panel methods, etc. It is not unusual for flutter calculations to be made using various aerodynamic methods (and with different assumptions) and for the results to be compared.

Combining the aerodynamic and structural models then gives

(96) $\begin{equation*} M \, \overbigddot{q} + (C + B) \, \overbigdot{q} + (K + A) \, q = 0 \end{equation*}$

where $A$ and $B$ are called aerodynamic influence matrices. Flutter is a dynamic instability, and its condition is determined by solving

(97) $\begin{equation*} \underbrace{K + A - \omega^2 M}_{\rm Real~part} \ + \ i \, \underbrace{\omega (C + B)) \Phi}_{\rm Imaginary~part} = \Re (\omega) + i \, \Im (\omega) = 0 \end{equation*}$

where $\omega$ is the complex eigenvalue. Notice that $\Phi$ typically represents the mode shape or eigenvector associated with the eigenvalue. In the context of flutter analysis, $\Phi$ represents the deformation pattern of the structure at the flutter condition. Each eigenvalue has an associated eigenvector that describes how different points in the structure move relative to each other during an oscillation.

The “real” part, $\Re (\omega )$ , represents the balance between stiffness (including aerodynamic stiffness) and inertial forces, which is associated with the natural frequency of the system. When $\Re (\omega ) > 0$ , it indicates oscillatory motion. If $\Re (\omega ) = 0$ , it means a static aeroelastic response. For example, divergence is a static instability analyzed by setting dynamic terms to zero, i.e.,

(98) $\begin{equation*} (K + A) q = 0 \end{equation*}$

which means $q \rightarrow \infty$ and the structure would theoretically deform infinitely under the aerodynamic loads. The critical divergence airspeed is found when the determinant of $(K + A) = 0$ . The “imaginary” part, $\Im (\omega )$ , represents the balance between damping (including aerodynamic damping) and the rate of change of displacement, which is associated with the damping of the system. If $\Im (\omega) > 0$ , it means positive damping (energy dissipation), while $\Im (\omega) < 0$ means negative damping (energy addition), which will extract energy from the flow and can lead to flutter.

The field of aeroelasticity also considers concepts such as limit cycle oscillations (LCO), where periodic self-sustained oscillations occur from nonlinearities in the aerodynamic or structural response. Engineers use sophisticated computational methods and experimental testing to predict and mitigate these aeroelastic phenomena, ensuring aircraft and spacecraft designs meet stringent safety, performance, and certification standards in diverse flight conditions.

Flutter Solution Process

The aeroelastic process must be conducted either by time-marching or by an iterative process. An iterative coupling process will typically proceed along these lines:

Define the geometry and grids for FEM and CFD. They will differ because the grids required to solve the aerodynamics have different requirements than the structural shape.
$\[\]$
Apply the material properties and boundary conditions in the FEM. Using isotropic material such as aluminum will have different properties to composites.
$\[\]$
Set up the CFD with the defined flow conditions or use any other selected aerodynamic model, e.g., a panel method.
$\[\]$
Compute the aerodynamic forces, $F_a$ , from the selected aerodynamic model.
Update the FEM model with aerodynamic forces $F_a$ .
Iterate between CFD and FEM to capture aeroelastic interactions, i.e., couple the aerodynamics and structural dynamics in the form
$\[\]$

(99) $\begin{equation*} M \, \overbigddot{q} + (C + B) \, \overbigdot{q} + (K + A) \, q = 0 \end{equation*}$

$\[\]$
Repeat the FEM analysis for the resulting structural response on the deformed structure.
$\[\]$
Conduct CFD simulations on the deformed structure to obtain new aerodynamic loads.
$\[\]$
Exchange data between FEM and CFD to formalize the coupled analysis.
$\[\]$
Analyze displacements, stresses, and aerodynamic forces.
$\[\]$
Evaluate the results against the stability criteria for flutter or divergence.

By integrating CFD and FEM, engineers can accurately predict and analyze aeroelastic phenomena, ensuring the safety and performance of aerospace structures. For example, this process has been applied to large aircraft such as the Boeing 747. These results provide valuable information on the stiffness of airframe components, helping to prevent aeroelastic and flutter issues before the first flight. However, flight testing is always essential to validate these calculations and ensure the aircraft is flutter-free across its operational flight envelope.

A model created by ZONA Technology shows aircraft deformations due to flutter. — A mathematical finite-element and aerodynamic model of a Boeing 747 shows the potential structural deformations that can occur because of flutter. (NASA/ZONA Technology.)

Flutter of a Cantilevered Wing

The understanding of the flutter mechanism is usually hidden in the relatively complex mathematics or inside some code that is being used as a “black box.” However, a classic exemplar to study for those being introduced to the field of aeroelasticity is a cantilever wing encastré at its root with a single degree of freedom in pitch (torsional motion), as shown in the figure below. Wings are well-damped aerodynamically in bending, so the torsional behavior is perhaps more instructive. The maximum airspeed of many aircraft is limited because of the onset of torsional flutter on the wings, which is a design compromise between wing strength, torsional rigidity, weight, and cost.

A cantilever wing encastré at its root with a single degree of freedom in pitch (torsional motion).

In this case $q$ = $\theta$ , so the equation of motion is

(100) $\begin{equation*} I \, \overbigddot{\theta} + C \, \overbigdot{\theta} + K \, \theta = M_a \end{equation*}$

where $I$ is the mass moment of inertia per unit span about the elastic axis, $C$ is the structural damping coefficient (usually not much), $K$ is the torsional stiffness (mostly dictated by the wing box structure), and $M_a$ is the aerodynamic moment. Assuming quasi-steady aerodynamics,^[14] the aerodynamic moment per unit span can be written as

(101) $\begin{equation*} M_a = \frac{1}{2} \varrho \, V^2 c^2 \, C_{M_\alpha} \theta \end{equation*}$

where $C_{M_\alpha}$ is the slope of the aerodynamic moment curve about the elastic axis. Notice that it can be assumed that any change in aerodynamic angle of attack, $\alpha$ , is proportional to the change in $\theta$ . Therefore, the equation of motion becomes

(102) $\begin{equation*} I \, \overbigddot{\theta} + C \, \overbigdot{\theta} + K \, \theta = \frac{1}{2} \varrho \, V_{\infty}^2 c^2 \, C_{M_\alpha} \theta \end{equation*}$

which is the equation to be solved for $\theta(t)$ .

Time-Marching Solution

There are different ways to find the solution to the torsional flutter problem. One approach is a numerical integration method, which can integrate the equations of motion by marching in time. First, a time step $\Delta t$ must be chosen that is small enough to capture the response accurately. Second, set initial values for $\theta(0)$ and $\overbigdot{\theta}(0)$ , which represent a disturbance to excite the aeroelastic system. Physically, such a disturbance could be a gust or a sudden flight control input. During flight testing, mechanical shakers excite the structure at different frequencies to examine possible aeroelastic problems or flutter.

Then, a numerical integration method, such as the Euler method (which is only first-order accurate) or a Runge-Kutta method (much more accurate but more computationally expensive), can be used to update the solution at each step, i.e.,

(103) $\begin{equation*} \overbigdot{\theta}_{n+1} = \overbigdot{\theta}_n + \overbigddot{\theta}_n \, \Delta t \quad \mbox{and} \quad \theta_{n+1} = \theta_n + \overbigdot{\theta}_n \, \Delta t \end{equation*}$

The values of $\overbigddot{\theta}_n$ can be computed based on the current values of $\theta_n$ and $\overbigdot{\theta}_n$ . The aerodynamic moment on the wing, $M_a$ , at each time step can be computed based on the current values of $\theta_n$ and $V_n$ ; usually $V_n$ = $V_{\infty}$ = constant. Continue marching until the desired end time or a specified condition, such as flutter onset. Suppose nothing happens after a long time period. In that case, the process can be repeated for other airspeeds or for variations in the different parameters, such as the center of gravity location, which affects the torsional inertial about the elastic axis (e.a.).

The figure below shows typical results that can be obtained from the time-marching process. Depending on the wing design, its aeroelastic response can be well-damped, showing a limit cycle oscillation, a classic form of diverging flutter, or rapid divergence. The process would be repeated for different airspeeds and values of air density (density altitude) to estimate the stability boundaries for the structure or aircraft.

Depending on the wing design, the aeroelastic response of the wing can be damped and show a limit cycle oscillation, flutter, or divergence.

Analytical Approach

A direct approach to finding the flutter speed can be used, which needs more mathematics. Assume a trial solution to the equation of motion of the form

(104) $\begin{equation*} \theta(t) = \theta_0 \, e^{i \omega t} \end{equation*}$

Differentiating gives

(105) $\begin{equation*} \overbigddot{\theta}(t) = -\omega^2 \, \theta_0 \, e^{i \omega t} \quad \mbox{and} \quad \overbigdot{\theta}(t) = i \, \omega \, \theta_0 \, e^{i \omega t} \end{equation*}$

Substituting into the equation of motion gives

(106) $\begin{equation*} I (-\omega^2 \, \theta_0 \, e^{i \omega t}) + C (i \omega \, \theta_0 \, e^{i \omega t}) + K (\theta_0 \, e^{i \omega t}) = \frac{1}{2} \varrho \, V^2 \, c^2 \, C_{M_\alpha} (\theta_0 \, e^{i \omega t}) \end{equation*}$

where $K$ is the torsional stiffness, i.e., $K = G \, J$ , where $G$ is the shear modulus and $J$ is the second moment of area of the section. Simplifying gives

(107) $\begin{equation*} -I \, \omega^2 + i \, C \omega + K = \frac{1}{2} \varrho \, V_{\infty}^2 \, c^2 \, C_{M_\alpha} \end{equation*}$

This latter equation can now be rearranged to form the characteristic equation, i.e.,

(108) $\begin{equation*} -I \, \omega^2 + i \, C \omega + K - \frac{1}{2} \varrho \, V_{\infty}^2 \, c^2 \, C_{M_\alpha} = 0 \end{equation*}$

The real and imaginary parts of the equation can be separated to determine if flutter will occur. For the real part, then

(109) $\begin{equation*} -I \omega^2 + K = \frac{1}{2} \varrho \, V_{\infty}^2 \, c^2 \, C_{M_\alpha} \end{equation*}$

and for the imaginary part, then

(110) $\begin{equation*} i \, C \, \omega = 0 \end{equation*}$

If it assumed that there is no structural structural damping, i.e., $C = 0$ , then

(111) $\begin{equation*} -I \, \omega^2 + K = \frac{1}{2} \varrho \, V_f^2 \, c^2 \, C_{M_\alpha} \end{equation*}$

where $V_f$ is the critical airspeed at which flutter will occur. Setting $\omega = 0$ allows the airspeed at which the system becomes unstable to be found, i.e.,

(112) $\begin{equation*} K = \frac{1}{2} \varrho \, V_f^2 \, c^2 C_{M_\alpha} \end{equation*}$

Rearranging gives the critical flutter speed as

(113) $\begin{equation*} V_f = \sqrt{\frac{2 K}{\varrho \, c^2 \, C_{M_\alpha}}} \end{equation*}$

This exemplar illustrates how aerodynamic forces interact with the structural dynamics of the wing, influencing its stability and behavior under various flight conditions. Aeroelastic analysis of such systems is crucial for designing aircraft wings that can withstand aerodynamic loads while maintaining stability and structural integrity. Notice from Eq. 113 that, in this case, the flutter speed is proportional to the square root of the ratio of the torsional stiffness to the aerodynamic effects. Making the wing stiffer will increase the flutter speed, but also increase structural weight and cost.

Flutter Occurences in Practice

In practice, aircraft structures can exhibit coupled responses between different “modes,” e.g., between wing torsion and bending, between the engine mount and the wing, or between the wing and a control surface. Two fatal crashes of the Lockheed Electra were attributed to a phenomenon called “pylon whirl flutter,” where the aerodynamic excitations from the propeller coupled with the natural frequency of the wing, which entered into a destructive resonant vibration and catastrophic flutter. Extensive structural modifications were made to the Electra’s design, including changes to the wings and the engine mounts to increase their stiffness. During its first flight tests, the Boeing Dash 80 (predecessor to the Boeing 707) experienced aileron flutter, after which modifications were made to the control system. This incident underscored the importance of flutter testing for jet airliners and led to improvements in the control surface design and testing protocols. The Boeing 747 encountered unexpected flutter on its horizontal stabilizers during its flight test program. The design was modified by adding additional structural reinforcements to ensure flutter-free operation.[/footnote]

The interaction between all structural modes must be represented on an aircraft to account for the overall coupled dynamic behavior of the structure. This process involves solving the coupled equations describing the entire structure and the aerodynamics, and finding the eigenvalues that represent the system’s natural frequencies and damping ratios, as outlined previously. The figure below shows an example of the information that can be extracted from this type of aeroelastic analysis, which shows the frequency and damping as a function of airspeed, or dynamic pressure should be used more appropriately. The likelihood of the onset of flutter increases if the structural frequencies begin to coalesce.

The frequency and damping of two structural modes as a function of airspeed, showing that flutter is likely to occur near frequency coalescence. — Frequency and damping of two structural modes as a function of airspeed show that the likelihood of flutter increases if the structural frequencies begin to coalesce.

As airspeed increases or the balance of aerodynamic to elastic forces changes, the likelihood of flutter increases when the damped frequency of structural modes begins to merge. The undamped zero airspeed frequencies are distinct at low speeds, and the system remains stable. As the airspeed increases, the aerodynamic forces acting on the structure also increase. These forces change the frequencies of the structure because they add an effective aerodynamic stiffness and damping to the system, which in some cases may cause the frequencies to change so that they start converging, a phenomenon known as frequency coalescence. This condition results in a continuous energy exchange between the aerodynamic forces and the elastic energy stored in structural vibrations and overcomes structural damping.

Flutter will likely arise if the net structural-aerodynamic damping of the aeroelastic system becomes negative. In the flutter condition, aerodynamic forces feed energy into the structural mode at a rate that exceeds the damping at a frequency that matches the structure’s natural frequency. If the energy input from the aerodynamic forces exceeds the inherent damping in the structure, the oscillations grow exponentially, leading to flutter. Flutter remains a critical consideration in aircraft design and testing. Historical incidents have highlighted the devastating potential of flutter, leading to significant advancements in our understanding and ability to prevent it.

Summary & Closure

This chapter has summarized various advanced computational methods for analyzing flight vehicles, emphasizing the critical role of computational fluid dynamics (CFD), mainly through the use of 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 essential for better understanding practical flight problems.

The extensive computational resources and algorithmic challenges posed by modern CFD methods, especially Navier-Stokes and Euler-based approaches, have limited their practical application in aircraft design. Modeling the flow around a complete aircraft is daunting, and validation of these models requires extensive work. This poses challenges for experimentalists in providing high-quality surface and off-surface flow field measurements. Proper validation is crucial for assigning the confidence levels necessary for using CFD in design.

The Finite Element Method (FEM) is employed to analyze structural integrity, using structured or unstructured grids to simulate stress distribution, deformation, and potential failure in mechanical components. Fine 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, optimizing designs, and ensuring adequate strength, durability, and margins of safety.

In aeroelastic analysis, engineers study phenomena such as flutter, divergence, and other dynamic instabilities that affect airplane structures. Mitigating these effects is crucial for designing wings, control surfaces, and other components that can withstand aerodynamic forces while maintaining stability and performance throughout the flight envelope. Aeroelasticity is vital in ensuring modern airplanes remain “flutter-free” and meet stringent safety and performance standards under all flight conditions.

5-Question Self-Assessment Quickquiz

For Further Thought or Discussion

To some, developing advanced CFD methods is considered the apotheosis of aerodynamicists and designers. Discuss why this perspective is misleading and hinders shorter-term research using other viable approaches.
Starting from the Euler equations in two-dimensional steady flow show that for an irrotational isentropic flow, the governing equation for the flow can be written as
$\left( \frac{\partial \phi}{\partial x} \right) + \left( \frac{ \partial \phi}{\partial y} \right) + \left( \frac{2}{\gamma - 1} \right) a^2 = \mbox{constant},$

where $\phi$ is a velocity potential. Explain why this equation is still challenging to solve.
If the flow is assumed to be incompressible, the vorticity transport equation can be derived directly by taking the curl of the Navier-Stokes momentum equation.
Laplace’s equation is the underlying governing equation for surface singularity (or panel) methods. Show that a solution to Laplace’s equation allows velocity potentials and velocity components to be added but not pressures.
Explain the types of computational grids that could be used for finite-difference solutions to the Navier-Stokes equations. What grids and resolutions would be preferable near the rotor instead of in the far wake? Justify your answer.
What is meant by a “chimera” grid?
Use a web search to find out about currently used grid generation methods for CFD applications to helicopter aerodynamics problems.
Explain the differences between a structured grid and an unstructured grid. Why might using unstructured grids be viewed as a beneficial approach in some types of problems found in helicopter aerodynamics?

Other Useful Online Resources

Visit the following websites to learn more about advanced computational methods for the aerospace field:

“Introduction to CFD Basics” by Drs. Bhaskararn and Collins from Cornell University
“What is CFD” – Introduction to Computational Fluid Dynamics.
Computational Fluid Dynamics (CFD) – History and Beginner’s Guide.
Getting started with CFD.
Finite Element Analysis Explained – Things you must know about FEM.
Understanding the Finite Element Method.
What is Finite Element Analysis? FEM explained for beginners.
Introduction to Aeroelasticity in NASTRAN.
Aeroelasticity – Introduction to Flutter.
Aerospace Structures – Aeroelasticity.

Moore's Law, observed by Gordon E. Moore in 1965, states that the number of transistors on a microchip doubles approximately every two years, leading to exponential increases in computing power and cost reductions. ↵
Time = money. In this regard, wall clock time = money plus compute time = more money. ↵
Published September 1933 ↵
Of "Taylor series" fame. ↵
"Numerical Solution of Convection-Diffusion Problems" by D. Brian Spalding (1972) ↵
"Numerical Heat Transfer and Fluid Flow" by Suhas V. Patankar (1980) ↵
Stokes assumed that $\lambda + (2/3) \mu = 0$ , although later research suggests it is not always accurate. ↵
By definition, this is where the local velocity approaches 99% of the external flow (or edge) velocity $U_e(x, t)$ . ↵
Jameson, A., "Transonic Potential Flow Calculation using Surface Transpiration Boundary Conditions," AIAA Journal, 13(7), 1975, pp. 826–833. ↵
In that case, There is no reason to use any "higher-order" distribution of vorticity as the numerical cost of the solution (order $N^3$ ) exceeds that of an increase in the number of panels (order $N^4$ ). ↵
"Modeling the Flow of Rarefied Gases at NASA," Forrest E. Lumpkin, NASA Johnson Space Center, 9/21/2012. ↵
A dimensionless constant typically between 0 and 0.5 for most engineering materials. ↵
Theodore Theodorsen, "General Theory of Aerodynamic Instability and the Mechanism of Flutter," NACA-TR-496 (1949). ↵
Meaning the the aerodynamics adjust instantaneously to changes in pitch angle with no time lag, even though this is physically incorrect for higher reduced frequencies. ↵

License

Icon for the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License

Introduction to Aerospace Flight Vehicles Copyright © 2022–2024 by J. Gordon Leishman is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, except where otherwise noted.