Multifidelity Comparison of Supersonic Wave Drag Prediction Methods Using Axisymmetric Bodies ^†

Abstract

Low-fidelity analytic and computational wave drag prediction methods assume linear aerodynamics and small perturbations to the flow. Hence, these methods are typically accurate for only very slender geometries. The present work assesses the accuracy of these methods relative to high-fidelity Euler, compressible computational-fluid-dynamics solutions for a set of axisymmetric geometries with varying radius-to-length ratios $(R/L)$. Grid-resolution studies are included for all computational results to ensure grid-resolved results. Results show that the low-fidelity analytic and computational methods match the Euler CFD predictions to around a single drag count ( ∼$1.0\times {10}^{-4}$) for geometries with $R/L\le 0.05$ and Mach numbers from 1.1 to 2.0. The difference in predicted wave drag rapidly increases, to over 30 drag counts in some cases, for geometries approaching $R/L\approx 0.1$, indicating that the slender-body assumption of linear supersonic theory is violated for larger radius-to-length ratios. All three methods considered predict that the wave drag coefficient is nearly independent of Mach number for the geometries included in this study. Results of the study can be used to validate other numerical models and estimate the error in low-fidelity analytic and computational methods for predicting wave drag of axisymmetric geometries, depending on radius-to-length ratios.

1. Introduction

Low-fidelity aerodynamic analysis utilizing linear supersonic theory offers fast and reasonably accurate methods for preliminary supersonic aircraft design. One such low-fidelity method is the widely known Harris wave drag code [1]. Modern versions of this code are readily available through open-source software, such as the NASA OpenVSP tool [2]. Note that OpenVSP Version 3.21.1 was used in this work. Previous publications include validation and comparison studies using OpenVSP and other similar wave drag prediction methods [3,4], but the pool of available cases is relatively small. The current work considers the accuracy of these linear wave drag prediction methods for simple geometries over a range of low supersonic Mach numbers. The results presented here offer a broader source of reference for the expected accuracy of these tools.

To characterize this accuracy, results are generated using three separate inviscid wave drag prediction methods, each representing a different level of fidelity. These methods include analytic solutions, the linear-theory-based OpenVSP wave drag tool, and Euler, compressible computational fluid dynamics (CFD) results. Simple axisymmetric bodies, defined by a radius-to-length ratio $(R/L)$, were chosen for this study because of their general relevance to preliminary supersonic aircraft design and the availability of analytic solutions relating their wave drag to geometric parameters. The current work includes grid convergence studies for each of these methods, including OpenVSP solution convergence with slicing parameters, and comparisons of converged results over the Mach number and $R/L$ design space.

Previous experimental and computational results offer some insight into trends presented in the current work. Results presented by Harris and Landrum [5] and Stivers [6] show that the wave drag coefficient of axisymmetric geometries is dependent on Mach number and the slenderness of the geometry. This Mach number dependence decreases as the slenderness of a geometry increases [7]. As will be noted in the Methods and Materials section that follows, the wave drag relationships derived from linear supersonic theory are largely independent of Mach number. This Mach number independence is also noted in works by Rallabhandi and Mavris [3] and Cheung [8]. For slender axisymmetric geometries, variations to equivalent area distributions due to changes in Mach number are relatively small. This results in wave drag results that are nearly independent of Mach number. These trends are seen in the results produced by all three methods presented in the current work.

The results and trends presented here offer a resource that can be used to understand the expected error in the wave drag predicted using these low-fidelity methods. Depending on the magnitude of this error, these low-fidelity methods may be deemed suitable or unsuitable for a specific stage of design or a specific aircraft configuration. As will be seen in the Results and Discussion section of this work, this accuracy appears to largely be dependent on the slenderness of the geometry, which conforms to the slender-body assumption used in the development of the linearized methods.

2. Methods and Materials

For completeness and convenient reference, this work begins with an overview of the methods and assumptions used in the development of each of the wave drag prediction methods outlined within. Following those developments, grid convergence studies pertaining to each of those methods are presented.

2.1. Development of the Linear Supersonic Wave Drag Equation

An estimate for the wave drag of a supersonic body can be found using an inviscid analysis of the momentum loss through a cylindrical control surface surrounding the body. The cylindrical control surface should have a radius, $R_1$, of several times the length of the supersonic body and the aft trailing disk should be located far enough downstream that the trailing flow leaving that surface can be considered two-dimensional. The disturbances traveling along Mach lines emanating from the body should cross the cylindrical walls upstream of the aft trailing disk. This control surface is shown in Figure 1 with respect to a generic wing-body configuration.

Ashley and Landahl [9] present the full mathematical development of the momentum analysis on this control surface. An abridged version of their work will be given here. The relationship resulting from this control surface analysis for the total pressure drag on the supersonic body is found to be

$$D=-{\rho }_{\infty }{U}_{\infty }^2\underset{S_2}{\int \int }{\varphi }_x{\varphi }_{r}d{S}_2+\frac{1}{2}{\rho }_{\infty }{U}_{\infty }^2\underset{S_3}{\int \int }({\varphi }_y^2+{\varphi }_z^2)d{S}_3$$

(1)

where the values of $\varphi$ are the supersonic velocity potential perturbations at a point in the flow. In Equation (1), the second integral represents the vortex, or induced, drag and the first integral represents the wave drag exclusive to supersonic flow. The solution for the vortex drag is the same as what is found for subsonic flow [9]. The focus of the current work will be restricted to estimates for the wave drag. From Equation (1), this gives

$$D_w=-{\rho }_{\infty }{U}_{\infty }^2\underset{S_2}{\int \int }{\varphi }_x{\varphi }_{r}d{S}_2$$

(2)

