The Lees–Dorodnitsyn (L–D) boundary layer equations for two-dimensional, non-reactive, laminar, hypersonic, boundary layer flows, and an assumption of an isentropic external flow are examined. They are applied to various geometries for which the Thin Shear Layer assumptions are valid. This study expands on previous work to develop a novel and robust methodology for computing high-temperature hypersonic flows using a uniform and compact computational stencil implemented through a computational tool, the Bulk-property Boundary Layer (BuBL) solver. In particular, we explore the impact of treating high-temperature effects present in hypersonic flows, namely, treating air as a thermally perfect gas with temperature-variable properties. The ability to solve these flows computationally using second-order finite difference methods is evaluated as are various models for viscosity, Prandtl number, and specific heat. The methodology for solving the external flow properties in the transformed L–D computational domain is also discussed. It is shown that the L–D equations evaluated using the “box” computational stencil are an effective means for evaluating laminar hypersonic boundary layer flows. Solutions for displacement and momentum thicknesses, skin friction, and Stanton number variations are obtained as a function of Prandtl number, specific heat model, and Mach number. Verification and validation measures are performed for the code. Excellent agreement is found in comparisons between BuBL and other computational fluid dynamics and experimental results, thus demonstrating the utility of the proposed methodology.

## I. INTRODUCTION

Modeling the aerothermodynamic behavior of high-speed flows is important for the development of hypersonic vehicles. Solvers that can quickly produce accurate solutions will enable effective design processes for hypersonic aircraft. To ensure utility, these solvers should be able to replicate wind tunnel flows with various working fluids and operating conditions to capture assorted flight conditions. They should also be able to model flight conditions that cannot be effectively recreated in wind tunnels. In this work, we are concerned with laminar, hypersonic boundary layers,^{1} as it is common to achieve extensive regions of laminar flow, even at extremely high speeds.^{2} To that end, we aim to model laminar, hypersonic, boundary layers without incurring unnecessary computational expense.

Since laminar solutions have been well studied in past decades, literature on their behavior is useful for grounding novel models for hypersonic boundary layers. These solutions maintain their relevance in modern hypersonics research^{3} for a number of reasons, perhaps most notably due to their role in understanding flow transition to turbulence. Laminar boundary layer solvers are frequently used to predict the mean flow profiles required for Linear Stability Theory (LST) and Parabolized Stability Equation solvers (PSE)^{4–8} as well as other transition solvers.^{9–11} In the same vein, the profiles generated by boundary layer solvers can be used to provide mean inflow boundary conditions for Direct Numerical Solvers (DNS).^{7,12} Accurately predicting the mean boundary layer profiles resulting from various flow parameters is also crucial in modeling output variables, such as wall heating, surface drag, and separation behavior, in addition to understanding a flow's turbulent transition regime. Zhao *et al.*^{13} underscored this importance while discussing the role of total temperature, Mach number, and Reynolds number as driving parameters for transitional boundary layer behavior—laminar boundary layers that progress to become turbulent. Kimmel^{14} also describes the role of the pressure gradient in the ultimate transition of the laminar boundary layer. Oz and Kara^{15} and Radenac *et al.*^{16} highlighted the need for efficient boundary layer solvers. As such, we introduce a methodology that would allow one to capture the impact of various flow parameters on the resulting laminar boundary layer.

Recent literature often takes fairly similar approaches to modeling laminar solutions, but their differences have large impacts on the applicability and accuracy of their results. Here, we will mostly focus on the combination of the Euler equation and differential boundary layer solvers, which is the methodology used in this paper. This methodology is preferable to integral methods which, despite their speed and simplicity, are generally less accurate.^{17} Differential Euler/boundary layer solvers use boundary layer equations derived from the Thin Shear Layer (TSL) approximations^{18} and the Euler equations for the external flow. Precedence for using the equations derived from the TSL assumptions is found throughout literature, especially for capturing flows with pressure gradients.^{19–21} Many times the effect of pressure gradients is neglected for simplicity.^{7,8,11,22} The Euler/boundary layer equation solution methodology is typically expected to be much less computationally expensive and complex than a full Navier–Stokes solver and to yield sufficiently accurate results.^{5,23}

A weakness of many existing methodologies is their failure to treat working fluid properties as temperature dependent. Because hypersonic flows are known to have large temperature gradients, the approximation of constant properties such as for viscosity, specific heat, and thermal conductivity fails to hold.^{24,25} Many past boundary layer solver methodologies neglect some of the effects of temperature variability^{5,26,27} or fail to model them sufficiently for our context by using curve fits or approximations to estimate flow properties.^{2,25,28–30} Solutions that do account for these properties tend to be implemented through more complex and computationally expensive means.^{4,31,32} Our methodology offers an alternative.

