THE DESIGN OF SHOCK-FREE TRANSONIC SLENDER BODIES

By Ron Buckmire A Thesis Submitted to the Graduate Faculty of Rensselaer Polytechnic Institute in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY Major Subject: Mathematics

Approved by the Examining Committee:

D.W. Schwendeman, Thesis Adviser

J.D. Cole, Thesis Adviser

A. Kapila, Member

Z. Rusak, Member

Rensselaer Polytechnic Institute Troy, New York June 1994 (For Graduation August 1994)

THE DESIGN OF SHOCK-FREE TRANSONIC SLENDER BODIES

By Ron Buckmire An Abstract of a Thesis Submitted to the Graduate Faculty of Rensselaer Polytechnic Institute in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY Major Subject: Mathematics The original of the complete thesis is on file in the Rensselaer Polytechnic Institute Library

Approved by the Examining Committee:

D.W. Schwendeman, Thesis Adviser

J.D. Cole, Thesis Adviser

A. Kapila, Member

Z. Rusak, Member

Rensselaer Polytechnic Institute Troy, New York June 1994 (For Graduation August 1994)

\begin{matrix}c \\ ◯ \end{matrix} Copyright 1994

by Ron Buckmire

All Rights Reserved

CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii ACKNOWLEDGEMENT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2. MATHEMATICAL THEORY . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1 Modeling the physical situation . . . . . . . . . . . . . . . . . . . . . 7 2.2 Transonic small-disturbance theory . . . . . . . . . . . . . . . . . . . 9 2.3 The hodograph transformation . . . . . . . . . . . . . . . . . . . . . . 12 2.4 The hodograph topology . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.4.1 Understanding the hodograph mapping . . . . . . . . . . . . . 21 2.5 Boundary value problem for R(w, ν) in the hodograph plane . . . . . 24

3. NUMERICAL PROCEDURE . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Change of variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Approximating the boundary value problem numerically . . . . . . . 27

3.2.1 Far-field boundary condition . . . . . . . . . . . . . . . . . . . 29 3.2.2 Grid refinement near free stream singularity at the origin … 31

3.3 Solving the hodograph boundary volume problem numerically …. 36 3.3.1 Solving in the elliptic region . . . . . . . . . . . . . . . . . . . 37 3.3.2 Solving in the hyperbolic region . . . . . . . . . . . . . . . . . 39

3.4 Computing the body of revolution from the numerical solution …. 40 3.5 Algorithm used in numerical solution of the hodograph BVP . . . . . 42

4. NUMERICAL RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1 Modeling the body of revolution . . . . . . . . . . . . . . . . . . . . . 45 4.2 Computed Shock-free Bodies . . . . . . . . . . . . . . . . . . . . . . . 48

4.3 Detailed features of shock-free solutions . . . . . . . . . . . . . . . . . 55 5. PHYSICAL PLANE CALCULATIONS . . . . . . . . . . . . . . . . . . . . 64

5.1 New numerical discretization scheme . . . . . . . . . . . . . . . . . . 68 5.2 Numerical solution of the physical plane boundary value problem . . 71 5.3 Numerical experiments in the physical plane . . . . . . . . . . . . . . 72

6. DISCUSSION AND CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . 80 LITERATURE CITED . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 A. The Dipole Strength Integral . . . . . . . . . . . . . . . . . . . . . . . . . . 84

A.1 Evaluating I_{1} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 A.2 Evaluating I_{2} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 A.3 Evaluating I_{3} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

B. Transonic Small-Disturbance Theory Asymptotics . . . . . . . . . . . . . . 90

B.1 Derivation of expression for pressure coefficient on the body . . . . . 93 C. Incompressible Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

LIST OF TABLES

Table 4.1 Values of β, α, ν_{∗} = S(x_{∗}; α) and ν^{∗} = S(x^{∗}; α) . . . . . . . . 48 Table 4.2 Values of α, β, K and M_{∞} for shock-free solutions . . . . . . . 48 Table 4.3 Values of K, D and a_{0}, a_{1} and a_{2} for shock-free solutions … 49 Table 4.4 Closure deviation values for computed shock-free bodies … 50

LIST OF FIGURES

Figure 2.1 Flow around a slender body of revolution . . . . . . . . . . . . 8 Figure 2.2 Boundary value problem in the physical plane . . . . . . . . . 12 Figure 2.3 Graph of a typical source distribution and its derivative …. 17 Figure 2.4 Graph of C(ν) showing different branches at ±∞ . . . . . . . 19 Figure 2.5 Diagram of the hodograph topology . . . . . . . . . . . . . . . 20 Figure 2.6 Domain in the physical plane . . . . . . . . . . . . . . . . . . 22 Figure 2.7 Domain in the hodograph plane . . . . . . . . . . . . . . . . . 23 Figure 3.1 Discretization of hodograph strip . . . . . . . . . . . . . . . . 27 Figure 3.2 Numerical stencil representing the partial differential equation 28 Figure 3.3 Grid geometry near origin . . . . . . . . . . . . . . . . . . . . 32 Figure 3.4 Interpolation method between coarse and fine grids . . . . . . 34 Figure 3.5 Plot showing R(0, ν) on successively-finer grids . . . . . . . . 36 Figure 4.1 α versus β . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Figure 4.2 Data for shock-free body number 1, SFB-1 . . . . . . . . . . . 51 Figure 4.3 Data for shock-free body number 2, SFB-2 . . . . . . . . . . . 52 Figure 4.4 Data for shock-free body number 3, SFB-3 . . . . . . . . . . . 53 Figure 4.5 Data for shock-free body number 4, SFB-4 . . . . . . . . . . . 54 Figure 4.6 Surface plot of R(u, ν) from numerical hodograph solution . . 55 Figure 4.7 Sonic line data used to compute the SFB-1 solution . . . . . . 56 Figure 4.8 Surface plots displaying effect of shocklessness on the jacobian 58 Figure 4.9 Characteristic curves associated with SFB-3 solution . . . . . 60 Figure 4.10 Body function F (x) data for the SFB-2 . . . . . . . . . . . . . 62 Figure 4.11 Source distribution S(x) data for the SFB-2 . . . . . . . . . . 62 Figure 4.12 S′(x) data for the SFB-2 . . . . . . . . . . . . . . . . . . . . . 63 Figure 4.13 Plot of G′(x) data for the SFB-2 . . . . . . . . . . . . . . . . 63 Figure 5.1 Exponentially-scaled numerical grid employed in (x, r̃) plane . 65 Figure 5.2 Γ = u(r̃_{0}) − S log r̃_{0} versus r̃_{0} (for S = 1, G = 5) . . . . . . . . 69 Figure 5.3 Contour plot of u = ϕ_{x}(x, r̃) for α = −.0925 and K = 3.919 . . 73 Figure 5.4 Physical plane G′(x) data for α = −.0925 and K = 3.919 … 73 Figure 5.5 Difference between input and output source distribution data

for SFB-2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

Figure 5.6 Contour plot of u = ϕ_{x}(x, r̃) for SFB-2 at K = 3.919 . . . . . 75 Figure 5.7 Plot of G′(x) for SFB-2 at K = 3.919 . . . . . . . . . . . . . . 75 Figure 5.8 Isobars around hodograph-designed body SFB-2 at off-design

conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

Figure 5.9 G′(x) associated with SFB-2 for various K above and below design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

LIST OF SYMBOLS

A^{∗} scaled cross-sectional area a local speed of sound a_{∞} free-stream sound speed C_{p} pressure coefficient C_{p_{body}} pressure coefficient on the body D dipole strength F (x) thin airfoil profile, r = δF (x) J Jacobian of the Hodograph transfor-

mation K transonic similarity parameter M local Mach number M_{∞} free-stream Mach number p pressure p_{∞} free-stream pressure r^{∗} inner asymptotic spatial variable r̃ outer asymptotic spatial variable R stream-tube area S(x) source distribution function U_{∞} nondimensional free-stream speed (x, r) cartesian coordinates in physical plane w perturbations from sonic speed α body symmetry factor β grid symmetry factor δ thickness ratio of slender body ϵ compressibility factor γ ratio of specific heats (u, ν) hodograph independent variables ρ density ρ_{∞} free-stream density ϕ outer disturbance potential function φ inner disturbance potential function Φ “exact” velocity potential function (x_{∗}, x^{∗}) extrema points of source distri-

bution (u_{∗}, u^{∗}) left and right boundaries of hodo-

graph plane (ν_{∗}, ν^{∗}) lower and upper boundaries of

hodograph plane

ACKNOWLEDGEMENT

I would like to thank my advisors, Professors Don W. Schwendeman and Julian D. Cole. Dr. Schwendeman’s role in the editing of this thesis is greatly appreciated. Thanks also go to Harriet Borton for assistance with layout and placement of the diagrams.

I would like to gratefully aknowledge the support of Air Force Office of Scientific Research Grant AFOSR 88-0037 and F49620-93-1-0022DEF.

ABSTRACT

The existence of shock-free flow around a slender body of revolution at near-sonic speeds is investigated using transonic small-disturbance theory. The governing partial differential equation, known as the Kármán-Guderley equation, and the boundary conditions are transformed into the hodograph plane. In the hodograph plane the spatial variables depend on the velocity components. The problem is solved numerically, employing an algorithm that involves a combination of finite-difference and iterative methods. A condition dependent upon the Jacobian of the transformation is developed to determine when a shock-free solution has been computed. A number of bodies possessing shock-free flows are calculated at different values of the transonic similarity parameter, K = (1 − M_{∞}^{2})/(δM_{∞})^{2}, where M_{∞} is the flow Mach number and δ is the body thickness. The body profiles computed are both fore-aft symmetric and fore-aft asymmetric. For moderate values of K (e.g. K = 3.5, corresponding to M_{∞} = 0.985 and δ = 0.1) there is little difficulty in finding shock-free solutions in the hodograph. Solutions have not been calculated for K less than about 3, corresponding to speeds very close to sonic.

Calculations in the physical plane of the flow around transonic slender bodies of revolution are performed at moderate values of K in order to confirm the shockfree nature of the body profiles obtained from the hodograph calculations and to explore off-design conditions. It is found that the flow field in the physical plane around the shock-free hodograph-designed bodies appears to be nearly shock-free, having at most a weak shock. Small perturbations to K or to the hodographdesigned body shape do not appear to create qualitatively different flow fields. The numerical evidence suggests that shock-free flows are isolated but that nearby flows possess, at most, weak shocks.

CHAPTER 1 INTRODUCTION

The problem of the existence of shock-free flows around slender axisymmetric bodies is the subject of this thesis. The question is an important one because it can be related to the efficiency of aircraft operating at high subsonic speeds. The transonic area rule given in Cole and Cook [4] implies that if shock-free flows can exist around bodies of revolution, then other bodies which have identical rates of change of crosssectional area can also possess flows with reduced drag due to their shock-free nature.

The specific problem posed is whether bodies of revolution can be constructed so that transonic flow is accelerated and then decelerated around the body in such a manner that despite the presence of a local supersonic region, no shock develops. The bodies possess local supersonic flow regions where the characteristic lines do not cross to form a shock. The accelerated flow over such bodies is able to decelerate smoothly at the back of the body and join the trailing transonic flow without forming a shock.

The idea that transonic shock-free flows exist is not a new one. Ringleb [21] was the first to postulate that shock-free transonic flow could exist, and demonstrated this using a two-dimensional nozzle flow. There have been many investigations of transonic flows around airfoils which are shock-free. After Pearcey [20] and Whitcomb and Clark [25] published experimental results that showed that shockfree flows could be achieved in practice, research interest was heightened. The first analytically designed shock-free transonic airfoil was a quasi-elliptical wing section produced by Nieuwland in [19]. Bauer, Garabedian and Korn in [1, 2, 3] built on the work of Nieuwland to produce numerical algorithms to systematically compute shock-free airfoils. The design method involved transforming the problem to the hodograph plane and using complex characteristics. In the hodograph plane, the dependence of the flow velocity components on the spatial co-ordinates is reversed. The algorithm used by Bauer et al. in [1, 2, 3] involves solving an inverse problem of first computing the shock-free flow in the hodograph and then determining the associated airfoil producing it. Wind tunnel experiments conducted by Kacprzynski et al. [10, 11, 12] verified that the computed airfoils of Bauer et al. do indeed possess shock-free supercritical flows. An airfoil is said to be supercritical if at some point the local flow speed over it exceeds the local speed of sound. The method of Sobieczky et al. in [23, 24] which involves looking at the physical plane characteristic curves in the supersonic zone and then trying to change flow parameters to disentangle them to produce a shock-free flow is a useful design tool. In a series of famous papers by Morawetz [16, 17, 18], she proved that if there exists shock-free flow around a two-dimensional airfoil, it is impossible to infinitesimally perturb the airfoil shape to produce another neighboring shock-free flow. The implication is that shock-free solutions of the equations governing transonic fluid flow are isolated.

The only known previous research that has been conducted concerning shockfree transonic flows in more than two dimensions is the work of Cole and Schwendeman [8], who computed a transonic shock-free flow around a fore-aft symmetric slender body of revolution. In this thesis the methods in [8] are expanded to calculate fore-aft asymmetric shock-free bodies. The algorithm used is influenced by ideas from the previous work around supercritical airfoils mentioned above. In some sense the work in [8] and in this thesis can be considered the first forays into the hodograph design of shock-free supercritical bodies.

The flow around a slender body traveling through a fluid at transonic speeds is modeled mathematically with the use of transonic small-disturbance theory as detailed in Cole and Cook [4]. A boundary value problem for the velocity perturbation potential of the flow is formulated using a triple-deck asymptotic analysis. Solution of the boundary value problem for the velocity potential function determines the outer flow field about the body. The approach of transforming the small-disturbance problem into the hodograph is applied. All continuous solutions in the hodograph plane are by definition shock-free flows in the physical plane. The body shape in the physical plane associated with the computed shock-free hodograph solution is then calculated. It is shown that the condition for whether the body can be regenerated from the hodograph solution is that the jacobian of the hodograph transformation is strictly less than zero at all points in the flow.

The results obtained in this thesis show that shock-free flows can be computed regularly. At least a half-dozen different shock-free solutions have been computed. A systematic algorithm to develop more such supercritical bodies is given. It should be clarified that the theorem Morawetz [16, 17, 18] proved was for supercritical airfoils, not supercritical bodies. There is no known corresponding theorem or proof for shock-free transonic flow around bodies, though the isolated nature of such solutions is widely believed to be true. The results in this thesis appear to suggest that if shock-free solutions are isolated, they appear to be surrounded by solutions which contain relatively weak shocks. This comment stems from the results that occur when the bodies which are believed to be shock-free are also tested by computing flows around them in the physical plane. The numerical results have features which may be due to the presence of weak shocks or which may just be unable to be resolved from the computed data. Regardless, it can certainly be said that for some range of off-design parameters no substantial shocks are observed around the hodograph-designed supercritical bodies.

Besides the inverse hodograph design methods of Cole and Schwendeman [8] and Bauer et al. [1, 2, 3], it is possible to attempt to design directly in the physical plane, as Sobieczky [23, 24] does. In the physical plane the algorithm used involves selecting a provisional body and then computing to see what the flow around such a body is. The body determines the flow. If the computed flow is shock-free then the initial body shape guess was a good one. A priori, the two algorithms look equivalent but it is well known that direct numerical calculation of the flow in the physical plane is more sensitive to perturbations in the body shape. In the hodograph, the computed body shape is not considered to be as sensitive to perturbations of the sonic bubble shape, so the hodograph method was selected. In the physical plane the criterion for a shock-free flow is that no discontinuities appear in the flow. In the hodograph plane the criterion is that the jacobian of the hodograph transformation does not become zero. This ensures that a shock-free flow computed in the hodograph plane can be transformed back to the physical plane.

1.1 Thesis outline

In Chapter 2 a mathematical description of the problem is discussed, involving the formulation of the boundary value problem in both the physical plane and in the hodograph plane. The governing equation is a mixed-type elliptic-hyperbolic partial differential equation. A detailed explanation of the hodograph mapping from the physical plane to the hodograph plane is given. Another good description of the hodograph mapping can be found in Kropinski [13]. A clear understanding of the boundary value problems in both the physical plane and the hodograph plane is necessary to facilitate comprehension of the numerical solution of the problem. For example, the slender body of revolution near the origin in the physical plane is mapped to infinity in the hodograph plane. In slender body theory the representation of the body in the physical plane is accomplished by a distribution of sources along the origin. The free stream flow far away from the body in the physical plane becomes an extremely singular point, located at the origin, in the hodograph. The flow-field behavior far away from the body is modeled by the flow due to a dipole situated at the physical plane origin. The entire semi-infinite physical plane is mapped into a thin infinitely long strip in the hodograph. The width of this hodograph strip is dependent upon the points where the source distribution function has its maxima and minima.

In Chapter 3 the boundary value problem in the hodograph is rewritten in scaled variables and then it is approximated numerically on a discrete grid. The boundary value problem is converted to a system of nonlinear algebraic equations defined on a series of overlapping numerical grids approximating the elliptic section of the hodograph plane. The overlapping grids have increasingly fine meshes, which is necessary to capture the rapid change in the solution near the singularity at the origin. The system of nonlinear equations is solved using a combination of point relaxation and Newton’s method. Using the newly-calculated solution in the elliptic region as initial data, a time-like marching scheme is used to compute the solution into the hyperbolic region. The jacobian of the hodograph transformation is checked during this marching scheme and if it is found to be negative everywhere a shock-free solution has been computed. The solution is not yet complete however. A presumed value was used for a constant which appears in the far-field boundary condition (the dipole strength of the far-field flow) and depends on the solution. A fixed-point iteration on this constant is conducted until the presumed and actual values are one and the same, within tolerance. Quantities of interest in the physical plane, such as the pressure coefficient on the body and the body itself, are computed from the complete hodograph solution.

In Chapter 4 the results obtained after using the algorithm discussed in Chapter 3 to produce a number of shock-free solutions are displayed. Contour plots of the isobars around the body, as well as graphs of the pressure coefficient on the body, are given to support the claim that no shock is present in the computed flows. Tables of the identifying values of parameters used to compute each particular shock-free body are also given. Surface plots, with contours, depicting the jacobian of shockless and shocked solutions in the hodograph are given in order to highlight the difference between the two situations. The characteristic curves in the hyperbolic region of the hodograph plane and those in the supersonic zone of the physical plane are displayed to show that in the computed shock-free solutions the characteristics do not intersect to form shocks.