A solution to Equation (2) is found using a lineal distribution of sources to represent an axisymmetric body in the flow [10]. This lineal distribution is a continuous line of sources with strengths varying in the axial direction. The perturbation velocity potential of this lineal source distribution is [9]

$$\varphi =-\frac{1}{2\pi }\underset{0}{\overset{x-Br}{\int }}\frac{f\left(x_1\right)d{x}_1}{\sqrt{{(x-x_1)}^2-B^2{r}^2}}$$

(3)

where $f\left(x_1\right)$ represents the strength of the lineal source distribution per unit axial length. The variable r in Equation (3) represents a radial distance from the axis of symmetry and is defined as $r=\sqrt{y^2+z^2}$. The value $\varphi$ is the velocity potential perturbation at a point as a result of the integration of the influence of the full lineal source distribution along its length. Ashley and Landahl [9] show the process of using Equation (3) to generate the derivatives for Equation (2). Substitution of those derivatives into Equation (2) produces the relationship

$$D_w=-\frac{{\rho }_{\infty }{U}_{\infty }^2}{4\pi }\underset{0}{\overset{l}{\int }}\underset{0}{\overset{l}{\int }}{f}^{\prime }\left(x_1\right)f^{\prime }\left(x_2\right)\ln \left|x_1-x_2\right|d{x}_{1}d{x}_2$$

(4)

which relates wave drag to the source strength of the lineal source distribution. The source strength, $f\left(x\right)$, is related to the rate of the area change of the axisymmetric body in the direction of the freestream [9], represented mathematically as $f\left(x\right)=S^{\prime }\left(x\right)$. This relationship assumes that the axisymmetric body is a slender body. A slender body is defined as a body with a span-wise dimension (radius or diameter) that is small compared to its length [9]. Additionally, disturbances to the flow generated by this body are assumed to be small [10,11]. Substituting this relationship for the source strength into Equation (4) yields

$$D_w=-\frac{{\rho }_{\infty }{U}_{\infty }^2}{4\pi }\underset{0}{\overset{l}{\int }}\underset{0}{\overset{l}{\int }}{S}^{″}\left(x_1\right)S^{″}\left(x_2\right)\ln \left|x_1-x_2\right|d{x}_{1}d{x}_2$$

(5)

where $S^{″}\left(x\right)$ is the second derivative of the area of the supersonic body with respect to axial location along the body.

Equation (5) is known as the linear supersonic wave drag integral, also referred to as the von Karman wave drag integral [12]. Its derivation can also be found in a number of additional sources, including work by Ward [13] and Lighthill [14]. Under the assumptions of linear supersonic theory, it represents the wave drag of a slender supersonic body and demonstrates that wave drag has a strong dependence on the area distribution of that body. The current work will examine a number of ways to produce solutions to Equation (5) for an estimate of the wave drag of an axisymmetric body. Additionally, the current work will present results that help further define what geometries could be considered slender.

2.2. Analytic Solutions of the Wave Drag Integral

Solving Equation (5) directly is challenging for two reasons. The first is due to the singularity encountered in the logarithmic term within the integral. The second is the need for a continuous representation of the area distribution, $S\left(x\right)$, of the supersonic body. A number of methods have been developed to simplify the integration, or to circumvent the issue of direct integration altogether. One such analytic approach utilizes a Fourier sine series to simplify the integration. Ashley and Landahl [9] provide the development of this method. Their work will be summarized here.

A change of variables for the axial location along a supersonic body of the form

$$x=\frac{l}{2}(1+\cos \theta )$$

(6)

allows for the area distributions in Equation (5) to be represented as an infinite Fourier sine series, defined as

$$S^{\prime }\left(x\right)=l\sum _{n=1}^{\infty }{A}_n\sin n\theta$$

(7)

Here, the $A_n$ values are the Fourier coefficients, l is the length of the body, and $\theta$ is referred to as the Glauert variable [9]. Note that the Fourier sine series in Equation (7) represents the first derivative of the area distribution, not the area distribution itself. Integration of Equation (7) yields

$$S\left(\theta \right)=\frac{l^2}{4}\left\{A_1\left(\pi -\theta +\frac{\sin 2\theta }{2}\right)+\sum _{n=2}^{\infty }{A}_n\left[\frac{\sin (n+1)\theta }{n+1}-\frac{\sin (n-1)\theta }{n-1}\right]\right\}$$

(8)

which represents the area distribution of an axisymmetric body as a function of the variable $\theta$, as defined in Equation (6).

The volume of a body represented by the Fourier sine series can be found by integrating Equation (8), which gives

$$V=\frac{\pi l^3}{8}\left(A_1-\frac{1}{2}{A}_2\right)$$

(9)

This result shows that the volume of the body is related to only the first and second Fourier coefficients, $A_1$ and $A_2$, and is independent of all other Fourier coefficients with $n>2$.

Using the Fourier sine series representation of the area distribution given in Equation (8), the double integration in Equation (5) can be solved to yield

$$D_w=\frac{\pi {\rho }_{\infty }{U}_{\infty }^2{l}^2}{8}\sum _{n=1}^{\infty }n{A}_n^2$$

(10)

which is the wave drag of a body as a function of the Fourier coefficients used in Equation (7). Nondimensionalizing Equation (10) with the freestream dynamic pressure, $\frac{1}{2}{\rho }_{\infty }{U}_{\infty }^2$, and an area based on the length of the body, $l^2$, yields a result for the wave drag coefficient as a function of the Fourier coefficients given by

$$C_{Dw}=\frac{\pi }{4}\sum _{n=1}^{\infty }n{A}_n^2$$

(11)

Equations (10) and (11) are simple yet powerful tools for studying the wave drag of axisymmetric supersonic bodies.