We offer a means to provide insight into the behavior of laminar hypersonic boundary layers by using the Lees–Dorodnitsyn (L–D) equations in coordination with models for the temperature-dependent quantities of air. In particular, we focus on one particular viscosity model and two models for specific heat. We also allow for flexibility in Prandtl number selection. These models are implemented through a new solver, the Bulk-property Boundary Layer (BuBL) solver, a MATLAB toolkit. BuBL is a robust, accurate, and efficient finite-difference solver capable of solving a variety of hypersonic boundary layer flows.

We organize this paper as follows: First, we introduce the L–D equations, and then, we detail the perfect gas model and thermally perfect gas extension used to define the working fluid's properties. Details and background of the computational approach are then given, followed by a discussion of results based on multiple comparison studies that establish BuBL as a useful tool to model laminar hypersonic boundary layer flow.

### A. Governing equations

The boundary layer equations used are referred by Anderson as the Lees–Dorodnitsyn (L–D) equations.^{24} These equations are also attributed to Levy *et al.*^{24,33} The L–D equations result from the Thin Shear Layer (TSL) approximations^{18} applied to the Navier–Stokes equations with an additional coordinate transformation. A detailed derivation of the L–D equations is provided by Anderson.^{24} Similar transformations are used in recent works.^{7,12,15,23,34,35} The L–D equations are valid for laminar, compressible boundary layers. In accordance with the TSL approximations, the flows included in this work are those whose boundary layers remain laminar, attached, and thin. The requirement of thinness means that while we can model flows over geometries with streamwise curvature by using their induced pressure gradients, we can expect a degradation of accuracy if the gradients become sufficiently adverse. In this paper, the highly nonlinear L–D equations are examined for hypersonic conditions. However, our implementation can also be used to solve other, slower speed, compressible flows.

The L–D transformations and the corresponding equations are as follows:

where $g=h/he,\u2009f\u2032=u/ue$, and the primes ($\u2032$) denote derivatives with respect to *η*. “*e*” subscripts describe flow quantities at the edge of the boundary layer, which is equivalent to the external flow quantity at the corresponding location, and

is the local value of the Chapman–Rubesin factor.^{24}

It should be noted that for the Lees–Dorodnitsyn transformation of the compressible boundary layer equations, the transformation from the *x* – *y* space to the $\xi \u2212\eta $ space requires integration. This integration is useful in a hypersonic context because it can accommodate the significant gradients in flow quantities that are characteristic of hypersonic flows. The transformation achieves this by scaling both the streamwise and normal coordinates by several external or edge (referring to the edge of the boundary layer) flow parameters: density, viscosity, and velocity. Boundary layer evolution information is also captured by the *ξ* parameter, which is possible due to its integral formulation.^{35} However, while this transformation has its advantages, this integral transformation creates additional computational complexity compared to other boundary layer methods. The L–D equations were chosen, despite being highly nonlinear and fifth-order, because they largely remove the need to increase the resolution of the boundary layer flow near the wall and also account for boundary layer growth.

For these equations, dimensionless profiles of the stream function $(\phi )$, velocity, shear stress, enthalpy, and heat transfer at various boundary layer locations were generated and are given by *f*, $f\u2032,\u2009f\u2033$, *g*, and $g\u2032$, respectively. With these variables, the system could be rewritten as a system of five first-order equations. This was done by reposing the momentum and energy equations in terms of the following five variables, which correspond to the five aforementioned dependent variables:

Here, we also note that for a calorically perfect gas, an ideal gas with a constant ratio of specific heats is given as follows:

Using these variables, the L–D equations can be equivalently rewritten as follows:

Note that Eq. (3) becomes a boundary condition for subsequent calculations.

This form is better suited for the numerical evaluation than the original. $F,U,S,H$, and *Q* are unknown quantities to be solved for at each location within each boundary layer profile. *C* and $C\u2032$, the Chapman–Rubesin factor and its derivative, in these equations are functions of *H*. As such, they are evaluated analytically as functions of the discretized variable *H* for the case of a calorically perfect gas. Later, we will discuss how for the thermally perfect gas case, they are evaluated as functions of temperature, *T*.

### B. Boundary conditions

The boundary layer equations are subject to a set of boundary conditions. Quantities at the leading edge are equivalent to the freestream values. There is a no-slip boundary condition at the wall. At the edge of the boundary layer, the pressure, velocity, and enthalpy are equivalent to their external/edge flow values. The following boundary conditions are used in all cases: *F _{wall}* = 0,

*U*= 0,

_{wall}*U*= 0, and

_{edge}*H*= 0. In the case of an isothermal (cold or hot) wall condition, the following fifth boundary condition was used:

_{edge}*T*=

_{wall}*T*. For this, the user would specify the wall temperature as a constant or as some function of

_{specified}*x*. For the case of an adiabatic wall, the fifth boundary condition is

*Q*= 0.

_{wall}## II. MODEL DEVELOPMENT

