Hypersonic boundary layer theory in the symmetry plane of blunt bodies

Abstract

Solving the three-dimensional boundary layer equations carries theoretical significance and practical applications, which also poses substantial challenges due to its inherent complexity. In this paper, the laminar boundary layer equations for the symmetry plane of three-dimensional bodies are derived in an orthogonal curvilinear coordinate system associated with the principal curvatures. The derivation of the boundary layer equations is based not only on the common symmetric properties of the flow, as given by Hirschel et al. (Three-Dimensional Attached Viscous Flow, 2014, Academic Press, pp. 183–187), but also incorporates the geometric symmetry properties of the body. The derived equations are more representative and simplified. Notably, these equations can degenerate to a form consistent with or equivalent to the commonly used boundary layer equations for special bodies such as flat plates, cones and spheres. Furthermore, for hypersonic flows, the crossflow velocity gradient at the boundary layer edge on the symmetry plane is derived based on Newtonian theory. Subsequently, this parameter can provide the necessary boundary condition needed for solving the boundary layer equations using existing methods. Finally, as examples, the equations developed in this paper are solved using the difference-differential method for several typical three-dimensional blunt shapes that appeared on hypersonic vehicles. They prove to be useful in the analysis and interpretation of boundary layer flow characteristics in the symmetry plane of blunt bodies.

1. Introduction

Symmetry is ubiquitous in this world, and almost all hypersonic flight vehicles, such as HTV-2, HIFiRE-5, X-43 and X-37B, also have bilateral symmetry. Conceivably, symmetry plane solutions are usually the most representative parts of complete solutions over an entire body and are frequently needed for practical purposes. Knowing these symmetry plane solutions, a clear picture of the complete flow may immediately follow. Although the contributions of modern computational fluid dynamics (CFD) have demonstrated significant progress over the past two decades, simulation of the symmetry plane flow characteristics usually requires computationally expensive three-dimensional calculations, which remain full of complexities, especially for the prediction of aerodynamic heating and skin friction in hypersonic flows (Candler et al. 2007; Mazaheri & Kleb 2007; Candler, Subbareddy & Brock 2015). Thus, solving the relatively simple boundary layer equations in a plane of symmetry is still attractive for theoretical analysis and can provide accurate results to verify other three-dimensional calculations. Patel & Baek (1987) have also reviewed the practical importance of plane-of-symmetry boundary layers.

Indeed, boundary layer theories have received considerable attention over the last hundred years or so. However, due to the complexities of fully three-dimensional flows, most boundary layer theoretical analyses have been devoted to two-dimensional or axially symmetric flows, such as flat plates or stagnation point flows. In the most general case of three-dimensional flow, the boundary layer cannot be analysed except by cumbersome numerical methods. Present available three-dimensional methods for hypersonic flows, such as HABP (Jain & Hayes 2004), MINIVER (Engel & Praharaj 1983), INCHES (Zoby & Simmonds 1985) and aeroheating (Hamilton, Dejarnette & Weilmuenster 1987), are engineering methods, and they are overly approximate and restrictive for general use. In this paper the boundary layer along the plane of symmetry for three-dimensional bodies is of considerable interest, where the flow can provide the boundary conditions to be satisfied by the flow on either side of it and, consequently, dictates the behaviour of the boundary layer over the entire body. The complexity of the plane-of-symmetry boundary layer is intermediate to that of the idealized two-dimensional boundary layer, on the one hand, and the fully three-dimensional boundary layer, on the other. It differs from the former in the presence of lateral flow convergence or divergence and from the latter in the absence of crossflow. Perhaps the best-known example thereof is the boundary layer on either the windward or the leeward side of a body of revolution, such as cones and spheroids (Cebeci & Su 1988). For these shapes, scholars have obtained specific boundary layer equations on a plane of symmetry based on their coordinate systems and solved the flow parameters for laminar or turbulent symmetric flows (Murdock 1972; Xu & Wang 1988). Even in the case of more generalized boundary layer equations, it is also not difficult to derive the boundary layer equations in symmetry planes under the corresponding coordinate system from fully three-dimensional boundary layer equations under the orthogonal curvilinear coordinate system (Hirschel, Cousteix & Kordulla 2014; Schlichting & Gersten 2016). A body-oriented coordinate system is generally used that has the advantage of being independent of the angle of attack and is easy to calculate even if the body is not defined analytically.

Deriving the boundary layer equations becomes manageable when a well-suited coordinate system is chosen. However, the challenge lies in closing the equations and achieving an exact solution. The system of equations by Schlichting & Gersten (2016) for the principal terms of the crossflow seem to be unclosed, as they contain terms with a longitudinal component of the peripheral velocity gradient, i.e. ${\partial w}/{\partial z}$, which must be taken into account in order to obtain an accurate description of the flow. Unfortunately, the crossflow leads to complexity in the boundary layer solutions at the symmetry plane. Currently, there are no general theoretical expressions for the crossflow velocity in the symmetry plane of three-dimensional bodies.

Previous studies have typically focused on specific shape flows, such as the symmetry plane of a spheroid (Wang 1970), a cone (Wu & Libby 1972) or a sharp cone (Murdock 1972), where the longitudinal velocity gradient can be approximately deduced based on these particular model shapes. For a universal symmetry plane, this crossflow has usually been simplified within the framework of ‘small crossflow’ or ‘axisymmetric analog’ approaches traditionally taken in integral-type formulations (Hecht 1976). For the small crossflow velocity assumption in the inviscid flow at the boundary layer edge (Cooke 1961), it is interpreted that the transverse curvature perpendicular to the symmetry plane is small. This condition imposes severe restrictions on the geometry of the body and the angles of attack. Adams (1972) evaluated the results of ‘effective cone’ calculations compared with a three-dimensional windward meridian analysis over a sharp cone. Agnone (1974) also developed a crossflow correction factor for the effective cone technique for the plane of symmetry of a body at the angle of attack. However, these methods involve variables related to the flow characteristics, and they may lead to underprediction of heat transfer in certain situations (Jain & Hayes 2004). Furthermore, except for relatively simple bodies such as wedges, cones, spheres and sphere–cones, it is also reasonably challenging to guess the correct form of pressure distribution on bodies of complex configuration and calculate a streamline path and metric coefficient for an equivalent axisymmetric body (Jain & Hayes 2004). Hence, the range of application of the above methods is limited in scope and they fail to yield unique results.

The present investigation aims at carrying out boundary layer studies in the symmetry plane of three-dimensional bodies and providing more general theoretical solutions for hypersonic boundary layer equations. No restrictions on the angle of attack or the ‘small crossflow’ assumption are imposed in deducing the equations, nor are there invoked any of the common assumptions such as similarity, conical flow and others. Consequently, the derived equations can provide a broader range of applications for the symmetry plane of three-dimensional bodies. The paper is organized as follows. Three-dimensional laminar boundary layer equations in orthogonal curvilinear coordinates are presented in § 2, and simplified equations in the symmetry plane are derived in § 3 based on a special orthogonal curvilinear coordinate system relevant to the principal curvatures. Transformation of the boundary layer equations using the Levy–Lees method is shown in § 4, and remarks on these equations are given in § 5. Based on Newtonian theory, the crossflow velocity gradient at the boundary layer edge in the symmetry plane for hypersonic flows is derived in § 6. Estimation of other boundary layer edge parameters and numerical solution of the equations are also included in § 6. In § 7, several typical blunt bodies are selected to show how the present equations are applicable to obtain the flow characteristics in the symmetry plane quickly. § 8 concludes the paper.

2. General three-dimensional laminar boundary layer equations under the body-oriented coordinate system

The choice of a coordinate system to describe the boundary layer equations over a three-dimensional body is of major importance and significantly facilitates analysing and solving these equations. In addition, the boundary layer approximation is impractical for implementation in any system other than a body-oriented coordinate system, which is independent of the angle of attack and takes into account the shape of the body surface under consideration. Here, a body-oriented orthogonal curvilinear coordinate system is chosen for the three-dimensional boundary layer equations. This coordinate system consists of two curvilinear coordinates, $x$ and $z$, locally parallel to the surface, and a third coordinate, $y$, locally normal to the surface, as illustrated in figure 1, where the Cartesian coordinate system ($X, Y, Z$) is used as the reference.

Figure 1. Orthogonal curvilinear coordinate system for a three-dimensional boundary layer. Here $(x, y, z)$ represents the body-oriented coordinate system and $(X, Y, Z)$ represents the Cartesian coordinate system.

The initial step in deriving boundary layer equations in a symmetry plane involves referencing the general three-dimensional boundary layer equations under the body-oriented coordinate system, which can be found in various references (Cebeci, Kaups & Moser 1976; Yousefi & Veron 2020). The laminar boundary layer equations for a steady compressible flow under a body-oriented orthogonal coordinate system can be expressed as follows: for continuity,

(2.1)\begin{equation} \frac{\partial}{\partial x}({{H_z}} \rho u)+\frac{\partial}{\partial z}({H_x} \rho w)+\frac{\partial}{\partial y}({H_x}{H_z} \rho v)=0; \end{equation}

for x-directional momentum,

(2.2a) \begin{align} \rho\left(\frac{u}{{H_x}} \frac{\partial u}{\partial x}+\frac{w}{{H_z}} \frac{\partial u}{\partial z}+v \frac{\partial u}{\partial y}+\frac{u w}{{H_x}{H_z}} \frac{\partial{H_x}}{\partial z}-\frac{w^2}{{H_x}{H_z}} \frac{\partial{H_z}}{\partial x}\right)={-}\frac{1}{{H_x}} \frac{\partial p}{\partial x}+\frac{\partial}{\partial y}\left(\mu \frac{\partial u}{\partial y}\right); \end{align}

for z-directional momentum,

(2.2b)\begin{align} \rho\left(\frac{u}{{H_x}} \frac{\partial w}{\partial x}+\frac{w}{{H_z}} \frac{\partial w}{\partial z}+v \frac{\partial w}{\partial y}-\frac{u^2}{{H_x}{H_z}} \frac{\partial{H_x}}{\partial z}+\frac{u w}{{H_x}{H_z}} \frac{\partial{H_z}}{\partial x}\right)={-}\frac{1}{{H_z}} \frac{\partial p}{\partial z}+\frac{\partial}{\partial y}\left(\mu \frac{\partial w}{\partial y}\right); \end{align}

for energy,

(2.3) \begin{align} \rho\left(\frac{u}{{H_x}} \frac{\partial h}{\partial x}+\frac{w}{{H_z}} \frac{\partial h}{\partial z}+v \frac{\partial h}{\partial y}\right)=\frac{u}{{H_x}} \frac{\partial p}{\partial x}+\frac{w}{{H_z}} \frac{\partial p}{\partial z}+\frac{\partial}{\partial y}\left(\frac{\mu}{P r} \frac{\partial h}{\partial y}\right)+\mu\left(\frac{\partial u}{\partial y}\right)^2+\mu\left(\frac{\partial w}{\partial y}\right)^2, \end{align}

where u, v and w are the velocity components in the x, y and z directions, respectively; $\rho$ is the density, p is the pressure, h is the static enthalpy, $\mu$ is the coefficient of viscosity and $Pr$ is the Prandtl number. Here, ${H_x}$ and ${H_z}$ denote the metric coefficients, which can be obtained from the relationships between (X, Y, Z) and (x, y, z). In accordance with conventional boundary layer theory, the metric coefficient in the normal direction ${H_y}$ is taken equal to unity, and ${H_x}$ and ${H_z}$ are assumed to be independent of y. Thus, ${H_x}$ and ${H_z}$ are given as