In order to investigate off-design behavior of bodies which were computed to be shock-free in the hodograph a code was developed to solve the original problem in the physical plane. It is difficult to attempt studies of off-design conditions in the hodograph plane. This is because the body shape is a result of hodograph computations of a given shock-free flow at a given particular speed. Off-design computations require being able to do the reverse: compute a flow at a given speed around a given body. The numerical approximation of the physical plane boundary value problem and conversion to a system of nonlinear algebraic equations is given in Chapter 5. Contour plots of the approximate isobars computed directly in the physical plane as well as the plots of the pressure coefficient for the hodograph-designed bodies are given in Chapter 5, to compare with the results displayed in Chapter 4. In the direct physical plane calculations, the flow around hodograph-designed bodies appears to contain relatively weak shocks. This is a dramatic reduction from the distinct shocks that are computed in the flow around bodies that were used as input into the hodograph algorithm.

The appendices provide details on some finer mathematical points that arise in Chapter 2. Appendix A gives the derivation for the volume integral defining the dipole strength of the far-field flow about the body. Appendix B details the asymptotics used in the formulation of the physical plane boundary value problem and contains the derivation of the expression for the pressure coefficient on the body.

CHAPTER 2 MATHEMATICAL THEORY

The mathematical modeling of a slender body of revolution moving through a fluid at transonic speeds is recounted in this chapter. The model of the physical situation is examined and the governing equations which describe this flow are detailed. A boundary value problem is formulated which, when solved determines the flow field everywhere. The hodograph plane and its properties are introduced and the transformation of the boundary value problem to this plane is described. The resulting boundary value problem in the hodograph plane is outlined.

2.1 Modeling the physical situation

Shock-free transonic flow around a slender body of revolution can be modeled by inviscid, compressible, steady and irrotational flow. The assumption of inviscid flow is reasonable since high-speed flow around a slender body of revolution can be expected to have relatively thin boundary layers and minimal viscous effects. The flow must be compressible since the flow is transonic. It is assumed that the flow is steady. To the order of approximation considered,the flow is also assumed to be irrotational.

An example of a typical flow past a slender body of revolution at subsonic speed is shown in Figure 2.1. By definition, a body of revolution is axisymmetric, so the variables x and r = \sqrt{y^{2} + z^{2}} can be used to describe the three-dimensional space about the body. The body is defined by r = δF (x) where −1 ≤ x ≤ 1. The flow at distances far from the body is uniform with speed U_{∞} and flow Mach number M_{∞} = U_{∞}/a_{∞}, where a_{∞} is the speed of sound in the free stream. If M_{∞} is close to one, a supersonic zone forms about the body as shown in Figure 2.1. For a typical body, a shock wave appears and it is the main goal of this study to construct body functions F (x) that possess a significant supersonic zone but are shock-free. This would be possible if a body could be constructed so that the characteristics in the supersonic zone emerging from the body do not cross to form a shock.

Figure 2.1: Flow around a slender body of revolution

The equations governing inviscid, compressible, steady, potential flow are

(a^{2} − Φ_{x}^{2})Φ_{xx} + (a^{2} − Φ_{r}^{2})Φ_{rr} + \frac{a^{2}}{r} Φ_{r} − 2Φ_{r}Φ_{x}Φ_{xr} = 0  (2.1)

and

\frac{a^{2}}{γ − 1} + \frac{Φ_{x}^{2} + Φ_{r}^{2}}{2} = \frac{a_{∞}^{2}}{γ − 1} + \frac{U_{∞}^{2}}{2} ,  (2.2)

where Φ(x, r) is the velocity potential and a(x, r) is the local speed of sound. Equation (2.1) is a statement of conservation of mass and (2.2) is Bernoulli’s integral. The constant γ is the ratio of specific heats. The boundary conditions for the velocity potential are

Φ(x, r) = U_{∞}x, as (x^{2} + r^{2}) → ∞, (2.3)

and

δF′(x) = \frac{Φ_{r}(x, r)}{Φ_{x}(x, r)},  on the body r = δF (x).  (2.4)

That is, flow far away from the body is uniform (2.3), and the flow along the surface of the body is tangent to the body (2.4). For flow past a slender body of revolution whose thickness δ is small and whose flow speed Mach number M_{∞} is close to unity, the full potential equation (2.1) can be approximated using transonic small-disturbance theory (TSDT), as discussed in the next section.

2.2 Transonic small-disturbance theory

Transonic small-disturbance theory shows that asymptotic expressions for flow around slender bodies of revolution with body thickness ratio δ can be written in terms of this small parameter. A brief description of the equations involved will be outlined in this section. For further details, see [4] and Appendix B. According to TSDT, a velocity potential Φ(x, r; M_{∞}, δ) exists and has an outer expansion which is valid as δ → 0 with x, r̃ = δM_{∞}r, and K = (1 − M_{∞}^{2})/(M_{∞}δ)^{2} all held fixed. The parameter K is referred to as the transonic similarity parameter. The outer expansion has the form

Φ = U_{∞}{x + δ^{2}ϕ(x, r̃; K) + …},  (2.5)

where the outer disturbance potential ϕ(x, r̃) satisfies the transonic small-disturbance equation (TSDE), also known as the Kármán-Guderley equation,

(K − (γ + 1)ϕ_{x})ϕ_{xx} + ϕ_{r̃r̃} + \frac{1}{r̃}ϕr̃ = 0.  (2.6)

An observation of the above equation reveals that it is of mixed elliptic-hyperbolic type. It is elliptic when ϕ_{x} < K/(γ + 1), corresponding to subsonic flow, and hyperbolic when ϕ_{x} > K/(γ + 1), corresponding to supersonic flow.

The near-field boundary condition for (2.6) comes from asymptotic matching with an inner expansion of Φ. This inner expansion is valid as δ → 0 with x, r^{∗} = r/δ, and K all held fixed and satisfies the tangent-flow boundary condition on the surface of the slender body (2.4). The inner expansion takes the form

Φ = U_{∞}{x + δ^{2} log(δ^{2}M_{∞})S(x) + δ^{2}φ(x, r^{∗}) + …}.  (2.7)

The “switch-back” term, δ^{2} log(δ^{2}M_{∞})S(x), is derived from an intermediate asymptotic matching as explained in [4] and Appendix B. The inner disturbance potential φ(x, r^{∗}) satisfies Laplace’s equation in cylindrical co-ordinates in each constant-x cross plane

∇^{∗^{2}}φ ≡ φ_{r^{∗}r^{∗}} + \frac{1}{r^{∗}} φ_{r^{∗}} = 0, |x| < 1. (2.8)

The general solution to (2.8) is

φ(x, r^{∗}) = S(x) log r^{∗} + G(x; K).  (2.9)

The function of integration S(x) is determined from the tangency boundary condition

δF′(x) = \frac{Φ_{r}(x, δF)}{Φ_{x}(x, δF)} = \frac{\{δ^{2} \frac{S(x)}{δF} + …\}}{\{1 + …\}} = \frac{δS(x)}{F (x)}  (2.10)

and from this one can see that

S(x) = F (x)F′(x) = \frac{1}{2π} \frac{dA^{∗}}{dx} ,  (2.11)

where A^{∗} = πF^{2}(x) is the scaled cross-sectional area of the body. The other function of integration, G(x; K), is an unknown function resulting from the solution to the partial differential equation (2.8) and is calculated by asymptotic matching to the solution of the outer problem after the perturbation potential ϕ is known. It is noted that the pressure coefficient C_{p}(x; K) on the surface of the body depends upon the function G′(x).

A derivation for C_{p_{body}} is given in Appendix B. The result is

C_{p_{body}}(x) = −δ^{2}{2S′(x) log(δ^{2}M_{∞}F (x)) + 2G′(x) + (F′(x))^{2}}.  (2.12)

Asymptotic matching of the inner and outer disturbance potentials as (r^{∗} → ∞ and r̃ → 0) shows that the boundary condition near the body for the outer disturbance potential ϕ(x, r̃) is a flow represented by a singularity distribution of sources along the axis of the body. The matching requires that

ϕ(x, r̃) = S(x) log r̃ + G(x; K),  as r̃ → 0, |x| < 1.  (2.13)

The far-field boundary condition for (2.6) is that velocity perturbations caused by the body of revolution attenuate at infinity. That is, (ϕ_{x}, ϕr̃) → 0 as (x^{2} + r̃^{2}) → ∞ so that the TSDE (2.6) becomes

Kϕ_{xx} + ϕ_{r̃r̃} + \frac{1}{r̃}ϕr̃ = 0,  (2.14)

to a first approximation. This is known as the Prandtl-Glauert equation. The relevant solution to (2.14) for a closed body is that given by a dipole. An analytic expression for the outer disturbance potential in the far-field can be written and the exact functional representation is

ϕ(x, r̃) = \frac{D}{4π} \frac{x}{(x^{2} + Kr̃^{2})^{3/2}} ,  (2.15)

where D is the dipole strength.

A relationship between the dipole strength D in (2.15), the outer perturbation potential ϕ(x, r̃) and the body shape F (x) can be derived by treating the non-linear terms in the TSDE (2.6) as a right hand side to the Prandtl-Glauert equation (2.14) and applying Green’s theorem. A brief derivation of an expression for the evaluation of D is given in Appendix A. The result is

D = D_{body} +D_{flow} = π ∫_{−1}^{+1} F^{2}(x) dx + π(γ +1) ∫_{−∞}^{∞} dx ∫_{0}^{∞} ϕ_{x}^{2}(x, r̃)r̃dr̃. (2.16)

Figure 2.2: Boundary value problem in the physical plane

In summary, the boundary value problem for the outer disturbance potential ϕ consists of the TSDE (2.6) defined for |x| < ∞, r̃ > 0 subject to the boundary condition (2.9) as r̃ → 0 for |x| < 1 (ϕ must be bounded as r̃ → 0, |x| > 1) and the far-field behavior given by (2.15) with (2.16). Figure 2.2 summarizes the physical plane boundary value problem for ϕ(x, r̃).

2.3 The hodograph transformation

In the previous section, the boundary value problem has been written in terms of the perturbation potential ϕ which is dependent upon independent spatial variables (x, r̃). In the hodograph plane, this dependence is reversed and the spatial variables are written as functions of the velocity components (ϕ_{x}, ϕr̃). The aim of this section is to discuss the transformation of the physical boundary value problem to the hodograph.

In order to transform the physical plane equations, they are written as a system of first-order PDEs. One choice of dependent variables to do this is

w = (γ + 1)ϕ_{x} − K,  ϑ = (γ + 1)ϕr̃,  (2.17)

where w measures perturbations from the sonic speed and ϑ measures the flow deflection from horizontal flow. The choice of variables (w, ϑ) follows that of [4], but later it will be changed for convenience in solving the equations numerically. When w is positive the flow is supersonic, when w is negative the flow is subsonic and when w equals zero the flow is sonic. Using these variables, (2.6) becomes

\{ \begin{matrix}ww_{x} = ϑr̃ + \frac{1}{r̃}ϑ \\ ϑ_{x} = wr̃ \end{matrix} \}  (2.18)

It is noticed that (2.18) can be simplified by choosing the variables

R = \frac{r̃^{2}}{2} ,  ν = r̃ϑ,  (2.19)

where R is a measure of the stream-tube area (note R ≥ 0) and ν is the radial mass flux. Using (2.19), (2.18) becomes

\{ \begin{matrix}ww_{x} = ν_{R} \\ 2Rw_{R} = ν_{x} \end{matrix} \}  (2.20)

The hodograph transformation is performed by using

\begin{matrix}w_{x} = \frac{1}{J} R_{ν}, & ν_{x} = −\frac{1}{J} R_{w}, \\ w_{R} = −\frac{1}{J} x_{ν}, & ν_{R} = \frac{1}{J} x_{w}, \end{matrix}  (2.21)

where

J = Jacobian of the transformation =\frac{∂(x, R)}{∂(w, ν)} = x_{w}R_{ν} − x_{ν}R_{w}.  (2.22)

The transformed system is

\{ \begin{matrix}wR_{ν} = x_{w} \\ 2Rx_{ν} = R_{w} \end{matrix} \}  (2.23)

A number of observations can be made about this new system of equations. The variables (x, R) have moved from being independent variables to being dependent variables. The factor \frac{1}{J} cancels out of the system (2.20) during the inversion process so that the jacobian does not appear in the transformed system (2.23).

Eliminating x from system (2.23), the partial differential equation that has to be solved for R(w, ν) is

(\frac{R_{w}}{2R} )_{w} − wR_{νν} = 0.  (2.24)

The above hodograph version (2.24) of the TSDE is hyperbolic in the supersonic region (w > 0) and elliptic in the subsonic region (w < 0). The behavior of the transformed partial differential equation is similar to Tricomi’s equation except that it is non-linear. They are both mixed-type partial differential equations. In the physical plane the region where the flow equation switches from being elliptic to hyperbolic depends on the solution to the equation and can not be predicted beforehand. In the hodograph plane, the equations are somewhat simpler in that the elliptic and hyperbolic regions are known from the onset of the calculation. In the physical plane, the determining factor in the equation (2.6) was the unknown term ϕ_{x}, but in the hodograph version of the equation (2.24) it is a known value of w, namely zero.

Eliminating x from the expression (2.22) for the Jacobian yields

J(w, ν) = wR_{ν}^{2} − \frac{R_{w}^{2}}{2R} .  (2.25)

Clearly, R is non-negative by definition in (2.19) so J < 0 in the subsonic region w < 0. In order for the mapping from the hodograph to the physical plane to be smooth, J must also be negative in the supersonic region w > 0, except possibly at isolated points. This is an essential requirement for a shock-free supersonic region. This is akin to checking that the flow characteristics in the physical plane are not intersecting to form a shock in the local supersonic region. The jacobian is the quantity that will be computed and checked to determine if a shock-free solution has been found. If it is true that J is negative everywhere then the hodograph mapping is order-reversing everywhere.

In the next section the topology of the hodograph plane, where the problem is solved, will be discussed.

2.4 The hodograph topology

The topology of the hodograph can be described by observing the result of applying the hodograph transformation to the boundary value problem in the physical plane as formulated in section 2.2 and illustrated in Figure 2.2.

In the physical plane, the far-field condition is that the disturbance caused by the body of revolution attenuates to zero far away from the body. The velocity perturbations (ϕ_{x}, ϕr̃) → 0, which correspond to w → −K and ν → 0 by substitution in (2.17). The entire far-field flow in the physical plane, (x^{2} +r̃^{2}) → ∞, maps to the single point (−K, 0) in the hodograph plane. The uniform free stream at infinity in the physical plane becomes a singularity in the hodograph, known as the free stream singularity. This is because all values of (x, R) emerge from the free-stream point w = −K, ν = 0. The flow behavior in the far-field of the physical plane is modeled by flow due to a dipole. Writing the dipole behavior in (2.15) in terms of the variables (w, ν) and (x, R) yields an implicit representation for the behavior of R(w, ν) near w = −K, ν = 0

\{ \begin{matrix}ν & = (γ + 1)r̃\frac{∂ϕ}{∂r̃} = & −(γ + 1)\frac{D}{4π} \frac{6KxR}{(x^{2} + 2KR)^{5/2}} \\ w + K & = (γ + 1)\frac{∂ϕ}{∂x} = & (γ + 1)\frac{D}{4π} \frac{2KR − 2x^{2}}{(x^{2} + 2KR)^{5/2}} \end{matrix} \}  (2.26)

The expression (2.16) for the dipole strength D in the physical plane involves an integral for D_{flow} which can also be written in hodograph variables. The result is

D_{flow} = \frac{π}{γ + 1} ∫_{−∞}^{∞} dx ∫_{0}^{∞}(w + K)^{2} dR = \frac{π}{γ + 1} ∫_{−∞}^{∞} dw ∫_{ν_{∗}}^{ν^{∗}} |J|(w + K)^{2} dν (2.27)

where ν ∈ (ν_{∗}, ν^{∗}) is an assumed bound for ν in the hodograph plane.

Near the body in the physical plane the perturbation potential is given by

ϕ(x, r̃) = S(x) log r̃ + G(x),  r̃ → 0, |x| < 1.

This can be differentiated to give

\{ \begin{matrix}ν & = (γ + 1)r̃\frac{∂ϕ}{∂r̃} = & (γ + 1)S(x) + … \\ w + K & = (γ + 1)\frac{∂ϕ}{∂x} = & (γ + 1)S′(x) log \sqrt{2R} + (γ + 1)G′(x) + … \end{matrix} \} (2.28)

The expressions in (2.28) show that as R → 0 (r̃ → 0), ν → (γ + 1)S(x) and w → ±∞ depending on the sign of S′(x). Thus, the body of revolution in the physical plane is mapped to w = ±∞ with (γ+1)min(S(x)) ≤ ν ≤ (γ+1)max(S(x)).

The definition of the source distribution function S(x) in (2.11) shows that it is related to the shape of the body F (x) and its derivative F′(x). It is also proportional to the rate of change of the cross-sectional area of the body. A typical source function is shown in Figure 2.3(a) where ν_{∗} ≡ (γ + 1)min(S(x)) and ν^{∗} ≡ (γ + 1)max(S(x)). These extrema occur at x = x^{∗} and x = x_{∗} where S′ = 0. In this study, it will be assumed that S′(x) has only two roots on x ∈ [−1, 1], where the body is defined.

(a) Source Distribution S(x)

(b) S′(x)

Figure 2.3: Graph of a typical source distribution and its derivative

In order to get a sense of what properties the hodograph dependent variables, R(w, ν) and x(w, ν) have near the body, expressions for them can be found by inverting the equations which represented the near-field boundary conditions in the physical plane (2.28). The dominant terms after inverting (2.28) are

\{ \begin{matrix}R(w, ν) = A(ν)e^{B(ν)w} + … \\ x(w, ν) = C(ν) + … \end{matrix} \} ,  (2.29)

where

A(ν) = \frac{1}{2} exp [B(ν) (K − (γ + 1)G′(C(ν)))]
B(ν) = \frac{2}{(γ + 1)S′(C(ν))}
C(ν) = inverse function of (γ + 1)S(x).