A key to this particular work and the capabilities of BuBL are the models used for air's flow properties. These models include Sutherland's law for viscosity, the Prandtl number selection, the perfect gas model, and Euler's equation for the inviscid edge of the boundary layer.

### A. Viscosity

The impact of variable viscosity tends to decrease temperatures and increase the velocity within the boundary layer.^{29} These trends are extensively highlighted for Falkner–Skan wedge flows^{36} in Pantokratoras.^{37} As such, noting the impact of viscosity-variability, we chose to use Sutherland's law for viscosity in this work. Sutherland's law is accurate in the range of several thousand degrees and is appropriate for many hypersonic applications.^{24} Sutherland's law is as follows:

where $Tref=288K,\u2009S=110K$, and $\mu ref=1.789\xd710\u22125kgm\u2009s$.

### B. Perfect gas model

In BuBL's implementation, air is represented as an ideal gas. Accordingly, the ideal gas law holds. In this study, we will begin by assuming air is not only an ideal gas but also calorically perfect. For calorically perfect gases, specific heat is constant and does not vary with temperature. It then follows that the equation relating enthalpy to temperature

simplifies to

In this case, the specific heat is

where *R* is the specific gas constant.

The assumption of a constant specific heat is only reasonable for flows of moderate temperatures.^{24} The ideal gas assumption is important because it is used as a state equation to relate density, temperature, pressure, and ultimately, enthalpy, velocity, and Mach number. Later, we discuss relaxing the calorically perfect gas assumption and treat air as only thermally perfect.

### C. Prandtl number

The Prandtl number is significant because it describes the ratio between viscous dissipation and heat conduction.^{38} For the purposes of this study, we model the Prandtl number as constant. This is a reasonable assumption because the Prandtl number generally varies a maximum of 20%–30% across a hypersonic boundary layer.^{24} Sufficient agreement with experimental data has been found when Prandtl number is treated as constant, as long as viscosity is allowed to vary.^{39} The expression for Prandtl number is given in (16). BuBL allows for an arbitrary Prandtl number to be specified. Values for air are usually limited to $0.7\u20130.75$ or approximated as 1.0

where *μ* is the viscosity in $(kgm\u2009s)$, *c _{p}* is the specific heat at constant pressure in $(Jkg\u2009K)$, and

*k*is the thermal conductivity of the gas in $(Wm\u2009K)$.

### D. External flow

The external or edge flow is not extensively considered in this work. We compute the external flow only as a means to provide inputs for the boundary layer flow calculations. It is assumed that the boundary layer does not substantially interact with the external flow—the boundary layer and the external flow have a one-way coupling. The properties of the external flow are allowed to vary in the streamwise direction. The external flow is also taken to be isentropic and inviscid. However, the external flow is treated as having a finite viscosity, calculated according to Sutherland's law to be consistent with the boundary layer.

BuBL calculates boundary layers for a given pressure gradient specified in real coordinates (the *x*-coordinate). This pressure gradient can be used to establish basic correlations or to correspond to specific external flow conditions or geometries using non-interacting inviscid theory. This external pressure gradient can be obtained from experimental results, Direct Numerical Simulation (DNS) results, local surface inclination methods,^{24} etc. From the pressure gradient, the quantities required at each computational grid-point (i.e., temperature, Mach number/velocity, corresponding location in the transformed *ξ*-coordinates) for the boundary layer calculations can be solved using isentropic relations, the *x*-*ξ* coordinate transformation equation, the Euler equation, and Sutherland's law. The Euler equation is provided in (17), and the complete calculations are provided by Onyeador,^{40}

This gives

in the L–D coordinate system.

The process of finding solutions for the external flow is rather straightforward using a basic Ordinary Differential Equation (ODE) solver when the gas is modeled to be calorically perfect. This methodology is reexamined later in Sec. III.

## III. THERMALLY PERFECT GAS EXTENSION

Relaxing the calorically perfect gas assumption can increase the accuracy of hypersonic calculations, particularly for high-temperature flows.^{4} Under the thermally perfect gas assumption, the ideal gas law still holds. However, the standard isentropic relations do not. This, in turn, complicates the solution methodology for the external flow. For a calorically perfect gas, we can no longer define specific heat to be constant. Instead, for these calculations, air is treated as a mixture of two molecules: diatomic oxygen (O_{2}) and diatomic nitrogen (N_{2}). The harmonic oscillator model described by Vincenti and Kruger^{41} is used to capture the activation of the vibrational energy modes that lead to variations in specific heat with increasing temperature. From this model, we arrive at the following:

The values for (19) are found in Table I. The resulting relationship for specific heat as a function of temperature is given in Fig. 1.

. | Mole fraction, m . | Specific gas constant, R ($kJkg\u2009K$) . | Θ_{vib} (K)
. |
---|---|---|---|

Nitrogen (N_{2}) | 0.789 | 0.2968 | 3390 |