(2.4a,b)\begin{align} {H_x}=\sqrt{\left(\frac{\partial X}{\partial x}\right)^2+\left(\frac{\partial Y}{\partial x}\right)^2+\left(\frac{\partial Z}{\partial x}\right)^2}, \quad {H_z}=\sqrt{\left(\frac{\partial X}{\partial z}\right)^2+\left(\frac{\partial Y}{\partial z}\right)^2+\left(\frac{\partial Z}{\partial z}\right)^2}. \end{align}

At the boundary layer edge, the pressure gradients in (2.2) must be matched to the flow outside the boundary layer, which are given by Euler equations for inviscid flow. Thus, (2.2a) and (2.2b) reduce to

(2.5) \begin{equation} \left.\begin{array}{c@{}} \displaystyle -\dfrac{1}{{H_x}} \dfrac{\partial p_e}{\partial x}=\rho_e\left(\dfrac{u_e}{{H_x}} \dfrac{\partial u_e}{\partial x}+\dfrac{w_e}{{H_z}} \dfrac{\partial u_e}{\partial z}+\dfrac{u_e w_e}{{H_x}{H_z}} \dfrac{\partial{H_x}}{\partial z}-\dfrac{w_e^2}{{H_x}{H_z}} \dfrac{\partial{H_z}}{\partial x}\right), \\ \displaystyle -\dfrac{1}{{H_z}} \dfrac{\partial p_e}{\partial z}=\rho_e\left(\dfrac{u_e}{{H_x}} \dfrac{\partial w_e}{\partial x}+\dfrac{w_e}{{H_z}} \dfrac{\partial w_e}{\partial z}+\dfrac{u_e w_e}{{H_x}{H_z}} \dfrac{\partial{H_z}}{\partial x}-\dfrac{u_e^2}{{H_x}{H_z}} \dfrac{\partial{H_x}}{\partial z}\right), \end{array}\right\} \end{equation}

where subscript e denotes the outer edge of the boundary layer.

The boundary conditions for the above equations are similar to those for the two-dimensional flow. That is, at the wall,

(2.6a)\begin{equation} y = 0, \quad u = v = w = 0, \quad T = T_w \text{ or } \left(\frac{\partial T}{\partial n}\right)_w = 0 \quad \text{(adiabatic wall),} \end{equation}

where subscript w denotes the wall. And at the boundary layer edge,

(2.6b)\begin{equation} y \rightarrow \infty, \quad u \rightarrow u_e, \quad w \rightarrow w_e, \quad T \rightarrow T_e. \end{equation}

Notably since the boundary layer thickness is not known a priori, the boundary condition at the edge of the boundary layer is given at large y, essentially y approaching infinity.

3. Boundary layer equations in the symmetry plane

We investigate the three-dimensional flow over three-dimensional bodies with a symmetry plane, and the vector of oncoming flow lies within that plane. The body-oriented orthogonal coordinate system ($x, y, z$) is applied in the present work. The schematic of the flow in the vicinity of the symmetry plane is shown in figure 2, where $R$ and $r$ are the two principal radii of curvatures of the surface, one in the symmetry plane and the other in its perpendicular plane, respectively. In this coordinate system, the plane of symmetry is taken as the plane of $z=0$, and the coordinate $z$ is aligned on the surface along the circumferential direction corresponding to the principal radius $r$. The coordinate $y$ is directed along the line normal to the body surface, and the coordinate $x$ is locally orthogonal to the other two coordinates. Differing from the streamline coordinate system, where the coordinate is formed by the inviscid streamlines and the metric coefficients must be calculated at every change of attitude of the body or Mach number for compressible flows, the proposed coordinates are related to the principal directions or the eigenvalues of the shape operator at a point in differential geometry. This coordinate system is independent of incoming flows and the metric coefficients are calculated only once.

Figure 2. Schematic geometry in the vicinity of the symmetry plane flows.

To obtain the metric coefficients appearing in the boundary layer equations, relationships between ($x$, $y$, $z$) and ($X, Y, Z$) are required. At any given surface point O on the symmetry plane in figure 2, its absolute coordinate values are independent of coordinate transformation for the metric coefficients. Thus, point O is still defined as the origin coordinates in the following analysing process. Points C and D are the centres of $r$ and $R$, respectively, i.e. $\textrm {CO}=r$ and $\textrm {DO}=R$ (figure 2). Points B and F are in the vicinity of the symmetry plane, and point E is the centre of the principal radii of point B in the $x$ direction. From the coordinate's illustration of these points, as shown in figure 2, the transformation between the proposed curvilinear coordinates ($x, y, z$) and the Cartesian coordinates ($X, Y, Z$) of point F are expressed as

(3.1a)$$\begin{gather} X=[R-r(1-\cos \varphi)] \sin \alpha, \quad \alpha=\frac{x}{R-r(1-\cos \varphi)}, \end{gather}$$

(3.1b)$$\begin{gather}Y=y, \end{gather}$$

(3.1c)$$\begin{gather}Z=r \sin \varphi, \quad \varphi=\frac{z}{r}, \end{gather}$$

where $\varphi$ is the angle between lines CO and CB, and $\alpha$ is the angle between EB and EF.

Because the plane of $z=0$ is the symmetry plane, the two principal radii of curvature of the surface R and r can be expanded with respect to the z coordinate as

(3.2a)$$\begin{gather} R(x, z)=R(x, 0)+\frac{1}{2} z^2 \frac{\partial^2 R}{\partial z^2}+\cdots, \end{gather}$$

(3.2b)$$\begin{gather}r(x, z)=r(x, 0)+\frac{1}{2} z^2 \frac{\partial^2 r}{\partial z^2}+\cdots. \end{gather}$$

Approximating these two equations to the first order in $z$ yields

(3.3)\begin{equation} \left.\begin{array}{c@{}} R(x, z) =R(x, 0), \\ r(x, z) =r(x, 0). \end{array}\right\} \end{equation}

Hence, $R$ and $r$ are a function of $x$ only.

To obtain the derivatives required for solving the metric coefficients, we let point F infinitely approach point O, i.e. $x \rightarrow 0$ and $z \rightarrow 0$, or $\varphi \rightarrow 0$ and $\alpha \rightarrow 0$. Thus, $\cos \alpha \approx 1$, $\cos \varphi \approx 1$, $\sin \alpha \approx {x}/{R}$ and $\sin \varphi \approx {z}/{R}$. By utilizing these equations and keeping the first truncation about z or x, the following equations can be obtained from (3.1) and (3.3):

(3.4)\begin{equation} \frac{\partial X}{\partial x}=1, \quad \frac{\partial Z}{\partial x}=0, \quad \frac{\partial X}{\partial z}=0, \quad \frac{\partial Z}{\partial z}=1, \quad \frac{\partial^2 X}{\partial x \partial z}={-}\frac{z}{R r}, \quad \frac{\partial^2 Z}{\partial x \partial z}=0. \end{equation}

Substituting these equations into (2.4a,b) yields the metric coefficients as

(3.5a)\begin{equation} {H_x}=1, \quad{H_z}=1. \end{equation}

Accordingly, the curvature parameters can be obtained as

(3.5b)\begin{equation} \frac{\partial{H_x}}{\partial z}={-}\frac{z}{R r}, \quad \frac{\partial{H_z}}{\partial x}=0. \end{equation}

So far, in the symmetry plane, all the metric coefficients and the curvature parameters needed in the boundary layer equations are obtained, consisting of lines of principal curvatures of the body surface.

Plane-of-symmetry flows are characterized by the fact that $w=0$ holds and all derivatives with respect to $z$ vanish, apart from the term ${\partial w}/{\partial z} \neq 0$. Then, the flow variables can be expanded into a Taylor series about the symmetry plane with respect to the $z$ coordinate as

(3.6a)$$\begin{gather} u(x, y, z)=u(x, y, 0)+\frac{1}{2} z^2 \frac{\partial^2 u}{\partial z^2}+\cdots, \end{gather}$$

(3.6b)$$\begin{gather}v(x, y, z)=v(x, y, 0)+\frac{1}{2} z^2 \frac{\partial^2 v}{\partial z^2}+\cdots, \end{gather}$$

(3.6c)$$\begin{gather}w(x, y, z)=z \frac{\partial w}{\partial z}+\frac{1}{6} z^3 \frac{\partial^3 w}{\partial z^3}+\cdots, \end{gather}$$

(3.6d)$$\begin{gather}\rho(x, y, z)=\rho(x, y, 0)+\frac{1}{2} z^2 \frac{\partial^2 \rho}{\partial z^2}+\cdots, \end{gather}$$

(3.6e)$$\begin{gather}p(x, z)=p(x, 0)+\frac{1}{2} z^2 \frac{\partial^2 p}{\partial z^2}+\cdots. \end{gather}$$

Substituting the above flow quantities of the first truncation and (3.5) into ((2.1)–(2.3)) yields the boundary layer equations for the symmetry plane. Terms with coefficients of z or higher orders of z are ignored in the resulting expressions. Additionally, $p(x, y)=p_e(x, y)$ is assumed throughout the boundary layer. Notably, the independent variables are functions of x and y only in the subsequent equations. The derived continuity, $x$-directional momentum, $z$-directional momentum and energy equations can then be written as

(3.7)\begin{gather} \frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial}{\partial z}\left(\rho z \frac{\partial w}{\partial z}\right)=0, \end{gather}

(3.8a)$$\begin{gather} \rho u \frac{\partial u}{\partial x}+\rho v \frac{\partial u}{\partial y}={-}\frac{\partial p_e}{\partial x}+\frac{\partial}{\partial y}\left(\mu \frac{\partial u}{\partial y}\right), \end{gather}$$

(3.8b)$$\begin{gather}\rho u \frac{\partial}{\partial x}\left(\frac{\partial w}{\partial z}\right)+\rho v \frac{\partial}{\partial y}\left(\frac{\partial w}{\partial z}\right)+\rho\left(\frac{\partial w}{\partial z}\right)^2={-}\frac{\partial^2 p_e}{\partial z^2}+\frac{\partial}{\partial y}\left[\mu \frac{\partial}{\partial y}\left(\frac{\partial w}{\partial z}\right)\right]-\rho \frac{u^2}{R r}, \end{gather}$$

(3.9)\begin{gather} \rho u \frac{\partial h}{\partial x}+\rho v \frac{\partial h}{\partial y}=u \frac{{\rm d} p_e}{{\rm d}\kern0.7pt x}+\frac{\partial}{\partial y}\left(\frac{\mu}{P r} \frac{\partial h}{\partial y}\right)+\mu\left(\frac{\partial u}{\partial y}\right)^2. \end{gather}

Accordingly, the Euler equations (2.6) reduce to

(3.10a)$$\begin{gather} -\frac{\partial p_e}{\partial x}=\rho_e u_e \frac{\partial u_e}{\partial x} , \end{gather}$$

(3.10b)$$\begin{gather}-\frac{\partial^2 p_e}{\partial z^2}=\rho_e\left[u_e \frac{\partial}{\partial x}\left(\frac{\partial w_e}{\partial z}\right)+\left(\frac{\partial w_e}{\partial z}\right)^2+\frac{u_e^2}{R r}\right]. \end{gather}$$