A minimum value of Equations (10) and (11) can be achieved when all coefficients unnecessary for defining a geometry are set to zero. To have a body with nonzero volume, either $A_1$ or $A_2$ must be defined, as seen in Equation (9). An area distribution with only $A_1$ defined will result in a nonzero base area, similar to a bullet. This geometry is often referred to as a von Kármán ogive [9]. An area distribution with only $A_2$ defined results in a start and end area of zero, visually represented as a geometry with a pointed nose and tail. This geometry is called a Sears–Haack body, which was developed independently by both Sears [15] and Haack [16]. This work will focus on the Sears–Haack geometry alone.

If the volume of the Sears–Haack body and its length are known values, Equation (9) can be solved for the coefficient $A_2$. This results in

$$A_2=-\frac{16V}{\pi l^3}$$

(12)

Using Equation (12) in Equation (8) yields

$$S\left(\theta \right)=\frac{4V}{\pi l}\left(\sin \theta -\frac{1}{3}\sin 3\theta \right)$$

(13)

which is the area distribution for the Sears–Haack body. Alternatively, the radius distribution of a Sears–Haack body can be represented as a function of the original axial location [4], or

$$r_s\left(x\right)=R{\left[4\frac{x}{l}\left(1-\frac{x}{l}\right)\right]}^{\frac{3}{4}}$$

(14)

which requires a maximum radius, R, to be defined. The value of R is related to the volume of the geometry via the relationship

$$V=\frac{3{\pi }^2}{16}{R}^{2}l$$

(15)

Using Equation (12) in Equation (10) yields a prediction for the wave drag of a Sears–Haack body based on linear theory:

$$D_w=\frac{64V^2}{\pi l^4}{\rho }_{\infty }{U}_{\infty }^2$$

(16)

These Sears–Haack equations provide a convenient point of comparison for designing a supersonic aircraft with a fixed length and volume. The Sears–Haack geometry and the related analytic relationships for wave drag will be used as the main point of comparison between the other methods studied in the current work.

2.3. Analytic Solutions at Higher Mach Numbers

The previous development of the analytic relationship for the wave drag of an axisymmetric body, Equation (5), did not result in an explicit dependence on Mach number. Instead, the Mach number influence comes from the area distribution itself, or, more specifically, $S^{″}\left(x\right)$. To utilize the previous analytic relationships at higher Mach numbers, an understanding of the equivalent body concept is necessary. At supersonic speeds, disturbances due to a body moving through the flow travel away the from the body along lines inclined at the Mach angle, commonly referred to as Mach lines. The Mach angle is related to the freestream Mach number through the relationship $\mu ={\sin }^{-1}\frac{1}{M}$. The flow properties related to the momentum loss, and thus wave drag, are propagated along these Mach lines. This indicates that it is not the area distribution normal to the longitudinal axis, but, rather, the area distribution inclined at the Mach angle that will dictate the source strength used in the linear wave drag equations [17].

The method of generating an equivalent body of revolution using planes inclined at the Mach angle, which will be referred to as Mach inclined planes, is outlined in a number of publications found in the linear supersonic literature. The work of Maglieri and Carlson [18] describes how the area related to the volume of an aircraft and the area related to the lift are considered for sonic boom analysis. The work by Jones [12] and Harris [1] each individually provide descriptions, including figures, of how the equivalent area needs to be considered when evaluating the wave drag of a supersonic body. The current work will summarize the description by Jones [12] with updated figures.

An equivalent body of revolution is generated using a series of parallel Mach planes to slice a supersonic body at locations along its longitudinal axis. The intersection of these parallel Mach planes with the supersonic body produces a cross-section of area that is dependent on the local geometry and the angle describing the Mach planes. The normal projection of these Mach-plane sliced-area cross-sections is then used to represent the normal, or noninclined, area cross-section of the equivalent body of revolution. A representation of the relationship between the original geometry and its equivalent body of revolution is shown in Figure 2. This process is repeated at incremental longitudinal locations until the area distribution of the equivalent body of revolution for the entire supersonic body has been generated. If the supersonic body is generating lift, an additional component of area that is dependent on the distribution of lift must be included. For the current work, only nonlifting bodies will be considered.

When the supersonic body is a nonaxisymmetric configuration, such as a wing-body combination, an equivalent body of revolution must be generated for each azimuthal angle around the body. An example of the Mach inclined planes at the same axial position, but at different azimuthal angles, is shown in Figure 3.

Slicing at these different azimuthal angles results in a different longitudinal distribution of area, as shown in Figure 4.

When estimating the wave drag of a nonaxisymmetric geometry, the wave drag of each equivalent body for different azimuthal angles is used to generate an integrated average over the azimuthal angle range. This integrated average can be applied to Equation (5), which becomes

$$D_w=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}-\frac{{\rho }_{\infty }{U}_{\infty }^2}{4\pi }\underset{0}{\overset{l}{\int }}\underset{0}{\overset{l}{\int }}{S}^{″}(x_1,\psi )S^{″}(x_2,\psi )\ln \left|x_1-x_2\right|d{x}_{1}d{x}_{2}d\psi$$

(17)

Note that the variable $\psi$, and not $\theta$, as is common in the literature, is used to represent the azimuthal angle in the current work. While this work will only study axisymmetric geometries, Equation (17) is used to capture the effects of azimuthal grid discretization on the predicted wave drag values.

2.4. A Numerical Technique for Solving the Wave Drag Integral and Its Modern Implementation in OpenVSP

The previous sections showed an analytic approach that can be used to represent arbitrary axisymmetric geometries, their area distributions, and a relationship for their wave drag. When the geometry is known a priori, a method capable of adequately representing more complex aircraft configurations, while also avoiding the introduction of numerical complications, is necessary.