Oxygen (O_{2}) | 0.211 | 0.2598 | 2270 |

. | Mole fraction, m . | Specific gas constant, R ($kJkg\u2009K$) . | Θ_{vib} (K)
. |
---|---|---|---|

Nitrogen (N_{2}) | 0.789 | 0.2968 | 3390 |

Oxygen (O_{2}) | 0.211 | 0.2598 | 2270 |

### A. Density ratio

Additionally, in the L–D equations (7)–(11), the following density ratio appears:

Previously, for a calorically perfect gas, this could be treated by equating the density ratio with *H*. However, because specific heat is not constant in the calorically perfect case, the density ratio is rewritten as a temperature ratio according to the ideal gas law. However, it can be simplified no further.

### B. Temperature (*T*) as an Unknown

In order to account for the more complex relationship between enthalpy and temperature due to the thermally perfect gas assumption, the computational method was altered to include temperature, *T*, as an unknown at each computational grid point. This is done because temperature cannot be found as an analytical function of enthalpy. As such, numerical methods must be used. The unknowns in the system then become *F*, *U*, *S*, *H*, *Q*, and *T*—introducing a sixth set of Eqs. (20) into the following system:

Furthermore, we must introduce a sixth boundary condition that provides the temperature at the wall. For the isothermal condition, this is straightforward and remains unchanged from before,

However, if the boundary condition is adiabatic, the *T _{w}* boundary condition is set by numerically calculating its value from that of

*H*. When calculating the residuals and the Jacobian required for the iterative solver,

_{w}*C*,

*Pr*, and $\rho e\rho $ are now treated as only functions of

*T*and not

*H*(as they were previously).

### C. External flow considerations

Under the thermally perfect gas model, the isentropic relations must be revisited, as the standard ratio relations no longer apply. Instead, the isentropic relations are rederived from the definition of entropy, *s*,

Because the flow is isentropic, *ds* integrates to zero. Therefore, (21) can be rewritten in differential form as

The relevant external flow quantities needed for the subsequent boundary layer calculations can then be computed as a function of $\xi (x)$ using the following three equations, in which all of the variables with subscript “*e*” are functions of *x*:

The adiabatic conservation of energy expression (26), the ideal gas law, $P=\rho RT$, and Sutherland's law for viscosity (12) are used to relate the quantities in (25) to *T _{e}*

With the character of the fluid and its parameters defined, the flowfield can then be simulated using computational methods.

## IV. COMPUTATIONAL METHOD

To solve the L–D equations, a computational approach is established. Previous work has delved into the myriad of methodologies for solving boundary layer equations. These involve the choice of numerical scheme, coupling or uncoupling of governing equations, and choice and implementation of linearization schemes, as well as treatment of non-similarity solutions.^{21,42} Blottner provides a review of many of these boundary layer codes and solutions.^{43} The methodology implemented through BuBL serves as an update to many of these dated approaches. This is approached with the goal of avoiding the pitfalls of large computational stencils and variable grid spacings.

In practice, as is generally the case with boundary layer equations, the L–D equations behave as a parabolic partial differential equation (PDE).^{30,44} Thus, a space-marching^{18} technique can be employed for their solution. For each *ξ*-location on the computational grid, the boundary layer is computed from the wall to the boundary layer edge, *η _{e}*. In theory, the edge of the boundary layer is located at infinity due to the asymptotic nature of the boundary layer. However, computationally, with the L–D boundary layer scaling, it was sufficient to use

*η*= 7.

_{e}To solve this fifth-order system, two discretization schemes were investigated: the “box” scheme as depicted by Keller and Cebeci,^{45} and a three-point backward difference adapted from Drela.^{18} Both computational schemes had potential disadvantages. The box scheme has second order truncation error in both *ξ* and *η*, as it is a central difference scheme in both dimensions. The box scheme is compact, which is advantageous numerically. However, because it is a centered difference implicit scheme, there is a possibility for oscillatory solution behavior similar to the Crank–Nicolson scheme.^{46} The three-point backward difference also theoretically yields a second-order truncation error in both dimensions. A key difference is that the three-point scheme achieves greater accuracy by approximating that first derivative in *ξ* using the two previous grid points, making the scheme explicit in the *ξ* direction. However, for the *i* = 2 station, a two-point backwards difference for the *ξ*-derivatives was employed (see Ref. 18). This two-point scheme has the potential to introduce error into the system as it has a first-order truncation error. We also must consider that “explicit difference procedures are not recommended” because of the stringent stability requirements when an adverse pressure gradient is present.^{47} The box scheme and the three-point difference schemes are shown in Fig. 2.

Ultimately, it was found that both schemes yielded similar performances for the many solution cases tested. Since neither model clearly outperformed the other, the box scheme was selected due to its compact nature. A similar computational scheme is employed by Su and Zhou.^{48}