Equations (3.7)–(3.9) are the boundary layer equations for the symmetry plane of a three-dimensional body within the proposed coordinate system. Compared with the general three-dimensional boundary layer equations, these equations are significantly simplified. While the x-directional momentum and energy equations are almost identical to those of planar boundary layers, it is important to note that the flow in the symmetry plane differs substantially from that of a planar boundary layer. Here, there may exist mass flow either into or out of the symmetry plane, a feature represented by the term ${\partial w}/{\partial z}$ in the equations. Notably, although figure 2 presents an illustrative diagram of a blunt body, along with the geometric relationships corresponding to the changes in coordinates, i.e. (3.1). These geometric relationships remain valid when dealing with concave surfaces, which can be understood as the inner surface of the blunt body depicted in figure 2. Meanwhile, the flow symmetry conditions in (3.6) are independent of whether the surface is blunt or concave. Thus, ((3.7)–(3.9)) are equally applicable to concave surfaces. Besides, the boundary layer equations provided above are derived from the general three-dimensional boundary layer equations for compressible flows. And for the specific symmetric problems addressed in this paper, both the flow's symmetric properties and the body's geometric symmetry properties are inherent, independent of the flow's Mach number. Hence, these equations are applicable across compressible flow regimes, including subsonic, supersonic and hypersonic flows, without differentiation. It is important to highlight that under hypersonic flow conditions, if chemical reactions are considered within the boundary layer, the diffusion of chemical species and energy changes resulting from chemical reactions must be incorporated into the aforementioned equations.

4. Transformation of the boundary layer equations

As in two-dimensional flows, expressing the boundary layer equations in transformed coordinates is more useful and convenient to yield instructive results. Several transformations can be used for this purpose. The one adopted herein, similar to the Levy–Lees transformation widely utilized in two-dimensional and axisymmetric flows, is

(4.1a)$$\begin{gather} \xi=\int_0^x \mu_e \rho_e u_e \,{{\rm d}\kern0.7pt x}, \end{gather}$$

(4.1b)$$\begin{gather}\eta=\frac{u_e}{\sqrt{2 \xi}} \int_0^y \rho\, {{\rm d} y}, \end{gather}$$

(4.1c)$$\begin{gather}z=z, \end{gather}$$

where $\rho _e$, $u_e$ and $\mu _e$ are the density, velocity and viscosity coefficients at the edge of the boundary layer, respectively. Since $\rho _e$, $u_e$ and $\mu _e$ are functions of x only in the symmetry plane, then $\xi =\xi (x)$. From the chain rule of differential calculus, the derivative transformations can be expressed as

(4.2a)$$\begin{gather} \frac{\partial}{\partial x}=\frac{\partial \xi}{\partial x} \frac{\partial}{\partial \xi}+\frac{\partial \eta}{\partial x} \frac{\partial}{\partial \eta}=\rho_e u_e \mu_e \frac{\partial}{\partial \xi}+\frac{\partial \eta}{\partial x} \frac{\partial}{\partial \eta}, \end{gather}$$

(4.2b)$$\begin{gather}\frac{\partial}{\partial y}=\frac{\partial \eta}{\partial y} \frac{\partial}{\partial \eta}=\frac{u_e}{\sqrt{2 \xi}} \rho \frac{\partial}{\partial \eta}. \end{gather}$$

To reduce the number of dependent variables appearing in the boundary layer equations, it is customary to introduce the streamfunctions $\psi$ and $\phi$, which identically satisfy the equation of continuity (3.7); they are defined as usual by

(4.3a)$$\begin{gather} \rho u=\frac{\partial \psi}{\partial y}, \end{gather}$$

(4.3b)$$\begin{gather}\rho z \frac{\partial w}{\partial z}=\frac{\partial \phi}{\partial y}, \end{gather}$$

(4.3c)$$\begin{gather}\rho v={-}\left(\frac{\partial \psi}{\partial x}+\frac{\partial \phi}{\partial z}\right). \end{gather}$$

As (3.6) indicate, near the symmetry plane, $w$ is dependent on the first-order and higher-order derivatives of $z$. In the following discussions, only the first truncation about $z$ is kept. Thus, the velocities $u$ and ${\partial w}/{\partial z}$ are functions of $x$ and $y$; hence, they are expressed as functions of $\xi$ and $\eta$, such that $u=u(\xi, \eta )$ and ${\partial w}/{\partial z}={\partial w(\xi, \eta )}/{\partial z}$. At the edge of the boundary layer, $u_e$ and ${\partial w_e}/{\partial z}$ are a function of $\xi$ only, such that $u_e=u_e(\xi )$ and ${\partial w_e}/{\partial z}={\partial w_e(\xi )}/{\partial z}$.

Let us define two functions of $\xi$ and $\eta$, $f(\xi, \eta )$ and $e(\xi, \eta )$, such that

(4.4a)$$\begin{gather} \frac{u}{u_e}=\frac{\partial f}{\partial \eta} \equiv f^{\prime}, \end{gather}$$

(4.4b)$$\begin{gather}\frac{\dfrac{\partial w}{\partial z}}{\dfrac{\partial w_e}{\partial z}}=\frac{\partial e}{\partial \eta} \equiv e^{\prime}. \end{gather}$$

where the prime denotes the partial derivative with respect to $\eta$. The relationship between $f$ and $\psi$ is the same as that in the two-dimensional boundary layer flows, such that

(4.5)\begin{equation} \psi=\sqrt{2 \xi} f(\xi, \eta). \end{equation}

From (4.3b), written in terms of the transformation given in (4.2b) and (4.4b), we have

(4.6a)\begin{equation} \frac{u_e}{\sqrt{2 \xi}} \rho \frac{\partial \phi}{\partial \eta}=\rho z \frac{\partial w}{\partial z}=\rho z e^{\prime} \frac{\partial w_e}{\partial z}\end{equation}

or

(4.6b)\begin{equation} \frac{\partial \phi}{\partial \eta}=\frac{\sqrt{2 \xi}}{u_e} z e^{\prime} \frac{\partial w_e}{\partial z}. \end{equation}

Integrating (4.6) with respect to $\eta$, we have

(4.7)\begin{equation} \phi=\frac{\sqrt{2 \xi}}{u_e} e(\xi, \eta) z{\frac{\partial w_e}{\partial z}}. \end{equation}

Substituting (4.5) and (4.7) into (4.3) yields

(4.8a)$$\begin{gather} \rho u=\frac{\partial \psi}{\partial y}=\rho u_e f^{\prime}, \end{gather}$$

(4.8b)$$\begin{gather}\rho w=\frac{\partial \phi}{\partial y}=\rho{\frac{\partial w_e}{\partial z}} z e^{\prime}, \end{gather}$$

(4.8c)$$\begin{gather}\rho v={-}\left(\frac{\partial \psi}{\partial x}+\frac{\partial \phi}{\partial z}\right)={-}\left[\mu_e \rho_e u_e\left(\frac{f}{\sqrt{2 \xi}}+\sqrt{2 \xi} \frac{\partial f}{\partial \xi}\right)+\frac{\partial \eta}{\partial x} \sqrt{2 \xi} f^{\prime}+\frac{\sqrt{2 \xi}}{u_e}{\frac{\partial w_e}{\partial z}} e\right]. \end{gather}$$

Defining the non-dimensional static enthalpy as

(4.9)\begin{equation} g=g(\xi, \eta)=\frac{h}{h_e}, \end{equation}

where $h_e$ is the static enthalpy at the boundary layer edge. Then, by substituting ((4.1)–(4.9)) into ((3.8) and (3.9)), the transformed boundary layer equations for the symmetry plane of three-dimensional bodies can be written as follows: for $\xi$-directional momentum,

(4.10a)\begin{equation} (C f^{\prime \prime})^{\prime}+ \left(f+\frac{2 \xi}{u_e} \frac{{\dfrac{\partial w_e}{\partial z}}}{u_e \mu_e \rho_e} e\right) f^{\prime \prime}=\frac{2 \xi}{u_e} \frac{\partial u_e}{\partial \xi}\left(f^{\prime 2}-\frac{\rho_e}{\rho}\right)+2 \xi f^{\prime} \frac{\partial f^{\prime}}{\partial \xi}-2 \xi \frac{\partial f}{\partial \xi} f^{\prime \prime}; \end{equation}

for $z$-directional momentum,

(4.10b)\begin{align} &\left(C e^{\prime \prime}\right)^{\prime}+\left[f+\frac{2 \xi}{u_e} \frac{{\dfrac{\partial w_e}{\partial z}}}{u_e \mu_e \rho_e} e\right] e^{\prime \prime} =\frac{2 \xi}{{\dfrac{\partial w_e}{\partial z}}} \frac{\partial\left(e^{\prime}{\dfrac{\partial w_e}{\partial z}}\right)}{\partial \xi} f^{\prime}\nonumber\\ &\quad +\frac{2 \xi}{u_e} \frac{{\dfrac{\partial w_e}{\partial z}}}{u_e \mu_e \rho_e}\left(e^{\prime 2}-\frac{\rho_e}{\rho}-\frac{u_e^2}{\left({\dfrac{\partial w_e}{\partial z}}\right)^2} \frac{1}{R r} \frac{\rho_e}{\rho}+\frac{u_e^2}{\left({\dfrac{\partial w_e}{\partial z}}\right)^2} \frac{1}{\rho R r} f^{\prime 2}\right)-2 \xi \frac{\partial f}{\partial \xi} e^{\prime \prime}; \end{align}

for energy,

(4.10c)\begin{align} \left(\frac{C}{P r} g^{\prime}\right)^{\prime}+\left(f+\frac{2 \xi}{u_e} \frac{{\dfrac{\partial w_e}{\partial z}}}{u_e \mu_e \rho_e} e\right) g^{\prime}=2 \xi\left(f^{\prime} \frac{\partial g}{\partial \xi}-g^{\prime} \frac{\partial f}{\partial \xi}+\frac{f^{\prime} g}{h_e} \frac{\partial h_e}{\partial \xi}+\frac{\rho_e u_e}{\rho h_e} f^{\prime} \frac{\partial u_e}{\partial \xi}\right)-C \frac{u_e^2}{h_e}\left(f^{\prime \prime}\right)^2, \end{align}

where

(4.11)\begin{equation} {Pr}=\frac{\mu_e c_p}{k}, \quad C=\frac{\rho \mu}{\rho_e \mu_e}, \end{equation}

$k$ is the thermal conductivity and $C$ is the ‘rho-mu’ ratio.

Accordingly, the transformed boundary layer conditions are as follows. At the wall,

(4.12a)\begin{equation} \eta=0,\ f=f^{\prime}=0,\ e=e^{\prime}=0,\ g=g_w \text{ (fixed wall temperature) or } g^{\prime}=0 \text{ (adiabatic wall),} \end{equation}

and at the boundary layer edge,

(4.12b)\begin{equation} \eta \rightarrow \infty, \quad f^{\prime}=1, \quad e^{\prime}=1, \quad g=1. \end{equation}

Equations (4.10) are the transformed compressible boundary layer equations for the symmetry plane. Their derivation and analysis draw partial reference from the planar boundary layer equations as studied by Anderson (2006). Apart from $f$ and $g$ that exist in the planar boundary layer equations, another variable $e$ is newly introduced in (4.10). The terms containing ${\partial w_e}/{\partial z}$ represent the convergence or divergence of the streamlines close to the symmetry plane. If ${\partial w_e}/{\partial z}=0$, these equations directly reduce to the same form as the planar boundary layer equations. Furthermore, (4.10) contain no further approximations or assumptions beyond those associated with the original boundary layer equations, namely ((3.7)–(3.9)).

5. Remarks on the equations

Both three-dimensional stagnation points and axisymmetric flows are regarded as special cases of more general symmetric flows, where all the above equations remain applicable. These equations can effectively degenerate into forms consistent with or equivalent to commonly used equations. We will discuss these two special cases in this section.

5.1. Stagnation point flows

Regarding the general three-dimensional stagnation point flow, the velocity at the boundary layer edge in the $x$ direction can be written as

(5.1)\begin{equation} u_e=\left(\frac{\partial u_e}{\partial x}\right)_s x, \end{equation}

where $({\partial u_e}/{\partial x})_s$ is the velocity gradient at the stagnation point external to the boundary layer. The subscript $s$ represents the stagnation point. Substituting (5.1) into (4.1a), we obtain