A method meeting this criteria was developed by Eminton and Lord in the 1950s [19] and remains widely used in modern linear supersonic theory-based wave drag analysis [1,3,4]. Tools such as the Harris wave drag code [1] and the more modern OpenVSP wave drag tool [4] are two of its better-known implementations. The method was designed to produce a solution to Equation (5) for complex aircraft configurations where the discrete Mach plane sliced-area distributions are known. Generating these discrete area distributions is generally a simple procedure that can be accomplished using modern geometry modeling and analysis tools. The Eminton and Lord method yields a continuous area distribution that both matches the specified discrete area points and is a minimizing solution to Equation (5) [20].

The current work will include a brief overview of the mathematics behind the Eminton and Lord method. Further details can be found in the original work by Eminton and Lord [19,20] and its implementation in OpenVSP by Waddington [4].

Some assumptions about the discrete area distributions are required when implementing the Eminton and Lord method. First, $S^{\prime }\left(x\right)$ is assumed to be a continuous function of x over the length of the aircraft. Second, it is assumed that the first derivative of the area at the start and end of the distribution is zero, or $S^{\prime }\left(0\right)=S^{\prime }\left(l\right)=0$. Under these assumptions, an approach using a Fourier sine series, similar to what was conducted for the analytic method, can be used to produce a solution to Equation (5) for more complex aircraft configurations. For area distributions matching these assumptions, a change of variable from x to $\theta$ of the form

$$\theta ={\cos }^{-1}\left(1-2x\right)$$

(18)

allows the first derivative of the area distribution to be expressed as an infinite Fourier sine series, written as

$$S^{\prime }\left(x\right)=\sum _{r=1}^{\infty }{a}_r\sin r\theta$$

(19)

The values for x represent a normalized distribution of locations along the aircraft with values between 0.0 and 1.0, corresponding to values of $\theta$ between $\pi$ and $0.0$, respectively. Eminton and Lord showed that to fit this Fourier series through n discrete points, the Fourier coefficients for $r\ge 2$ can be determined using

$$a_r=\frac{1}{r}\sum _{i=1}^n{\lambda }_i\left(\frac{\sin \left(r-1\right){\theta }_i}{r-1}-\frac{\sin \left(r+1\right){\theta }_i}{r+1}\right)$$

(20)

where the values of ${\lambda }_i$ are solved using the n discrete area values in a system of n linear equations [19,20]. The value of $a_0$ is set by the nose area, or

$$a_0=S_N$$

(21)

and $a_1$ is set by the difference in nose and base area, or

$$a_1=\frac{4}{\pi }\left(S_B-S_N\right)$$

(22)

where the nose and base areas are the first and last values of the discrete area distribution, respectively. Direct integration of Equation (19) yields

$$S\left(x\right)=a_0+\frac{1}{4}{a}_1\left(\theta -\sin \theta \cos \theta \right)+\sum _{r=2}^{\infty }\frac{1}{4}\left(\frac{\sin \left(r-1\right)\theta }{r-1}-\frac{\sin \left(r+1\right)\theta }{r+1}\right)a_r$$

(23)

which is an expression for a continuous area distribution as a function of the Fourier coefficients, $a_r$. From Equation (23), the derivatives necessary for evaluating Equation (5) can be obtained. The complete algorithm is straightforward to implement and runs in a fraction of a second. The more complicated component of the algorithm is producing the necessary discrete Mach inclined area distributions of an aircraft. At the time of publishing the current work, the NASA OpenVSP software [2] has the most modern and convenient implementation of this wave drag prediction method [4]. OpenVSP utilizes the Eminton and Lord method to solve for the wave drag of a complete aircraft configuration via Equation (17). The OpenVSP implementation of this tool is the foundation for the linear wave drag results presented in this work. OpenVSP can be operated as a standalone executable, including a user-friendly GUI, or alternatively run through its Python API for easy automation and design space studies.

The wave drag estimate produced by the Eminton and Lord [19] method, and, thus, the OpenVSP linear wave drag tool, is of the form

$$\frac{D_w}{q}=I=-\frac{1}{2\pi }\underset{0}{\overset{l}{\int }}\underset{0}{\overset{l}{\int }}{S}^{″}\left(x_1\right)S^{″}\left(x_2\right)\ln \left|x_1-x_2\right|d{x}_{1}d{x}_2$$

(24)

The x variables within the integrals can be nondimensionalized using the length of the body, yielding

$$\frac{D_w}{q}=-\frac{1}{l^2}\frac{1}{2\pi }\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1}{\int }}{S}^{″}\left(\frac{x_1}{l}\right)S^{″}\left(\frac{x_2}{l}\right)\ln \left|\frac{x_1}{l}-\frac{x_2}{l}\right|d\left(\frac{x_1}{l}\right)d\left(\frac{x_2}{l}\right)$$

(25)

The nondimensionalization of the wave drag can be completed by dividing both sides of Equation (25) by a reference area. To be consistent with the analytic relationship in Equation (11) above, a reference area defined by the length of the body squared, $l^2$, is used. The result is

$$C_{Dw}=-\frac{1}{l^4}\frac{1}{2\pi }\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1}{\int }}{S}^{″}\left(\frac{x_1}{l}\right)S^{″}\left(\frac{x_2}{l}\right)\ln \left|\frac{x_1}{l}-\frac{x_2}{l}\right|d\left(\frac{x_1}{l}\right)d\left(\frac{x_2}{l}\right)$$

(26)

Equations (11) and (26) are the relationships for $C_{Dw}$ that are used through the remainder of this work.

2.5. Euler CFD Solutions