### A. Newton–Raphson solver

Our methodology involves finding the solution of the finite difference method using a Newton–Raphson solver as described by Drela.^{18} For the Jacobian of the Newton solver, the analytical derivatives of the discretized equations are employed. This results in a banded matrix system whose sparse structure can be exploited to quickly and reliably solve the matrix system using a native, MATLAB linear equation solver. At each solution *ξ*-location, the initial guess for the Newton solver was supplied by the previous boundary layer profile. A solution at a location was determined to be converged when the solution residuals fell below a user-specified tolerance. Values of $1\xd710\u221211$ or less are used in this paper.

### B. External flow method

The external flow calculations are responsible for transforming input values into the L–D space and for resolving the edge boundary condition. BuBL is written such that the external boundary layer can be computed using a finite-difference method with arbitrary resolution. A centered difference scheme is used in which (22)–(26) are solved simultaneously as a two-point boundary value problem. This was done so that the external flow could be taken as “exact,” minimizing the additional error introduced by the numerical calculation of the external flow. This procedure was selected such that it would not significantly contribute to the computational expense of the solution. A schematic of the resulting computational grid is shown in Fig. 3.

Importantly, increasing the resolution of the edge of the boundary layer for the preprocessing calculations did not significantly increase the runtime of the resulting code. In one case, even with an edge resolution that was 200 times that of the rest of the boundary layer stencil, the runtime increased by less than 5%. This is an insignificant increase given that runtimes are on the order of seconds for even well-resolved solutions.

## V. RESULTS AND DISCUSSION

With the implementation and methodology established, we can analyze the accuracy of the resulting computations and the impact of the various fluid properties used.

### A. Quantities of interest

In this study, there are several quantities of interest that can be used to describe the size and behavior of the boundary layer. The first two quantities to be discussed are the displacement thickness, $\delta *$, and the momentum thickness,^{33} *θ*. Expressions for these integral quantities are given in (27) and (28). The displacement thickness and the momentum thickness describe the amount of freestream inertia and momentum, respectively, displaced by the boundary layer,

where

The local skin friction coefficient, which is a non-dimensional quantity for surface shear stress, is also of interest. It is given below in L–D quantities,

Finally, the Stanton Number *C _{H}* is a dimensionless number characterizing the ratio of the heat transfer at the wall to the thermal capacity of the fluid. The complete expression to calculate

*C*using L–D quantities is given in (31). To find the most accurate

_{H}*C*for the isothermal wall case, the full L–D boundary layer equations for the corresponding adiabatic wall case must be solved to find

_{H}*h*, the adiabatic wall temperature. The Stanton number can also be reasonably approximated as (32) which is useful for the adiabatic case

_{aw}^{24}

where *s*, the Reynolds analogy factor, equals $Pr23$ for a laminar flow.^{24}

### B. Verification and general results

We aim to verify that the BuBL implementation of this methodology is sound. Verification of laminar, hypersonic flowfields is defined by Roy *et al.*^{39} as the determination of whether the model implementation accurately represents the conceptual model. Verification can be achieved by examining the results of a grid convergence study. We also choose to examine the predicted values for a number of sample cases. The results of these cases not only provide physical insight but also assist with code verification since the results are found to be consistent with the mathematical framework.

#### 1. Grid convergence study

A study of the convergence behavior of the computational stencil was completed using the derived quantities of displacement thickness and skin friction to verify the expected performance. The convergence study compared the root mean squared error of the quantities of interest. The error in boundary layer displacement thickness and wall skin friction compared to a highly resolved solution decreased as the computational grid was refined, as shown in Fig. 4. The comparison is given as a function of the grid lengthscale, which was defined as the inverse of the square root of the number of degrees of freedom (DoFs). This convergence study shows that the stencil achieved the expected second-order rate of convergence that is predicted by the truncation error of the finite difference scheme.^{43}

#### 2. Prandtl number

To establish the importance of the Prandtl number for laminar, hypersonic flows, we can examine the impact of varying the Prandtl number. Table II gives the resulting ratio of the thicknesses of the velocity and thermal boundary layers for various *Pr* values 1 m along an adiabatic flat plate with $T\u221e=500K$ and $P\u221e=1000Pa$. Most of the *Pr* values included are commonly used in literature, and *Pr* = 1.15 is included for demonstrative purposes. *δ _{u}* is the thickness of the velocity boundary layer, and

*δ*is the thickness of the thermal boundary layer. These thicknesses are taken as the value of

_{h}*y*when the respective quantity within the boundary layer reaches 99.9% of its asymptotic edge value. $||\delta u/\delta h||$ is the ratio of the two values normalized to be consistent with this ratio being equal to one when the Prandtl number equals 1. At Mach 5, we observe the expected trends based on the definition of Prandtl number: the $\delta u/\delta h$ ratio increases as