(5.2)\begin{equation} \xi=\int_0^x \mu_e \rho_e{\left(\frac{\partial u_e}{\partial x}\right)_s} x\, {{\rm d}\kern0.7pt x}=\frac{1}{2} \mu_e \rho_e{\left(\frac{\partial u_e}{\partial x}\right)_s} x^2. \end{equation}

Consider a self-similar solution to the governing boundary layer equations for the stagnation point case, $f$, $e$ and $g$ are assumed to be functions of $\eta$ only; hence, ${\partial f}/{\partial \xi }={\partial f^{\prime }}/{\partial \xi }={\partial e}/{\partial \xi }={\partial e^{\prime }}/{\partial \xi }={\partial g}/{\partial \xi }=0$ in (4.10). Substituting (5.1) and (5.2) into (4.10), and neglecting the terms with higher orders of $x$, the final equations for the general three-dimensional stagnation point flow are written as follows: for $\xi$-directional momentum,

(5.3a)\begin{equation} (C f^{\prime \prime})^{\prime}+\left(f+\left(\frac{{\dfrac{\partial w_e}{\partial z}}}{\dfrac{\partial u_e}{\partial x}}\right)_s e\right) f^{\prime \prime}=f^{\prime 2}-\frac{\rho_e}{\rho}; \end{equation}

for $z$-directional momentum,

(5.3b)\begin{equation} (C e^{\prime \prime})^{\prime}+\left(f+\left(\frac{{\dfrac{\partial w_e}{\partial z}}}{\dfrac{\partial u_e}{\partial x}}\right)_s e\right) e^{\prime \prime}=\left(\frac{{\dfrac{\partial w_e}{\partial z}}}{\dfrac{\partial u_e}{\partial x}}\right)_s\left(f^{\prime 2}-\frac{\rho_e}{\rho}\right); \end{equation}

for energy,

(5.3c)\begin{equation} \left(\frac{C}{P r} g^{\prime}\right)^{\prime}+\left(f+\left(\frac{{\dfrac{\partial w_e}{\partial z}}}{\dfrac{\partial u_e}{\partial x}}\right)_s e\right) g^{\prime}=0. \end{equation}

Equations (5.3) are totally identical to the general three-dimensional stagnation point flow derived by other authors (Libby 1967). It can also be used directly as a starting point for solving the boundary layer equations over the whole body with a difference-differential method.

Consider special cases: if $(({{\partial w_e}/{\partial z}})/({{\partial u_e}/{\partial x}}))_s=0,$ (5.3) reduce directly to represent the planar stagnation point flow, where the equations are as follows: for $\xi$-directional momentum,

(5.4a)\begin{equation} (C f^{\prime \prime})^{\prime}+f f^{\prime \prime}=f^{\prime 2}-\frac{\rho_e}{\rho}; \end{equation}

for energy,

(5.4b)\begin{equation} \left(\frac{C}{P r} g^{\prime}\right)^{\prime}+f g^{\prime}=0. \end{equation}

Since the flows in the z direction are uninteresting in this case, (5.3b) is inessential and is dropped here directly. The stagnation flow for the two-dimensional cylinder is the most special case. Equations (5.4) for the two-dimensional cylinder stagnation flow are totally the same as those deduced by other authors (Anderson 2006).

If $(({{\partial w_e}/{\partial z}})/({{\partial u_e}/{\partial x}}))_s=1$, (5.3) reduce to represent the axisymmetric stagnation flow, with a sphere stagnation point flow as the most special case. Then, the equations are given as follows: for $\xi$-directional momentum,

(5.5a)\begin{equation} (C f^{\prime \prime})^{\prime}+(f+e) f^{\prime \prime}=f^{\prime 2}-\frac{\rho_e}{\rho};\end{equation}

for $z$-directional momentum,

(5.5b)\begin{equation} (C e^{\prime \prime})^{\prime}+(f+e) e^{\prime \prime}=f^{\prime 2}-\frac{\rho_e}{\rho}; \end{equation}

for energy,

(5.5c)\begin{equation} \left(\frac{C}{P r} g^{\prime}\right)^{\prime}+(f+e) g^{\prime}=0. \end{equation}

Further discussions regarding (5.5) will be conducted in the next section since it is also a special case of axisymmetric flow.

5.2. Axisymmetric flows

The aforementioned cases of two-dimensional and axisymmetric stagnation point flow can be regarded as specific cases of more generalized symmetry plane flow. In the case of axisymmetric flow where the incoming flow aligns with its axis, the situation becomes slightly more complicated. In the commonly used axisymmetric coordinate system applied for axisymmetric boundary layers, as detailed by Schlichting & Gersten (2016), for example, the velocity and its gradients in the circumference direction are zero. Consequently, the boundary layer equations over an axisymmetric body are similar to those for two-dimensional flows. However, although $w_e$ is also zero in our proposed coordinates, the extra variable ${\partial w_e}/{\partial z}$, which is the gradient of the velocity at the outer edge of the boundary layer along the z direction, is not necessarily zero and is required in advance for solving (4.10). In principle, deriving this variable entails solving (3.10b) with respect to the pressure distribution on the surface. However, for an axisymmetric flow, the situation can be simplified, and this variable can be directly deduced from the kinematic property of the axisymmetric flow, with no requirement of the pressure distribution, as illustrated in figure 3.

Figure 3. Schematic geometry of the axisymmetric boundary layers.

At the stagnation point, the velocity gradients can be calculated from Newtonian theory. Consider an arbitrary point D on the surface of the axisymmetric body, apart from the stagnation point. The velocity at the boundary layer edge of this point is $u_{e, \textrm {D}}$. Imagine a line drawn tangent to the surface at point D and it intersects with the axis of symmetry at point A, making an angle $\theta$ between them, as shown in figure 3. Here CD represents the principal radius of curvature of the surface at point D, i.e. $r$ in the present study. Due to the axisymmetric property, AC also coincides with the symmetry line of the physical body.

To derive ${\partial w_e}/{\partial z}$, the velocity near point D at the $z$ coordinate is required. Assume a nearby point $\textrm {D}^{\prime }$ along the $z$ coordinate, and the angle between lines AD and $\textrm {A} \textrm {D}^{\prime }$ is $\omega$. As $\textrm {D}^{\prime }$ infinitely approaches D, i.e. $\omega \rightarrow$ 0; and from the property of the axisymmetric flow, it can be obtained that the magnitude of the velocity at a nearby point $\textrm {D}^{\prime }$ along the $z$ coordinate is approximately equal to that at point $\textrm {D}, u_{e, \textrm {D}^{\prime }} = u_{e, \textrm {D}}$, which is also consistent with (3.6a). Accordingly, the velocity component in the $z$ direction at point $\textrm {D}^{\prime }$ can be expressed as

(5.6)\begin{equation} w_{e, {\rm D}^{\prime}}={u_{e, {\rm D}^{\prime}} \tan \omega} =u_{e, {\rm D}} \tan \omega . \end{equation}

In addition, the coordinate of the nearby point $\textrm {D}^{\prime }$ is

(5.7)\begin{equation} z=L_{\rm A D} \tan \omega, \end{equation}

where $L_{\textrm {AD}}$ represents the length between points A and D, and $L_{\textrm {AD}}={r}/{\tan \theta }$. Then, substituting (5.7) into (5.6) yields

(5.8a)\begin{equation} w_{e, {\rm D}^{\prime}} = u_{e, {\rm D}} \tan \omega=z \frac{u_{e, {\rm D}} \tan \theta}{r} \end{equation}

or

(5.8b)\begin{equation} \left({\frac{\partial w_e}{\partial z}}\right)_{\rm D}=\frac{u_{e, {\rm D}} \tan \theta}{r}. \end{equation}

With this relation, there is no substantial difficulty in solving the boundary layer equations for axisymmetric flows. Since point D is an arbitrary point, the subscript D is dropped in the following discussions.

Moreover, at the wall, $w={\partial w}/{\partial z}=u=0$, i.e. (5.8) hold both on the boundary layer edge and on the wall at any point. Since the boundary layer thickness is assumed to be much less than the radii of the surface, (5.8) should be valid throughout the boundary layer, that is,

(5.9)\begin{equation} \frac{\partial w(y)}{\partial z}=\frac{u(y) \tan \theta}{r}. \end{equation}

This equation indicates that ${\partial w(y)}/{\partial z}$ is directly dependent on $u(y)$ and the geometric factors throughout the boundary layer. And ${\partial w(y)}/{\partial z}$ has the same profile as $u(y)$, that is,

(5.10a)\begin{equation} e^{\prime}=f^{\prime} \end{equation}

and

(5.10b)\begin{equation} e=f. \end{equation}

By substituting (5.10) into (4.10), we have, for $\xi$-directional momentum,

(5.11a)\begin{align} (C f^{\prime \prime})^{\prime}+\left(1+\frac{2 \xi}{\mu_e \rho_e u_e} \frac{\tan \theta}{r}\right) f f^{\prime \prime}=\frac{2 \xi}{u_e} \frac{\partial u_e}{\partial \xi}\left(f^{\prime 2}-\frac{\rho_e}{\rho}\right)+2 \xi f^{\prime} \frac{\partial f^{\prime}}{\partial \xi}-2 \xi \frac{\partial f}{\partial \xi} f^{\prime \prime}; \end{align}

for energy,

(5.11b)\begin{align} \left(\frac{C}{P r} g^{\prime}\right)^{\prime}+\left(1+\frac{2 \xi}{\mu_e \rho_e u_e} \frac{\tan \theta}{r}\right) f g^{\prime}=2 \xi\left(f^{\prime} \frac{\partial g}{\partial \xi}-g^{\prime} \frac{\partial f}{\partial \xi}+\frac{f^{\prime} g}{h_e} \frac{\partial h_e}{\partial \xi}+\frac{\rho_e u_e}{\rho h_e} f^{\prime} \frac{\partial u_e}{\partial \xi}\right)-C \frac{u_e^2}{h_e}\left(f^{\prime \prime}\right)^2. \end{align}

These equations are the expressions of the axisymmetric boundary layer flow in the proposed coordinate system.

Consider special cases: for a flat plate case or a sharp cone, $u_e=$ const, $T_e=$ const, $p_e=$ const, $f$ and $g$ are assumed to be functions of $\eta$ only.

For a flat plate, we have $\tan \theta =0$ or $r \rightarrow \infty$, (5.11) reduce to

(5.12a)$$\begin{gather} (C f^{\prime \prime})^{\prime}+f f^{\prime \prime}=0, \end{gather}$$

(5.12b)$$\begin{gather}\left(\frac{C}{P r} g^{\prime}\right)^{\prime}+f g^{\prime}+C \frac{u_e^2}{h_e}\left(f^{\prime \prime}\right)^2=0. \end{gather}$$

The above equations for a flat plate are the same as those deduced by Anderson (2006).

For a sharp cone, we have

(5.13)\begin{equation} \xi=\int_0^x \mu_e \rho_e u_e\, {{\rm d}\kern0.7pt x}=\mu_e \rho_e u_e x \end{equation}

and

(5.14)\begin{equation} \tan \theta=\frac{r}{x}. \end{equation}

Under these conditions, (5.11) for a sharp cone reduce to

(5.15a)$$\begin{gather} (C f^{\prime \prime})^{\prime}+3 f f^{\prime \prime}=0, \end{gather}$$

(5.15b)$$\begin{gather}\left(\frac{C}{P r} g^{\prime}\right)^{\prime}+3 f g^{\prime}+C \frac{u_e^2}{h_e}\left(f^{\prime \prime}\right)^2=0. \end{gather}$$

Moreover, if the transformation (4.1b) is replaced by