The CFD results presented in the current work were generated using the NASA Cart3D tool suite [21], which in recent years has been commonly used in preliminary supersonic aircraft design and sonic boom analysis [22,23]. Cart3D solves the Euler equations with a nonlinear, finite-volume, second order, MPI parallel flow solver operating on a nonconformal Cartesian volume mesh utilizing cut-cell technology. As an inviscid CFD Euler solver, the only inputs related to flow properties are the freestream Mach number, the aerodynamic angles, and the ratio of specific heats, $\gamma$. In this work, the Mach number is varied according to each study, the aerodynamic angles are kept as zero, and $\gamma =1.4$.

When combined with its fully-integrated adjoint solver, Cart3D is capable of output-based mesh adaptation to minimize discretization error in the volume mesh, and it is also capable of conducting gradient-based design shape-optimization with user-defined parameterizations of the models undergoing design. In this study, the van Leer limiter option was used to provide smooth and aggressive limiting with low dissipation. In order to capture both wave drag estimates and off-body equivalent area information, the output-based volume mesh adaptation was implemented using a pressure functional along a line sensor a small distance (approximately 3.2 times the model radius) below the axisymmetric model. A generic functional, appropriate for mesh adaptation or design, was defined as

$$J={\int }_{s_{min}}^{s_{max}}{\left(\frac{p\left(s\right)-p^{*}\left(s\right)}{p_{\infty }}\right)}^{2}ds$$

(27)

with line sensor coordinate $s\in [s_{min},s_{max}]$, the static off-body pressure $p\left(s\right)$ interpolated from the solution, a target static pressure $p^{*}\left(s\right)$, and freestream static pressure $p_{\infty }$. The target was set as $p^{*}\left(s\right)=p_{\infty }$ so that adjoint sensitivities would track disturbances to the freestream flow impacting the line sensor throughout the mesh adaptation process. For smaller models, five cycles of mesh adaptation were conducted over four hours of wall-clock time on 140 CPUs, resulting in a final volume mesh with approximately 56 million hexes. Larger models required up to seven cycles of mesh adaptation in over five hours of wall-clock time on the same computing hardware, resulting in approximately 63 million hexes. Even though the discretization error was minimized with respect to the off-body functional, sufficient mesh resolution was obtained along the model to provide confidence in the inviscid drag estimate.

At supersonic speeds, the components that make up the total inviscid drag are induced drag (drag resulting from the generation of lift) and wave drag (drag resulting from the generation of shock waves). Because the current work studies axisymmetric supersonic bodies at zero angle of attack, the induced drag and the wave drag due to lift are assumed to be zero. In theory, this allows for the inviscid drag predicted by Cart3D to serve as an adequate reference point for the wave drag due to volume predicted by the analytic and linear methods. An example of a Cart3D adapted volume mesh, solution, and convergence history are shown in Figure 5 for a Mach 1.4 case where the model has a 0.02 radius-to-length ratio ($R/L$). The nonzero inviscid drag coefficient result is expected due to wave drag. Five cycles of mesh adaptation are visible in the residual convergence history and the inviscid drag coefficient convergence history.

The scope of the current work originally included a comparison study between the equivalent area distributions predicted from off-body pressure measurements to the area distributions generated with linear theory Mach plane slicing. Theoretically, these area distributions should be the same for axisymmetric bodies that are not generating lift. Unfortunately, that study has been omitted from the current work and is reserved for future studies. Regardless, the full process that was used to generate both the Cart3D surface and volume meshes has been discussed here for completeness.

A grid resolution study for the surface mesh model was conducted using the same Cart3D output-based mesh adaptation approach described above. The following section begins with a summary of that study to show that the inviscid drag values were produced using well-resolved surface meshes, thus minimizing the error of the $C_{Dw}$ results that will be used for comparisons in later studies.

2.6. Generating Sears–Haack Geometries and Mach Plane Sliced-Area Distributions

With an established analytic solution available, the Sears–Haack body is an ideal reference for comparisons between different wave drag prediction methods. The current work will look at converged solutions for the wave drag coefficient, $C_{Dw}$, generated using the analytic, OpenVSP-based linear, and Cart3D CFD approaches outlined in the previous sections. A description of the approach used to generate the geometries and the $C_{Dw}$ values for each of these methods will now be provided. From these results, the converged solutions will be used to generate contours of the $C_{Dw}$ predicted by each method over the design space, as well as the difference in $C_{Dw}$ between each set of results.

The Sears–Haack geometries were generated using a Cart3D compatible unstructured mesh format with triangle elements. The individual Sears–Haack geometries were defined using a length, L, maximum radius, R, number of azimuthal vertex locations, $n_t$, and number of axial vertex locations, $n_a$. The distribution of azimuthal and axial vertex locations on the Sears–Haack geometry and the resulting panel distribution is demonstrated in Figure 6. Note that the image on the right of Figure 6 is displaying only half of the Sears–Haack geometry.

The radius distribution of each geometry was determined using Equation (14) and a specified value for R. The length of all cases studied was fixed at 10 m and the maximum radius had a range between 0.1 and 1.0 m, resulting in $R/L$ values between 0.01 and 0.1. The number of azimuthal vertices was varied between 10 and 210 for the linear method meshes and between 50 and 800 for the Cart3D meshes. The number of axial vertices was varied between 10 and 210 for the linear method meshes and between 100 and 1600 for the Cart3D meshes.

In addition to inputs defining the Sears–Haack geometries themselves, the OpenVSP linear wave drag tool requires values for the Mach number, reference area, number of azimuthal angles, $n_T$, and number of axial slices, $n_A$. An example of the OpenVSP GUI and wave drag tool is shown in Figure 7. A single Mach inclined plane is shown slicing through the Sears–Haack geometry at the undertrack azimuthal angle, $\psi =180.0^{\circ }$. Within the OpenVSP wave drag tool, azimuthal angles are defined as a rotation around the x-axis and axial locations are defined along the x-axis.