*Pr*increases. Rerunning the computations for $Ma\u221e=15$, we see similar trends for the constant Prandtl number cases. However, in general, the thicknesses of the boundary layers are larger than the $Ma\u221e=5$, as would be anticipated given the greater amounts of heating in higher Mach number flows, due to kinetic energy dissipation, and resulting lower densities within the boundary layer.

^{24}For both Mach number flows, we also evaluate the skin friction and the Stanton number. Again, as expected wall heat transfer and wall shear stress decrease with Prandtl number.

^{37}This variation is more pronounced for the lower Mach number case which is consistent with the canonical studies of Van Driest.

^{49}For the following studies, the Prandtl number is set to a common value of 0.72.

. | $Ma\u221e=5$ . | $Ma\u221e=15$ . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Pr
. | 0.71 . | 0.72 . | 0.75 . | 1 . | 1.15 . | 0.71 . | 0.72 . | 0.75 . | 1 . | 1.15 . |

$\delta u\xd710\u22122\u2009(m)$ | 1.71 | 1.72 | 1.72 | 1.77 | 1.79 | 3.73 | 3.73 | 3.75 | 3.88 | 3.97 |

$\delta h\xd710\u22122\u2009(m)$ | 2.00 | 2.00 | 1.98 | 1.94 | 1.91 | 4.00 | 4.01 | 4.01 | 4.07 | 4.10 |

$||\delta u/\delta h||$ | 0.939 | 0.939 | 0.954 | 1 | 1.031 | 0.975 | 0.975 | 0.980 | 1 | 1.014 |

c _{f} | 6.715 | 6.710 | 6.696 | 6.601 | 6.558 | 2.509 | 2.507 | 2.500 | 2.454 | 2.435 |

C _{H} | 4.218 | 4.176 | 4.056 | 3.301 | 2.987 | 1.577 | 1.560 | 1.514 | 1.227 | 1.109 |

. | $Ma\u221e=5$ . | $Ma\u221e=15$ . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Pr
. | 0.71 . | 0.72 . | 0.75 . | 1 . | 1.15 . | 0.71 . | 0.72 . | 0.75 . | 1 . | 1.15 . |

$\delta u\xd710\u22122\u2009(m)$ | 1.71 | 1.72 | 1.72 | 1.77 | 1.79 | 3.73 | 3.73 | 3.75 | 3.88 | 3.97 |

$\delta h\xd710\u22122\u2009(m)$ | 2.00 | 2.00 | 1.98 | 1.94 | 1.91 | 4.00 | 4.01 | 4.01 | 4.07 | 4.10 |

$||\delta u/\delta h||$ | 0.939 | 0.939 | 0.954 | 1 | 1.031 | 0.975 | 0.975 | 0.980 | 1 | 1.014 |

c _{f} | 6.715 | 6.710 | 6.696 | 6.601 | 6.558 | 2.509 | 2.507 | 2.500 | 2.454 | 2.435 |

C _{H} | 4.218 | 4.176 | 4.056 | 3.301 | 2.987 | 1.577 | 1.560 | 1.514 | 1.227 | 1.109 |

#### 3. Ideal gas model

We can observe the impact of relaxing the assumption that the working fluid, air in this case, is a calorically perfect gas. As an example, a case for a $Ma\u221e=10$ flow with an adiabatic wall condition and adverse pressure gradient over a 1 m geometry was simulated. Distributions of the resulting Stanton numbers and displacement thicknesses for the thermally and calorically perfect cases are given in Fig. 5. We note, here, we expect from Hayes and Probstein^{50} that *C _{H}* at

*x*= 0 is a singularity point since $CH\u221dqw$, and

*q*scales approximately with $x\u22121/2$ for sufficiently thin flows.

_{w}^{24}

When the specific heat is a function of temperature, the heat transfer in the boundary layer increases compared to the fixed specific heat case. The increase is more pronounced as the boundary layer grows. This is likely a result of the moderated temperatures within the boundary layer that result from the activation of the vibrational modes in the thermally perfect gas case.^{41} For the sample case, neglecting the leading edge singularity, there is, on average, a 16.5% difference between the Stanton number distributions (Fig. 5).

Additionally, we note that the use of the thermally perfect gas assumption results in a significant difference in the displacement thickness of the boundary layer. Temperature drives the growth of the boundary layer,^{24} so with the reduced temperatures that stem from variable specific heats, the boundary layer would not grow to quite the same size. For the example case, after 1 m, there is a 17.4% difference in the two solutions for the boundary layer displacement thickness (Fig. 5).

Because of the nontrivial impact of the thermally perfect gas assumption, for all remaining calculations, we use this assumption over the calorically perfect assumption unless otherwise stated.

### C. Comparison studies