The form of the function C(ν) changes depending on which side of the hodograph strip it is evaluated on. This is obvious by looking at Figure 2.3(a) and recalling that C(ν) is essentially the inverse of the function plotted. It can be defined as

C(ν) = \{ \begin{matrix}C_{∗}(ν) & if w → −∞ and ν_{∗} ≤ ν ≤ ν^{∗} \\ C^{∗}(ν) & if w → +∞ and ν_{∗} ≤ ν ≤ ν^{∗} \end{matrix}

where C^{∗}(ν) and C_{∗}(ν) are plotted in Figure 2.4.

(a) Graph of C_{∗}(ν)

(b) Graph of C^{∗}(ν)

Figure 2.4: Graph of C(ν) showing different branches at ±∞ By looking at (2.29) it is seen that R decays exponentially to zero as w → ±∞

with ν fixed, depending on the sign of B(ν), while x approaches a constant, C(ν).

Looking at the plot of S′(x) in Figure 2.3(b) it is seen that for −1 < x < x^{∗} and x_{∗} < x < 1, S′(x) and thus B(ν) are positive so that these segments of the body are mapped to w = −∞, ν_{∗} < ν < ν^{∗} in the hodograph. The section of the body for x^{∗} < x < x_{∗} where S′(x) and B(ν) are negative lies at w = +∞, ν_{∗} < ν < ν^{∗}. The points x = x^{∗} and x = x_{∗} on the body are mapped to lines ν = (γ + 1)S(x^{∗}) = ν^{∗} and ν = (γ+1)S(x_{∗}) = ν_{∗} along which w takes all values from −∞ to +∞. Thus the domain in the physical plane is mapped to an infinite rectangular strip, ν_{∗} < ν < ν^{∗}, |w| < ∞ in the hodograph plane.

Figure 2.5: Diagram of the hodograph topology

The two extrema x_{∗} and x^{∗} are important features of the body since they determine the bounds ν_{∗} and ν^{∗} of the hodograph plane (see Figure 2.5). They also determine the exact spot on a transonic slender body that a shock will form. In a paper by Cole and Malmuth [5] it is proved by applying transonic small-disturbance theory that if a shock wave forms on a slender body of revolution, it must do so at x^{∗} or x_{∗}. Also, the supersonic zone as displayed in Figure 2.1 is anchored to the body at these two points. If the direction of the flow is from nose to tail, the shock is expected to form at x_{∗}, the extrema point nearer to the tail. Later it will be observed that the greatest numerical difficulties arise near these extrema points, because they are the site of the activity which leads to shock discontinuities.

Another important feature of the problem is the property that R(w, ν) decreases exponentially to zero in the hodograph strip as one approaches the body of revolution at w = ±∞. Knowledge of this property of the solution implies that in practice one does not need to evaluate R on an infinite strip in order to approximate the near-field boundary condition of R → 0 as one approaches the body, since R will be exponentially small at moderate values of w. In fact, because R decays exponentially one would think that the domain could be truncated without difficulty. However, since G(x; K) and other data on the body itself are desired, and the body is situated at w ± ∞ in the hodograph, care must be taken when truncating the domain. There is important information to be gleaned from the exponential “tail” of the solution at w = ±∞.

2.4.1 Understanding the hodograph mapping

Further insight into the nature of the hodograph mapping can be derived by comparing Figure 2.6 and Figure 2.7. This is accomplished by tracing a path along a contour at r̃ = 0 from points A to F in the physical plane and observing how it is mapped to the hodograph plane. Point A in Figure 2.6 is in the far-field where the flow is uniform. This point gets mapped to the free stream singularity at w = −K, ν = 0 in the hodograph (point A′). Moving from A to the nose at point B corresponds to moving from the free stream singularity along the branch cut until one reaches the point B′ at the leftmost edge of the hodograph strip at w = −∞. The points B, C, D and E are all on the body of revolution itself, which is mapped to the perimeter of the hodograph plane. As one moves from the nose to the special point x^{∗} (point C in Figure 2.6) this corresponds to moving up along the leftmost

Figure 2.6: Domain in the physical plane

boundary to the top of the strip at point C′. Moving from x^{∗} to x_{∗} in the physical plane corresponds to moving along the rightmost edge of the hodograph strip at w = +∞, from C′ to D′. The points x^{∗} and x_{∗} correspond to the entire upper and lower horizontal boundaries of the hodograph strip itself. As one moves from the second extrema x_{∗} at point D to the tail of the body in the physical plane at point E this corresponds to moving up from the bottom left corner of the hodograph region to just below the branch cut along ν = 0. The tail (B) and the nose (E) are mapped to two points which are very close to each other (B′ and E′) in the hodograph plane, but exist on different sides of the branch cut. Proceeding from the tail back to the far-field far away from the body (point F) in the physical plane corresponds to moving in along ν = 0^{−} from the leftmost edge of the hodograph at w = −∞ to the free stream singularity at w = −K (point F′). A second contour in the physical plane at r̃ = r̃_{c} > 0 is mapped to a slightly curved contour just inside the rectangular perimeter contour in the hodograph. This second contour r̃_{c} is mapped to the dotted contour shown in Figure 2.7. Contours in the physical

Figure 2.7: Domain in the hodograph plane

plane very far away from the body in the physical plane in Figure 2.6 are mapped to the small contours near the free stream singularity shown in Figure 2.7.

2.5 Boundary value problem for R(w, ν) in the hodograph plane

The boundary value problem in the hodograph plane is summarized in this section. The partial differential equation to be solved in the hodograph strip is

(\frac{R_{w}}{2R} )_{w} − wR_{νν} = 0

The near-field boundary condition is represented by the exponential behavior of R as it approaches the body in the hodograph plane at w = ±∞

\begin{aligned}R(w, ν) & = & A(ν)e^{B(ν)w}, as w → ±∞ \\ R(w, ν_{∗}) & = & 0 \\ R(w, ν^{∗}) & = & 0 \end{aligned}

where A(ν) and B(ν) are given in Section 2.4. The far-field boundary condition is given implicitly by noting the behavior of (w, ν) near the origin is

\begin{aligned}ν & = & −(γ + 1)\frac{D}{4π} \frac{6KxR}{(x^{2} + 2KR)^{5/2}} \\ w + K & = & (γ + 1)\frac{D}{4π} \frac{2KR − 2x^{2}}{(x^{2} + 2KR)^{5/2}} . \end{aligned}

CHAPTER 3 NUMERICAL PROCEDURE

This chapter details how the boundary value problem summarized in Section 2.5 is solved numerically. A numerical grid is overlaid on the domain of interest, the hodograph plane, and the governing partial differential equation is written as a system of discrete equations. In the region where the equations are elliptic Newton’s method and successive over-relaxation are used. In the hyperbolic region, a timelike marching scheme is used. A method of iteration is used in which the equations representing the boundary value problem are repeatedly solved to obtain a body possessing a shock-free flow, i.e. a numerical solution in the hodograph with J < 0 everywhere.

3.1 Change of variables

It is found that a change of variables is useful when solving the boundary value problem numerically. The change of variables moves the free stream singularity to the origin of the hodograph plane and shifts the sonic line to K/(γ + 1). From now on, the problem shall be described in terms of the new variables

u′ = \frac{w + K}{γ + 1} ,  ν′ = \frac{ν}{γ + 1},  R′ = KR.  (3.1)

Writing the TSDE (2.24) in the new variables and dropping the primes gives the new version of the hodograph PDE to be solved

(\frac{R_{u}}{2R})_{u} + (1 − ϵu)R_{νν} = 0,  where ϵ = \frac{γ + 1}{K}.  (3.2)

A parameter ϵ measuring the the compressibility of the flow is introduced. It is the solution of the incompressible version of (3.2) obtained by setting ϵ = 0 that is used as an initial guess during the iterative solution method of the compressible boundary value problem. There are a number of quantitative differences between the problem formulation in these variables and the earlier version introduced in chapter 2. It is noted that the free stream singularity is now located at u = ν = 0 and the sonic line becomes u = 1/ϵ. In the limit ϵ → 0 (K → ∞) the sonic line moves off to infinity and the problem becomes subsonic everywhere. The solution behavior in this incompressible limit is used to start the numerical iteration procedure to obtain a shock-free solution with ϵ > 0.

The boundary conditions in the new variables are given below. The near-field boundary condition is

\begin{aligned}R(u, ν) & = & Ā(ν)e^{B̄(ν)u}, & as u → ±∞, ν_{∗} < ν < ν^{∗} \\ R(u, ν^{∗}) & = & 0, & |u| < ∞, ν^{∗} = S(x^{∗}) \\ R(u, ν_{∗}) & = & 0, & |u| < ∞, ν_{∗} = S(x_{∗}). \end{aligned}  (3.3)

The functions Ā(ν) and B̄(ν) are dependent on the source distribution function S(x) similar to A(ν) and B(ν) described in section 2.4

Ā(ν) = \frac{1}{2K} exp [−B̄(ν)G′(C̄(ν))]
B̄(ν) = \frac{2}{S′(C̄(ν))}
C̄(ν) = inverse function of S(x).

Notice that the ν′ introduced in (3.1) and now written without the prime means that the factor (γ + 1) will be missing from a number of expressions introduced in Chapter 2 but which also appear in this chapter. The far-field boundary condition for R(u, ν) as (u, ν) → 0 is calculated implicitly from an assumed dipole flow, similar to that shown in section 2.4. The equations representing the dipole flow near the free stream singularity are

\begin{aligned}u & = & \frac{D}{4π} \frac{2R − 2x^{2}}{(x^{2} + 2R)^{5/2}} \\ ν & = & −\frac{D}{4π} \frac{6xR}{(x^{2} + 2R)^{5/2}} \end{aligned}  (3.4)

where

D = π ∫_{−1}^{1} F^{2}(x) dx + ϵπ ∫_{−∞}^{∞} ∫_{ν_{∗}}^{ν^{∗}} |J|u^{2} dνdu,  (3.5)

and

J = −(1 − ϵu)R_{ν}^{2} − \frac{R_{u}^{2}}{2R}.  (3.6)

3.2 Approximating the boundary value problem numerically

The boundary value problem given in the previous section is approximated using finite differences. The infinite hodograph strip is replaced by a truncated computational domain u_{∗} ≤ u ≤ u^{∗}, ν_{∗} ≤ ν ≤ ν^{∗}. The function R(u, ν) is evaluated on a numerical grid, which is illustrated in Figure 3.1. The region near the origin is omitted and dealt with separately using a grid refinement procedure outlined in subsection 3.2.2. The discrete approximation of R(u, ν) is R_{i,j} ≈ R(u_{i}, ν_{j}) where u_{i} and ν_{j} are discrete points in the computational domain given by

\begin{matrix}u_{i} = iΔu, & i = N_{∗} to N^{∗}, \\ ν_{j} = jΔν, & j = M_{∗} to M^{∗}. \end{matrix}  (3.7)

where u_{∗} = N_{∗}Δu, u^{∗} = N^{∗}Δu, ν_{∗} = N_{∗}Δν and ν^{∗} = N^{∗}Δν. For flexibility, the grid spacings Δu and Δν are not necessarily equal.

Figure 3.1: Discretization of hodograph strip

On the numerical grid, finite differences are used to produce the following approximation to the TSDE (3.2):

\frac{1}{Δu^{2}} (\frac{R_{i+1,j} − R_{i,j}}{R_{i+1,j} + R_{i,j}} − \frac{R_{i,j} − R_{i−1,j}}{R_{i,j} + R_{i−1,j}} ) + (1 − ϵu_{i})\frac{(R_{i,j+1} − 2R_{i,j} + R_{i,j−1})}{Δν^{2}} = 0. (3.8)

Figure 3.2: Numerical stencil representing the partial differential equation

This discrete form of the partial differential equation is valid in the interior of the numerical domain, N_{∗} + 1 ≤ i ≤ N^{∗} − 1, and M_{∗} + 1 ≤ j ≤ M^{∗} − 1. The first term on the left hand side of (3.8) results from a difference in “fluxes” {R_{u}/2R} centered around u_{i+1/2} and u_{i−1/2} respectively. For example,

\frac{R_{u}}{2R}|_{u=u_{i+1/2}} = \frac{R_{i+1,j} − R_{i,j}}{R_{i+1,j} + R_{i,j}} .

The 5-point stencil illustrating (3.8) is given in Figure 3.2. To numerically approximate (3.8) at the point marked by a ◯, all five points marked by an × are employed.

The boundary conditions are also discretized. The near-field boundary condition (3.3) is applied on the perimeter of the numerical domain. It is written as a combination of Dirichlet conditions

R_{i,j} = 0,  when j = M_{∗} or j = M^{∗}, for N_{∗} ≤ i ≤ N^{∗},  (3.9)

and Neumann conditions

\frac{2}{Δν} \frac{R_{i+1,j} − R_{i,j}}{R_{i+1,j} + R_{i,j}} = B̄_{j} when i = N_{∗} or i = N^{∗} − 1 for j = M_{∗} ≤ j ≤ M^{∗},

(3.10)

where B̄_{j} = B̄(ν_{j}) = \frac{2}{S′(x_{j})} with ν_{j} = S(x_{j}).

The far-field boundary condition is applied close to the origin in the hodograph strip. The discrete condition on R_{i,j} is that it conforms to the behavior of the flow about a dipole of strength D. Near the origin the numerical grid is refined by the addition of a number of increasingly fine grids in order to better capture the singular behavior of R. The far-field condition is applied on the innermost boundary of the finest grid. More details of the grid refinement structure near the origin are given in subsection 3.2.2. The dipole strength D is computed numerically by evaluating (3.5) on all the grids

D = D_{body} + ϵπ ∑_{grids^{all}} |J_{i+1/2,j+1/2}|u_{i+1/2}^{2} ΔνΔu  (3.11)

where the Jacobian J is computed at cell centers by using

\begin{aligned}J_{i+1/2,j+1/2} & = & \frac{1}{2Δu^{2}} \frac{(R_{i+1,j} + R_{i,j} − R_{i+1,j−1} − R_{i,j−1})^{2}}{(R_{i+1,j} + R_{i+1,j−1} + R_{i,j} + R_{i,j−1})} \\ & & +(1 − ϵu_{i})\frac{(R_{i,j+1} + R_{i,j} − R_{i−1,j} − R_{i−1,j−1})^{2}}{4Δν^{2}} . \end{aligned}  (3.12)

More details of exactly how the far-field boundary condition is dealt with during the numerical solution of the problem are given in the next subsection.

3.2.1 Far-field boundary condition

The far-field boundary condition is evaluated on the perimeter of a rectangular contour enclosing the free stream singularity at the origin. The standard procedure of excising a small region containing a singular point from the rest of the numerical domain is performed. Other methods of handling the singularity, such as “subtracting off the singularity” are not convenient because the equations are non-linear. In practical terms, a Dirichlet condition is evaluated on the perimeter of the cut-out region. For known values of {u^{(k)}, ν^{(k)}} along the rectangular contour {x^{(k)}, R^{(k)}} is obtained by solving the system

u^{(k)} = \frac{D}{4π} \frac{2R^{(k)} − 2x^{(k)}^{2}}{(x^{(k)}^{2} + 2R^{(k)})^{5/2}}
ν^{(k)} = −\frac{D}{4π} \frac{6x^{(k)}R^{(k)}}{(x^{(k)}^{2} + 2R^{(k)})^{5/2}} .

The system is solved using Newton’s method. It is written as finding the root of a vector function f of length 2

f(z) = 0,  where z = \binom{x}{R}  (3.13)

and

f = ( \begin{matrix}\frac{−6xR}{(x^{2} + 2R)^{5/2}} + \frac{4π}{D} u_{i} \\ \frac{2R − 2x^{2}}{(x^{2} + 2R)^{5/2}} − \frac{4π}{D} ν_{j} \end{matrix} )  (3.14)

For notational convenience the superscript (k) is dropped. The Jacobian of the system is calculated and is

J = \frac{∂f(z)}{∂z} = ( \begin{matrix}\frac{12R^{2} − 24x^{2}R}{(x^{2} + 2R)^{7/2}} & \frac{6x^{3} − 18xR}{(x^{2} + 2R)^{7/2}} \\ \frac{6x^{3} − 18xR}{(x^{2} + 2R)^{7/2}} & \frac{12x^{2} − 6R}{(x^{2} + 2R)^{7/2}} \end{matrix} )  (3.15)

The system is solved using the iterative algorithm

z^{n+1} − z^{n} = J^{−1}(z^{n})f(z^{n}),  n = 0, 1, …  (3.16)
z^{0} given.

Due to the presence in (3.15) of terms which have different orders of magnitude depending on position, the method of solution is not trivial. A continuation procedure is used to obtain R and x along the inner boundary contour. The value of R and x previously calculated at some point on the inner boundary usually suffices as an initial value z^{(0)}. On ν = 0 the solution to (3.13) and (3.14) can be found explicitly, so these values are used as starting points for the continuation procedure. In some cases, values at a neighboring point not on the boundary can be used instead. Values of R are systematically calculated on the perimeter of the cut-out region by solving the system above repeatedly and then treated as Dirichlet data in the implementation of the far-field boundary condition.

3.2.2 Grid refinement near free stream singularity at the origin

A uniform computational grid is not well suited to capturing the behavior of R as it changes rapidly near the singularity at the origin. The grid spacings Δu and Δν on the uniform grid are chosen as small as possible given the limits of the computer, but this is still too large to adequately represent the solution there. To rectify this problem it is decided to superimpose a series of increasingly fine grids on the region of the computational domain near the origin. The rationale for this procedure, which is known as grid refinement, is that the presence of more numerical cells close to the origin will be able to increase the accuracy and efficacy of resolving R near the free stream singularity. A diagram of the grid refinement structure near the origin, showing k − 1 grids superimposed on the standard numerical grid, is given in Figure 3.3(a) (k = 3).

At all interior points of all grids the difference equation (3.8) is solved numerically with grid spacing (\frac{Δu}{2^{ℓ−1}} , \frac{Δν}{2^{ℓ−1}} ), for ℓ = 1, 2, …, k. The value of R_{i,j} at points on the outer perimeter of grid ℓ are obtained via interpolation from the next coarser grid ℓ − 1 (ℓ > 1). At the points on the inner perimeter of grid ℓ, R_{i,j} are obtained via interpolation from the next finer grid ℓ + 1 (ℓ < k − 1). On grid 1, R_{i,j} on the outer perimeter are determined by applying the near-field boundary condition (3.3). On grid k, the value of R_{i,j} on the inner perimeter are obtained by applying the far-field boundary condition and solving the system (3.14) using (3.16).

(a) Figure showing refinement near origin

(b) Details of grid at origin

Figure 3.3: Grid geometry near origin

For example, in Figure 3.3(b) there are two superimposed grids on the main grid, thus k = 3. Grids 2 and 3 have 12 by 12 grid points each. The refinement grids are assumed to be symmetric about the u-axis, so only the refinement region ν ≥ 0 in Figure 3.3(b) will be discussed. The points marked by ◯ represent points on the outer perimeter of grid 2 where data is extrapolated from the coarser grid 1. The □ points represent points on grid 2 which are extrapolated from the finer grid 3. The four innermost points on grid 3 marked by × represent the points where values are determined by solving the implicit system (3.4). There is no extrapolation along the dotted line u < 0, ν = 0 because the branch cut in the hodograph is situated there.

The chosen method of interpolation is a third-order tensor product interpolation. The interpolation method must be greater than second-order accurate in order to maintain the consistency of the numerical method used to solve the secondorder PDE. Since the interpolation chosen is third-order accurate, the numerical discretization of the solution R_{i,j} will be only first-order accurate on the boundary, but the overall finite difference scheme is still second-order accurate. As can be seen in Figure 3.3, the number of points involved in the interpolation is very small compared to the number of points in the entire numerical domain. Near a point (u_{0}, ν_{0}) on the perimeter, the form of the interpolation is

R(u, ν) = ∑_{i=−1}^{1} ∑_{j=−1}^{1} R_{i,j}p_{i}(u)q_{j}(ν)  (3.17)

where p_{i}(u) and q_{j}(ν) are Lagrange interpolating polynomials

p_{i}(u) = ∏_{k≠i^{k=−1}}^{1} \frac{u − u_{i}}{u_{k} − u_{i}}  and  q_{j}(ν) = ∏_{k≠j^{k=−1}}^{1} \frac{ν − ν_{j}}{ν_{k} − ν_{j}} .  (3.18)

Figure 3.4 shows the interpolation stencil and illustrates how the interpolation between coarse and fine grids proceeds. The interpolation method uses a local coordinate system with R(u, ν) as the origin and involves the closest nine points within one grid spacing.

Figure 3.4: Interpolation method between coarse and fine grids For example, in Figure 3.4 the interpolated value at the point marked by an

× can be calculated using the expression

R(u, \frac{Δν}{2} ) = ∑_{i=−1}^{1} ∑_{j=−1}^{1} R_{i,j}p_{i}(0)q_{j}(\frac{Δν}{2} )  (3.19)
= R_{−1,−1}p_{−1}(0)q_{−1}(\frac{Δν}{2} ) + R_{−1,0}p_{−1}(0)q_{0}(\frac{Δν}{2} ) + R_{−1,1}p_{−1}(0)q_{1}(\frac{Δν}{2} ) +

R_{0,−1}p_{0}(0)q_{−1}(\frac{Δν}{2} ) + R_{0,0}p_{0}(0)q_{0}(\frac{Δν}{2} ) + R_{0,1}p_{0}(0)q_{1}(\frac{Δν}{2} ) +

R_{1,−1}p_{1}(0)q_{−1}(\frac{Δν}{2} ) + R_{1,0}p_{1}(0)q_{0}(\frac{Δν}{2} ) + R_{1,1}p_{1}(0)q_{1}(\frac{Δν}{2} )  (3.20)

where the interpolating polynomials are

p_{−1}(u) = \frac{u(u − Δu)}{2Δu^{2}},  q_{−1}(ν) = \frac{ν(ν − Δν)}{2Δν^{2}},
p_{0}(u) = \frac{(u^{2} − (Δu)^{2})}{Δu^{2}}, q_{0}(ν) = \frac{(ν^{2} − (Δν)^{2})}{Δν^{2}},
p_{1}(u) = \frac{u(u + Δu)}{2Δu^{2}},  q_{1}(ν) = \frac{ν(ν + Δν)}{2Δν^{2}}.

(3.21)

Evaluating p_{i}(u) at u = 0 means that it will only have a non-zero value when i = 1, as can be seen from (3.21). The expression for R(u, \frac{Δν}{2} ) is simplified from (3.20) to

R(u, \frac{Δν}{2} ) = \frac{3}{8}R_{0,−1} + \frac{3}{4}R_{0,0} − \frac{1}{8}R_{0,1}.  (3.22)

The efficacy of the grid refinement procedure near the origin can be evaluated. The incompressible boundary value problem (ϵ = 0 in 3.2) is solved repeatedly, each time with successively more grid refinement at the origin. That is, the problem is solved with numerical grids which have an increasing number of overlapping grids and then the behavior of R(u, ν) as (u, ν) → 0 is examined. The asymptotic behavior of R(u, 0) along ν = 0 as u → 0 and along u = 0 as ν → 0 can be found explicitly from equation (3.4). This produces the expressions

R(u, 0) = \frac{1}{2}(\frac{D}{4πu})^{(2/3)},  u → 0  (3.23)

and

R(0, ν) = (3)^{(−3/2)} \frac{D}{2πν} , ν → 0. (3.24)

Figure 3.5 is a plot showing how the data approaches the limiting function as more grids are superimposed upon the region close to the origin. The data obtained from each grid is denoted by a different symbol. Since the refinement grids overlap, there is often more than one symbol at a particular value of ν. For example, at ν ≈ 0.04 there are four symbols. What this shows is that data obtained from a grid with five overlapping grids ( denoted by the symbol △) is more accurate than data obtained from a grid with four overlapping grids (denoted by the symbol ▽), which in turn is more accurate than the initial data obtained from the grid with three overlapping grids (denoted by ◯). The closer the data points are to the curve representing the explicit asymptotic solution given in (3.24), the more accurate they are. The data from the solution with just two overlapping grids (denoted by ×) is exactly on the curve because for that grid arrangement that is the closest grid point to the origin,

Figure 3.5: Plot showing R(0, ν) on successively-finer grids

and the value there was set explicitly to agree with the asymptotic expression in (3.24).

In the following section, details of how the boundary value problem which has been transformed into a set of discrete equation is solved numerically will be discussed.

3.3 Solving the hodograph boundary volume problem numerically

The boundary value problem for R(u, ν) has been restated as a series of discrete equations which need to be solved numerically. Due to the mixed-type elliptichyperbolic nature of the governing partial differential equation (3.2), the solution method is dependent on whether the PDE is elliptic or hyperbolic at any particular point. To avoid this complication, it is decided to solve the problem in the elliptic region first and then to use this solution as an initial condition to compute the solution in the hyperbolic region. The solution method in the hyperbolic region involves an explicit time-like marching scheme. Details will be given in subsection 3.3.2. The solution method in the elliptic region involves solving a slightly different boundary value problem then that described in section 3.2. The entire problem has been reduced to this: given an S′(x), compute R everywhere so that J is always negative. The details of solving the elliptic problem will be given below in subsection 3.3.1.

3.3.1 Solving in the elliptic region

The line of separation between the elliptic and hyperbolic regions corresponds to the boundary of the sonic bubble in the physical plane and is referred to as the sonic line. In the hodograph plane, the sonic line is the rightmost boundary of the elliptic region and exists along u = u^{∗} = N^{∗}Δν = 1/ϵ. Given provisional data for R(u^{∗}, ν) along the sonic line, the boundary value problem given in section 3.2 is solved using this new Dirichlet data along the new boundary. The idea is to solve for R(u, ν) in the elliptic region first and then later solve the problem in the hyperbolic region. The sonic line data is represented by a curve with free parameters to allow flexibility in selecting the shape. The chosen function is

Q(z) = P (z) exp [\frac{σz^{2}}{1 − z^{2}} ] ,  |z| ≤ 1  (3.25)

where

z(ν) = \frac{2ν − ν_{∗} − ν^{∗}}{ν_{∗} − ν^{∗}} , and P (z) = a_{0} + a_{1}z + a_{2}z^{2} + …, and σ < 0.

The functional form of P (z; {a_{i}}) is influenced by Ā(ν), since the near-field behavior at u = u^{∗} according to (3.3) is R = Āe^{B̄u^{∗}} ≈ Ā(ν). The initial data used is from the corresponding incompressible solution along u = u^{∗}. Then the function P (z) is chosen to conform to this data, but with parameters that can be perturbed to give the flexibility needed to alter the sonic line data when looking for shock-free solutions. The new boundary condition along the far right edge of the numerical domain is

R_{i,j} = Q(z(ν_{j})), where i = N^{∗}, j = M_{∗} to M^{∗}. (3.26)

Once the problem has been discretized on the elliptic numerical domain and equations have been identified for R_{i,j} at every single point the question of what method to use to solve the equations arises. Newton’s method is chosen. This is the same method that was used to solve the implicit equations of the far-field boundary conditions. In that situation the function f had two components and was dependent on two variables. In the case of the solution R_{i,j} in the hodograph plane there are (N^{∗} + N_{∗} + 1)(M_{∗} + M^{∗} + 1) components of the function F(Z), with

Z_{k} = R_{i,j}, where k = (i + N_{∗})(M_{∗} + M^{∗} − 1) + j + M_{∗},

with −N_{∗} + 1 ≤ i ≤ N^{∗} − 1 and M_{∗} ≤ j ≤ M^{∗} − 1. It is noted that the problem of solving the boundary value problem numerically for R_{i,j} has been reduced to the solution of F(Z) = 0 using Newton’s method. It is found that the solution method is very sensitive to the initial condition, F(Z^{(0)}), which makes it difficult to use this method to produce a converged solution. Successive over-relaxation (SOR) is first performed on R_{i,j} throughout the domain to partially solve the problem and then Newton’s method is applied. The relaxation steps are conducted by sweeping from i = N^{∗} − 1 to N_{∗} + 1 and calculating the residual (i.e. the difference from zero) involved in computing the left hand side of (3.8), the finite difference approximation to the governing PDE. The relaxation is repeated until the maximum value of the residual on grid 1 is below a given tolerance. At that point the relaxation procedure is conducted on the next finer grid, until the residual on this new grid falls below the desired tolerance. This procedure is complicated by the interpolation that occurs between the overlapping section of the two grids. Whenever interpolation occurs as one moves from grid to grid the points involved have a larger error than surrounding points, which changes the maximum value of the residual for the entire grid.

Remember, the points involved in the interpolation end up being only first-order accurate. Thus, the relaxation procedure will oscillate between the two grids until the residuals on both grids fall below a certain tolerance. At that point another grid is added and the process will repeat. This is known as grid relaxation. SOR is being conducted on a grid by grid basis. After the residual has been driven down below tolerance on grid k, the method will bootstrap its way up to grid 1 by performing SOR on the intermediary grids k − 1 through 2 successively. Once the method is back on grid 1 it drives the residual below tolerance there and starts to move down the series of grids again. This process is repeated until all the grids consecutively give a residual error below the desired tolerance, at which time the residuals are so small that the values of R_{i,j} are very close to the solution values. Using these values of R_{i,j} to set Z^{(0)}, Newton’s method is then used to obtain R_{i,j} on all the grids at once. It takes no more than five Newton iterations at this point to bring the maximum residual on all grids to below 10^{−8}. Thus R has now been computed everywhere in the elliptic region N_{∗} ≤ i ≤ N^{∗∗}, M_{∗} ≤ j ≤ M^{∗} to within tolerance.

3.3.2 Solving in the hyperbolic region

After R has been computed everywhere in the elliptic region it is computed in the hyperbolic region using an explicit time-like marching scheme. The difference scheme (3.8) can be rearranged to produce

R_{i+1,j} = \frac{1 + T_{i,j}}{1 − T_{i,j}} R_{i,j},  (3.27)

where

T_{i,j} = \frac{R_{i,j} − R_{i−1,j}}{R_{i,j} + R_{i−1,j}} + (ϵu_{i} − 1)(R_{i,j+1} − 2R_{i,j} + R_{i,j−1}).  (3.28)

Starting from the sonic line at u = 1/ϵ and using R_{i,j} and R_{i−1,j} from the previously computed solution in the elliptic region, the two-step scheme (3.27) and (3.28) is utilized to obtain R_{i+1,j} in the hyperbolic region. The marching scheme is used to compute R_{i,j} until i = N^{∗}. Using (3.12), the jacobian is computed along with R_{i,j} in the hyperbolic region. In order for the hodograph transformation to remain valid the jacobian must be negative everywhere. This criterion must be met in order for a shock-free solution to be calculated. If the jacobian is determined to be positive at some point then the sonic line data is perturbed slightly by changing the parameters in the P (z) function and the boundary value problem is solved again in the elliptic region. Using this new solution in the elliptic region, the solution is extended into the hyperbolic region by applying (3.27) and (3.28). The jacobian is re-calculated. This procedure is repeated until a P (z; {a_{i}}) function for the sonic line data is found which after the elliptic boundary value problem is solved produces a solution with a jacobian that remains negative wherever it is calculated in the hyperbolic region.

Once the numerical solution is known everywhere in both the hyperbolic and elliptic regions, and it possesses a jacobian which is negative everywhere, the strength D of the dipole modeling the far-field flow is computed. The integrals in (3.11) are evaluated numerically to determine whether the computed dipole strength is equivalent to the inputted value used in computing the far-field boundary conditions. If the two numbers are different a type of fixed-point iteration is performed, by which the computed value of D is used as the correct dipole strength and the solution in the entire hodograph plane is re-computed. Since the dipole affects the flow mainly in the far-field (near the free stream singularity at the origin in the hodograph) the solution as a whole is relatively insensitive to changes in dipole strength and the method converges to a single value of D after a few iterations.

3.4 Computing the body of revolution from the numerical solution

Now all aspects of the problem in the hodograph have been computed. The goal of this thesis is to discover transonic slender bodies of revolution that are shockfree in the physical plane. Thus R(u, ν) must be used to compute x(u, ν) using (2.23) and the two hodograph dependent variables must be inverted to obtain information about functions of interest in the physical plane, scilicet F (x) and S(x).

If x is defined at the center of the numerical cells of the grid then the first order system of PDEs (2.23), written in new variables, can be approximated by finite differences

\begin{aligned}x_{i+1/2,j+1/2} − x_{i−1/2,j+1/2} & = & (ϵu_{i} − 1)(R_{i,j+1} − R_{i,j}) \\ x_{i+1/2,j+1/2} − x_{i−1/2,j−1/2} & = & \frac{(R_{i,j+1} − R_{i,j})}{(R_{i,j+1} + R_{i,j})}. \end{aligned}  (3.29)

Using (3.29), one value of x must be given as an initial point to start the calculation. The chosen value was x_{N_{∗}+1/2,1/2} = −1, representing the nose of the body of revolution. In order to recover F (x) and S(x) from the hodograph data it is recalled that information about the body went into the formula for the near-field boundary condition on the perimeter of the hodograph strip. In fact, at large values of |u_{i}|, ν_{j+1/2} ≈ S(x_{i+1/2,j+1/2}). This is the discrete version of part of the near-field boundary condition (3.3). Once S(x) is known F (x) can be competed using numerical integration of

F^{2}(x) = 2 ∫_{−1}^{x} S(x̄) dx̄.  (3.30)

Information about the body at infinity in the hodograph comes from the formula for the asymptotic behavior of R in the hodograph (3.3), which is written below in the new scaled variables

R ∼ \frac{K}{2} exp [\frac{2}{S′(x)}(u − G′(x))] , as u → ±∞.

By applying this formula on the two far edges of the numerical grid representing |u_{i}| large (i = N_{∗} or i = N^{∗} − 1), expressions for discrete values S′, G′ and F′ can be written

S′_{j} = \frac{2Δν}{log(R_{i+1,j}/R_{i,j})}, G′_{j} = \frac{S′_{j}}{2} log [\frac{2R_{i,j}}{K} (\frac{R_{i,j}}{R_{i+1,j}} )^{i}] , F′_{j} = \frac{S′_{j}}{F_{j}} .  (3.31)

The values of F′_{j}, S′_{j} and G′_{j} on the body are all needed to compute the pressure coefficient C_{p} on the body in accordance with (2.12).

The numerical procedure by which the boundary value problem formulated in Chapter 2 is scaled slightly and then converted to a system of discrete equations has been given in this chapter. The details of how multiple grids are used to represent the numerical domain, the method of solving the discrete equations and how to generate the necessary data to recover a shock-free body of revolution are also given. An algorithm, or list of actions executed in sequence during the numerical calculations, is given in section 3.5.

3.5 Algorithm used in numerical solution of the hodograph BVP

1 Choose body function F (x)

2 Calculate S(x) and S′(x)

3 Find roots of S′(x): x^{∗} and x_{∗}

4 Discretize the hodograph strip as shown in Figure 3.1

5 Set on-the-body boundary conditions at u = ±∞

6 Choose P (z; {a_{i}}) function to replace boundary condition at u = +∞ with given sonic line data at u = 1/ϵ

7 Assume a value for D, the dipole strength: D_{new}

8 Set far-field boundary conditions near the free stream singularity at the origin by solving the given implicit equations describing the solution behavior

9 Choose an initial condition for R everywhere in the elliptic region

10 Solve the boundary value problem numerically for R in the elliptic region 11 Use time-marching scheme to extend solution into hyperbolic region and check

jacobian is negative everywhere its computed

12 If the computed jacobian is greater than zero at any point, go to item 6 and iterate

13 Integrate to obtain a value for D_{num} from the numerically generated flow field

14 If |D_{num} − D_{new}| greater than tolerance, go to item 7 and iterate

15 Convert variables back into the physical plane: compute F(x) and S(x) from numerical solution. Also G′, S′, F′ are all needed to compute C_{p_{body}}

In the next chapter, a survey of the results obtained using this algorithm are given.

CHAPTER 4 NUMERICAL RESULTS

In this chapter, the results of applying the numerical procedures outlined in the previous chapter to calculate specific shock-free bodies will be given. The exact details of a number of computed shock-free transonic slender bodies will be given.

Each particular shock-free body can be identified uniquely by a number of different quantities. For a shock-free solution to have been computed there are a number of functions and parameters whose values were known explicitly in order to initiate the search: the equation of the body of revolution F (x), the dipole strength D, the value of the transonic similarity parameter K and the functional form of the sonic line data P (z; {a_{i}}). Knowing the values of these parameters and repeating the numerical procedures outlined in Chapter 3 will, according to TSDT, lead to computation of bodies possessing flows that are shock free.

Only one of the four pieces of data necessary to compute and identify a particular shock-free solution of the boundary value problem remains invariant after being inputted initially into the numerical method: the value of K. The other variables are altered during the iteration process described in section 3.3. The data representing the body shape will not likely be identical to the inputted function F (x), though it will probably be similar. Iteration on D and P (z; {a_{i}}) are vital parts of the algorithm listed in Section 3.5. The reason why K is invariant is because the transonic similarity parameter appears directly in the governing PDE through the value of ϵ in equation (3.2). This value is linked to the decision of where the sonic line will be in the hodograph plane (u = u^{∗}). Thus the value of K is a fixed feature of the boundary value problem.

In practice the choice of where the sonic line will be in the hodograph is made by setting u^{∗} = 2ν^{∗}. This fixes the value of K for the problem being solved. Since u^{∗} = 1/ϵ, an expression for K in terms of the dimensions of the hodograph strip can be written

K = (γ + 1)u^{∗} = 2(γ + 1)ν^{∗},  where ν^{∗} = S(x^{∗}).  (4.1)

So, K depends on the shape of the particular body being solved, since the value of ν^{∗} is determined by the roots of S′(x; α). How the form of the initial body function is chosen is given in the next section.

4.1 Modeling the body of revolution

The function representing the input shape of the slender body of revolution is chosen so that the body can be either fore-aft symmetric or fore-aft asymmetric, depending on the value of a parameter α, referred to here as the body symmetry parameter. A body is said to be fore-aft symmetric if its functional representation is an even function, i.e. F (x; α) = F (−x; α). The functional representation of F (x; α) is chosen to be

F (x; α) = (1 − x^{2})(1 + αx) = 1 + αx − x^{2} − αx^{3},  |α| small, |x| ≤ 1. (4.2)

The parameter α is taken to be small because the initial idea was to try and obtain solutions that are perturbations of the first computed shock-free body given by Cole and Schwendeman in [8]. That body was a fore-aft symmetric body, which corresponds to α = 0. Once a form of F (x) is known the source distribution function S(x) and its derivative S′(x) can be obtained. (Recall that S(x) = F (x)F′(x) from (2.11).) Explicit expressions for S(x) and S′(x) are

S(x; α) = α + (α^{2} − 2)x − 6αx^{2} + 2(1 − 2α^{2})x^{3} + 5αx^{4} + 3α^{2}x^{5} (4.3)
S′(x; α) = α^{2} − 2 − 12αx + 6(1 − 2α^{2})x^{2} + 20αx^{3} + 15α^{2}x^{4}  (4.4)

In the case of a fore-aft symmetric body, α = 0, which implies that x_{∗} = −x^{∗} and ν^{∗} = −ν_{∗}. The hodograph strip will be symmetric about the u-axis (ν = 0). This solution simplifies the numerical calculation recounted in Chapter 3 because the solution of the entire problem can be obtained with computations performed solely in the upper half-plane. Using the methods described in this thesis both fore-aft symmetric and fore-aft asymmetric bodies can be computed.

In practice, it is decided to follow the convention of discretizing the hodograph strip by choosing the grid spacings Δν = Δu = h, so that each numerical cell is square as opposed to being rectangular. For convenience, the hodograph strip has been discretized so that grid lines fall at ν = ν_{∗}, ν = 0 and ν = ν^{∗}. There is a branch cut all along the negative half of the u-axis, so to simplify the computational procedure a grid line is desired along ν = 0. Two Dirichlet boundary conditions occur along the extreme upper and lower edges of the hodograph strip, namely along ν = ν_{∗} and along ν = ν^{∗}, respectively, so grid lines are desired there also. This results in a side problem to be solved in order to find α such that S(x^{∗}; α)/S(x_{∗}; α) is a rational number, where x^{∗} and x_{∗} are determined by

\begin{aligned}S′(x^{∗}; α) & = & 0 \\ S′(x_{∗}; α) & = & 0. \end{aligned}  (4.5)

The rational number will depend on the desired ratio of grid lines. The side problem can be summarized as

\{ \begin{matrix}S′(x^{∗}; α) & = & 0 \\ S′(x_{∗}; α) & = & 0 \\ S(x^{∗}; α) + βS(x_{∗}; α) & = & 0 \end{matrix} \} .  (4.6)

A grid symmetry parameter, β = \frac{M^{∗}}{M_{∗}}, has been introduced that represents the chosen ratio of grid lines in the upper and lower half-planes. The system (4.6) can be solved simply, using a version of the bisection method. If β = 1, then the number of grid lines in the upper and lower half-planes is the same, so the strip is symmetric. The third equation in the system then implies that ν^{∗} = S(x^{∗}; α) = −S(x_{∗}; α) = −ν_{∗}. As shown previously, that corresponds to the fore-aft symmetric problem, where α = 0 which implies x^{∗} = −\frac{1}{\sqrt{3}} and x_{∗} = \frac{1}{\sqrt{3}}. The asymmetric values (α ≠ 0) will be perturbations of these numbers.

When discretizing the hodograph domain, the procedure is to select the ratio of asymmetry in the hodograph region, which implies a value of β. For example, a typical choice is 40 grid lines in the positive ν half-plane compared to 36 in the negative ν half-plane. This corresponds to a value of β = \frac{40}{36} = \frac{20}{18}. The system (4.6)

Figure 4.1: α versus β

is solved and the values of x^{∗}, x_{∗} and α are computed for this particular slender body of revolution. These values set the dimensions ν^{∗} and ν_{∗} of the hodograph strip to be used in the solution of the problem. A graph of the body symmetry factor versus the grid symmetry factor, α versus β, can be seen in Figure 4.1. This figure displays the range of possible values for the two symmetry parameters. A table of choices of β with associated values of α, ν_{∗} and ν^{∗} is listed in Table 4.1. These choices correspond to actual input parameters used to compute shock-free solutions employing the methods described in this thesis. Further details about the computed shock-free solutions are give in the next section.

\begin{matrix}β & α & ν_{∗} & ν^{∗} \\ 20/17 & -0.14560 & -0.72073 & 0.84792 \\ 20/18 & -0.09251 & -0.73484 & 0.81648 \\ 20/19 & -0.04456 & -0.75139 & 0.79093 \\ 20/20 & 0.00000 & -0.76980 & 0.76980 \\ 18/20 & 0.09251 & -0.81648 & 0.73484 \end{matrix}

Table 4.1: Values of β, α, ν_{∗} = S(x_{∗}; α) and ν^{∗} = S(x^{∗}; α)

4.2 Computed Shock-free Bodies

In this section, detailed parameters will be given that are associated with actual computed shock-free bodies as well as plots of the bodies themselves. Some of these parameters are shown in Table 4.2 and Table 4.3. Table 4.2 contains the parameters β, which represents the degree of asymmetry of the solution and K, which represents the speed of the computed shock-free flow. The design Mach number M_{∞} for a body of standard 10% thickness is also given.

Recall from Chapter 3 that after the parameters in Table 4.1 and Table 4.2 are known, they are used to set up the numerical grid. Then the sonic line data is adjusted and the equations solved until a shock-free solution results. Afterwards a fixed-point iteration on the dipole strength D is employed. In Table 4.3 the value

Name  β K = (1 − M_{∞}^{2})/(δM_{∞})^{2} δ  M_{∞}
SFB-1 20/17  4.06999  0.1 0.98025
SFB-2 20/18  3.91912  0.1 0.98096
SFB-3 20/19  3.79648  0.1 0.98154
SFB-4 18/20  3.52721  0.1 0.98282

Table 4.2: Values of α, β, K and M_{∞} for shock-free solutions

of the dipole strength and the value of the parameters a_{0}, a_{1}, a_{2} in P (z) describing the equation of the given sonic line are given for each computed shock-free body.

Name  K  D  {a_{0}, a_{1}, a_{2}}
SFB-1 4.06999 3.43760 .023, .000, .000
SFB-2 3.91912 3.23787 .011, -.002, .013
SFB-3 3.79648 3.53389 .030, -.001, .005
SFB-4 3.52721 3.15588 .016, .000, .000

Table 4.3: Values of K, D and a_{0}, a_{1} and a_{2} for shock-free solutions

Knowledge of the parameters in these two tables is enough to describe a particular shock-free solution. In the next few pages, relevant data for the four transonic shock-free slender bodies of revolution, SFB-1 through SFB-4 will be plotted.

The reason why contours of the horizontal velocity are plotted is that they indicate lines of constant pressure, i.e. they are approximate isobars. The importance of the region near the points x^{∗} and x_{∗} is dramatically illustrated. All isobars go through these points, which means that the pressure takes on all values there. Generating the contours is an easy task. In the hodograph plane a vertical line at u = u_{const} with the values R(u_{const}, ν) and x(u_{const}, ν) becomes a contour of u = u_{const} in the physical plane. In the plots of the approximate isobars of u = ϕ_{x} = u_{const}, the dashed contours indicate subsonic flow when u_{const} < K/(γ + 1), and solid contours indicate supersonic flow when u_{const} > K/(γ +1). The sonic bubble is represented by a dotted line between the two different types of lines. Note the substantial size of the sonic bubble possessed by each of the shock-free bodies. A close up of the contours confirms the fact that the contours do not intersect. The appearance of “holes” in the isobar plots near to the x_{∗} and x^{∗} points is merely a visual effect due to the fact that the last data point is not plotted. The sharp change in the curves of C_{p} versus x near the point x_{∗} and x^{∗} is not considered to be a shock, only a region of rapid change, because the computed flow in the hodograph plane has J < 0 everywhere, and thus the flow must be smooth by definition. The convention of plotting values of the negative pressure coefficient is followed throughout this chapter.

A difficult feature to discern from the plot of each particular body is that they

Name Δ F |_{x=1}
SFB-1 0.02423
SFB-2 0.04758
SFB-3 0.02457
SFB-4 0.01414

Table 4.4: Closure deviation values for computed shock-free bodies

are not closed. In other words, the value at the tail is not exactly zero, which is different from the value at the nose, which is zero. The deviations are small (between 1% and 5%), and are tabulated in Table 4.4. It is worthwhile to mention that the computed airfoils of Bauer, Garabedian and Korn [1] have this property also. The shock-free body of Cole and Schwendeman [8] is closed by definition, because it is a fore-aft symmetric body. The closure deviation in the body data at the tail is most visible on SFB-2, which also explain the slight kink in the C_{p} curve at x = 1.

(a) Isobars around SFB-1 (b) Pressure Coefficient on SFB-1

Figure 4.2: Data for shock-free body number 1, SFB-1

(a) Isobars around SFB-2 (b) Pressure Coefficient on SFB-2

Figure 4.3: Data for shock-free body number 2, SFB-2

(a) Isobars around SFB-3 (b) Pressure Coefficient on SFB-3

Figure 4.4: Data for shock-free body number 3, SFB-3

(a) Isobars around SFB-4 (b) Pressure Coefficient on SFB-4

Figure 4.5: Data for shock-free body number 4, SFB-4 4.3 Detailed features of shock-free solutions

In this section more details of the shock-free solutions will be given. Examples of four computed shock-free bodies, together with isobar contours, have been shown in the previous pages. In Figure 4.6 a surface plot with contours illustrates an actual computed solution R(u, ν) to the hodograph version of the transonic smalldisturbance equation given in Chapter 3. Notice the singular behaviour at the origin and the exceedingly rapid decay as one moves away from that point. In fact Figure 4.6 indicates that the solution R(u, ν) decreases exponentially as u → +∞. This behavior is in accordance with the near-field boundary condition (3.3) given in Chapter 3. So, despite the fact that this boundary condition was replaced by the

Figure 4.6: Surface plot of R(u, ν) from numerical hodograph solution Dirichlet data along the sonic line at u_{∗} in practice, it is observed that the computed solution satisfies the exponential decay condition as desired.

Other features of the computed solution (this one corresponds to SFB-3) can be observed by looking at the contours in Figure 4.6. For example, along ν = 0 for u < 0 the solid line indicates the branch cut which was described in the mathematical theory given in Chapter 2. The shape of the contours illustrates the fact the far-field resembles that of a dipole.

Figure 4.7: Sonic line data used to compute the SFB-1 solution

In Figure 4.7 the dashed lines are examples of alternate sonic line data that were used to try and generate the shock-free solution corresponding to SFB-1 are shown. The solid line is the actual sonic line data used to generate SFB-1. The initial curve comes from the solution to the incompressible version of the transonic smalldisturbance equation. As described in Chapter 3, the algorithm involves setting the sonic line data, solving the entire hodograph boundary value problem and checking the jacobian. If the jacobian is non-negative somewhere, the sonic line data is altered and the procedure repeated until J < 0 in the entire hyperbolic region. The precise sonic line data used to generate the shock-free solution corresponding to SFB-1 is given by the function

R(u = u^{∗}, ν) = 0.023 exp [\frac{−z^{2}}{1 − z^{2}} ] , where z(ν) = \frac{2ν − ν_{∗} − ν^{∗}}{ν_{∗} − ν^{∗}} .

The determination that a shock-free solution has been computed is derived from the behavior of the jacobian of the hodograph transformation. The jacobian must be negative everywhere it is computed. How the jacobian changes as the sonic line data is massaged to produce a shock-free solution is shown in Figure 4.8. The importance of the extrema of the source distribution function, which in the hodograph are mapped to the singular lines ν = ν_{∗} and ν = ν^{∗}, is demonstrated. It is a very subtle effect. Along these singular boundaries of the hodograph plane the Jacobian appears to “bunch up” and become positive. As soon as that happens the solution is no longer single-valued and information about the hodograph solution can not be translated into happenings in the physical plane. In other words, it no longer corresponds to a solution which is smooth. We refer to these non-smooth solutions as shocked. The “bunching” phenomenon can be observed by comparing the figures in Figure 4.8. By looking at the surface plots and the contours below one can see that Figure 4.8(a) and Figure 4.8(b) have a ridge along the left boundary, which in the contour plots correspond to a concentration of contour lines. In Figure 4.8(c), which is displaying the jacobian of the shock-free solution, there is no such ridge and the contour plot is noticeably different from the one in Figure 4.8(a), which displays the jacobian for a shocked solution.

There are other indicators of a shock-free solution apart from looking at the jacobian. The calculation of the jacobian is the most important one, because if it is non-negative then computations in the hodograph are meaningless, because the solution can no longer be translated to the physical plane.

(a) J(u, ν) for shocked solution  (b) J(u, ν) for lightly-shocked solution  (c) J(u, ν) for shock-free solution

Figure 4.8: Surface plots displaying effect of shocklessness on the jacobian

Another indicator of a shock-free solution is that like characteristics do not cross each other. The characteristics of the governing partial differential equation, the hodograph version of the transonic small-disturbance equation (3.2), can be computed and plotted in the hodograph plane. This is done be re-writing (3.2) as

R_{uu} + 2R(1 − ϵu)R_{νν} − \frac{R_{u}^{2}}{R} = 0  (4.7)

so that the characteristics are given by

\frac{dν}{du} = ±\sqrt{2R(ϵu − 1)}.  (4.8)

The equation for the characteristics in the hodograph (4.8) is integrated numerically using the trapezoidal rule and plotted. Knowing R everywhere on the grid, the right-hand side of (4.8) becomes a function of (u, ν), using a suitable interpolation scheme. This implicit equation is then solved using the bisection method to obtain a set of (u, ν) values representing the hodograph characteristic curves. The results of doing this for the SFB-4 body are shown in Figure 4.9(a). One can calculate the physical plane characteristic curves by using the (u, ν) values representing the hodograph characteristics to calculate corresponding (x, R) values and use these values to plot the characteristics in the physical plane. In Figure 4.9(b) the physical plane characteristics are the image of the hodograph characteristics plotted in Figure 4.9(a). In both figures it is clear that characteristics of a similar type do not cross each other. The only exception is at x^{∗} and x_{∗} in Figure 4.9(b) where like characteristics appear to be very close to one another. An extreme closeup reveals that they do not, in fact, intersect. However, once again these plots illustrate the importance of the extrema points.

Unfortunately, the hodograph data near the extrema points is difficult to obtain accurately enough to definitively resolve the details of shock formation in the physical plane. This can be seen by looking at the plots of F (x), S(x), S′(x)

(a) Hodograph plane characteristics (b) Physical plane characteristics

Figure 4.9: Characteristic curves associated with SFB-3 solution and G′(x) associated with the hodograph-designed body SFB-2 in Figure 4.10, Figure 4.11, Figure 4.12 and Figure 4.13 respectively. Since the expression for C_{p_{body}} involves the functions plotted in the above figures, this is the reason for the discontinuities seen near the x^{∗} and x_{∗} points for all the plots of C_{p} given for SFB-1, SFB-2, SFB-3 and SFB-4 earlier in the chapter. However, there are some features which are worth commenting on. Notice the flat nature of the pressure coefficient plot between the extrema points and the sudden drop near both the nose and the tail of the body. From the placement of the actual data points obtained from the hodograph it is easy to see why no definitive statement about the shock-free nature of the flow can be made. The data is too coarse to resolve the question of whether a jump in C_{p} is occurring at x = x^{∗}. It is clear that even if a shock is present, it is not a very substantial one.

The claim that shock-free solutions have indeed been computed is supported by all the figures displayed in this chapter. It is supported by the lack of intersections in the plots of the characteristics in Figure 4.9(a) and Figure 4.9(b), as well as the degree of smoothness in the jacobian surface plot of Figure 4.8(c) and all the isobar plots.

The computations were performed for the most part on an IBM 3090-200S though at times a Sun SparcStation-2 and later, a Sun SparcStation-10 were employed. Typical values of N_{∗}, N^{∗}, M_{∗} and M^{∗} were N_{∗} = 20, N^{∗} = 18, M_{∗} = 180 and M^{∗} = 180, resulting in a grid with about 16 000 data points. On a Sparc-10 the maximum residual is reduced to below tolerance (1.0x10^{−8}) in around 4 hours.

So far in this chapter the results and significant details of numerical solution of the compressible boundary value problem were given to show that shock-free flows can be computed in the hodograph plane. In the next chapter, information will be given about solving the problem in the physical plane and verifying that the hodograph-designed bodies also possess shock-free flows in the physical plane.

Figure 4.10: Body function F (x) data for the SFB-2 Figure 4.11: Source distribution S(x) data for the SFB-2

Figure 4.12: S′(x) data for the SFB-2 Figure 4.13: Plot of G′(x) data for the SFB-2

CHAPTER 5 PHYSICAL PLANE CALCULATIONS

Another means of determining shock-free transonic slender bodies of revolution is to directly solve the boundary value problem for ϕ(x, r̃), the perturbation potential in the physical plane. The ability to compute the flow around a given slender body of revolution at a given speed is useful in confirming that shock-free flow exists around the hodograph-designed bodies of Chapter 4. More importantly, “off-design” computations of the flow around these bodies can be conducted.

The boundary value problem in the physical plane for the velocity perturbation potential ϕ(x, r̃) is solved numerically. The salient details of the physical plane boundary value problem from Chapter 2 are summarized below. The governing partial differential equation to be solved (2.6) is

(K − (γ + 1)ϕ_{x})ϕ_{xx} + ϕ_{r̃r̃} + \frac{1}{r̃}ϕr̃ = 0,  |x| < ∞, r̃ > 0.  (5.1)

The near-field boundary condition (2.13) is

ϕ(x, r̃) → S(x) log r̃ + G(x), as r̃ → 0, |x| ≤ 1
ϕ(x, r̃) bounded,  for r̃ = 0, |x| > 1.

(5.2)

The far-field boundary condition (2.15) is

ϕ(x, r̃) → \frac{D}{4π} \frac{x}{(x^{2} + Kr̃^{2})^{3/2}} ,  as (x^{2} + r̃^{2})^{1/2} → ∞  (5.3)

where

D = D_{body} + D_{flow} = π ∫_{−1}^{+1} F^{2}(x) dx + π(γ + 1) ∫_{−∞}^{∞} dx ∫_{0}^{∞} ϕ_{x}^{2}(x, r̃)r̃dr̃. (5.4)

In order to solve the problem numerically, it must be approximated discretely similar to the procedure outlined in Chapter 3 for the hodograph version. A rectangular numerical grid is overlaid on a truncated version of the domain of interest of the boundary value problem, the upper half-plane. A grid is defined as

\begin{matrix}x_{i} = iΔx, & −n ≤ i ≤ n + 1, \\ r̃_{j} = r̃_{j−1} + Δr̃_{j}, & 1 ≤ j ≤ m,  r̃_{0} given. \end{matrix}  (5.5)

The numerical grid is scaled exponentially so that there are more grid points as the

Figure 5.1: Exponentially-scaled numerical grid employed in (x, r̃) plane

x-axis is approached at r̃ = 0. This can be seen from the diagram of a typical grid given in Figure 5.1. The governing PDE (5.1) is discretized using a conservative finite difference scheme with ϕ_{i,j} ≈ ϕ(x_{i+1/2}, r̃_{j+1/2}). That is, ϕ_{i,j} is defined at cellcenters. (5.1) is re-written in conservative form as

\frac{∂P}{∂x} + \frac{∂Q}{∂r̃} = 0,  (5.6)

where

P = r̃(K − \frac{(γ + 1)}{2} \frac{∂ϕ}{∂x}) \frac{∂ϕ}{∂x} and Q = r̃\frac{∂ϕ}{∂r̃}. (5.7)

The discretized form of (5.6) is

\frac{1}{Δx} (P_{i+1/2,j} − P_{i−1/2,j}) + \frac{2}{Δr̃_{j} + Δr̃_{j+1}} (Q_{i,j+1/2} − Q_{i,j−1/2}) = 0  (5.8)

where

P_{i+1/2,j} = \{ \begin{matrix}r̃_{j} (K − \frac{γ + 1}{2} \frac{ϕ_{i+1,j} − ϕ_{i,j}}{Δx} )(\frac{ϕ_{i+1,j} − ϕ_{i,j}}{Δx} ) & if elliptic \\ r̃_{j} (K − \frac{γ + 1}{2} \frac{ϕ_{i,j} − ϕ_{i−1,j}}{Δx} )(\frac{ϕ_{i,j} − ϕ_{i−1,j}}{Δx} ) & if hyperbolic \end{matrix} (5.9)

and

Q_{i,j+1/2} = \frac{ϕ_{i,j+1} − ϕ_{i,j}}{log r̃_{j+1} − log r̃_{j}} .  (5.10)

The test for whether the flow is elliptic or hyperbolic is

if (\frac{ϕ_{i+1,j} − ϕ_{i−1,j}}{Δx} ) > \frac{K}{γ + 1} ⇒ elliptic flow,
if (\frac{ϕ_{i+1,j} − ϕ_{i−1,j}}{Δx} ) < \frac{K}{γ + 1} ⇒ hyperbolic flow.

(5.11)

The two separate formulas given in (5.9) for P are an attempt to recognize the nature of the mixed-type PDE (5.1). The discretization scheme for Q given in (5.10) is used in order to accurately discretize this flux near r̃ = 0. Note that farther from r̃ = 0 the scheme becomes equivalent to standard centered differencing. Further discussion of the discretization of Q will be postponed to section 5.1.

It is commonly known that elliptic equations and hyperbolic equations require different types of discretizations schemes. Stable hyperbolic schemes must reflect the fact that information carried by the equation is directional. In (5.9), the direction of the flow is the same direction information travels in the hyperbolic scheme. In the discretization of elliptic PDEs the property that information at each point in the domain affects every other point in the domain must be represented in the difference scheme. The difference scheme given in (5.9) recognizes the dual nature of the Kármán-Guderley equation and alters the scheme when the flow is elliptic or hyperbolic. This technique is known as the Murman-Cole switching method and was introduced in [7]. It is the classical solution to the problem of numerically discretizing mixed-type elliptic-hyperbolic partial differential equations. Other methods of discretization have been developed, most notably the use of artificial viscosity, so that instead of a sudden binary switch between the different discretizations, a phased-in approach is used. The work of Engquist and Osher in [9] is another possible method.

Two observations about the chosen method of numerical solution of (5.1) can be made. The first is that it is the discretization of the conservative form of the TSDE that is being conducted. The conservative form was chosen because the work of Lax and Wendroff [15] suggests that possible shock formation will be captured accurately. The second observation is that the local truncation error of the discretizations in (5.9) are different. The hyperbolic discretization is first order O(Δx), while the elliptic discretization method is second order O(Δx^{2}). These discretization methods were chosen because it is believed that only a small fraction of the points in the numerical grid will be hyperbolic so a more complex discretization method was not necessary. The shock-jump conditions used to model the flow across shock waves are first-order accurate, which sets the accuracy of shock capturing. Other higher-order methods of hyperbolic discretizations are also given in [7].

The boundary conditions for ϕ(x, r̃) are also approximated numerically. The near-field boundary condition is treated by imposing the condition that the flux Q = (r̃ϕr̃)|_{r̃=0} is equal to the source distribution function S(x). This is discretized as

\frac{ϕ_{i,j+1} − ϕ_{i,j}}{log r̃_{j+1} − log r̃_{j}} = S(x_{i}), at j = 0, −n_{b} ≤ i ≤ n_{b} − 1. (5.12)

The body of revolution exists between −n_{b} ≤ i ≤ n_{b}, which corrresponds to |x| ≤ 1. The far-field boundary condition is that ϕ_{i,j} is prescribed along the outer perimeter of the numerical grid by

ϕ_{i,j} = \frac{D}{4π} \frac{x_{i}}{(x_{i}^{2} + Kr̃_{j}^{2})^{3/2},}  (5.13)

where

D ≈ D_{body} + π(γ + 1) ∑_{i} ∑_{j} r̃_{j} (\frac{ϕ_{i,j} − ϕ_{i−1,j}}{Δx} )^{2} ΔxΔr̃_{j}.  (5.14)

The integral for D is approximated numerically and evaluated using a two-dimensional version of the Midpoint rule.

5.1 New numerical discretization scheme

The discretization scheme for Q involving the r̃-derivatives uses a new scheme discovered while trying to accurately calculate the perturbation potential near the body at r̃ = 0, namely ϕ_{i,0}. In order to calculate the function G(x) in (5.2) the perturbation potential ϕ needs to be evaluated as r̃ → 0. However, at r̃ = 0 the Kármán-Guderley equation (5.1) is singular. Many differential equations written in cylindrical or spherical co-ordinates have this property. Often the boundary condition is simply that the solution is bounded at the singular point. However, in the case of the boundary value problem for ϕ(x, r̃) the boundary condition comes from an inner solution (described in Appendix B) which imposes the condition that ϕ blows up along r̃ = 0, |x| ≤ 1. The important features of the singular behavior of ϕ can be investigated by solving the associated boundary value problem for u(r̃) shown below:

\frac{d}{dr̃} (r̃\frac{du}{dr̃}) = 0,  (5.15)
r̃\frac{du}{dr̃}|_{r̃=0} = S,  (5.16)
u(1) = G.  (5.17)

The exact solution to the above boundary value problem is

u(r̃) = S log r̃ + G, 0 ≤ r̃ ≤ 1. (5.18)

This solution has a singularity along r̃ = 0 (u → −∞ as r̃ → 0 for S ≠ 0). An investigation of the singular nature of ϕ as r̃ → 0 can be made by numerically solving the boundary value problem for u(r̃) using standard centered differencing for (5.15) on a uniform grid given by r̃_{j} = jΔr̃ where u(r̃_{j}) = u_{j}:

\frac{1}{Δr̃} \{(r̃_{j+1/2} \frac{u_{j+1} − u_{j}}{r̃_{j+1} − r̃_{j}} ) − (r̃_{j−1/2} \frac{u_{j} − u_{j−1}}{r̃_{j} − r̃_{j−1}} )\} = 0,  (5.19)

with

r̃_{1/2} \frac{u_{1} − u_{0}}{r̃_{1} − r̃_{0}} = S,  (5.20)

and

u_{m} = G,  (5.21)

where S = 1, G = 5 and Δr̃ = (r̃_{max} − r̃_{0})/m. Note that (5.19) can also be

Figure 5.2: Γ = u(r̃_{0}) − S log r̃_{0} versus r̃_{0} (for S = 1, G = 5)

represented as

\frac{1}{Δr̃} \{Q̂_{j+1/2} − Q̂_{j−1/2}\} = 0, where Q̂_{j} = r̃_{j} \frac{u_{j+1} − u_{j}}{r̃_{j+1} − r̃_{j}} .

When the system defined by (5.19),(5.20) and (5.21) is solved repeatedly with m = 100 and r̃_{max} = 1 using decreasing values of r̃_{0} the graph in Figure 5.2 is the result. This indicates that the closer to r̃ = 0 the boundary condition (5.17) is applied, the less accurate the numerical solution becomes. In other words, the value of Γ(r̃) = u(r̃) − S log r̃ diverges from the expected value of G as r̃_{0} → 0.

However, by using an alternative differencing similar to that of Q introduced in the previous section (5.15) becomes

\frac{1}{Δr̃} \{(\frac{u_{j+1} − u_{j}}{log r̃_{j+1} − log r̃_{j}} ) − (\frac{u_{j} − u_{j−1}}{log r̃_{j} − log r̃_{j−1}} )\} = 0.  (5.22)

By observation, (5.22) is solved exactly by u_{j} = S log r̃_{j} + G, the numerical representation of the exact solution (5.18). Another way to show the difference between the usual discretization method (5.19) and the special discretization (5.22) for the sample boundary value problem (5.15)-(5.17) is to look at the respective truncation errors, τ and τ̂:

τ_{j} = −\frac{Δr̃^{2}}{24} (r̃_{j}u′′′_{j})′  (5.23)
τ̂_{j} = −\frac{Δr̃^{2}}{24} [r̃_{j}u′′′_{j} − (\frac{2u′_{j}}{r̃_{j}} )]′ ,  (5.24)

where ′ signifies an r̃ derivative. If u(r̃) has the form given in (5.18), τ ≡ 0 while

τ̂_{j} = \frac{1}{6} \frac{Δr̃^{2}}{r̃_{j}^{3}} .

If r̃_{j} ∼ O(Δr̃) then clearly τ̂_{j} ∼ O(\frac{1}{Δr̃}), which does not go to zero as Δr̃ → 0. Thus the new method will be exact, while the other method will not even be consistent. This explains the behavior of Γ displayed in Figure 5.2 which is the result of using the same Δr̃ and repeatedly solving the boundary value problem for u(r̃) with decreasing values of r̃_{0}. As soon as r̃_{0} becomes less than the Δr̃ used on the grid, the truncation error of the usual scheme begins to grow, causing the diversion of Γ(r̃_{0}) from the expected value of u(r̃_{0}) − S log r̃_{0} = G. For the alternate scheme, Γ(r̃_{0}) = G regardless of how small r̃_{0} is, as the dotted line in Figure 5.2 illustrates.

It should be noted that the principle behind the alternate scheme is similar to the method used to develop the Scharfetter-Gummel scheme introduced in [22], which occurs in semiconductor design. This new differencing scheme can be applied to equations which resemble the Laplace equation in spherical and cylindrical coordinates which have similar singularities at r̃ = 0.

5.2 Numerical solution of the physical plane boundary value problem

The boundary value problem for ϕ(x, r̃) in the physical plane consisting of (5.1), (5.2) and (5.3) has been reduced to a set of simultaneous discrete equations for ϕ_{i,j} which must be solved numerically. The method chosen is successive line relaxation because of its simplicity of implementation. Generally, the iteration method converges slowly, requiring thousands of iterations. Other methods, such as Newton’s method or using multiple grids could be used to increase the speed of convergence but this was not considered necessary for the calculation at hand.

Let the solution ϕ_{i,j} along a vertical column of the computational mesh be given by

X_{i} = ( ϕ_{i,1} ; ⋮ ; ⋮ ; ⋮ ; ϕ_{i,m} )

where X_{i} is a m-dimensional vector. The system of discrete equations can be written as F(X_{i}) = 0 and one iteration of Newton’s method is performed by solving the system

J_{i}dX_{i} = −F

for dX_{i}, where J_{i} = \frac{∂}{∂ϕ_{i,k}} F_{i}(ϕ_{i,1}, ϕ_{i,k}, …, ϕ_{i,m}) for 1 ≤ k ≤ m. This is equivalent to doing one line relaxation step on X_{i} by updating X_{i}^{new} = X_{i}^{old} + dX_{i}. The line relaxation process is repeated as the method sweeps from i = −n + 1 to n − 1. The direction of the sweep corresponds to the direction of the fluid flow across the body, from nose to tail. In the far-field, ϕ is updated using the expression (5.3) for the assumed dipole flow. The value of the dipole strength D, which depends on the solution ϕ(x, r̃), is computed by evaluating the integrals in (5.4) approximately using the Midpoint Rule with the current value of ϕ_{i,j} every one hundred relaxation sweeps.

5.3 Numerical experiments in the physical plane

Once it has been established that a computer code has been developed which accurately solves the boundary value problem for ϕ(x, r̃), investigations can be conducted of the shock-free results obtained from the hodograph calculations outlined in Chapter 4. One aspect of the usefulness of the code is that off-design conditions can also be investigated. The computer code acts as a “numerical wind tunnel” in which experiments can be conducted to investigate flow around slender bodies at transonic speeds.

One such experiment involves comparing the computed flow about the initial shape, which is shocked, and the designed final shape, which is shock-free. Referring to Chapter 4 which contains details on a number of bodies possessing shock-free flows computed in the hodograph, input values of α and K to use in the physical plane calculation are obtained. The physical plane computations are used a means of validating the hodograph results and investigating off-design conditions.

For example, computing around a body generated by substituting α = −.09251 in F (x; α) at K = 3.919 in the physical plane results in a solution whose isobars are displayed in Figure 5.3. By setting m = n = 100 in (5.5) and n_{b} = 50 in (5.12) the numerical grid used to produce this solution has 201 points in the x-direction with 100 points in the r̃-direction, exponentially scaled, similar to Figure 5.1. Running on a Sparcstation-10, the maximum residual on all ∼20 000 points on the numerical grid falls below 10^{−4} after ∼8000 sweeps in about 25 minutes. A distinct shock is clearly visible by the coalescence of the isobars. Figure 5.4 shows the plot of G′(x) versus x for this K = 3.919 solution. The typical graph to show the shock-free nature of the flow would be of the pressure coefficient C_{p_{body}} versus x. However since C_{p_{body}} is essentially the same as G′(x), in this chapter graphs of G′(x) are Figure 5.3: Contour plot of u = ϕ_{x}(x, r̃) for α = −.0925 and K = 3.919 Figure 5.4: Physical plane G′(x) data for α = −.0925 and K = 3.919 given. If one examines the expression for C_{p_{body}}(x) given in (2.12) and derived in Appendix B, it is clear that the function which acts as a source for jumps is the function G′(x). The other functions F (x), F′(x) and S′(x) in (2.12) are all relatively smooth functions, since F is smooth. The jump in the G′ data at x ≈ 0.38 makes it clear the flow around this particular body possesses a shock.

Figure 5.5: Difference between input and output source distribution data for SFB-2

The outputted source distribution data that is computed from the discrete hodograph variables x_{i+1/2,j+1/2} and R_{i,j} is not the same as the inputted function S(x; α). The iterative solution algorithm described in section 3.5 has resulted in data which are perturbations of the inputted function S(x; α). Figure 5.5 illustrates the difference between input and output source distributions used in the hodograph. This output source distribution is the source distribution data for the shock-free hodograph-designed body SFB-2. Figure 5.3 and Figure 5.4 illustrate that the computed flow around a body with a source distribution function represented by Figure 5.6: Contour plot of u = ϕ_{x}(x, r̃) for SFB-2 at K = 3.919

Figure 5.7: Plot of G′(x) for SFB-2 at K = 3.919 the dotted line in Figure 5.5 clearly possesses a distinct shock. Using a smoothed version of the SFB-2 source distribution data as input (shown as the solid line in Figure 5.5), the physical plane code is run again at the same value of K used to produce Figure 5.3. This value of K, 3.919, is the design value for SFB-2. The result is a contour plot shown in Figure 5.6 which is dramatically different from the previous contour plot in Figure 5.3. It is clear that the distinct shock of Figure 5.3 has disappeared, replaced by a possible weak shock near the extrema point x_{∗} in Figure 5.6. It is not a surprise that the isobars indicate a shock-free flow, since SFB-2 (and the others SFB-1, SFB-3 and SFB-4) were designed to be shock-free at their own particular value of K.

The G′(x) data for flow around the hodograph-designed body SFB-2 in the physical plane is shown in Figure 5.7. It is qualitatively different from the G′(x) data shown in Figure 5.4. But G(x) in the physical plane calculation is a difficult expression to evaluate on the body since it is the difference between two large quantities |ϕ_{i,0} − S(x_{i}) log r̃_{0}| as r̃_{0} → 0. It was during the search to improve the accuracy of this computation that the new discretization scheme given in section 5.1 was found. It should not be hard to see the connection between computing Γ(r̃_{0}) accurately in Section 5.1 and computing G(x) on the body. Using the new discretization improves confidence in the calculation of G(x), with the proviso that the data for S(x) must be reliable. The computation of G(x) is extremely sensitive to the nature of S(x). However, due to the geometry of the hodograph plane, S(x) is difficult to recover accurately near x^{∗} and x_{∗}, the extrema of S′(x). As noted in chapter 2, these points are mapped to the singular horizontal lines ν = ν^{∗} and ν = ν_{∗} in the hodograph. Unfortunately, it is exactly at these extrema where information about the source distribution needs to be the most accurate. It is shown by Cole and Malmuth in [5] that shocks if they do form, will attach themselves to the body at precisely these points. As can be seen in Figure 4.8 from Chapter 4, it is along these singular boundaries of the hodograph that the Jacobian surface plots exceed the shock-inducing value of zero. Using the data available one is not able to say with certainty whether the slight dips in G′(x) apparent in Figure 5.7 indicate the occurrence of weak shocks or the result of inaccurate or missing data in S(x) obtained from the hodograph near the extrema x^{∗} and x_{∗}. However the flat shape of the G′(x) curve over more than half the body is reminiscent of the numerical results for shocked supercritical flow over bodies given by Krupp and Murman in [14]. It is clear that the pressure distribution has been dramatically altered by using the hodograph-designed source distribution as input.

Investigations into how the computed flow around hodograph-designed bodies changes for off-design conditions can be conducted. According to the hodograph data the flow is only shock-free at the design value of K because it is only at this particular value of K that the jacobian remains negative everywhere in the hodograph plane. In Figure 5.8(a) and Figure 5.8(b) contour plots can be seen for values of K 10 percent above and below the design K of 3.919. When K is decreased the likelihood of shock development increases because the disturbances caused by the body are greater. Also, the volume of the supersonic bubble is larger. However, no substantial difference in the flow is observed, as the similarity in the contour plots indicates. In fact, Figure 5.9 depicts plots of the pressure coefficient (actually G′) along the body at various values of K other than the design K for SFB-2. The fact that a shock only forms after a significant variation in K (about 25% smaller than the design K) supports the claim that indeed shock-free transonic flows exist around bodies of revolution at speeds different from the design parameters. The results of exploring the off-design behavior of the hodograph-designed shock-free body SFB-2 given in Figure 5.8(b) and Figure 5.9 support the claim that transonic flows around slender bodies of revolution with little or no shock waves present can be computed with regularity.

(a) Contour plot at off-design K=4.311 (ΔK = +10%)

(b) Contour plot at off-design K=3.527 (ΔK = -10%)

Figure 5.8: Isobars around hodograph-designed body SFB-2 at off- design conditions 79

Figure 5.9: G′(x) associated with SFB-2 for various K above and below design

CHAPTER 6 DISCUSSION AND CONCLUSIONS

Numerical computations have been made in the hodograph plane to compute solutions of the transonic small-disturbance equation which are everywhere continuous. These solutions correspond to flows around shock-free transonic slender bodies. The condition for a shock-free solution in the hodograph plane is that the jacobian of the hodograph transformation is strictly less than zero in the entire plane. The algorithm involves setting the initial shape of the sonic line and then computing the corresponding subsonic and supersonic flows. A number of shock-free bodies, both fore-aft symmetric and fore-aft asymmetric, were calculated using this method in the hodograph plane.

Numerical calculations are also carried out in the physical plane to verify the shock-free nature of the hodograph-designed bodies and to investigate off-design conditions. In the physical plane, the bodies appear to have small weak shocks attached to them. These shocks do not substantially increase in size or strength when the physical plane computations are conducted at values of the transonic similarity parameter K = (1 − M_{∞}^{2})/(M_{∞}δ)^{2}, which differ from the hodograph design value by up to 20%.

These results show that a systematic algorithm to design (nearly) shock-free transonic slender bodies has been developed. The fact that the flows remain nearly shock-free for a range of the parameters M_{∞} and δ illustrates the practical nature of the design method introduced in this thesis.

Future work will involve improving some technical numerical features of the calculations in order to more accurately resolve the apparent small weak shocks apparent at the extrema points. Other methods of verifying the shock-free nature of the flow around the hodograph-designed bodies by solving the full potential equation instead of the transonic small-disturbance approximation can be considered. The idea of using the hodograph-designed bodies as starting points for shock-free design in the physical plane is another avenue for future research.

LITERATURE CITED

[1] Bauer F., P. Garabedian and D. Korn, Supercritical Wing Sections, Lecture Notes in Economics and Mathematical Systems, No. 66, Springer-Verlag, New York, 1972.

[2] Bauer F., P. Garabedian and D. Korn, Supercritical Wing Sections II, Lecture Notes in Economics and Mathematical Systems, No. 108, Springer-Verlag, New York, 1975.

[3] Bauer F., P. Garabedian and D. Korn, Supercritical Wing Sections III, Lecture Notes in Economics and Mathematical Systems, No. 150, Springer-Verlag, New York, 1977.

[4] Cole J.D., and L.P. Cook, Transonic Aerodynamics, North Holland, 1986.

[5] Cole J.D., and N. Malmuth, Shock Wave Location on a Slender Transonic Body of Revolution, Mechanics Research Communications, 16(6), (1989), pp. 353-357.

[6] Cole J.D., and A.F. Messiter, Expansion Procedures and Similarity Laws for Transonic Flow, Z.A.M.P., 8, (1957), pp. 1-25.

[7] Cole J.D., and E.M. Murman, Calculation of Plane Steady Transonic Flows, AIAA Journal, Volume 9 Number 1, (January 1971), pp. 114-121

[8] Cole J.D., and D.W. Schwendeman, Hodograph Design of Shock-free Transonic Slender Bodies, 3rd International Conference on Hyperbolic Problems, Editors Bjorn Engqvist and Bertil Gustafsson, Uppsala, Sweden, June 11-15 1990.

[9] Engquist, B. and S. Osher, Upwind Difference Schemes for Systems of Conservation Laws: potential flow equations, Scientific and Technical Aerospace Reports, Volume 19 (March 1981), N81-246424.

[10] Kacprzynski, J.J., L. H. Ohman, P.R. Garabedian, and D.G. Korn, Analysis of the Flow Past a Shockless Lifting Airfoil in Design and Off-design Conditions , N.R.C. of Canada Aeronautical Report LR-554, Ottawa (1971).

[11] Kacprzynski, J.J., Wind Tunnel Test of the Shockless Lifting Airfoil No. 1, Project Report 5x5/0062, National Research Council of Canada, Ottawa (1972).

[12] Kacprzynski, J.J., Wind Tunnel Test of the Shockless Lifting Airfoil No. 2, Labaratory Technical Report LTR-HA-5x5/0067, National Research Council of Canada, Ottawa (1973).

[13] Kropinski, M.C., A Study of Optimal Critical Airfoils, Ph.D. thesis, Rensselaer Polytechnic Institute, Troy, NY (1993).

[14] Krupp J.A., and E.M. Murman, Computation of Transonic Flows past Lifting Airfoils and Slender Bodies, AIAA Journal, Volume 10 Number 7, (July 1972), pp. 880-886

[15] Lax P.D., and B. Wendroff, Systems of Conservation Laws, Commun. Pure Appl. Math., 13, pp. 217-237 (1960).

[16] Morawetz C.S., On the non-existence of continuous transonic flows past airfoils at high Reynolds numbers I., Commun. Pure Appl. Math., 9, pp. 45-68 (1956).

[17] Morawetz C.S., On the non-existence of continuous transonic flows past airfoils at high Reynolds numbers II., Commun. Pure Appl. Math., 10, pp. 107-131 (1957).

[18] Morawetz C.S., On the non-existence of continuous transonic flows past airfoils at high Reynolds numbers III., Commun. Pure Appl. Math., 11, pp. 129-44 (1958).

[19] Nieuwland G.Y., Transonic potential flow around a family of quasi-elliptical sections , Technical Report T. 172 (1967).

[20] Pearcey H.H., The aerodynamic design of section shapes for swept wings, Advances in Aeronautical Sciences, Vol. 3 (1962), London.

[21] Ringleb, F., Exacte Lösungen der Differentialgleichunsen eineradiabatischen Gasströmung, ZAMM 20 , pp. 185-198 (1941).

[22] Scharfetter, D.L. and H.K. Gummel, Large-Signal Analysis of a Silicon Read Diode Oscillator, IEEE Transactions on Electron Devices, Vol. ED-16, Number 1, (January 1969) pp. 64-84.

[23] H. Sobieczky and A.R. Seebass, Supercritical Airfoil and Wing Design, Ann. Rev. Fluid Mech. (1984), 16, pp. 337-363.

[24] H. Sobieczky, K.-Y. Fung, A.R. Seebass and N.J. Yu, A New Method for Designing Shock-free Transonic Configurations, AIAA J., 17, (1979) pp. 722-729.

[25] Whitcomb R.T., and L. Clark, An airfoil shape for efficient flight at supercritical Mach numbers, NASA Tech. Man. X-1109 (1965).

APPENDIX A The Dipole Strength Integral

The value of the strength of the dipole D is obtained from re-writing the Kármán-Guderley equation (2.6) as an inhomogeneous form of the Prandtl-Glauert equation (2.14) and then applying Green’s Theorem. To wit, the TSDE is re-written as

Kϕ_{xx} + ϕ_{r̃r̃} + \frac{1}{r̃}ϕr̃ = (γ + 1)ϕ_{x}ϕ_{xx}.  (A.1)

This can be written as

ϕ_{xx} + \frac{1}{K} ϕ_{r̃r̃} + \frac{1}{Kr̃}ϕr̃ = \frac{(γ + 1)}{2K} \frac{∂ϕ_{x}^{2}}{∂x} .  (A.2)

If one lets r̂ = \sqrt{K}r̃ then

\frac{∂}{∂r̂} = \frac{dr̃}{dr̂} \frac{∂}{∂r̃} = \frac{1}{\sqrt{K}} \frac{∂}{∂r̃},  and  \frac{∂^{2}}{∂r̂^{2}} = \frac{1}{K} \frac{∂^{2}}{∂r̃^{2}} .

Thus the previous equation can be rewritten as

ˆ∇^{2}ϕ ≡ \frac{1}{r̂}ϕr̂ + ϕ_{r̂r̂} + ϕ_{xx} = \frac{ϵ}{2} \frac{∂ϕ_{x}^{2}}{∂x} .

Recall that the fundamental solution for the Laplace equation in three-dimensional cartesian co-ordinates is

G(x, y, z; ξ, η, ζ) = −\frac{1}{4πσ} ,  where σ = \sqrt{(x − ξ)^{2} + (y − η)^{2} + (z − ζ)^{2}}. (A.3)

By definition, the fundamental solution vanishes at infinity and solves the equation

ˆ∇^{2}G = δ(σ) = δ(x, y, z; ξ, η, ζ).

Applying this information to Green’s formula

∫_{V} (G ˆ∇^{2}ϕ − ϕ ˆ∇^{2}G) dV = ∫_{∂V} (G \frac{∂ϕ}{∂n} − ϕ\frac{∂G}{∂n}) dA  (A.4)

produces

\begin{matrix}∫∫∫_{V} (G(x − ξ, y − η, z − ζ) ˆ∇^{2}ϕ(ξ, η, ζ) − ϕ(ξ, η, ζ) ˆ∇^{2}G(x, y, z; ξ, η, ζ)) dξdηdζ \\ = ∫_{∂V} (G \frac{∂ϕ}{∂n} − ϕ\frac{∂G}{∂n}) dA \end{matrix}

(A.5)

and

\begin{matrix}∫∫∫_{V} (G(x − ξ, y − η, z − ζ)\frac{ϵ}{2} \frac{∂ϕ_{ξ}^{2}}{∂ξ} ) dξdηdζ − ϕ(x, y, z) \\ = −∫_{A_{1}}(G \frac{∂ϕ}{∂n} − ϕ\frac{∂G}{∂n}) dA − ∫_{A_{2}}(G \frac{∂ϕ}{∂n} − ϕ\frac{∂G}{∂n}) dA. \end{matrix}  (A.6)

The above can be re-written so that

ϕ(x, y, z) = I_{1} + I_{2} − I_{3}  (A.7)

where

I_{1} = ∫_{A_{1}}(G \frac{∂ϕ}{∂n} − ϕ\frac{∂G}{∂n}) dA,  (A.8)
I_{2} = ∫_{A_{2}}(G \frac{∂ϕ}{∂n} − ϕ\frac{∂G}{∂n}) dA,  (A.9)

and

I_{3} = ∫∫∫_{V} \frac{∂G}{∂ξ} (x, y, z, ξ, η, ζ)\frac{ϵ}{2}ϕ_{ξ}^{2} dξdηdζ.  (A.10)

The two areas A_{2} and A_{1} are the external surface of an expanding sphere, and the extenal surface of a contracting cylinder, respectively.

A.1 Evaluating I_{1}

The knowledge that the integral is in cylindrical co-ordinates should be used so that the expressions can be made more specific.

I_{1} = ∫_{−∞}^{∞} ∫_{0}^{2π}(G \frac{∂ϕ}{∂ρ̂} − ϕ\frac{∂G}{∂ρ̂})ρ̂dθdξ  (A.11)

The Green’s function and ϕ can be written in the (ξ, ρ̂) variables

G(ξ, ρ̂) = \frac{−1}{4π}[(x − ξ)^{2} + r̂^{2} + ρ̂^{2} − 2r̂ρ̂cos θ]^{−1/2}.

In Chapter 2 it is given that as r̃ → 0, the perturbation potential ϕ takes the form ϕ(x, r̃) = S(x) log r̃ + G(x), which can be written as

ϕ(ξ, ρ̂) = S(ξ) log(\frac{ρ̂}{\sqrt{K}}) + G(ξ) + …

The derivatives of the functions ϕ and G will also be needed

\frac{∂G}{∂ξ} = −\frac{x − ξ}{4π} [(x − ξ)^{2} + r̂^{2} + ρ̂^{2} − 2r̂ρ̂cos θ]^{−3/2}  (A.12)
\frac{∂G}{∂ρ̂} = \frac{ρ̂ − r̂cos θ}{4π}[(x − ξ)^{2} + r̂^{2} + ρ̂^{2} − 2r̂ρ̂cos θ]^{−3/2}  (A.13)
\frac{∂ϕ}{∂ξ} = S′(ξ) log(\frac{ρ̂}{\sqrt{K}}) + G′(ξ) + …  (A.14)
\frac{∂ϕ}{∂ρ̂} = \frac{S(ξ)}{ρ̂} + …  (A.15)

The integral I_{1} is evaluated on the surface of a cylinder of radius ρ̂_{0}. By substituting the above derivatives into the expression (A.11) for I_{1} the integral becomes

I_{1} = ∫_{−∞}^{∞} ∫_{0}^{2π} ρ̂_{0}dθdξ \{(\frac{S(ξ)}{ρ̂_{0}}) (\frac{−1}{4π}[(x − ξ)^{2} + r̂^{2} + ρ̂_{0}^{2} − 2r̂ρ̂_{0} cos θ]^{−1/2})

− (S(ξ) log(\frac{ρ̂_{0}}{\sqrt{K}}) + G(ξ)) (\frac{ρ̂_{0} − r̂cos θ}{4π}[(x − ξ)^{2} + r̂^{2} + ρ̂_{0}^{2} − 2r̂ρ̂_{0} cos θ]^{−3/2})\}

(A.16)

Taking the limit as ρ̂_{0} → 0 of I_{1} produces

I_{1} = ∫_{−∞}^{∞} ∫_{0}^{2π} −\frac{S(ξ)}{4π}[(x − ξ)^{2} + r̂^{2}]^{−1/2} dθdξ

− lim_{ρ̂_{0}→0} \frac{1}{4π} ∫_{−∞}^{∞} ∫_{0}^{2π} (S(ξ) log(\frac{ρ̂_{0}}{\sqrt{K}}) + G(ξ)) \frac{(ρ̂_{0} − r̂cos θ)ρ̂_{0} dθdξ}{[(x − ξ)^{2} + r̂^{2} + ρ̂_{0}^{2} − 2r̂ρ̂_{0} cos θ]^{3/2}} .

(A.17)

This simplifies in the limit to

I_{1} = −\frac{1}{2} ∫_{−1}^{1} \frac{S(ξ)dξ}{\sqrt{(x − ξ)^{2} + r̂^{2}}} .  (A.18)

A.2 Evaluating I_{2}

Since I_{2} involves integrating around an expanding sphere of radius R the integral can be re-written more specifically as

I_{2} = lim_{R→∞} ∫_{−π}^{π} ∫_{0}^{2π} \{G(R⃗, R⃗_{0})\frac{∂ϕ}{∂R} − ϕ(R⃗, R⃗_{0})\frac{∂G}{∂R}\} R^{2} sin ϑdθdφ  (A.19)

Using the fact that the Green’s function is in spherical co-ordinates

G(R⃗, R⃗_{0}) = \frac{−1}{4π} \frac{1}{|R⃗ − R⃗_{0}|}  (A.20)

where R⃗ = xx⃗ + yy⃗ + zz⃗ and R⃗_{0} = ξx⃗ + νy⃗ + ζz⃗. The derivative of G can be calculated also

\frac{∂G}{∂R} = G_{R} = \frac{1}{4πR} \frac{R⃗ · (R⃗ − R⃗_{0})}{|R⃗ − R⃗_{0}|^{3}}  (A.21)

If one assumes ϕ has the form of a dipole, then

ϕ ∼ \frac{D}{|R⃗ − R⃗_{0}|^{3}} ,  D constant  (A.22)

and

ϕ_{R} ∼ \frac{3DR⃗ · (R⃗ − R⃗_{0})}{R|R⃗ − R⃗_{0}|^{5}}.  (A.23)

Using these values for ϕ, G and their derivatives,

I_{2} = −\frac{3D}{4π} ∫_{−π}^{π} ∫_{0}^{2π} \frac{RR⃗ · (R⃗ − R⃗_{0})}{|R⃗ − R⃗_{0}|^{6}} sin φ dθdφ

(A.24)

−\frac{D}{4π} ∫_{−π}^{π} ∫_{0}^{2π} \frac{RR⃗ · (R⃗ − R⃗_{0})}{|R⃗ − R⃗_{0}|^{6}} sin φ dθdφ.  (A.25)

It can be seen that lim_{R→∞} I_{2} = 0.

A.3 Evaluating I_{3}

The third integral involves the volume integral around the contracting cylinder of radio ρ̂_{0}. Using this fact, I_{3} can be re-written in cylindrical co-ordinates as

I_{3} = ∫_{−∞}^{∞} ∫_{−∞}^{∞} ∫_{−∞}^{∞} \frac{∂G}{∂ξ} \frac{ϵ}{2}ϕ_{ξ}^{2} dξdηdζ = 2π ∫_{−∞}^{∞} ∫_{0}^{∞} G_{ξ}(ξ, ρ̂)\frac{ϵ}{2}ϕ_{ξ}^{2}(ξ, ρ̂) ρ̂dρ̂dξ (A.26)

I_{3} = ∫_{−∞}^{∞}∫_{0}^{∞} −\frac{ϵ}{4} \frac{x − ξ}{[(x − ξ)^{2} + r̂^{2}]^{3/2}} ϕ_{ξ}^{2}(ξ, ρ̂) ρ̂dρ̂dξ.  (A.27)

If one looks at the expression \frac{x − ξ}{[(x − ξ)^{2} + r̂^{2}]^{3/2}}

\frac{x − ξ}{[(x − ξ)^{2} + r̂^{2}]^{3/2}} = \frac{x − ξ}{(x^{2} + r̂^{2})^{3/2}} \frac{1}{[1 − \frac{2ξx}{(x^{2} + r̂^{2})} + \frac{ξ^{2}}{(x^{2} + r̂^{2})}]^{3/2}}  (A.28)

as (x, r̂) → ∞, (A.27) can be re-written

I_{3} = \frac{x}{(x^{2} + r̂^{2})^{3/2}} ∫_{−∞}^{∞}∫_{0}^{∞} −\frac{ϵ}{4} ϕ_{ξ}^{2}(ξ, ρ̂) ρ̂dξdρ̂.  (A.29)

Using this same expansion technique as (x, r̂) → ∞ for the integral I_{1}, (A.18) can be re-written as

I_{1} = −\frac{1}{2} \frac{x}{(x^{2} + r̂^{2})^{3/2}} ∫_{−1}^{1} ξS(ξ)dξ − \frac{1}{2} ∫_{−1}^{1} \frac{S(ξ) dξ}{(x^{2} + r̂^{2})^{1/2}} ,  (A.30)

since

\begin{aligned}\frac{1}{[(x − ξ)^{2} + r̂^{2}]^{1/2}} & = & \frac{1}{(x^{2} + r̂^{2})^{1/2}} \frac{1}{[1 − \frac{2ξx}{(x^{2} + r̂^{2})} + \frac{ξ^{2}}{(x^{2} + r̂^{2})}]^{1/2}} \\ & = & \frac{1}{(x^{2} + r̂^{2})^{1/2}} [1 + \frac{ξx}{(x^{2} + r̂^{2})} + …] . \end{aligned}  (A.31)

But ∫_{−1}^{1} S(ξ)dξ = \frac{F^{2}(ξ)}{2}|_{−1}^{1} = 0, since F (1) = F (−1) = 0. Remembering (A.7), an expression for ϕ(x, r̂) can be written as

ϕ = \frac{x}{(x^{2} + r̂^{2})^{3/2}} \{−\frac{1}{2} ∫_{−1}^{1} ξS(ξ)dξ + \frac{ϵ}{4} ∫_{−∞}^{∞}∫_{0}^{∞} ϕ_{ξ}^{2}(ξ, ρ̂) ρ̂dξdρ̂\},  (A.32)

since

∫_{−1}^{1} ξS(ξ)dξ = ξ \frac{F^{2}(ξ)}{2}|_{−1}^{1} − ∫_{−1}^{1} \frac{F^{2}(ξ)}{2}dξ.

So a final integral expression for ϕ can be written

ϕ(x, r̂) = \frac{1}{4} \frac{x}{(x^{2} + r̂^{2})^{3/2}} \{∫_{−1}^{1} F^{2}(ξ)dξ + ϵ ∫_{−∞}^{∞}∫_{0}^{∞} ϕ_{ξ}^{2}(ξ, ρ̂) ρ̂dξdρ̂\} (A.33)

which, when compared to the expression for a dipole of strength D below

ϕ(x, r̂) = \frac{D}{4π} \frac{x}{(x^{2} + r̂^{2})^{3/2}} ,

gives an expression for D. Namely

D = π {∫_{−1}^{1} F^{2}(ξ)dξ + ϵ ∫_{−∞}^{∞}∫_{0}^{∞} ϕ_{ξ}^{2}(ξ, ρ̂) ρ̂dξdρ̂} .  (A.34)

Note ρ̂ = \sqrt{K}ρ̃, so in (ξ, ρ̃) variables similar to (x, r̃) variables of Chapter 2 this equation becomes

D = π {∫_{−1}^{1} F^{2}(ξ)dξ + (γ + 1) ∫_{−∞}^{∞}∫_{0}^{∞} ϕ_{ξ}^{2}(ξ, ρ̃) ρ̃dξdρ̃} .

which is identical to (2.16) with (x, r̃) ⟺ (ξ, ρ̂).

APPENDIX B Transonic Small-Disturbance Theory Asymptotics

Starting from the Full Potential Equation (2.1) an assumption is made that the potential Φ(x, r; M_{∞}, δ) can be written as an asymptotic expansion, and the variable r scaled accordingly. This can be accomplished by introducing the unknown functions μ(δ), ν(δ), A(M_{∞}) and B(M_{∞}). The expression (2.5) for the outer potential ϕ can be written more generally

Φ = U_{∞} \{x + \frac{ν(δ)}{A(M_{∞})}ϕ(x, r̃; K) …\},  where r̃ = μ(δ)B(M_{∞})r.  (B.1)

The expressions in (B.1) will be used in a limit process type expansion with M_{∞} → 1 and δ → 0 simultaneously with (r̂, K) remaining fixed. The derivatives of Φ with respect to x and r will be needed and are given below

\frac{∂Φ}{∂x} = U_{∞}\{1 + \frac{ν(δ)}{A(M_{∞})}ϕ_{x}(x, r̃)\}
\frac{∂^{2}Φ}{∂x^{2}} = U_{∞} \frac{ν(δ)}{A(M_{∞})}ϕ_{xx}
\frac{∂Φ}{∂r̃} = U_{∞} \frac{ν(δ)}{A(M_{∞})}μ(δ)B(M_{∞})ϕr̃(x, r̃)
\frac{∂Φ}{∂r̃^{2}} = U_{∞} \frac{ν(δ)}{A(M_{∞})}(μ(δ)B(M_{∞}))^{2}ϕ_{r̃r̃}(x, r̃)
\frac{∂^{2}Φ}{∂x∂r̃} = U_{∞} \frac{ν(δ)}{A(M_{∞})}(μ(δ)B(M_{∞}))ϕ_{xr̃}(x, r̃)  (B.2)

(Of course \frac{∂}{∂r} = μ(δ)B(M_{∞})\frac{∂}{∂r̃}.) Bernoulli’s integral (2.2) can be re-written (by dividing by U_{∞}^{2}) as

(\frac{a}{U_{∞}} )^{2} = \frac{1}{M_{∞}^{2}} + \frac{γ − 1}{2} \{1 − (\frac{Φ_{x}}{U_{∞}} )^{2} − (\frac{Φ_{r}}{U_{∞}} )^{2}\} .  (B.3)

The full potential equation (2.1) can be re-written (by dividing by U_{∞}^{3}) as

[(\frac{a}{U_{∞}} )^{2} − (\frac{Φ_{x}}{U_{∞}} )^{2}] \frac{Φ_{xx}}{U_{∞}} + [(\frac{a}{U_{∞}} )^{2} − (\frac{Φ_{r}}{U_{∞}} )^{2}] \frac{Φ_{rr}}{U_{∞}}

+ (\frac{a}{U_{∞}} )^{2} \frac{Φ_{r}}{rU_{∞}} − 2 (\frac{Φ_{r}}{U_{∞}} ) (\frac{Φ_{x}}{U_{∞}} ) (\frac{Φ_{xr}}{U_{∞}} ) = 0.

(B.4)

Using (B.2) and (B.3), (B.4) becomes

\begin{matrix}\frac{ν}{A}ϕ_{xx} [\frac{1}{M_{∞}^{2}} − 1 + \frac{(γ + 1)}{2} \{\frac{ν^{2}}{A^{2}} ϕ_{x}^{2} − \frac{2ν}{A} ϕ_{x}\} − \frac{(γ − 1)}{2} \{\frac{ν^{2}μ^{2}B^{2}}{A^{2}} ϕ2 ; r̃ϕ_{x}\}] \\ +\frac{ν}{A}μ^{2}B^{2}ϕ_{r̃r̃} [\frac{1}{M_{∞}^{2}} + \frac{(γ − 1)}{2} \{\frac{ν^{2}}{A^{2}} ϕ_{x}^{2} − \frac{2ν}{A} ϕ_{x}\} − \frac{(γ + 1)}{2} \frac{ν^{2}μ^{2}B^{2}}{A^{2}} ϕ2 ; r̃] \\ +\frac{ν}{A}μ^{2}B^{2} \frac{ϕr̃}{r̃} [\frac{1}{M_{∞}^{2}} + \frac{(γ − 1)}{2} \{\frac{−ν^{2}}{A^{2}} ϕ_{x}^{2} − \frac{2ν}{A} ϕ_{x} − \frac{ν^{2}μ^{2}B^{2}}{A^{2}} ϕ2 ; r̃\}] \\ −2\frac{ν^{2}μ^{2}B^{2}}{A^{2}} ϕr̃ϕ_{xr̃} (1 + \frac{ν}{A}ϕ_{x}) = 0. \end{matrix} (B.5)

Equating all the dominant terms in (B.5) produces

\frac{1}{M_{∞}^{2}} − 1 = K \frac{ν}{A},  (B.6)
\frac{ν}{A} = \frac{(μB)^{2}}{M_{∞}^{2}} .  (B.7)

Then at O(\frac{ν^{2}}{A^{2}} ) the equation (B.5) becomes

Kϕ_{xx} − (γ + 1)ϕ_{x}ϕ_{xx} + ϕ_{r̃r̃} + \frac{1}{r̃}ϕr̃ = 0,

which is clearly a re-written version of the TSDE given in (2.6):

(K − (γ + 1)ϕ_{x})ϕ_{xx} + ϕ_{r̃r̃} + \frac{1}{r̃}ϕr̃ = 0.

Looking at a general form for the inner expansion produces the expression

Φ = U_{∞}{x + σ(δ)φ(x, r^{∗}) + …}  (B.8)

where r^{∗} = \frac{r}{δ} and φ(x, r^{∗}) = S(x) log r^{∗} + G(x). The boundary condition on the body is

δF′(x) = \frac{Φ_{r}(x, r)}{Φ_{x}(x, r)}, at r = δF (x) ⟺ r^{∗} = F (x). (B.9)

This implies that

\frac{Φ_{r}}{Φ_{x}} = \frac{σ(δ)\frac{S(x)}{r^{∗}} \frac{1}{δ}}{1 + σ(δ)[S′(x) log r^{∗} + G′(x)]} = δF′(x).

This can be simplified as

σ(δ)\frac{S(x)}{F (x)} = δ^{2}F′(x)  (B.10)

which determines the expressions σ(δ) = δ^{2} and S = F F′. If the inner expansion (B.8) is matched with the outer expansion (B.1)

\frac{ν(δ)}{A(M_{∞})} = σ(δ) = δ^{2} ⇒ ν(δ) = δ^{2}, A(M_{∞}) = 1

and

(μ(δ)B(M_{∞}))^{2} = δ^{2}M_{∞}^{2} ⇒ μ(δ)B(M_{∞}) = δM_{∞}.

Thus all previously unknown functions are now known

ν(δ) = δ^{2}  (B.11)
μ(δ) = δ  (B.12)
A(M_{∞}) = 1  (B.13)
B(M_{∞}) = M_{∞}  (B.14)
K = \frac{1 − M_{∞}^{2}}{(M_{∞}δ)^{2}}  (B.15)
r̃ = δM_{∞}r  (B.16)
r^{∗} = \frac{r̃}{δ^{2}M_{∞}} = \frac{r}{δ}  (B.17)

We would like the inner potential to match with the outer potential. This principle can be written as

ϕ = S(x) log r̃ + G(x) as r̃ → 0

must equal

φ = S(x) log r^{∗} + G(x) as r^{∗} → ∞

But ϕ = S(x) log r̃ + G(x) = S(x) log r^{∗} + G(x) + S(x) log(δ^{2}M_{∞}) so for the two expressions of the first terms to be identically equal a “switch-back” term needs to be added or subtracted from one of the expansions. In order to match to higher order, intermediate regions may be needed. The switch-back term is added to the inner expansion. Thus the correct inner expansion for the perturbation potential Φ involving φ(x, r^{∗}) is

Φ = U_{∞}{x + δ^{2}S(x) log(δ^{2}M_{∞}) + δ^{2}φ(x, r^{∗}) + …}.  (B.18)

B.1 Derivation of expression for pressure coefficient on the body

The expression for the pressure coefficient is

C_{p} = \frac{p − p_{∞}}{ρ_{∞}U_{∞}^{2}/2} = \frac{p_{∞}}{ρ_{∞}U_{∞}^{2}/2}(\frac{p}{p_{∞}} − 1).  (B.19)

An expression for p/p_{∞} can be derived from (2.2) and the universal gas law p/ρ^{γ} = constant. That expression is

\frac{p}{p_{∞}} = [1 + \frac{γ − 1}{2}M_{∞}^{2} \{1 − (\frac{Φ_{x}}{U_{∞}} )^{2} − (\frac{Φr̃}{U_{∞}} )^{2}\}]^{γ/(γ−1)}  (B.20)

To calculate the pressure on the body we use the expansion for the inner disturbance potential to substitute into (B.20). Assuming Φ has the form given in (B.8)

\frac{∂Φ}{∂x} = U_{∞} \{1 + δ^{2}S′(x) log(δ^{2}M_{∞}) + δ^{2}φ_{x}(x, r^{∗}) + …\}

(B.21)

\frac{∂Φ}{∂r̃} = U_{∞} \{δφ_{r^{∗}}(x, r^{∗}) + …\}  (B.22)

Substituting the above expressions in (B.20) produces

\frac{p}{p_{∞}} = 1 + \frac{γM_{∞}^{2}}{2} \{−2δ^{2}S′(x) log(δ^{2}M_{∞}) − 2δ^{2}φ_{x} − δ^{2}φ_{r^{∗}}^{2} + …\}  (B.23)

where

C_{p}(x, r^{∗}) = −2δ^{2}φ_{x} − δ^{2}φ_{r^{∗}}^{2} − 2δ^{2}S′(x) log(δ^{2}M_{∞})  (B.24)

But the pressure coefficient on the body can be calculated by substituting the fact that φ = S(x) log r^{∗} and evaluate on r^{∗} = F (x) C_{p}(x, r^{∗} = F (x)) = −2δ^{2}S′(x) log(δ^{2}M_{∞})−2δ^{2}G′(x)−2δ^{2}S′(x) log F (x)−δ^{2}(F′(x))^{2}

(B.25)

This is equivalent to the expression given in (2.12).

APPENDIX C Incompressible Details

A number of insights into the solution of the elliptic part of the compressible problem can be gained by investigating the incompressible version of the problem. The incompressible problem is a simplified version of the full problem, but many of the important features remain intact. The domain of interest is the same, though in the incompressible case there is no hyperbolic region. The governing equation is elliptic everywhere. It is obtained by setting ϵ = 0 in (3.2) to produce

(\frac{R_{u}}{2R})_{u} + R_{νν} = 0.  (C.1)

Some of the initially unknown details of the compressible problem which are only determined after solving the boundary value problem are now known explicitly. Expressions for the dipole strength D and the function G(x) can be written assuming that F (x; α) has the form defined in 4.2. In the case of D the D_{flow} component is now zero, so that

D = D_{body} = ∫_{−1}^{1} F^{2}(x) dx.  (C.2)

For G(x) the explicit expression is

G(x) = −\frac{1}{2} ∫_{−1}^{+1} S′(ξ) sgn(x − ξ) log 2|x − ξ| dξ.  (C.3)

Using the definition of S′(x; α) from (4.4) the corresponding expression for G′(x; α) can be written by evaluating the integral in (C.3):

G′(x; α) = S′(x; α) log[4(1 − x^{2})] − \frac{3}{4}(4 − 3α^{2}) − 22αx + \frac{3}{4}(12 − 29α^{2})x^{2}
+\frac{110}{3}αx^{3} + \frac{625}{4}α^{2}x^{4}.  (C.4)

Also, using the definition of F (x; α) in (4.2) one can evaluate the integral in (C.2) to obtain an exact value for D_{body}, namely

D_{body} = \frac{16π}{105} (7 + α^{2}).  (C.5)