Because the geometries in this study are axisymmetric, the Mach plane sliced-area distributions generated at different azimuthal angles are nearly identical. This results in $C_{Dw}$ approximations that are largely independent of the number of azimuthal angles used. The current work used a fixed number of 30 azimuthal angles. The number of axial slices was varied between 10 and 210 and the range of Mach numbers used was between 1.0 and 2.0. The reference area was defined as the length of the body squared, as was conducted in developing Equations (11) and (26).

Table 1. Range of Mach number and geometric values.

M [min, max]	$\mathit{R}/\mathit{L}$ [min, max]	L [m]	${\mathit{A}}_{\mathit{ref}}$ [sq-m]
[1.0, 2.0]	[0.01, 0.1]	10	100

gives a summary of the Mach number range and general geometric values used in the current work.

Table 2. Number of vertices and slices used in the OpenVSP linear wave drag and Cart3D Euler CFD methods.

Method	${\mathit{n}}_{\mathit{a}}$ [min, max]	${\mathit{n}}_{\mathit{t}}$ [min, max]	${\mathit{n}}_{\mathit{A}}$ [min, max]	${\mathit{n}}_{\mathit{T}}$
Linear	[10, 210]	[10, 210]	[10, 210]	30
Euler CFD	[100, 1600]	[50, 800]	-	-

gives the range of the number of vertices and slices used in the respective methods.

While these axisymmetric geometries are much simpler than complete supersonic configurations, their area distributions are useful as minimum wave drag solutions. To place the geometries used in the current study in context with actual aircraft, the approximate fuselage radius-to-length ratios of a variety of supersonic aircraft are shown in Figure 8. These values were approximated using dimensions pulled from different technical and promotional sources [24,25,26,27,28].

The aircraft designed for commercial and long-distance supersonic cruise tend to have radius-to-length ratios toward the lower end of the range in Table 1, while the military and experimental aircraft occupy the middle to high end. Many of these aircraft were designed to operate at Mach numbers higher than the range studied in the current work, but, as the results will show, the accuracy of the different methods relative to each other is far less dependent on Mach number than on the radius-to-length ratio of the geometry. Figure 8 demonstrates that the results of the current work can be a relevant reference to supersonic aircraft designers who want to understand if these lower-order estimates are useful in their work.

2.7. Euler CFD Convergence Study

To obtain an accurate comparison between the Cart3D Euler CFD and linear theory results, grid-resolved CFD solutions are necessary. A study was completed to determine the axial and azimuthal vertex distribution necessary to produce converged Cart3D $C_{Dw}$ results. This study looked at the corner-cases, or combinations of max and min Mach and $R/L$ values, and one interior point of the Mach and $R/L$ range used in this work. To verify convergence, Richardson extrapolation, as outlined in Phillips et al. [29], was used to estimate the order of convergence and a converged value of $C_{Dw}$ for the set of cases ran at each Mach and $R/L$ combination.

In all cases, an apparent order of convergence of at least 2.0 was obtained. Additionally, the error, defined as the absolute value of the difference between the Cart3D value and the Richardson extrapolation predicted value, was calculated and is shown in Figure 9. The error is given relative to the ratio of the maximum triangular edge length to the total length of the Sears–Haack geometries. In every case, the error was reduced to at most $1.0\times {10}^{-7}$, well below the total predicted $C_{Dw}$, when the finest mesh resolution was used. The percent error of the solutions generated using the finest mesh resolution relative to the value predicted using Richardson extrapolation was on the order of $1.0\times {10}^{-5}\%$. Based on these results, the finest mesh resolution, using a distribution of 800 azimuthal vertices and 1600 axial vertices, was deemed adequate for producing the Cart3D $C_{Dw}$ values for the remainder of this work. This fine mesh resolution resulted in a ratio of maximum triangular edge length to total length of between $6.4\times {10}^{-4}$ and $1.3\times {10}^{-3}$ for the $R/L=0.01$ and $R/L=0.1$ cases, respectively.

2.8. Linear Tool Convergence Study

To ensure accurate $C_{Dw}$ values were obtained using the OpenVSP based linear method, a study was conducted to determine the distribution of axial and azimuthal vertices and the number of axial slices necessary to obtain a converged solution. This study looked at the lower bound, midpoint, and upper bound of the Mach number range. These Mach numbers are 1.0, 1.5, and 2.0, respectively. Similarly, the lower bound, upper bound, and three midpoint values of $R/L$ were used to define the geometries. As mentioned previously, the number of axial vertices, azimuthal vertices, and axial slices were each varied from 10 to 210. As each parameter was varied in their respective studies, the number of axial vertices, azimuthal vertices, or axial slices was kept constant at $n=100$. For example, as the effect of changing the number of axial slices was studied, the number of axial and azimuthal vertices was held constant at $n_a=100$ and $n_t=100$, respectively. This was performed to minimize the combined effects of changing multiple parameters at the same time.

The value of $C_{Dw}$ predicted by the OpenVSP linear wave drag tool was compared to the Cart3D Euler CFD solution to provide an error representing the percentage difference relative to the CFD value. The exact definition of this error is

$$Erro{r}_{C_{Dw}}\left[\%\right]=\frac{C_{Dw,Linear}-C_{Dw,Cart3D}}{C_{Dw,Cart3D}}\times 100$$

(28)