(5.16)\begin{equation} \eta=\sqrt{3} \frac{u_e}{\sqrt{2 \xi}} \int_0^y \rho\, {{\rm d} y} \end{equation}

then the boundary layer equations for a cone will be totally the same as that for a flat plate, i.e. (5.12). The derivation process is basically the same as that provided in § 4 and will not be repeated in this section. Hence, it is apparent that if the parameters at the boundary layer edge and at the wall are the same for a flat plate and a sharp cone, the variations in flow properties throughout the boundary layer will exhibit the same profiles in their transformed plane. The only difference is a factor of $\sqrt {3}$ between (4.1b) and (5.16), which indicates that the boundary layer thickness of a flat plate is $\sqrt {3}$ times greater than that of a cone. Subsequently, the relationship of heat flux for a cone and a flat plate can be easily obtained as

(5.17)\begin{equation} q_{cone}=\sqrt{3} q_{{flat\ plate }}. \end{equation}

Then, let us return to further simplify (5.5) for the axisymmetric stagnation flow, where (5.10) also hold. Thus, (5.5) can be further reduced to the following: for $\xi$-directional momentum,

(5.18a)\begin{equation} (C f^{\prime \prime})^{\prime}+2 f f^{\prime \prime}=f^{\prime 2}-\frac{\rho_e}{\rho}; \end{equation}

for energy,

(5.18b)\begin{equation} \left(\frac{C}{P r} g^{\prime}\right)^{\prime}+2 f g^{\prime}=0. \end{equation}

In summary, (5.4) for a cylinder stagnation point, (5.12) for a flat plate, (5.15) for a sharp cone and (5.18) for a sphere stagnation point are ordinary differential equations (ODEs) and are self-similar. Besides, they are all deduced in the same coordinate.

6. Solutions of the hypersonic laminar boundary layer equations

From §§ 2–4, no assumptions are made beyond the scope of general boundary layer theory. The boundary layer equations (4.10) given above apply to compressible flows in the symmetry plane of three-dimensional bodies; they are also equally applicable to subsonic and supersonic flows with no distinction made for such cases. They can be applied directly to hypersonic flows if the high temperature effects are neglected.

In this section we subsequently attempt to find exact solutions to the aforementioned equations. Of course, solution of the boundary layer equations, whether for planar equations (Anderson 2006) or for (4.10), both require the properties at the boundary layer edge in advance. While numerous studies have investigated methods for determining outer boundary layer parameters (Griffith, Adams & Majors 1981; Liechty et al. 2003), these do not encompass the newly introduced variable ${\partial w_e}/{\partial z}$ at the boundary layer edge that appears in (4.10), which also necessitates prior deduction. Unfortunately, even modern CFD poses significant uncertainty in deducing this value, where both the numerator ${w_e}$ and the denominator $z$ are zero in the symmetry plane, but ${\partial w_e}/{\partial z}$ does not need to be zero. The numerical solution of ${\partial w_e}/{\partial z}$ is sensitive to factors such as grid resolution and numerical schemes, similar to heat flux calculation that remains a challenge. The dissipation that is required to control error across the bow shock and in the stagnation region can modify the boundary layer and result in poor heat transfer rate predictions (Candler et al. 2015). Unlike the situation discussed in § 5.2, where ${\partial w_e}/{\partial z}$ can be derived from the kinematic property of axisymmetric flow with no requirement of pressure distribution, there currently exist no theoretical formulas to derive this value for the symmetry plane of general three-dimensional blunt bodies.

In the context of hypersonic flows, the Newtonian flow model provides a convenient approximation for obtaining the pressure distribution. We also try to deduce ${\partial ^2 p_e}/{\partial z^2}$ based on the Newtonian theory and subsequently obtain ${\partial w_e}/{\partial z}$ by solving (3.10b). However, Newtonian theory is only applicable under hypersonic conditions, and we have not found a more general method to deduce ${\partial w_e}/{\partial z}$ under subsonic or supersonic conditions. Therefore, our interest here is in applying the results of this boundary layer theory to discussions focused on hypersonic flows.

6.1. Approximation of the crossflow velocity gradient at the boundary layer edge for hypersonic flows

Newtonian theory for prediction of the pressure distribution at a point depends only on the knowledge of its local surface inclination relative to the free stream. However, the second-order pressure gradient in the $z$ direction in (3.10b) will be affected by its surrounding flow field or the shape geometry. An equivalent cone is formed to consider the effect of surface geometry on the crossflow velocity. Point D is considered to be on the symmetry plane. Due to the symmetry property, the surface near the symmetry plane is approximated as an equivalent cone here, as shown in figure 4. Moreover, the semiangle of the cone is defined as

(6.1)\begin{equation} \theta=\arctan \frac{{\rm d} r}{{\rm d}\kern0.7pt x} \quad\text{ or }\quad \tan \theta=\frac{{\rm d} r}{{\rm d}\kern0.7pt x}, \end{equation}

where AC is the axis and the length of CD represents the principal radius of curvature of the surface at point D. Line AD is tangent to the surface at point D. Unlike the tangent-cone method commonly used to estimate the local surface pressure, where $\theta$ represents the local deflection angle of the surface, the parameter $\theta$ in this study is defined by (6.1). This definition allows for the consideration of the influence of surface geometry variation on the second-order pressure gradient or the crossflow velocity gradient. For axisymmetric flows, a sharp cone, for example, AC coincides with the axis line of the actual body. For most cases, AC will continuously vary along $x$ as the surface geometry changes, and it is specifically determined by (6.1).

Figure 4. Illustration of the tangent-cone method for the approximation of the crossflow velocity gradient.

We assume a nearby point $\textrm {D}^{\prime \prime }$ near the symmetry plane, which also comes from point D rotated by a small angle $\varphi$ around the axis of the cone, i.e. $\varphi$ represents the angle between line BD and $\textrm {BD}^{\prime \prime }$. The direction of the normal vector to the surface at point $\textrm {D}^{\prime \prime }$ is consistent with the direction of the vector $\boldsymbol {\textrm {CD}}^{\prime \prime }$. Thus, in Cartesian coordinates $(X, Y, Z)$, the unit normal vectors to the surface at points D and $\textrm {D}^{\prime \prime }$ are expressed as

(6.2a,b) \begin{equation} \boldsymbol{n}_{\rm D}=\left(\begin{array}{@{}c@{}} -\sin \theta \\ \cos \theta \\ 0 \end{array}\right), \quad \boldsymbol{n}_{{\rm D}^{\prime \prime}}=\left(\begin{array}{@{}c@{}} -\sin \theta \\ \cos \theta \cos \varphi \\ \cos \theta \sin \varphi \end{array}\right). \end{equation}

To calculate the vector dot product of the free-stream velocity and the normal vector to the surface, it is necessary to transform both to the same coordinate system. Here, we transform both to the coordinate system $(x, y, z)$. Transformation between the coordinates $(x, y, z)$ and the coordinates $(X, Y, Z)$ is expressed as

(6.3) \begin{equation} \left(\begin{array}{@{}c@{}} x \\ y \\ z \end{array}\right)=\left(\begin{array}{@{}ccc@{}} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{@{}c@{}} X \\ Y \\ Z \end{array}\right). \end{equation}

Thus, the unit normal vector to the surface at point $\textrm {D}^{\prime \prime }$ in the coordinates $(x, y, z)$ is expressed as

(6.4) \begin{align} \boldsymbol{n}_{{\rm D}^{\prime \prime}}=\left(\begin{array}{@{}ccc@{}} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array}\right)\left(\begin{array}{@{}c@{}} -\sin \theta \\ \cos \theta \cos \varphi \\ \cos \theta \sin \varphi \end{array}\right)=\left(\begin{array}{@{}c@{}} -\cos \theta \sin \theta+\sin \theta \cos \theta \cos \varphi \\ \sin ^2 \theta+\cos ^2 \theta \cos \varphi \\ \cos \theta \sin \varphi \end{array}\right). \end{align}

Consider the free-stream velocity as a vector and $\beta$ is the angle between the tangent to the surface and the free stream. Then, the unit velocity vector in the coordinates $(x, y, z)$ is expressed as

(6.5) \begin{equation} \boldsymbol{n}_{V_{\infty}}=\left(\begin{array}{@{}c@{}} \cos \beta \\ -\sin \beta \\ 0 \end{array}\right). \end{equation}

From Newtonian theory, the pressure coefficient $C_p$ at point $\textrm {D}^{\prime \prime }$ is given by

(6.6)\begin{equation} C_p=\frac{p-p_{\infty}}{\dfrac{1}{2} \rho_{\infty} V_{\infty}^2}=2(\boldsymbol{n}_{{\rm D}^{\prime \prime}} \boldsymbol{\cdot} \boldsymbol{n}_{V_{\infty}})^2, \end{equation}

where the subscript $\infty$ represents the free-stream flow parameter. Substitution of (6.4) and (6.5) into (6.6) yields

(6.7)\begin{align} C_p=2(-\cos \beta \cos \theta \sin \theta+\cos \beta \sin \theta \cos \theta \cos \varphi-\sin \beta \sin ^2 \theta-\sin \beta \cos ^2 \theta \cos \varphi)^2. \end{align}

Differentiating the above equation with respect to $\varphi$ yields

(6.8)\begin{equation} \left.\begin{array}{c} \dfrac{\partial C_p}{\partial \varphi}=4(-\cos \beta \cos \theta \sin \theta+\cos \beta \sin \theta \cos \theta \cos \varphi-\sin \beta \sin ^2 \theta-\sin \beta \cos ^2 \theta \cos \varphi) \\ (-\cos \beta \sin \theta \cos \theta \sin \varphi+\sin \beta \cos ^2 \theta \sin \varphi). \end{array}\right\} \end{equation}

As $\textrm {D}^{\prime \prime }$ infinitely approaches $\textrm {D}, \varphi \rightarrow 0$ and $\omega \rightarrow 0$. Thus, $\cos \varphi \approx 1, \cos \omega \approx 1$ and $\sin \varphi \approx \varphi$. From figure 4 we can obtain the following expression:

(6.9)\begin{equation} \varphi \approx \frac{L_{{\rm D D}^{\prime \prime}}}{L_{{\rm B D}}} = \frac{L_{{\rm A D}} \omega}{r \cos \theta}. \end{equation}

Substitution of (5.7) into (6.9) yields

(6.10)\begin{equation} \varphi \approx \frac{z}{r \cos \theta}. \end{equation}

Thus, (6.8) can be simplified as

(6.11a)\begin{equation} \frac{\partial C_p}{\partial z} r \cos \theta=4(-\sin \beta)(-\cos \beta \sin \theta \cos \theta+\sin \beta \cos ^2 \theta) \frac{z}{r \cos \theta} \end{equation}

or

(6.11b)\begin{equation} \frac{\partial C_p}{\partial z}=4 \sin ^2 \beta\left(\frac{\tan \theta}{\tan \beta}-1\right) \frac{z}{r^2}. \end{equation}

Differentiating the above equation with respect to $z$ yields

(6.12a)\begin{equation} \frac{\partial^2 C_p}{\partial z^2}=4 \sin ^2 \beta\left(\frac{\tan \theta}{\tan \beta}-1\right) \frac{1}{r^2} \end{equation}

or

(6.12b)\begin{equation} \frac{\partial^2 p}{\partial z^2}=2 \rho_{\infty} V_{\infty}^2 \sin ^2 \beta\left(\frac{\tan \theta}{\tan \beta}-1\right) \frac{1}{r^2}. \end{equation}

In the symmetry plane, i.e. $\varphi =0$, the pressure coefficient is given by