To showcase the utility of BuBL and the use of the L–D equations for solving laminar, hypersonic boundary layers, we performed extensive verification and validation. A selection of those results can be found in Ref. 40. Here, we choose to show further, selected validation studies that utilize computational and experimental results found by other researchers. These validation exercises are informed by the framework described by Roy *et al.*^{39}

#### 1. Gallo et al.

Gallo *et al.* performed experiments to measure boundary layer phenomena for hypersonic flows with adverse pressure gradients.^{21} They achieved this adverse pressure gradient using a curved, 2D compression surface in a hypersonic wind tunnel running nominally at Mach 10.5. They also include results from their numerical solution of the L–D boundary layer equations. This study is examined here due to its inclusion of boundary layer profiles. Boundary layer quantity profiles are the primary output of BuBL and other calculations are derived from these profile quantities. Therefore, comparing against these profiles is critical for solution validation.

The Gallo procedure differs from our simulations in a few ways. The original intention of the Gallo study was to underscore the importance of capturing the nonsimilar properties of boundary layers, i.e., accounting for boundary layer history in the solution of the boundary layer profiles. They achieve their goal by solving the boundary layer equations using the total enthalpy form of the L–D equations,^{50} a different, but largely equivalent form of the L–D equations. Additionally, their experimental geometry has a small, rounded nose region that results in a bow shock. BuBL's current implementation assumes an infinitely sharp leading edge, and as a shock-fitting solver, cannot capture the impact of this shock. Gallo *et al.* also use an alternative methodology to solve for leading edge effects. However, as we will show that these leading edge effects can be reasonably neglected. We will also show that differences in flow modeling do not significantly impact the resulting solutions. As such, we expect substantial but not complete agreement with the Gallo *et al.* results.

To model the Gallo experiments, we dimensionalize their calculated external pressure gradient and find a corresponding curve-fit to use as a boundary condition for the BuBL solver, $Pe(x)=Pspecified$. The specific pressure gradient is omitted here for brevity but can be found in Gallo *et al.*^{21} We also calculate the freestream flow conditions felt by the leading edge of the geometry and use an isothermal boundary condition as provided by Gallo *et al.*

The results for the predicted velocity profiles at three different streamwise locations are given in Fig. 6. There is excellent agreement between the results predicted by BuBL and the computational fluid dynamics (CFD) from the study. Our solution overlaps with theirs in several areas. Since the Gallo *et al.* data are presented without quantified uncertainty, we do not compare directly to these points. Rather, we consider the difference in the computed boundary layer velocities calculated at the y location of the experimental data points as done by Gnoffo *et al.*^{51} It is found that the average differences between the profiles were 2.8%, 2.7%, and 0.8% for $X/L=0.23,0.55$, and 0.81, respectively. Worth noting is that, due to the adverse pressure gradient, this flow approaches separation. However, for the length of the given geometry and flow conditions, separation is not reached.

Comparisons are conducted for the predicted displacement and momentum thicknesses over the geometry (Fig. 7). Again, we see acceptable agreement with both the Gallo *et al.* CFD and experiments. Due to both of these boundary layer quantities being derived quantities, it is perhaps not surprising that there is a greater discrepancy in these quantities than what was seen in the velocity profiles. Although the shapes of the thickness distributions are generally similar, there is some deviation in the solutions, particularly in the displacement thickness. Because our treatment of the pressure boundary condition and the leading edge presumably differed from Gallo, it is therefore reasonable that the thickness distributions from the two solvers would also have some dissimilarities. Gallo *et al.* also used different flow models and had a less stringent threshold for the residuals of their iterative solver.

Despite those factors and all other experimental uncertainty, we also find that the BuBL results are able to predict experimental results for the two boundary layer thicknesses with similar accuracy to the Gallo CFD. Over the domain for which the Gallo results are calculated, the displacement thickness distributions differed by no more than 3% at any location and had an average difference of less than 2%. The momentum thickness differed by no more than approximately 13% with an average relative difference of less than 6%.

#### 2. Holden and Chadwick

Next we compare the results of Holden and Chadwick.^{52} We include this study because it allows us to validate our code against another CFD code GASP (General Aerodynamic Simulation Program). GASP is a finite-volume Navier–Stokes solver capable of a variety of chemical modeling options. Details of the implementation of GASP are described in the original study. Holden and Chadwick also include experimental data with approximations of associated uncertainty. This study is for a curved compression ramp with a nominal Mach number of 12, although the measured value was slightly lower.

The original purpose of this study is to capture transitional boundary layers. However, we can compare with their Run 39 for which the flow remains laminar for a significant portion of its run due to its relatively low Reynolds number of $3.211\xd7105$ per foot. In this experiment, transition begins around *x* = 24 in. As such, we expect our results to diverge significantly from the experimental data around that point and can only compare to the GASP laminar results from there on.