Results from this convergence study are shown in Figure 10, Figure 11 and Figure 12 for Mach 1.0, 1.5, and 2.0, respectively. The error in $C_{Dw}$ associated with changes in the number of axial vertices is shown relative to the resulting ratio of axial vertex spacing to total body length. The error in $C_{Dw}$ associated with changes in the number of azimuthal vertices is shown relative to the resulting ratio of azimuthal vertex spacing to maximum body circumference. In most cases, clear convergence is achieved when the axial and azimuthal vertex spacing ratio is at or below $1.0\times {10}^{-2}$, which requires a distribution of at least 100 vertices in those directions. This begins to differ at higher Mach numbers and the largest values of $R/L$, as seen in Figure 11 and Figure 12.

In general, these results show that an axial vertex spacing ratio of around $6.0\times {10}^{-3}$ and an azimuthal vertex spacing ratio of around $1.0\times {10}^{-2}$ provide an adequately converged solution over the range of Mach numbers and $R/L$ values. These vertex spacing ratios were achieved using a number of axial and azimuthal vertices of $n_a=150$ and $n_t=100$, respectively. Additionally, the number of axial slices that provided converged solutions was around $n_A=100$. An exception to this is the solutions with $R/L=0.1$ geometries, which do not appear to converge with increasing axial slices for $M=1.5$ and $M=2.0$. A similar result is seen for the $R/L=0.07$ results at $M=2.0$. As will be seen in the Results and Discussion section that follows, solutions with the largest $R/L$ values, and $R/L=0.1$ in particular, tend to have the largest error relative to the Euler CFD results. Because of this, the convergence of their linear solutions was placed at a lower priority, and a value of $n_A=100$ for the number of axial slices was used for all Mach numbers in the current work.

2.9. Matching Area Distributions at Higher Mach Numbers with a Fourier Sine Series

At Mach numbers greater than 1.0, the area distributions resulting from Mach inclined plane slicing are no longer the same area distribution representing the original Sears–Haack geometry. This result is demonstrated in Figure 13 for the case of $R/L=0.1$ and $M=2.0$.

It should be noted that Figure 13 demonstrates the case with the most extreme differences in area distributions studied in the current work. These differences are less severe for lower Mach numbers and smaller $R/L$ values. The differences between these area distributions indicate that the Fourier coefficients representing the original geometry may not be adequate for predicting $C_{Dw}$ at higher Mach numbers. To address this, the Fourier series area distribution of Equation (8) is used to match the Mach plane sliced-area distribution generated using the OpenVSP linear wave drag tools. Once the necessary Fourier coefficients of Equation (8) are known, Equation (11) can be used to generate the analytic estimate for $C_{Dw}$.

To determine the number of Fourier coefficients necessary to adequately match a Mach plane sliced-area distribution with one generated using Equation (8), a study that minimized the root mean squared error between the two area distributions was completed. This study looked at a subset of the range of Mach numbers and $R/L$ values described previously. A gradient-based optimization routine was used to minimize the root mean square error (RMSE) between the area distribution produced using Equation (8) and the area distribution generated through the Mach plane slicing. Each case used an initial condition of zero for each of the Fourier coefficients, rather than initializing to a previous solution. The resulting $C_{Dw}$ and RMSE values for the cases at $M=2.0$ are shown in Figure 14.

These results show that convergence of the area distributions is achieved with around 20 Fourier coefficients. Some variation in the RMSE is seen beyond 20 Fourier coefficients, which is likely due to the initialization of each case being independent of the previous solution. Figure 14 also shows that the $C_{Dw}$ predicted using Equation (8) converged to the solution produced using the OpenVSP linear tools at around 16 coefficients, in most cases. The one exception to this is the $R/L=0.1$ case, which struggles with convergence as additional coefficients are used. It was decided that 20 coefficients would still be acceptable for this work. Lower priority was placed on the $C_{Dw}$ convergence for $R/L=0.1$ because convergence of the area distribution RMSE was clearly demonstrated and the $C_{Dw}$ results for this case will prove to have higher error in general.

3. Results and Discussion

The following study provides a comparison of converged $C_{Dw}$ results generated using the analytic, linear, and Cart3D Euler CFD methods. This study examined 110 cases over 11 different Mach numbers and 10 different $R/L$ values, ranging from 1.0 to 2.0 and 0.01 to 0.1 respectively. Contours of the resulting $C_{Dw}$ predictions are shown in Figure 15.

These results indicate that $C_{Dw}$ is largely independent of Mach number. The wave drag coefficient increases with increasing $R/L$ but little change is seen with an increase in Mach number. This trend is seen in both the linear and analytic results and confirmed by the Cart3D results for $R/L$ values below 0.05. Some Mach number dependence is seen in the Cart3D results for $R/L>0.05$. This indicates that Mach plane slicing is less impactful for these slender axisymmetric geometries where the longitudinal changes in cross-sectional area are inherently small. As mentioned in the Introduction, this trend in linear wave drag estimates for slender geometries was noted in the work by Rallabhandi and Mavris [3]. Based on these results, a good estimate for $C_{Dw}$ of a Sears–Haack body, and potentially other axisymmetric bodies, could be obtained using Equation (11) and the original Fourier coefficients used to define the geometry. This may not be the case for more complex geometries.

To complete this study, the difference in $C_{Dw}$ between the linear, analytic, and Cart3D solutions was computed at each Mach and $R/L$ combination. Contours of this difference in terms of drag counts are shown in Figure 16.

These results show that the difference between solutions is largely dependent on the $R/L$ values. The largest differences, upwards of 30 or more drag counts, correspond to the largest $R/L$ values. This is not entirely surprising because for $R/L=0.1$ the Sears–Haack geometry has a diameter of $20\%$ its total length, which is a very nonslender body. The assumption of a slender supersonic body was loosely defined in the development of the linear wave drag equation and these results appear to show that assumption being violated at the largest $R/L$ values. Evidence of this is seen by comparing example Cart3D solutions for the $R/L=0.1$ and $R/L=0.01$ models at Mach 1.4 in Figure 17. The larger radius model demonstrates much stronger and distributed shock and expansion features (note the $10\times$ increase in pressure-coefficient contour levels between the small radius solution in Figure 17a and the large radius solution in Figure 17b) compared to the distinct and weaker features expected in slender-body assumptions in the linear wave drag derivation.