(6.13)\begin{equation} C_p|_{\varphi=0}=2 \sin ^2 \beta=\frac{p_{\varphi=0}-p_{\infty}}{\dfrac{1}{2} \rho_{\infty} V_{\infty}^2}. \end{equation}

Substitution of this equation into (6.12) yields

(6.14)\begin{equation} \left.\frac{\partial^2 p}{\partial z^2}\right|_{\varphi=0}=\frac{2(p_{\varphi=0} -p_{\infty})}{r^2}\left(\frac{\dfrac{{\rm d} r}{{\rm d}\kern0.7pt x}}{\tan \beta}-1\right). \end{equation}

Given our focus on the symmetry plane flow, the subscript $\varphi =0$ is dropped in the following discussions. Equation (6.14) indicates that the contribution to ${\partial ^2 p}/{\partial z^2}$ can be divided into two parts. The first part is the contribution of the kinetic factor, which is characterized by a constant, i.e. 1 in (6.14). This part can also be seen as the contribution of an equivalent cylinder with a radius of $r$ to ${\partial ^2 p}/{\partial z^2}$. The second part is the contribution of the geometric factor, which is characterized by ${\textrm {d} r}/{\tan \beta }$, and it can also be seen as the influence of the local geometric factor or the deviation of the geometry from an equivalent cylinder with a radius of $r$ on ${\partial ^2 p}/{\partial z^2}$. Furthermore, ${\textrm {d} r}/{\textrm {d}{\kern0.8pt}x}$ represents the bluntness effects of the blunt body on ${\partial ^2 p}/{\partial z^2}$, where $\tan \beta$ represents the equivalent axisymmetric effects or axisymmetric effects in simple terms.

Consider special cases: for a hemisphere or a cylinder, ${\textrm {d} r}/{\textrm {d}{\kern0.8pt}x}=0$ or $({{\textrm {d} r}/{\textrm {d}{\kern0.8pt}x}})/{\tan \beta }=0,$ (6.14) degenerates to a formula similar to that at the stagnation point, i.e. ${\partial ^2 p}/{\partial z^2}=-{2(p_e-p_{\infty })}/{r^2}$; for a sharp cone with zero angle of attack, $\theta =\beta$ or $({{\textrm {d} r}/{\textrm {d}{\kern0.8pt}x}})/{\tan \beta }=1$, namely, the contributions of geometric and kinetic factors offset each other out. As a result, ${\partial ^2 p}/{\partial z^2}=0$.

Associating (6.14) with (3.10b) yields

(6.15)\begin{equation} {-}2(p_e-p_{\infty})\left(\frac{\dfrac{{\rm d} r}{{\rm d}\kern0.7pt x}}{\tan \beta}-1\right) \frac{1}{r^2}=\rho_e\left[u_e \frac{\partial}{\partial x}\left({\frac{\partial w_e}{\partial z}}\right)+\left({\frac{\partial w_e}{\partial z}}\right)^2+\frac{u_e^2}{R r}\right]. \end{equation}

Similarly, the kinetic factor and the geometric factor in (6.15) will affect ${{\partial w_e}/{\partial z}}$, just as they affect ${\partial ^2 p}/{\partial z^2}$.

In this study, Newtonian theory is used to obtain the pressure distribution, which in turn facilitates the calculation of the crossflow velocity gradient ${\partial w_e}/{\partial z}$. However, if more accurate pressure distributions are obtained through alternative methods, they can be applied directly to (3.10b), resulting in a more precise crossflow velocity gradient. Notably, the range of applicability of (6.14) and (6.15) is consistent with Newtonian theory, which is applicable only under hypersonic conditions. Anderson (2006) extensively discussed the historical background, various improved forms and the scope of application of Newtonian theory in his book. In the Newtonian model of fluid flow, the particles in the free stream impact only on the frontal area of the body; they cannot curl around the body and impact on the back surface. Hence, the application of Newtonian theory appears most suitable for bodies turning convex surface parts toward hypersonic oncoming gas. Moreover, it is worth noting that, similar to various methods for solving wall pressure, the precision and range of applicability of these methods are not entirely definitive. Nevertheless, Newtonian theory is straightforward and easy to apply in hypersonic flows.

6.2. Boundary layer edge parameters of blunt bodies under hypersonic conditions

Until now, this paper has articulated the new contributions made in the realm of existing boundary layer theories. We have selected several common blunt bodies encountered in hypersonic flight for developing our solving program and validating our theory. Another reason for choosing blunt bodies is the existence of efficient approximate methods for solving the boundary layer edge parameters for such shapes. Notably, for non-blunt bodies, considerations must be given to the applicability of Newtonian theory and the selection of appropriate methods for solving the boundary layer edge parameters.

Parameters at the boundary layer edge, i.e. the values of ${u_e}$, ${\rho _e}$, ${h_e}$, ${\mu _e}$ and ${\partial w_e}/{\partial z}$, must be determined to solve (4.10). While ${\partial w_e}/{\partial z}$ is deduced from (6.15), other boundary layer edge parameters are obtained using the simplest and, therefore, the most commonly employed method (Griffith et al. 1981). From the exact normal shock wave theory, ${p_s}$ behind a normal shock is then determined as

(6.16)\begin{equation} {\frac{p_{s}}{p_{\infty}}=\left[\frac{(\gamma+1)^{2}{M_{\infty}}^{2}}{4 \gamma{M_{\infty}}^{2}-2(\gamma-1)}\right]^{\gamma/(\gamma-1)} \left[\frac{1-\gamma+2\gamma{M_{\infty}}^{2}}{\gamma+1}\right],} \end{equation}

where $\gamma$ is the ratio of specific heats and $p_s$ is the total pressure behind a normal shock wave. Subsequently, $p_e$ over the body of interest is obtained from the modified Newtonian method firstly as

(6.17) \begin{equation} {C_p=\frac{p_e-p_\infty}{\dfrac12\rho_\infty v_\infty^2}=\frac{p_s-p_\infty}{\dfrac12\rho_\infty v_\infty^2}\cos^2\theta ,} \end{equation}

where $C_p$ is the pressure coefficient. The subscript $\infty$ denotes the free-stream flow parameter. Subsequently, the local flow conditions along the boundary layer edge, i.e. $u_e$, $\rho _e$, $T_e$, are deduced by isentropically expanding from these stagnation conditions to the known local pressure $p_e$, thus treating the boundary layer edge as an isentropic surface. The isentropic relations are given as

(6.18a,b)\begin{align} {\left(\frac{\rho_{e}}{\rho_{\infty}}\right)^{\gamma}= \frac{p_{e}}{p_{\infty}}=\left(\frac{1+\dfrac{\gamma-1}{2} M_{\infty}^{2}}{1+\dfrac{\gamma-1}{2}M_{e}^{2}} \right)^{{\gamma}/({\gamma-1})}=\left(\frac{T_{e}}{T_{\infty}} \right)^{{\gamma}/({\gamma-1})},\quad u_{e}=M_{e}\sqrt{\gamma\frac{p_{e}}{\rho_{e}}}.} \end{align}

The viscosity $\mu$ and $\mu _e$ are calculated according to Sutherland's law in this study. And for air at standard conditions, ${Pr}$ equals 0.71 in this study.

To complete the system, the following two thermodynamic state equations are required:

(6.19a)\begin{equation} T=T(p,\rho),\quad h=h(p,\rho). \end{equation}

Here $T$ is the temperature. A calorically perfect gas is applied in the present study, then

(6.19b)\begin{equation} p=\rho RT,\quad h=c_pT, \end{equation}

where $R$ is the specific gas constant of the gas and $c_p$ is the specific heat capacity at constant pressure; $h_e$ is calculated from (6.19b). At this point, all the necessary outer boundary layer parameters required for solving (4.10) have been obtained.

Notably, the above method does not consider the effect of the entropy layer on the boundary layer edge parameters. This method is satisfactory for the forward partition of a blunt body. Nevertheless, alternative existing methods for taking the outer edge properties can replace the methods outlined in this section, such as extracting the boundary layer edge quantities from the computational fluid dynamic solutions (Hamilton & Greene 1994; Liechty et al. 2003). Solving the Navier–Stokes equations is more time consuming. It can yield more accurate parameters such as $u_e$, $\rho _e$, $h_e$ and $\mu _e$. However, even modern CFD poses significant uncertainty in deducing ${\partial w_e}/{\partial z}$, as already described earlier. Thus, we mainly opt for ((6.15)–(6.19)) to provide a rapid estimation of the boundary layer conditions. Anderson (2006) also discussed the entropy effects on aerodynamic heating under hypersonic flows in Chapter 6 and provided some methods to mitigate their influence.

6.3. Numerical solution of the boundary layer equations

Compared with the planar boundary layer equations (Anderson 2006), exact solutions of (4.10) do not add inherent difficulty if the parameters at the boundary layer edge are known. Commonly used methods for solving these general, non-similar planar boundary layer equations, such as local similarity, difference-differential approach and finite-difference solutions (Anderson 2006), are all applicable for solving (4.10). Furthermore, the existing literature in this field is quite expansive. Finite-difference solutions discretize the equations into finite grids and solve them numerically, requiring appropriate grid resolution and numerical schemes similar to those used for solving the Navier–Stokes equations. Instead, we choose the difference-differential approach, which can simplify (4.10) to ODEs for solution in some cases and is also an exact solution for general non-similar boundary layers. The Crank–Nicolson type differential method is similar to that used by other authors (Blottner & Ellis 1973; Wang 1975; Popinski & Baker 1976) is adopted here.

Schematic for finite-difference solution of the boundary layer is shown in figure 5. Assume that we wish to calculate the boundary layer profiles at the station denoted by $i$. The derivatives in the direction along the surface are replaced with finite-difference relations. The first-order derivatives with respect to $\xi$ in (4.10) are replaced by the two-point backward finite-difference quotients as

(6.20)\begin{equation} {\frac{\partial f}{\partial\xi}=\frac{f_i-f_{i-1}}{\Delta \xi_i} \text{ at a given } j.} \end{equation}

Figure 5. Schematic for finite-difference solution of the boundary layer.

Identical expressions are used for other parameters as those used for the derivative with respect to $\xi$, such as ${\partial f'}/{\partial \xi }$, ${\partial e'}/{\partial \xi }$ and ${\partial g}/{\partial \xi }$. At the stagnation point, i.e. $i = 1$ in figure 5, the velocity gradients at the stagnation point in the $x$ and $z$ directions are calculated from Newtonian theory, that are

(6.21a,b)\begin{equation} {\left(\frac{{\rm d}u_e}{{\rm d}\kern0.7pt x}=\frac{1}{R}\sqrt{\frac{2(p_e-p_\infty)}{p_e}}\right)_s, \quad \left(\frac{{\rm d}w_e}{{\rm d}z}= \frac{1}{r}\sqrt{\frac{2(p_e-p_\infty)}{p_e}}\right)_s,} \end{equation}

where the subscript $s$ represents the stagnation point. Afterward, we proceed to solve the ordinary differential boundary layer equations (5.3) to obtain the boundary layer profiles at $i = 1$. Choosing the stagnation point as the starting point, the difference-differential method is then implemented to move to the next downstream location. In other words, $f_{i-1}$ in (6.20) are known numbers. The only unknown number in (6.20) is $f_i$. Recall that ${\textrm {d}u_e}/{\textrm {d}\xi }$ and ${\partial w_e}/{\partial z}$ are also known numbers, obtained from the known velocity variation at the boundary layer edge. Thus, when the difference expressions such as (6.20) are substituted into (4.10), the only derivatives that appear are $\eta$. Hence, (4.10) are ODEs that can be integrated across the boundary layer at station $i$. A more detailed description of the solution method in this section can be found in the book by Anderson (2006).