Using a curve fit for the pressure measurements taken over the ramp, BuBL's results closely predict the experimental data over its laminar region, accounting for the estimated 5% experimental uncertainty (Fig. 8).^{52} We also find close agreement with the GASP Navier Stokes solver with an average difference of 11%–12% between the two when neglecting the leading edge. The difference between our solution and that of the CFD solver used by Holden and Chadwick can be explained by BuBL's inexact treatment of the leading edge as well as the approximate nature of boundary layer equations as compared to full Navier–Stokes simulators like GASP.

#### 3. Holden et al.

Finally, we compare a NASA Experimental Database^{53} which gives the Holden *et al.*^{54} for a $12\xb0$ ramp geometry tested in the CUBRC facility. This study is included not only because of its status as a standard for CFD verification but also because it is particularly illuminative of the impact of the ideal gas model used. It has been established that the temperature dependence in specific heat becomes significant at high temperatures (see Sec. III), and this significance is evident in this particular work.

Although this set of experimental results is for a geometry with a compression corner, we can effectively model the laminar region of this flow as a flat plate. In BuBL, this entails specifying a constant external boundary condition. Doing so is acceptable because the ramp leading up to the compression corner itself is a flat region where the flow is known to be laminar.

Holden *et al.* provide the initial conditions for the flow, and we assume the pressure gradient in the external flow is zero. The Mach number is 11.7 and the surface boundary condition is isothermal. Leading edge effects are neglected by assuming the plate is sharp. In Fig. 9, we compare the results generated by BuBL using both the calorically perfect and thermally perfect gas assumptions. Experimental results through and past the compression corner are also shown as a reference.

It can be determined that the flow is laminar up to the recirculation region at the compression corner, where the heating drops somewhat and then spikes significantly as the transition to turbulence occurs. This transition region occurs about 15 in. downstream of the leading edge. We estimate the first five experimental data points to be in the laminar region.

In Fig. 9, we see excellent agreement with the laminar data points, especially with the four downstream, laminar data points. In the plot, error bars are included to capture a conservative estimate of experimental uncertainty as reported by Holden *et al.*^{54} and Gnoffo *et al.*^{51} The discrepancy in the first data point is expected due to BuBL's lack of treatment of the leading edge. Of particular note is the difference in BuBL's results for the constant- and variable-specific heat cases. This experiment has relatively high enthalpies. As such, we can anticipate the activation of vibrational modes of the gas resulting in notably different heating properties. As the variable specific heat case is more true to real gases, it is not surprising that its results more closely align with the Holden *et al.* experiments while the constant specific heat case tends to underpredict the heating.

## VI. CONCLUDING REMARKS

In summary, we were able to show that the Lees–Dorodnitsyn boundary layer equations are well suited for modeling high-temperature, laminar, hypersonic boundary layers. Using the second-order box scheme, finite difference computational stencil, BuBL can recover convergent solutions that can be adapted for working fluids of various properties with relatively low computational expense.

Through this study, we are also able to capture insight into the physics of hypersonic flows. Using the L–D equations and the proposed methodology, the impact of the gas model assumptions and the Prandtl number can be demonstrated. However, this is not a comprehensive study of the Prandtl number. More research would be required to determine a definitive framework for the Prandtl number selection. However, it would also be possible within the BuBL framework to allow the Prandtl number to vary as a function of the flow conditions or as a function of temperature, similar to specific heat. A temperature-dependent Prandtl number would also require further study.

Another aspect that future research must address is the impact of the model input uncertainty, as described by Roy and Oberkampf,^{55} on BuBL's accuracy. Uncertainty is reported in some of the measurement quantities featured in this paper; however, there is also uncertainty in the measurements of the nominal values of the experimental operating conditions. For example, in the latter two Holden experimental test cases,^{53,54} the freestream quantities are provided with percentage measures of uncertainty. Therefore, it is possible that disagreement with the experimental findings is due to incorrect values being input to the CFD code—the CFD simulation would not accurately represent the experimental conditions it is being compared to. Such model input uncertainty is assured to bias the results of the CFD. Gnoffo *et al.* do propose a method to quantify this uncertainty.^{51} However, we do expect the impact of this uncertainty to vary with the exact CFD implementation. In order to remedy this issue and determine the full implications of model input uncertainty, a full sensitivity analysis would need to be explored to determine to what extent resulting solutions fluctuate with varying input parameters.

In the face of experimental and model uncertainty, BuBL is still able to be verified and validated against other CFD and experimental results. As is, the proposed framework can be used for laminar boundary layer profile prediction and for the prediction of other quantities of interest which has important implications for modeling laminar hypersonic boundary layers and for future turbulent transition studies.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Chelsea Nneka Onyeador:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Alexander James Hodge:** Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Resources (supporting); Software (supporting); Validation (supporting); Writing – review & editing (supporting). **Wesley Harris:** Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Resources (supporting); Software (supporting); Validation (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

_{∞}≤ 9

*Aerodynamics of Viscous Fluids*(unpublished).