The results of Figure 16a,b show very little dependence on Mach number. An exception to this is seen at the lowest Mach number cases. There appears to be a jump in the $C_{Dw}$ difference as the Mach number approaches 1.0. It should be noted that the data used to generate these contours only contain points for $M=1.1$ and $M=1.0$, meaning that it is not entirely clear from the current results how quickly the difference in $C_{Dw}$ is actually changing in that region. The contour plots in Figure 16 were generated using the Python-based Matplotlib contour function [30], which uses a marching squares method for generating the contour lines. A more sophisticated surface fitting approach, or simply increasing the number of design space samples, could improve the $C_{Dw}$ comparison in this Mach number region.

Figure 16c shows that the linear and analytic solutions match well over most of the Mach and $R/L$ design space. For geometries below an $R/L=0.05$, the difference in $C_{Dw}$ is less than approximately $2\%$ of a drag count for all Mach numbers. The increase in the $C_{Dw}$ difference at higher $R/L$ values is likely attributable to the same difficulties outlined in the linear and analytic method convergence studies discussed previously. Recall that for $R/L\approx 0.1$, converged solutions were difficult to obtain using both the analytic and linear methods.

In general, the more slender the axisymmetric body, and for Mach numbers greater than 1.0, the closer the linear and analytic $C_{Dw}$ solutions are to the Cart3D Euler CFD values. For solutions with $R/L$ larger than about 0.05, the difference in $C_{Dw}$ relative to the Cart3D solutions becomes greater than one drag count. All of this considered, the $C_{Dw}$ solutions produced by the linear and analytic methods still show trends similar to the Cart3D solutions at $R/L$ greater than 0.05, as seen in Figure 15. This similarity in trends implies that the lower-fidelity tools are still useful for the preliminary design analysis they were primarily intended for, even for less-slender bodies.

4. Conclusions

The accuracy of linear theory for axisymmetric bodies was considered here. A summary of the development of the wave drag integral, Equation (5), and the linear and analytic equations fundamental to solving it were included. Additionally, the equations representing the geometry of the Sears–Haack body were outlined. Finally, the concept of the equivalent body of revolution was discussed in addition to some background on the Cart3D Euler CFD tool. A comparison of different methods for predicting the wave drag of an axisymmetric body was also presented. The Sears–Haack geometry was the axisymmetric body of choice for these studies. Wave drag results produced using the OpenVSP-based linear method and an analytic approach using a Fourier sine series were compared to solutions generated using the Cart3D Euler CFD tool.

The convergence of the analytic, linear, and Euler CFD $C_{Dw}$ solutions was also considered. In general, good convergence was achieved using each of the different methods for all but the least-slender geometries defined in the design space. For the OpenVSP linear wave drag tool, grid-converged solutions were obtained using geometries defined with a distribution of 150 axial vertices and 100 tangential vertices, which resulted in axial and azimuthal vertex spacing ratios of $6.0\times {10}^{-3}$ and $1.0\times {10}^{-2}$, respectively. Additionally, $n_A=100$ axial slice locations were deemed adequate for producing converged solutions using the OpenVSP linear method. For the Cart3D solutions, grid-converged solutions were obtained using geometries defined with a distribution of 1600 axial vertices and 800 tangential vertices which resulted in a ratio of maximum triangular edge length to total length of between $6.4\times {10}^{-4}$ and $1.3\times {10}^{-3}$ for the $R/L=0.01$ and $R/L=0.1$ cases, respectively.

Converged $C_{Dw}$ solutions from the linear, analytic, and Euler CFD methods were generated over the design space of Mach number and $R/L$ combinations ranging in values from 1.0 to 2.0 and 0.01 to 0.1, respectively. The $C_{Dw}$ results were then used to produce contour plots of the difference in predicted $C_{Dw}$ between the linear, analytic, and Euler CFD solutions. These results showed that, for Mach numbers greater than 1.0, the more slender the axisymmetric body, the more accurate the solution would be relative to the Euler CFD solutions.

For bodies with $R/L\le 0.05$, the linear tool results showed a difference in $C_{Dw}$ relative to the Euler CFD results, ranging linearly from a single drag count ($1.0\times {10}^{-4}$) for bodies with $R/L\approx 0.05$ to around $0.001$ drag counts ($1.0\times {10}^{-7}$) for bodies with $R/L\approx 0.01$. Similar trends were seen for the analytic results, with a slight increase in error over the whole design space. At $R/L>0.05$, the difference in $C_{Dw}$ of both the linear and analytic methods begins to increase more rapidly with increasing $R/L$, and at $R/L\approx 0.1$, the difference is upwards of 30 or more drag counts ($30\times {10}^{-4}$) relative to the Euler CFD results. These large differences at the higher $R/L$ results seem to indicate that the slender-body assumption used in developing the linear wave drag equation is being violated. Even at these large $R/L$ cases, the linear and analytic methods appear to capture the same trends in $C_{Dw}$ that the Euler CFD results show over the design space.

Additional work could expand these types of studies to more complex supersonic configurations. These configurations might include axisymmetric geometries with distributions of area more representative of an aircraft’s fuselage or even wing-body combinations at a zero-lift flight condition. Future work could also revisit the comparison of off-body pressure-derived equivalent area distributions to the Mach plane sliced-area distributions and the potential for predicting the wave drag of axisymmetric bodies using these off-body equivalent area distributions.