After substituting the derivatives with respect to $\xi$ into (4.10) they are reformulated into a system of linear algebraic equations in terms of the single independent variable $\eta$. The successive over-relaxation iterative method is then adopted to solve these equations. To accelerate convergence, a relaxation factor is applied in the iterative progress. Once (4.10) has been solved for all grid points, the transformed streamfunction values at discretized points within the boundary layer are known. These values can be used to calculate detailed parameters such as velocity, enthalpy and temperature at each grid. Consequently, they can be employed to compute other boundary layer characteristics of interest, such as skin friction and heat transfer.

7. Examples and discussions

7.1. Comparison with experimental results

7.1.1. Stagnation point heat flux

Considerable attention is often given to the stagnation point of blunt bodies, where the local heat flux is usually the greatest. Additionally, the subsequent calculations for the blunt bodies in this paper also choose the stagnation point as the starting point, requiring the boundary layer information at this initial point. Thus, we firstly compare the stagnation point heat flux of spheres obtained by solving (4.10) with both the experimental data, and the results estimated using the analytical approach proposed by Fay & Riddell (1958). Because of its simplicity, the Fay and Riddell approach has been used extensively for re-entry flight design as well as ground test facility applications, where the expression is shown in (7.1) as below:

(7.1)\begin{equation} {q_s=0.763Pr^{{-}0.6}(\rho_e\mu_e)^{0.5}\sqrt{\frac{{\rm d}u_e}{{\rm d}\kern0.7pt x}}(h_0-h_w).} \end{equation}

Equation (7.1) is its simplified form and is valid for equilibrium and a non-reacting gas.

The heat flux values corresponding to the ‘present method’ column in table 1 are obtained by solving (5.18). Equation (5.18) are ODEs and are self-similar. The boundary layer edge parameters are obtained based on ((6.15)–(6.19)). It can be observed from table 1 that the stagnation points heat flux obtained by the present method aligns very well with both experimental values and the results from the commonly used Fay and Riddell approach.

Table 1. Heat flux at the stagnation point of spheres in air between the present theoretical results and the experimental data and the Fay–Riddell results. Experimental data of case 1 are from Irimpan & Menezes (2023), cases 2 and 3 are from Irimpan et al. (2015), case 4 are from Kaattari (1978).

7.1.2. Heat flux distribution on the sphere

The second validation case involves a sphere, and the results from the present method by solving (4.10) are also compared with the experimental data. Two sets of experimental data are considered. The experiments were conducted in the Langley 31-inch Mach tunnel (Hamilton, Millman & Greendyke 1992). This study primarily employs ((6.15)–(6.19)) to obtain the boundary layer edge parameters. As a comparison, we also obtain the boundary layer edge parameters of these cases by solving the Navier–Stokes equations, which are solved by an in-house code named the advanced research tool in science and technology CFD (ARTIST-CFD). The ARTIST-CFD is a parallel, finite volume code that solves two-dimensional, axisymmetric and three-dimensional Navier–Stokes equations based on multi-block structured meshes. Roe upwind scheme with second-order MUSCL reconstruction and minmod limiter is applied for inviscid fluxes, and viscous fluxes are discretized by second-order central difference scheme. For details about the ARTIST-CFD, we refer the reader to Wang et al. (2023); Guo, Wang & Li (2024). The boundary layer edge in the computations is defined where the local stagnation enthalpy is equal to 99.5 % of the free-stream stagnation enthalpy. However, due to the difficulty in solving ${\partial w_e}/{\partial z}$, it is still obtained through (6.15).

Figure 6(a) first presents a comparison of the boundary layer edge parameters, i.e. $T_e$ and $u_e$, derived from ((6.16)–(6.19)) with those obtained from CFD. In the front portion of the model, the differences between the two methods are minor, while the main discrepancies appear in the rear part of the model. As mentioned in § 6.2, ((6.16)–(6.19)) provide a quick and straightforward method to solve for the boundary layer edge parameters but do not account for the entropy effect. Computational fluid dynamics methods, although inevitably introducing additional complexity, can provide more reliable boundary layer edge parameters, except for ${\partial w_e}/{\partial z}$. After obtaining the boundary layer edge parameters using the two methods described above, we used them as inputs to solve the boundary layer equation (4.10), resulting in the respective boundary layer parameter profiles. These were then compared with the boundary layer parameter distributions obtained from CFD, as shown in figures 6(b) and 6(c), which respectively present the temperature and velocity profiles. When the boundary layer edge parameters are consistent, the theoretical solutions from this study match quite well with the CFD-derived boundary layer parameter distributions at both $\theta =0^{\circ }$ and $\theta =41.7^{\circ }$ positions. Here $\theta$ is the angle between the radius vector and the flight direction, as illustrated in figure 6(b). Since the velocity at the edge of the boundary layer at the stagnation line is zero, as also illustrated in figure 6(a), only the velocity profiles at $\theta =41.7^{\circ }$ are provided.

Figure 6. (a) Temperature and velocity profiles at the boundary layer edge. (b) Temperature profiles in the boundary layer. (c) Velocity profiles in the boundary layer. Two methods are used to obtain the boundary layer edge parameters; one involves deriving them through ((6.16)–(6.19)) (dash dot), while the other involves solving the Navier–Stokes equations (solid). The results are obtained under the conditions of $R = r = 50.8$ mm, $T_w = 300$ K, $\gamma = 1.4$, $M_{\infty } = 9.74$, $p_{\infty } = 65.2$ Pa, $T_{\infty } = 53.17$ K and zero angle of attack.

Furthermore, the heat flux distribution on the sphere between the present theoretical results and the experimental data are shown in figure 7, where two boundary layer edge parameters are utilized. One involves deriving them through ((6.16)–(6.19)) (dash dot), while the other involves using CFD (solid). As shown in figure 7, in the forward region of the model, particularly where $\theta$ is smaller than $70^{\circ }$, the results by the present method exhibit a close agreement with experimental data (i.e. within $\pm 10\,\%$) for the two boundary layer parameters. However, as $\theta$ further increases, disparities between theoretical and experimental results become more pronounced while using ((6.16)–(6.19)) to deduce the boundary layer edge parameters. This discrepancy primarily arises because the simplified method employed in ((6.16)–(6.19)) neglects the influence of the entropy layer. This issue becomes prominent at the rear of blunt bodies, where the boundary layer thickness increases and interacts with the entropy layer, resulting in vorticity interaction. Consequently, the entropy layer poses analytical challenges when calculating the above boundary layer parameters on the surface. This is also evident in figure 6(a), where the discrepancy between the boundary layer edge parameters derived from ((6.16)–(6.19)) and those obtained via CFD increases progressively towards the rear part of the model. Certainly, while choosing alternative methods, such as extracting the boundary layer edge quantities from CFD, these discrepancies in the rear part of the sphere are subject to improvement (solid line in figure 7).

Figure 7. Heat flux distribution on the sphere in air between the present theoretical results and the experimental data under the conditions of $R = r = 50.8$ mm, $T_w = 300$ K, $\gamma = 1.4$ and zero angle of attack. Two methods are used to obtain the boundary layer edge parameters, one involves deriving them through ((6.16)–(6.19)) (dash dot), while the other involves solving the Navier–Stokes equations (solid). Other conditions are as follows: Case 1: $M_{\infty } = 9.74$, $p_{\infty } = 65.2$ Pa, $T_{\infty } = 53.17$ K; Case 2: $M_{\infty } = 9.55$, $p_{\infty } = 29.8$ Pa, $T_{\infty } = 54$ K.

7.1.3. Heat flux distribution on the spherical blunt cone

The third comparison case involves a $15^{\circ }$ semiapex spherical blunt cone. The experimental tests were conducted in air at the Ames 3.5-foot hypersonic wind tunnel (Cleary 1969). The tests were made at a Mach number of 10.6 and a total temperature of 1111 K. The authors conducted tests at different Reynolds numbers, angles of attack and nose radius. We selected a specific set of conditions with a nose radius of $R = 9.5$ mm, a total model length of $L = 541.53$ mm and an incoming unit Reynolds number of $3.9 \times 10^6\textrm{ m}^{-1}$ for comparison. The heat flux results at three different angles of attack were compared. The boundary layer edge parameters were obtained both through ((6.16)–(6.19)) and CFD, similar to the previous sphere case. Notably, ${\partial w_e}/{\partial z}$ was still obtained through (6.15).

The comparative results are depicted in figure 8. Notably, the $X$ in figure 8 is based on values in the Cartesian coordinate system, with the coordinate origin at the stagnation point of the model. For the heat flux results (solid) obtained by solving (4.10) with boundary layer edge parameters derived from CFD, there is relatively good agreement with experimental data across all three conditions as shown in figure 8. When using boundary layer edge parameters quickly obtained from ((6.16)–(6.19)), excellent agreement is observed between the theoretical (dash dot) and experimental heat flux values for the model's forward region at $0^{\circ }$ angle of attack, as well as for the entire model at $15^{\circ }$ and $20^{\circ }$ angles of attack. However, at the rear section of the model with a $0^{\circ }$ angle of attack, disparities between theoretical (dash dot) and experimental results become more pronounced. This discrepancy is attributed to the same reason for the increasing disparity in heat flux observed at the rear section of the model in figure 7. At the rear section of the model with a $0^{\circ }$ angle of attack, the boundary layer is thicker, making it more susceptible to the entropy effects, thereby affecting the accuracy of solving ((6.16)–(6.19)) to obtain boundary layer edge parameters. Similar to the cases shown in figure 7, CFD simulations provide more accurate boundary layer edge parameters, thereby enhancing the accuracy of heat flux calculation at the rear portion of the model under $0^{\circ }$ angle of attack.

Figure 8. Heat flux distribution along the windward meridian of the $15^{\circ }$ semiapex spherical blunt cone in air, comparing the present theoretical results with the experimental data under three different angles of attack. The test conditions are as follows: $R = r = 9.5$ mm, $T_w = 294.4$ K, $\gamma = 1.4$. $M_{\infty } = 10.6$, $p_{\infty } = 132.1$ Pa, $T_{\infty } = 47.33$ K, $L = 541.53$ mm. Two methods are used to obtain the boundary layer edge parameters, one involves deriving them through ((6.16)–(6.19)) (dash dot), while the other involves solving the Navier–Stokes equations (solid).

Based on the above analysis, it can be concluded that solution accuracy of the boundary layer equations (4.10) primarily relies on the precision of the boundary layer edge parameters. With reliable boundary layer edge parameters, solving (4.10) can yield reliable results. Besides, the impact of entropy on the boundary layer edge parameters is itself an important topic, with dedicated research being conducted in this area. This paper does not delve into the discussion of how entropy affects the accuracy of solving ((6.16)–(6.19)) to obtain boundary layer edge parameters. Additionally, for most regions of the blunt body considered in this study, ((6.16)–(6.19)) can still produce satisfactory results. Thus, ((6.16)–(6.19)) are still employed in the following examples.

7.2. Leading edge of a saucer-shaped body

One example to be examined here has a very simple shape, namely, a semi-circular disk with a semi-circular edge. Similar shapes can be found in many hypersonic vehicles, such as X51 and some other wave riders. In these shapes, the most severe aerodynamic heating occurs at the leading edge, and existing methods can only predict it approximately. Thus, this paper aims at providing more precise results using the proposed method. The schematic geometry of the semi-circular disk is shown in figure 9, where $R$ and $r$ remain constant along the whole symmetry plane. In the case where $R / r \rightarrow \infty$, the body can be treated as a cylinder with a radius of $r$ and if $R / r=1$, the body is a sphere.

Figure 9. Schematic geometry of a semi-circular disk body.

Figure 10 shows the theoretical heat flux at the stagnation point and at the symmetry plane using the non-dimensional form of $q / q_{{s}, {sphere}}$, where $q_{s}$ represents the heat transfer rate at the stagnation point and the subscript ‘sphere’ denotes the case under $R / r=1$. As shown in figure 10(a), the heat flux rate at the stagnation point $(\theta =0)$ decreases rapidly as $R / r$ increases, and the final result converges to that of a cylinder where $q_{{s}, {cylinder}} / q_{{s}, {sphere}}=0.7274$. This value can be interpreted as the difference between the leading coefficient in the approximate expression for the stagnation point heat transfer rate, i.e. 0.57 for a cylinder and 0.763 for a sphere, as shown by Anderson (2006).

Figure 10. Theoretical heat flux under the conditions of $M_{\infty }=7$, $p_{\infty }=1197$ Pa, $T_{\infty }=226.5$ K, $T_w=293$ K, $\gamma =1.4$, $r=5$ mm and zero angle of attack. (a) At the stagnation point. (b) At the symmetry plane. The boundary layer edge parameters are obtained based on ((6.16)–(6.19)).

Figure 10(b) shows $q / q_{{s}, {sphere}}$ as a function of $\theta$ at the symmetry plane, including the results from Lees’ approximation (Lees 1956) under $R / r=1$. Lees’ approximation is a method used to estimate aerodynamic heating in hypersonic flows. It involves an analytical approximation developed by Lees (1956), which provides a simplified approach to calculate heat flux at the surface of a body moving at hypersonic speeds. Despite its simplicity, Lees’ approximation has been widely used in the field of hypersonic aerodynamics due to its ease of implementation and reasonable accuracy for certain flow conditions. At $\theta$ greater than $60^{\circ }$, the change in $q / q_{{s}, {sphere}}$ with $R / r$ is relatively small. At $\theta$ smaller than $60^{\circ }, q / q_{{s}, {sphere}}$ decreases with an increasing $R / r$ and converges to the heat flux distribution in the symmetry plane of a cylinder. At $R / r=10$, the maximum deviation of $q / q_{{s},{sphere}}$ from that for the cylinder heat flux distribution is within $4\,\%$. Besides, the results presented in this paper also align well with Lees’ approximation under the condition $R / r=1$, which is acceptable for aerodynamic heating calculations.

7.3. Leading edge of an ellipsoid-cone-shaped body

Blunt shapes are of particular interest for various applications, including aerodynamic decelerators for aerobraking systems or as re-entry vehicles. Among these configurations, the spheroid has garnered attention as an orbital, ballistic-entry vehicle. Thus, a spheroid, defined as an ellipsoid of revolution about the major axis of an ellipse, is examined here. The equation of a spheroid with $X$ as the symmetry axis is given as

(7.2)\begin{equation} \frac{X^2}{a^2}+\frac{Y^2+Z^2}{c^2}=1 . \end{equation}

The semi-axis $c$ is the equatorial radius of the spheroid and $a$ is the distance from centre to pole along the symmetry axis. The model is an oblate spheroid if $c > a$ and a prolate spheroid if $c < a$. The case of $a = c$ corresponds to a sphere. The schematic geometry of the spheroid is shown in figure 11. A full mathematical description of the analytical expression for the principal curvature is given by Weatherburn (1927).

Figure 11. Schematic geometry of a spheroid body.

In figure 12 the crossflow velocity gradient and heat flux at the symmetry plane are plotted as functions of $\theta$ for $c / a=0.6$, 0.8, 1, 1.3, 1.5 and 2. Both are displayed using the non-dimensional form and normalized by their value at the stagnation point. It is important to note that at the stagnation point for the model illustrated in figure 11, the flow diverges, indicated by $({\partial w_e}/{\partial z})_s$ being positive in the current coordinate system. Therefore, values of ${({\partial w_e}/{\partial z})}/{({\partial w_e}/{\partial z})_s}$ larger than zero indicate that the flows are diverging from the symmetry plane, and there are mass flows diverging from the symmetry plane for all cases, as evidenced in figure 12(a). Comparatively, the prolate spheroid exhibits a smaller crossflow velocity gradient than the sphere, whereas the oblate spheroid displays a larger crossflow velocity gradient. Generally, a larger crossflow velocity gradient corresponds to greater flow divergence from the symmetry plane, thinner boundary layer thickness or higher temperature gradients near the wall. This trend is also reflected in figure 12(b), where the heat flux values along the symmetry plane gradually increase from $c / a=0.6$ to $c / a=2$.

Figure 12. Crossflow velocity gradient (a) and heat flux (b) distributions at the symmetry plane under the conditions of $M_{\infty }=7, p_{\infty }=1197\ \mathrm {Pa}, T_{\infty }=226.5\ \mathrm {K}, T_w=293\ \mathrm {K}, \gamma =1.4, a=50\ \mathrm { mm}$ and zero angle of attack. The boundary layer edge parameters are obtained based on ((6.16)–(6.19)).

Furthermore, the heat flux distribution on the oblate spheroid exhibits greater complexity, where the heat flux on some aft section of the body is even larger than that at the stagnation point, as observed for $c / a=2$ in figure 12(b). On a general three-dimensional blunt body, the longitudinal and circumferential compression both contribute to heat fluxes. In the case of $c / a=2$, the crossflow velocity gradient is larger than that at the stagnation point. Consequently, the accompanying increase in circumferential compression on some portion of the boundary layers may dominate the compression effect, resulting in higher heat flux values. These findings also reflect that the impact of three dimensionality on heat fluxes is much more complicated than the heat fluxes analysed in two-dimensional cases where only longitudinal compression exists.

Another cone-shaped body examined in this study is an elliptic cone, which shares similarities with vehicles like the HIFiRE-5 (Juliano, Jewell & Kimmel 2019) and HyTRV (Xiang et al. 2022). The schematic geometry of the elliptic cone considered here is shown in figure 13, where $r = 20$ mm, $L = 100$ mm and the half-angle in the minor axis $(X\unicode{x2013} Z)$ plane is 7$\deg$. The ellipticity $c/a$ is set to 1, 1.5, 2 and 3. Notably, the case of $c/a = 1$ corresponds to a 7-deg sphere–cone configuration.

Figure 13. Schematic geometry of an elliptic cone body (not to scale).

The crossflow velocity gradient and heat flux at the $\phi =0^{\circ }$ symmetry plane exhibit similar trends to those observed for a spheroid body, as shown in figure 14. Hence, we will not discuss this in detail here. Instead, more attention is drawn to the parameters at the $\phi =90^{\circ }$ symmetry plane, as shown in figure 15. The values of ${({\partial w_e}/{\partial z})}/{({\partial w_e}/{\partial z})_s}$ under $c / a=1$ are positive, indicating that the flows are diverging from the symmetry plane. However, for $c / a$ that are larger than 1, the values of ${({\partial w_e}/{\partial z})}/{({\partial w_e}/{\partial z})_s}$ can become negative at certain sections in the aft body of the cone, indicating converging flow to the symmetry plane; and the larger the $c / a$ is, the further forward ${({\partial w_e}/{\partial z})}/{({\partial w_e}/{\partial z})_s}$ becomes negative.

Figure 14. Crossflow velocity gradient (a) and heat flux (b) distributions at the $\phi =0^{\circ }$ symmetry plane under the conditions of $M_{\infty }=7$, $p_{\infty }=1197$ Pa, $T_{\infty }=226.5$ K, $T_w=293$ K, $\gamma =1.4$, $r=20$ mm and zero angle of attack. The boundary layer edge parameters are obtained based on ((6.16)–(6.19)).

Figure 15. Crossflow velocity gradient (a) and heat flux (b) distributions at the $\phi =90^{\circ }$ symmetry plane under the conditions of $M_{\infty }=7$, $p_{\infty }=1197$ Pa, $T_{\infty }=226.5$ K, $T_w=293$ K, $\gamma =1.4$, $r=20$ mm and zero angle of attack. The boundary layer edge parameters are obtained based on ((6.16)–(6.19)).

In the case of converging flow, the shock will be pushed farther away, resulting in a thicker boundary layer thickness and reduced heat flux, as illustrated in figure 15(b). Nevertheless, the non-dimensional heat fluxes at the $\phi =90^{\circ }$ symmetry plane are relatively close to each other for all cases. Notably, we do not consider the presence of separated flows in this analysis.

8. Concluding remarks

The general boundary layer equations are derived specifically for the symmetry plane of three-dimensional bodies. These equations are established on the basis of a body-oriented orthogonal curvilinear coordinate system relevant to the principal curvatures. Derivation of the boundary layer equations relies on the symmetric properties of the flow and the geometry, both of which are independent of the flow state. Thus, the equations derived in this paper based on compressible laminar flows are also applicable to incompressible flows. In addition, a turbulence model such as the Baldwin–Lomax model can be incorporated into the current laminar boundary layer equations ((3.7)–(3.9)) to address turbulence effects. Moreover, ((3.7)–(3.9)) can also be extended to account for chemically reacting viscous flows by adding the species continuity equations and modifying the energy equation to incorporate diffusion effects. These extensions align with those presented by Anderson (2006), where he augmented the two-dimensional laminar boundary layer equations with turbulence models and expanded the equations to accommodate high temperature reacting flows.

In this paper the choice of a coordinate system relevant to the principal curvatures is a good strategy, which can avoid the singularity appearing in the metric coefficient and geodesic curvatures in other coordinates (Cebeci, Chen & Kaups 1990). This simplifies the solution of the equations to some extent. Besides, the calculation of metric coefficients is required only once in the proposed coordinate, which is also independent of the flow characteristics.

With no restrictions on the magnitude of the crossflow, the crossflow velocity gradient at the boundary layer edge on the symmetry plane is derived based on Newtonian theory, where the three-dimensional effects are considered. This theoretical expression has not been previously seen in other studies. Other boundary layer edge parameters are solved based on the normal shock and isentropic relations, which is the simplest and, therefore, the most commonly employed method for blunt bodies. Notably, the entropy-layer effects on the boundary layer edge conditions are not considered in this method. Subsequently, the difference-differential method is applied to obtain an exact solution of our derived boundary layer equations. Comparisons between the theoretical and experimental heat flux data regarding the stagnation point heat flux and heat flux distribution on both a sphere and a sphere-cone body are displayed. The flow characteristics at the leading edge of saucer-shaped bodies and ellipsoid-cone-shaped bodies are finally presented and discussed based on the developed boundary layer theory. Besides, theoretical expression of the crossflow velocity gradient has the same range of applicability with Newtonian theory, i.e. applicable only under hypersonic conditions. However, other existing methods for taking the outer edge properties are all applicable to the present boundary layer equations.

Given the known boundary layer edge parameters, we can degrade (4.10) to an ODE equation using the difference-differential method to solve it, thus minimizing the influence of grid resolution and numerical schemes. Hence, its solution accuracy heavily relies on the precision of the boundary layer edge parameters. One reliable approach to solving these parameters is through fully viscous three-dimensional calculations (Anderson 2006). Notably, the difference-differential approach chosen in this paper for solving the boundary layer equations requires a starting point of known information for computation, i.e. the information of station 1 in figure 5. Initially, the three-dimensional stagnation point boundary layer equations are solved, and the obtained results are used as the starting point for downstream calculations. However, for non-blunt bodies, the selection of the starting point needs additional consideration if employing the difference-differential approach proposed in this paper.

Notably, theoretical framework developed herein for hypersonic flows shares similarities with the discussions on viscous boundary layers outlined in Anderson (2006, Chapter 6). Nonetheless, the present theory can provide a much more comprehensive range of applications for the symmetry plane of three-dimensional bodies. This theory completely covers the planar boundary layer equations and the two typical cases examined in (Anderson 2006, Chapter 6).