Evaluating the Hypersonic Leading-Edge Phenomena at High Reynolds and Mach Numbers

Abstract

Computational Fluid Dynamics (CFD) solutions have played an important role in the design and evaluation of complex problems where analytical solutions are not available. Among many practical applications, hypersonic flows have been an area of intense research because of the important challenges found in this flow regime. The numerical study conducted herein, focuses on solving the hypersonic flat plate problem under realistic conditions, at high Reynolds and Mach numbers. The numerical scheme implemented in this study solves the two-dimensional unsteady Navier Stokes Equations, using a novel technique called Integro-Differential Scheme (IDS) that combines the traditional finite volume and the finite difference methods. Moreover, this scheme is built on the premise of reducing the numerical errors through the implementation of a consistent averaging procedure. Unlike other numerical approaches, where free molecular effects are considered, this study enforces no-slip and fixed temperature as boundary conditions. The IDS approach accurately predicted the physics in the vicinity of the hypersonic leading edge at such realistic conditions. Even though there are slight discrepancies between the numerical solution and the available experimental data, the IDS solution revealed some interesting details about the flow field that was previously undiscovered.

1. Introduction

The flow over a flat plate is a classic yet fundamental fluid dynamic problem. Although the flow boundaries appear to be simple, the resulting flow field depends greatly on the prescribed free stream conditions. Of course, the free stream conditions are mainly defined by the Mach and Reynolds numbers as well as the ratio of specific heats. It is the ranges at which the free stream conditions are set that dictate the physics of the resulting flow field over the flat plate, and the complexities associated with it. Herein lie the many technical challenges of predicting the flat plate flow field. For example, at low Reynolds number and for subsonic Mach numbers at constant specific heats ratio of 1.4, the resulting flat plate only encourages the growth of a laminar boundary layer. Simulating such flow fields is relatively simple. As the Reynolds number increases, the boundary layer transitions to turbulent and the flow field becomes more challenging to simulate numerically. In the cases where the Reynolds number gets in the order of several million, the Mach number gets into the Hypersonic range, and the ratio of specific heat gets closer to 1.2, making the flow field interactions get complicated, and numerical simulations become unpredictable.

This chapter is concerned with the flow field over a flat plate at hypersonic conditions and at high Reynolds number. Understanding the flow field dynamics at these conditions will provide aerospace designers valuable insights into the complex interactions found in space vehicles such as rockets, space shuttles as well as military applications, for instance, hypersonic and long-range missiles. Under these conditions, the major aerodynamic concerns are aerodynamic heating and shockwave boundary layer interactions. In addition, the flow field may consist of two flow regimes; one mainly governed by the kinetic flow theory and another governed by the continuum flow theory [1]. In the case of the flat plate, especially near the tip, it is speculated that the displacement thickness increases rather drastically, causing the flow to move upward, initiating a compression shock wave and the formation of a strong interaction region. A weak interaction region follows this region. The resulting flow field becomes even more complicated because of the complex dynamics associated with the two regions. This shock is called a bow shock, due to its characteristic curvature. The region between the surface and the shock wave is called the shock layer [2], refer to Figure 1. Further, the shock layer is divided into two sublayers, each dominated by either inviscid or viscous effect. The sub-layer closest to the plate surface is known as the boundary layer, and the outer sublayer is the so-called entropy layer. Typically, the boundary layer undergoes an important transition; usually from a laminar to a turbulent boundary layer.

Figure 1.

Sketch of the flat plate problem.

Two regions can also characterize the flow along the plate; one near the leading edge and another further away. In the leading-edge region, the viscous-inviscid interactions are very strong, and they affect both sublayers: the inviscid entropy and the viscous boundary sublayers. Further, this strong interaction results in the merging of the entropy and boundary layer. In contrast, further away from the leading edge, the inviscid-viscous interaction is weak, and the two sub-layers remain separated. The two zones that are mainly characterized by the inviscid-viscous interactions are referred to as the strong and weak interaction regions, respectively. The flow phenomena in the strong and weak interaction regions at the leading edge of the hypersonic plate problem are of paramount importance to this analysis. Because of the inherent complexity of the flow physics, analytical models are scarce, and reliable analyses can only be obtained exclusively by experiments and numerical simulations.

CFD emerged as a valuable tool for these types of flow studies. Nevertheless, the CFD tool must be capable of resolving sharp gradients while negotiating systems of partial differential equations of varying types. In other words, not only are the grids expected to be extremely fine to fully resolve the sharp gradients manifesting in these regions, the CFD schemes are also expected to remain computationally stable, accurate and timely.

Many numerical solutions have been proposed. Blottner [3] solved the boundary layer problem with finite chemical reactions using finite differences. In this study, 11 chemical species and 20 reactions were considered. Another numerical study was carried out by [4], where the full time dependent Navier Stokes Equations (NSE) were solved using particle-in-cell and fluid-in-cell computing methods. In addition, the study revealed that pressure gradient is appreciable near the leading edge. Unlike reference [3], the boundary conditions used in [4] were velocity slip and temperature jump at the surface plate. These types of boundary conditions are widely applied in rarefied hypersonic flows, which have been an active area of research. These types of flows are found near the leading-edge and experimental results suggest that strong interaction theory overpredicts the surface pressure [4]. These discrepancies are attributed to the transition between continuum and free molecular flow. An extensive comparison was presented by [5], where Direct Simulation Monte Carlo (DSMC) results were compared to the NSE solution in order to evaluate the accuracy of the NSE in this regime. They concluded that including the slip conditions improved the predicted values on the surface properties. However, knowledge about the Knudsen layer is required to properly define the slip conditions at the surface. Although the Knudsen number is small in the freestream, its value is considerably high close the surface where density gradients are large [5]. Under these scenarios, the continuum hypothesis of the NSE falls apart and the accuracy of the technique is no longer ensured. Furthermore, the numerical implementation of such boundary conditions requires further simplifications and assumptions. For example, the tangential and energy accommodation factors affect the CFD solutions. Tangential accommodation values of 0.5 seem to provide accurate results near the leading edge, whereas values between 0.75 and 1.0 yield the best agreement further along the plate [6]. The same author in another publication [7] claimed that the difficulty of defining these slip conditions is in determining the correct values for the coefficients mentioned above and other empirical terms required for the implementation.

The major objective of this book chapter is to numerically solve the hypersonic flow over a flat plate problem with a novel numerical method called Integro-Differential Scheme (IDS) [8]. This study ignores the slip and jump boundary conditions introduced by [7] and directly prescribes the boundary conditions applicable to the continuum flow theory. Based on the literature review presented above, the authors of this chapter suggested that the slip and jump conditions are more appropriate for use with the Burnet equations. Further, many literature reviews suggested that the hypersonic leading-edge phenomena are best explained through the use of the transition regime, which intersects continuum and free molecular flow theories, where Burnet equations are appropriate. Although, technical evidence exists to support this hypothesis [5] this demonstration is not the focus of this book chapter.

2. The governing equation

Numerical solutions of fluid dynamic problems are governed by conservation laws. These laws can be expressed mathematically either in the differential or the integral form. In the case of compressible fluid flows, these coupled laws form a closed system of partial differential equations that is called the system of Navier-Stokes Equations (NSE). Herein, the conservation of mass, momentum and energy principles in the integral form are of interest to this study, and they are expressed as follows:

$$\frac{\partial }{\mathrm{\partial t}}\underset{\mathrm{V}}{∰}\mathrm{\rho dv}+\underset{\mathrm{S}}{∯}\left({\left(\mathrm{\rho u}\right)}_{\mathrm{i}}{\mathrm{ds}}_{\mathrm{i}}\right)=0$$E1

$$\frac{\partial }{\mathrm{\partial t}}\underset{\mathrm{V}}{∰}{\left(\mathrm{\rho u}\right)}_{\mathrm{k}}\mathrm{dv}+\underset{\mathrm{S}}{∯}\left({\left(\mathrm{\rho u}\right)}_{\mathrm{i}}{\mathrm{ds}}_{\mathrm{i}}\right){\mathrm{u}}_{\mathrm{k}}=-\underset{\mathrm{S}}{∯}\left({\mathrm{pds}}_{\mathrm{i}}\right)+\underset{\mathrm{S}}{∯}{\mathrm{\tau }}_{\mathrm{ik}}{\mathrm{ds}}_{\mathrm{i}}$$E2

$$\frac{\partial }{\mathrm{\partial t}}\underset{\mathrm{V}}{∰}\left(\mathrm{\rho e}\right)\mathrm{dv}+\underset{\mathrm{S}}{∯}\left({\left(\mathrm{\rho eu}\right)}_{\mathrm{i}}{\mathrm{ds}}_{\mathrm{i}}\right)=-\underset{\mathrm{S}}{∯}\left({\mathrm{pu}}_{\mathrm{i}}{\mathrm{ds}}_{\mathrm{i}}\right)+\underset{\mathrm{S}}{∯}\left(\left({\mathrm{\tau }}_{\mathrm{ik}}{\mathrm{u}}_{\mathrm{k}}\right){\mathrm{ds}}_{\mathrm{i}}\right)+\underset{\mathrm{S}}{∯}\left({\mathrm{q}}_{\mathrm{i}}{\mathrm{ds}}_{\mathrm{i}}\right)$$E3

In Eqs. (1)–(3) the symbols: $\rho ,$ $u,$ $t$ represent the density, the velocity components of an elementary control fluid element, and time, respectively. In addition, the symbols $e,$ $p,$ ${\tau }_{\mathit{ij}}$ and $q_i$ in Eqs. (1)–(3) represent the internal energy, pressure, the stress tensor and the heat flux associated with an elementary control volume, respectively. Internal energy, pressure, stress tensor and heat flux are defined by Eqs. (4)–(7)

$$\mathrm{e}={\mathrm{C}}_{\mathrm{v}}\mathrm{T}+\frac{1}{2}{\mathrm{u}}_{\mathrm{k}}{\mathrm{u}}_{\mathrm{k}}$$E4

$$\mathrm{p}=\mathrm{\rho RT}$$E5

$${\mathrm{\tau }}_{\mathrm{ik}}=\mathrm{\mu }\left(\frac{\partial {\mathrm{u}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{k}}}+\frac{\partial {\mathrm{u}}_{\mathrm{k}}}{\partial {\mathrm{x}}_{\mathrm{i}}}-\frac{2}{3}{\mathrm{\delta }}_{\mathrm{ik}}\frac{\partial {\mathrm{u}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right)$$E6

$${\mathrm{q}}_{\mathrm{k}}=-\mathrm{k}\frac{\mathrm{\partial T}}{\partial {\mathrm{x}}_{\mathrm{k}}}$$E7

In Eq. (5), $R$ is the gas constant. The symbols $\mu$ and $k$ represent the viscous and thermal properties of the fluid of interest. For air, the viscosity of the fluid is evaluated using Sutherland’s law and the thermal conductivity expression,

$$\mathrm{k}=\left(4.76\times {10}^{-6}\right){\mathrm{T}}^{3/2}/\left(\mathrm{T}+112.0\right)$$E8

In the case of 3D aerothermodynamics, the NSE (1–8) represent a closed system of five equation relative to five unknowns. These unknowns are called Primitive Variables (PV), and are defined in the vector form as follows:

$$\mathrm{PV}≔{\left[\mathrm{\rho }\mathrm{u}\mathrm{v}\begin{array}{ll}\mathrm{w}& \mathrm{T}\end{array}\right]}^T$$E9

The goal of any numerical solution to the NSE is to determine the primitive variables at every grid point. However, obtaining a unique solution to the NSE (1–8) requires the prescription of initial and boundary conditions.

Of course, the full set of the NSE (1–8) does not readily lend itself to analytical solutions. It is only in recent decades, with the advent of modern computers that the non-existence of analytical solutions to the NSE ceased to be a limitation to our understanding of the physics underpinning flow fields. Modern computers also gave birth to the many modern numerical methods capable of solving the NSE. Among these methods are the Conservation-Element Solution Element (CESE) method, Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), Discontinuous Galerkin Methods (DGM) and the Integro-Differential Scheme (IDS). These computational methods have all been applied to the task of solving the NSE, and have all provided varying degrees of success when it comes to elucidating the details associated with the various flow field physics of interest to the engineers at realistic Reynolds and Mach numbers. This report highlights the IDS procedure of solving the NSE.

3. The IDS fluid model

Consider the Integro-Differential Model (IDM) as it is applied to the computational solution to the NSE (1–8). In general, the IDS solution of a given fluid dynamic problem is built on an interconnecting set of spatial and temporal fluid cells. In the Cartesian system of coordinates, a typical fluid cell is nothing more than a carefully chosen elementary rectangular prism, defined by the dimension; dx, dy and dz. It is the application of a specified fluid cell in relationship to the NSE equations that determines whether it becomes a spatial or a temporal cell. Nevertheless, consider the fluid cell illustrated in Figure 2 where its implementation in the NSE is irrelevant at this time.

Figure 2.

Spatial cell with notation at surface nodes.

3.1. IDS cell properties

The Cartesian cell defined in Figure 2, has the following properties:

1. The rectangular prism has 6 elementary surfaces and each surface is defined through the use of 3 directional normals. Further, each normal is defined in either a positive or a negative direction in relationship to its respective axis. This definition is not unique to the IDM, but it is used here for illustrative purposes only

2. Each of the six elementary surfaces, ds, of the fluid cell, is considered a vector and is defined as follows: ${\mathrm{\delta s}}_{\mathrm{i}}^{\pm }={\mathrm{n}}_{\mathrm{i}}^{\pm }{\mathrm{d\xi }}_{\mathrm{j}}{\mathrm{d\eta }}_{\mathrm{k}}$, where the indices i, j, and k vary from 1 to 3, representing the x, y or z coordinate direction. In addition, the area is considered a vector, having both direction and magnitude. Likewise, the volume of the elementary cell can be defined by $\mathrm{\delta v}=\text{dxdydz}$ for uniform grids. Note, this type of evaluation also works well for non-uniform and unstructured grids.

3. The fluid cell defined in Figure 2 also allows for mass, momentum and energy fluxes to traverse its surfaces. At any given instance, the net spatial fluxes traversing a given surface are defined by a combination of their inviscid and viscous counterparts. The inviscid and viscous fluxes on the cell surfaces are defined by the symbols: E, E_vis, F, F_vis, G, and G_vis, representing the inviscid and viscous fluxes in the x, y and z directions, respectively.

4. In accordance with the IDM, the flux values are approximated from their edge quantities using their arithmetic averages, and are assumed to be located at the center of the respective surfaces. Consequently, all quantities evaluated on any of the cell surfaces are labeled as averaged quantities. The fluxes of interest on the cell surfaces are the average inviscid and viscous spatial fluxes on the ${\mathrm{n}}_{\mathrm{x}}^{\pm }$ surfaces, and they are defined by the symbols: ${\mathrm{E}}_{\mathrm{Avg}}^{\text{Surf},\mathrm{nx}\pm }$ and ${\mathrm{E}}_{\mathrm{Vis},\mathrm{Avg}}^{\text{Surf},\mathrm{nx}\pm }$. Likewise, the average inviscid and viscous fluxes on the ${\mathrm{n}}_{\mathrm{y}}^{\pm }$ and ${\mathrm{n}}_{\mathrm{z}}^{\pm }$ surfaces are defined by the symbols: ${\mathrm{F}}_{\mathrm{Avg}}^{\text{Surf},\mathrm{nx}\pm }$, ${\mathrm{F}}_{\mathrm{Vis},\mathrm{Avg}}^{\text{Surf},\mathrm{nx}\pm }$, ${\mathrm{G}}_{\mathrm{Avg}}^{\text{Surf},\mathrm{nz}\pm }$ and ${\mathrm{G}}_{\mathrm{Vis},\mathrm{Avg}}^{\text{Surf},\mathrm{nz}\pm }$, respectively. Having established the fact that the average cell flux quantities can only be defined on one of its elementary surfaces, the expressions needed for the evaluation of the spatial fluxes can now be summarily expressed as,

   $$\mathrm{E}={\mathrm{\delta s}}_{\mathrm{i}}^{\pm }{\left\{\begin{array}{l}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2+\mathrm{p}\\ \mathrm{\rho uv}\\ \mathrm{\rho uw}\\ \left(\mathrm{\rho }{\mathrm{e}}_{\mathrm{T}}+\mathrm{p}\right)\mathrm{u}\end{array}\right\}}_{{\mathrm{Avg}}^{\mathrm{S}}},{\mathrm{E}}_{\mathrm{Vis}}={\mathrm{\delta s}}_{\mathrm{i}}^{\pm }{\left\{\begin{array}{l}0\\ {\mathrm{\tau }}_{\mathrm{xx}}\\ {\mathrm{\tau }}_{\mathrm{xy}}\\ {\mathrm{\tau }}_{\mathrm{xz}}\\ {\mathrm{u\tau }}_{\mathrm{xx}}+{\mathrm{v\tau }}_{\mathrm{xy}}+{\mathrm{w\tau }}_{\mathrm{xz}}-{\mathrm{q}}_{\mathrm{x}}\end{array}\right\}}_{{\mathrm{Avg}}^{\mathrm{S}}}$$E10

   where the subscripts in the right-hand terms represent the location and type of operations used in evaluating the required average quantities. Similarly, for the inviscid and viscous fluxes, F, F_vis, G, and G_vis.

5. In a similar manner, the average quantities within the cell volume, such as the time fluxes, ${\mathrm{U}}_{\mathrm{Avg}}^{\text{Cell}}$, and the rate of change of the time fluxes, ${\left(\partial U/\partial t\right)}_{\mathit{Avg}}^{\mathit{\text{Cell}}},$ are defined as

$${\mathrm{U}}_{\mathrm{Avg}}^{\text{Cell}}={\left\{\begin{array}{l}\mathrm{\rho }\\ \mathrm{\rho u}\\ \mathrm{\rho v}\\ \mathrm{\rho w}\\ \mathrm{\rho }{\mathrm{e}}_{\mathrm{T}}\end{array}\right\}}_{{\mathrm{Avg}}^{\mathrm{v}}},{\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}^{\text{Cell}}=\frac{\partial }{\mathrm{\partial t}}{\left\{\begin{array}{l}\mathrm{\rho }\\ \mathrm{\rho u}\\ \mathrm{\rho v}\\ \mathrm{\rho w}\\ \mathrm{\rho }{\mathrm{e}}_{\mathrm{T}}\end{array}\right\}}_{{\mathrm{Avg}}^{\mathrm{v}}}$$E11

3.2. Computing the cell properties

At this point, the concept of evaluating the volume and surface areas of an elementary fluid cell is fully formulated. However, the computation of the average time and spatial fluxes within and on the surface of an IDS cell is still not uniquely defined. For example, how are the primitive variables, which are defined at a point and the average time fluxes within the cell or the spatial fluxes on an elementary surface related? How is the flow field defined in relationship to the IDS fluid cell concept? To answer these questions and others, consider Figure 2 once more. Assume the red dots in Figure 2 represent the physical grid points in the flow field of interest, and at each of these points, the primitive variables are known. A typical elementary fluid cell is then built around eight such points, with each point separated by no more than one grid point. Using these assumptions, the average time fluxes within the cell can be computed from the arithmetic mean, as:

$${\mathrm{U}}_{\mathrm{Avg}}^{\text{Cell}}=\frac{1}{8}\left[{\left\{\begin{array}{l}\mathrm{\rho }\\ \mathrm{u}\\ \mathrm{v}\\ \mathrm{w}\\ \mathrm{T}\end{array}\right\}}^1+{\left\{\begin{array}{l}\mathrm{\rho }\\ \mathrm{u}\\ \mathrm{v}\\ \mathrm{w}\\ \mathrm{T}\end{array}\right\}}^{1\text{'}}+\dots +{\left\{\begin{array}{l}\mathrm{\rho }\\ \mathrm{u}\\ \mathrm{v}\\ \mathrm{w}\\ \mathrm{T}\end{array}\right\}}^{4\text{'}}\right]$$E12

Similarly, cell surface quantities, such as, the ${\mathrm{n}}_{\mathrm{x}}^{+}$ surface fluxes defined by ${\mathrm{E}}_{\mathrm{Avg}}^{\text{Surf},\mathrm{nx}+}$ can be computed as follows:

$${\mathrm{E}}_{\mathrm{Avg}}^{\text{Surf},\mathrm{nx}+}=\frac{{\mathrm{ds}}_{\mathrm{nx}}^{+}}{4}\left[{\left\{\begin{array}{l}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2+\mathrm{p}\\ \mathrm{\rho uv}\\ \mathrm{\rho uw}\\ \mathrm{u}\left(\mathrm{\rho e}+\mathrm{p}\right)\end{array}\right\}}_{\mathrm{nx}+}^3+{\left\{\begin{array}{l}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2+\mathrm{p}\\ \mathrm{\rho uv}\\ \mathrm{\rho uw}\\ \mathrm{u}\left(\mathrm{\rho e}+\mathrm{p}\right)\end{array}\right\}}_{\mathrm{nx}+}^{3\text{'}}+{\left\{\begin{array}{l}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2+\mathrm{p}\\ \mathrm{\rho uv}\\ \mathrm{\rho uw}\\ \mathrm{u}\left(\mathrm{\rho e}+\mathrm{p}\right)\end{array}\right\}}_{\mathrm{nx}+}^4+{\left\{\begin{array}{l}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2+\mathrm{p}\\ \mathrm{\rho uv}\\ \mathrm{\rho uw}\\ \mathrm{u}\left(\mathrm{\rho e}+\mathrm{p}\right)\end{array}\right\}}_{\mathrm{nx}+}^{4\text{'}}\right]$$E13

In a similar manner, the other inviscid flux quantities at the ${\mathrm{n}}_{\mathrm{x}}^{+}$, ${\mathrm{n}}_{\mathrm{y}}^{\pm }$, and ${\mathrm{n}}_{\mathrm{z}}^{\pm }$ surfaces can be found. Unfortunately, computing the viscous flux quantities; ${\mathrm{E}}_{\mathrm{Vis},\mathrm{Avg}}^{\text{Surf},\mathrm{nx}\pm }$, ${\mathrm{F}}_{\mathrm{Vis},\mathrm{Avg}}^{\text{Surf},\mathrm{nx}\pm }$, and ${\mathrm{G}}_{\mathrm{Vis},\mathrm{Avg}}^{\text{Surf},\mathrm{nz}\pm }$, are somewhat complicated and greater care is required. Refer to Ref. [8] for details. In summary, the IDM allows for the grid points, and the primitive variables allocated at those points to be used in uniquely formulating the elementary fluid cells and completely defining their flow field characteristics.

3.3. The IDS control volume and its properties

A typical IDS control volume in the 3D Cartesian system of coordinates is illustrated in Figure 3.

Figure 3.

Illustration of control volume.

As can be observed, the IDS control volume consists of eight cells. The properties of the IDS control volume is as follows:

1. It is of interest to note that the centers of the eight cells form the vertex of an overlapping cell. This overlapping cell is called the temporal cell. As such, the control volume consists of eight neighboring cells that allow for the formation of a temporal cell. In other words, analogous to the manner in which the eight grid points formed the vortex of a given cell, so too, the center of eight neighboring cells formed the vortex of a temporal cell.

2. Also of interest to note is the fact that at the vertex of the temporal cell, the rate of change of the time fluxes are known. Consequently, at the center of the temporal cell, the average time rate of change of the time fluxes are computed as,

   $${\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}^{\mathrm{CV}}=\frac{1}{8}\left[{\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}^{{\text{Cell}}_1}+{\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}^{{\text{Cell}}_2}+..\dots +{\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}^{{\text{Cell}}_8}\right]$$E14

3. Similarly, at the center of an IDS control volume, the average time fluxes are defined by the arithmetic averages

   $${\mathrm{U}}_{\mathrm{Avg}}^{\mathrm{CV}}=\frac{1}{8}\left[{\mathrm{U}}_{\mathrm{Avg}}^{{\text{Cell}}_1}+{\mathrm{U}}_{\mathrm{Avg}}^{{\text{Cell}}_2}+..\dots +{\mathrm{U}}_{\mathrm{Avg}}^{{\text{Cell}}_8}\right]$$E15

4. The IDS solution can only be advanced from within the temporal cells. Within a given temporal, the average temporal fluxes are updated as follows:

   $${\mathrm{U}}_{\text{Updated},\mathrm{Avg}}^{\mathrm{CV}}={\mathrm{U}}_{\mathrm{Avg}}^{\mathrm{CV}}+{\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}^{\mathrm{CV}}\mathrm{\delta t}$$E16

   where the symbol $\mathit{\delta t}$ is the time increment, and methods for its computations are defined later in this report.

5. In the 3D Cartesian system of coordinates, the grids can be developed such that the center of the IDS control volume always overlaps with the center of a grid point. If this were the case, then the updated primitive variables can be computed from the expression

   $${\left\{\begin{array}{l}\mathrm{\rho }\\ \mathrm{u}\\ \mathrm{v}\\ \mathrm{w}\\ \mathrm{T}\end{array}\right\}}_{\text{Updated}}={\left\{\begin{array}{l}{\mathrm{U}}_1\\ {\mathrm{U}}_2/{\mathrm{U}}_1\\ {\mathrm{U}}_3/{\mathrm{U}}_1\\ {\mathrm{U}}_4/{\mathrm{U}}_1\\ \left(\frac{{\mathrm{U}}_5}{{\mathrm{U}}_1}-\frac{1}{2}\left(\frac{{\mathrm{U}}_2^2}{{\mathrm{U}}_1^2}+\frac{{\mathrm{U}}_3^2}{{\mathrm{U}}_1^2}+\frac{{\mathrm{U}}_4^2}{{\mathrm{U}}_1^2}\right)\right)\mathrm{\gamma }\left(\mathrm{\gamma }-1\right){\mathrm{M}}_{\infty }^2\end{array}\right\}}_{\text{Updated},\mathrm{Avg}}^{\mathrm{CV}}$$E17

It is with this IDM in mind, where the spatial cell, the temporal cell and the control volume concepts are paramount, that the NSE are transformed into their Integro-Differential counterparts. In addition to this, Eqs. (14)–(16) highlight the main differences between the IDS and other CFD schemes, where the solution vector and the time rate are computed using local values, whereas the IDS computes the right-hand side terms from Eq. (16) through the use of Eqs. (14) and (15). From the computational perspective, the IDS performs more floating points operations per node and therefore, the method is computationally expensive. Roughly speaking, for 2D flows the method performs 8 times more floating-point operations. Nevertheless, the calculation of the viscous stresses and heat fluxes, Eq. (10), demands the evaluation of the stresses and heat fluxes at the faces of the control volumes, where the mean value theorem is used and hence, an extra averaged is required. Off course, for 3D fluid flows these operation increases because of the spanwise component where the averaging procedure is also implemented.

4. The Integro-differential scheme

Consider the IDM described in the preceding section for the special 3D Cartesian system of coordinates. If the NSE (1–3) were directly applied to a fluid element, such as the one illustrated in Figure 2, the analytical solution arising from this process, especially when the mean value principles are invoked, will yield the following transformational equations:

$${\left(\frac{\mathrm{\partial \rho }}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}=-\left[\sum _{\mathrm{m}=1}^6{\left(\text{ρuδs}\right)}_{\mathrm{m}}^{\pm }\right]/\mathrm{\delta v}$$E18

$${\left(\frac{\partial \left({\mathrm{\rho u}}_{\mathrm{k}}\right)}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}=-\left[\sum _{\mathrm{m}=1}^6{\left(\left({\mathrm{\rho u}}_{\mathrm{i}}+{\mathrm{p\delta }}_{\mathrm{ik}}\right)\mathrm{\delta s}\right)}_{\mathrm{m}}^{\pm }{\mathrm{u}}_{\mathrm{k}}-\sum _{\mathrm{m}=1}^6{\left({\mathrm{\tau }}_{\mathrm{ik}}\mathrm{\delta s}\right)}_{\mathrm{m}}^{\pm }{\mathrm{u}}_{\mathrm{k}}\right]/\mathrm{\delta v}$$E19

$${\left(\frac{\mathrm{\partial }\left(\mathrm{\rho e}\right)}{\mathrm{\partial t}}\right)}_{\mathrm{Avg}}=-\left[\sum _{\mathrm{m}=1}^6{\left(\left(\mathrm{\rho e}+\mathrm{p}\right){\mathrm{u}}_{\mathrm{i}}\mathrm{\delta s}\right)}_{\mathrm{m}}^{\pm }-\sum _{\mathrm{m}=1}^6{\left({\mathrm{\tau }}_{\mathrm{ik}}{\mathrm{u}}_{\mathrm{k}}\mathrm{\delta s}\right)}_{\mathrm{m}}^{\pm }-\sum _{\mathrm{m}=1}^6{\left({\mathrm{q}}_{\mathrm{i}}\mathrm{\delta s}\right)}_{\mathrm{m}}^{\pm }\right]/\mathrm{\delta v}$$E20

where the index m, m = 1, 6, defines the surfaces with positive and negative normal, respectively, and the indices i and k, go from 1 to 3, defining the coordinate directions. The technical challenges in computing Eqs. (14)–(16) lie in the careful and consistent manner in which the NSE auxiliary/closure Eqs. (4)–(8) are evaluated as they are applied to a fluid cell. It is worthwhile to repeat that Eqs. (4)–(8) must be evaluated in accordance with the requirements of each cell.

In the special case of 3D Cartesian systems with the spatial flow field domain defined on uniform cells, the viscous and inviscid fluxes can be expressed in vector form, as:

$${\left(\frac{\mathrm{\partial U}}{\mathrm{\partial t}}\right)}_{\mathit{Avg}}+\frac{{\mathrm{\Delta }}_{\mathrm{nx}}^{\pm }\left(\mathrm{E}-{\mathrm{E}}_{\mathrm{vis}}\right)}{\mathrm{\Delta \nu }}+\frac{{\mathrm{\Delta }}_{\mathrm{ny}}^{\pm }\left(\mathrm{F}-{\mathrm{F}}_{\mathrm{vis}}\right)}{\mathrm{\Delta \nu }}+\frac{{\mathrm{\Delta }}_{\mathrm{nz}}^{\pm }\left(\mathrm{G}-{\mathrm{G}}_{\mathrm{vis}}\right)}{\mathrm{\Delta \nu }}=0$$E21

where the U, E, E_vis, F, F_vis, G and G_vis vectors were defined earlier. The subscripts in Eq. (17) define the location and type of operations used in evaluating respective average quantities. In addition, the difference operators; ${\mathrm{\Delta }}_{\mathrm{nx}}^{\pm }$, ${\mathrm{\Delta }}_{\mathrm{ny}}^{\pm }$, and ${\mathrm{\Delta }}_{\mathrm{nz}}^{\pm }$ represent the difference in the surface fluxes across each cell, such that, the surface information of each cell are independently computed.

In summary, the integral form of the NSE equations (1–3) were analytically solved using the mean value theorem over an elementary control volume. The resulting solution was expressed in the form an Integro-differential formulation as described by Eq. (17). Finally, as in all explicit schemes, Eq. (16) can be used to compute the update solution vector U, such that, where the time fluxes and rate of change of the time fluxes vectors were defined in Eqs. (14) and (15). The symbol $\mathrm{\delta t}$ is the time increment and is computed by using the Courant-Friedrichs-Lewy (CFL) criterion. In this book chapter, the CFL criterion is computed with the aid of the expression [2],

$$\mathrm{\delta t}=C{\left[\frac{\left|u\right|}{\mathrm{\Delta }x}+\frac{\left|v\right|}{\mathrm{\Delta }y}+\frac{\left|w\right|}{\mathrm{\Delta }y}+a\sqrt{\frac{1}{\mathrm{\Delta }{x}^2}+\frac{1}{\mathrm{\Delta }{y}^2}+\frac{1}{\mathrm{\Delta }{z}^2}}+\frac{2}{{\text{Re}}_L}\nu \left(\frac{1}{\mathrm{\Delta }{x}^2}+\frac{1}{\mathrm{\Delta }{y}^2}+\frac{1}{\mathrm{\Delta }{z}^2}\right)\right]}^{-1}$$E22

where $a$ is the local speed of sound, C is the Courant number, and $\gamma$ is the specific heat ratio, and $\mathrm{\nu }$ is computed from the expression

$$\mathrm{\nu }=\text{max}\left[\left(\frac{4}{3}\mathrm{\mu },\left(\mathrm{\gamma \mu }/\text{Pr}\right)\right)/\mathrm{\rho }\right]$$

The typical values used for C in this analysis range as follows: $0.5\le C\le 0.8$.

The IDS offers an important numerical advantage given that it is a very stable and accurate method. Further, the IDS provides a significant reduction in both spatial and temporal numerical dispersion through its use of the mean value theorem in computing the finite volume quantities. In the IDS approach, the time and spatial fluxes are appropriately approximate. In addition, the method has been shown to be consistent and has a minimum discretization error of order p = 2. A major advantage of the IDS method is that computations involving the compressible NSE are very stable, and no numerical oscillations are typically detected when the grid is fully refined. Another advantage of the IDS method is that it is suitable for solving both steady and unsteady flows at realistic Reynolds and Mach numbers. Experiences have shown that it is quite a satisfactory method for solving high Reynolds number flows, where the viscous regions are very thin, and shock boundary layer interactions are significant. For these flows, once the mesh is highly refined the IDS is able to resolve the flow physics within the viscous regions with a great order of accuracy.

5. The hypersonic flow over a flat plate

Consider the hypersonic flow over a 1 meter long flat plate at zero angle of attack. The freestream Mach number is set at 8.6, the Reynolds number set to 3.4757 × 10⁶ (based on the plate length), the Prandtl number to 0.70 and the specific heat ratio, $\gamma$, to 1.4. These conditions are similar to those presented by [9] except that their experiments considered a slightly longer plate. Nevertheless, in comparison to this effort, the Reynolds number used in [9] was in perfect agreement, but the Mach number was slightly greater (approximately 0.93%). The boundary conditions were set as follows: free stream values were assigned to the inflow and the far field boundaries, the interior flow primitive variables were extrapolated to the exit plane, and three separate sets of conditions were assigned to the base of the domain. Symmetric boundary conditions were assigned to the leading and trailing edge gaps, and a combination of no-slip and fixed wall temperature assigned to the solid wall. It is of interest to mention here, in this effort the dimensionless temperature is set to 1.0 at the wall, whereas it was set to a value of 0.828 in the experimental study [9]

In efforts to obtain a grid independent solution, five sets of grids of sizes ranging from 1001 by 1001 nodes to 5001 by 16001 nodes in the streamwise and vertical direction, respectively, were studied. Since the gradients of the flow field parameters in the direction normal to the wall are sharper than those in the direction parallel to the wall, substantially finer grids were placed in the vertical direction [2]. In addition to the five sets of grids described above, an extra case, termed the modified grid, was also considered. This was done in efforts to more thoroughly evaluate the IDS capabilities in predicting the viscous dissipation effects inside the boundary layer. In the case of the modified grid, the height of the domain was reduced by half, resulting in an equivalent grid size of 6001 by 32001 nodes. This reduction in height effectively reduced the cell height by a factor of 2, resulting in a finer set of grids without a substantial increase in computational load.

The IDS flow field solutions resulting from the grid independence study are summarized in the Figures 4 and 5. Note, that in Figures 4 and 5, the modified grid is represented by the grid size of 6001 by 16001*. Note, Figure 4 illustrates the behavior of the streamwise velocity component in the y-direction at a location of x = 0.5 m from the leading edge. A careful observation of the data demonstrates that the height of the boundary layer is approximately 0.0088 of non-dimensional units, representing 8.8 mm based upon $0.99U_{\infty }$. In addition, no evidence of shock is shown in Figure 4. Similarly, Figure 5 depicts the temperature profile at 0.5 m from the leading-edge. However, in this case, the effect of viscous dissipation within the boundary layer is clearly demonstrated [10]. As noted in Figure 5, the temperature increases from the outer edge of the boundary layer toward the wall, reaching two peaks. The outer peak indicates the presence of a shock, while the inner peak indicates the effects of boundary layer dissipation due to viscous friction. Similar trends were also found in [11, 12].

Figure 4.

U velocity profile at 0.5*L.

Figure 5.

Temperature profile at 0.5*L.

The grid study data suggest that grid independence was obtained for a grid size of 5001 by 16001. A closer observation of Figure 5, reveals the existence of the two expected sub-layers; namely the entropy layer and the viscous boundary layer. As supported by Figure 5, although the height of the shock wave was fully resolved with mesh sizes; 4001 by 8001; 5001 by 16001 and 10001 by 16001*, the dissipation effects were clearly not. Herein, the conclusion is the boundary layer needs an extremely finer set of grids to resolve its physics when compared to the mesh size needed to resolve the entropy layer and the shock wave.

In this chapter, the normal Mach number is measured in the direction of the pressure gradient, and computed with the aid of Eq. (23), which is fully described in [13]. Since pressure and density are the two variables that produce the greatest change as they traverse discontinues in the flow field, they are also very efficient in detecting shocks. In some cases, however, false indications may occur, so a small degree of filtering is required [14]. The filtering criteria proposed by [14] with a threshold of ϵ = 0.007 is used in Eq. (23).

$${\mathit{Ma}}_n=\frac{\mathit{Ma}\cdot \nabla p}{\left|\nabla p\right|}=1$$E23

Consider the IDS flow field distribution of the normal Mach number computed from expression (23) illustrated in Figure 6. Using a filtering threshold of ϵ = 0.007, the bow shock, the viscous boundary layer, entropy layer and other flow field features are extracted. It can be observed that the bow shock starts slightly ahead of the leading-edge tip and it certainly displays the characteristic curvature. This characteristic was also reported in [1]. In addition, Figure 6 accurately predicts the growth of the boundary layer. More importantly, the similarities of the predicted features illustrated in Figure 1 closely matches the IDS computed results illustrated in Figure 6. Of greater significance is the fact that the two interaction zones; namely the strong and weak inviscid-viscous interactions zones merged, and all are vividly computed.

Figure 6.

Normal Mach number.

Now, consider the Q-criterion introduced by [15] and computed with the aid of expression,

$$Q=\frac{1}{2}\left({‖\mathrm{\Omega }‖}^2-{‖S‖}^2\right)$$E24

Note that the symbols $‖\mathrm{\Omega }‖$ and $‖S‖$ in Eq. (24) represent the Euclidean norm of the vorticity and rate of strain tensor, respectively. Further, Eq. (24) is an effective tool for the extraction of the vector field topology that represents the local balance between the rate of shear strain and vorticity. When Eq. (24) is applied to the IDS flow field solution over the flat plate, the regions with dominant rates of strains and vorticity are revealed. The Q-criterion results are documented in Figure 7. Again, as observed the bow shock, the viscous boundary layer, entropy layer and other flow field features defined by the vortical structures are effectively captured.

Figure 7.

Q-criterion contour plot.

In efforts to closely observe the flow physics in the weak interaction regime, the variation of the Q-criterion in the y-direction is plotted and illustrated in Figure 8. Note, Figure 8 provides quantitative information about Q-criterion at the location x = 0.27 along the plate and merely complements the information already presented in Figure 7. Nevertheless, Figure 8 reveals that the strain dominates over vorticity very near the wall; as the Q-criterion becomes negative as it approaches the wall. The Q-criterion behavior observed in this analysis is typical within viscous sub layers where the shear stress is laminar [16]. A second layer, the so-called turbulent layer where the swirling motions are common causes the Q-criterion to turn positive. In this region, the viscous effects contribute to large increases in entropy, and consequently vorticity [17]. Thus, this layer is characterized by positive values of Q-criterion. Moving deeper into the flow field, vertically above the wall, the inviscid region is revealed. In this region, the rate of strain dominates over the rate of rotation, and the Q-criterion gets deeper in the negative direction, only to be reversed as the shock wave is penetrated.

Figure 8.

Q-criterion profile at 0.27 from the leading edge.

A close-up investigation of the flow physics at the hypersonic leading edge was conducted. To this end, an extra case was analyzed where the length and height of the domain were reduced, while the dimensionless parameter, such as Reynolds and Mach numbers were kept constant. Boundary conditions and the freestream values of the primitive variables were the same from the previous analyses. However, unlike the previous solution, where the full plate of length 1 m was studied, the following results do not consider trailing edge. The focus of this analysis is to analyze the flow physics near the leading edge thoroughly. The length and height were selected as 0.1 and 0.02 m, respectively and the leading-edge gap was defined as 0.1 dimensionless units. Figure 9 illustrates IDS prediction in the form of the Q-Criterion at the leading edge. These results indicate that the rate of strain is the dominant effect at the leading-edge followed by large rotational motions that cause a delay in the growth of the boundary layer. Figure 9 also shows that the bow shock wave is formed ahead of the plate.

Figure 9.

Q-criterion contour plot at the leading edge.

It is of interest to note that at the tip of the plate, the region with the greatest rate of strain within the flowfield is observed, albeit a small region. Immediately following this region there is a similarly small region with the greatest rate of rotation, refer to Figure 10. The two regions with the greatest rate of strain and rate of rotation are located at the leading edge, and it is from these two regions that the shock wave and the boundary layer respectively emanate. At this point, the technical details associated with the growth of these sub-layers is not fully resolved, and as such, further analysis is warranted.

Figure 10.

Horizontal profile of Q-criterion at y = 2 × 10⁻⁵.

5.1. Validating the IDS hypersonic leading-edge solution

The experimental data developed by [9] was used to validate the accuracy of the IDS solution to the Hypersonic flat plate problem. In Figure 11, the temperature and u-velocity profiles at 80 mm from the leading edge for both cases: experimental and IDS, respectively, are compared. Note that Figure 11 shows that the IDS solution under-predicts the maximum temperature inside the boundary layer according to the experimental data. Figure 11 also shows the comparative behavior of the horizontal component of the velocity vector for the experimental study and the IDS solution.

Figure 11.

Experimental and IDS solution.

Obviously, there are differences in the experimental and IDS solutions. The reasons for these differences were analyzed. First, it turns out that the required freestream conditions used during the experiment were not public and could not be easily reproduced. Reference [9] described that the study was carried out in a divergence nozzle, and as such, the flat plate experienced a somewhat favorable pressure gradient during the course of the experiment. Second, the freestream conditions were computed using a computational package called STUBE [18], that used the piston theory to estimate the pressure in the reservoir; a huge approximation. Other discrepancies were found between the experimental data and the STUBE predictions [9]. Moreover, [9] also compared the experimental data with the numerical solution from a commercial CFD package called CFD-FASTRAN. The CFD simulations used the freestream conditions calculated by STUBE, and even then, the freestream velocity was scaled to match the measured external velocity [9]. Nevertheless, even with these uncertainties, the data presented in Figure 5 is considered to be reasonably close and the IDS code is considered validated.

Figure 12 illustrates the vertical profiles of the density and the y-component of the velocity vector at 80 mm from the leading edge, respectively. However, unlike in the case of the temperature and U velocity profiles, there are no equivalent experimental data available for direct comparison. However, as can be observed in Figure 12 the three most important phenomena, namely the boundary layer, the entropy layer and the shock wave along with their respective flow filed characteristics are distinctly captured.

Figure 12.

Density and vertical velocity profile.

5.2. Approximate boundary layer analysis

The boundary layer thickness in the absence of adverse pressure gradient can be computed as [19]:

$$\frac{\delta }{x}=\left[5.0+\left(0.2+0.9\frac{T_w}{T_{\mathit{aw}}}\right)\left(\gamma -1\right)M_e^2\right]\sqrt{\frac{C_w}{{\text{Re}}_x}}$$E25

Where C_w represent the Chapman-Rubesin parameter, Re_x is the Reynolds number at the measurement point and T_aw, the adiabatic wall temperature, is given by:

$$T_{\mathit{aw}}=T_e\left(1+\sqrt{\text{Pr}}\left(\frac{\gamma -1}{2}\right)M_e^2\right)$$E26

Equations (25) and (26) were used to compute the theoretical boundary layer, whose value is 1.87 mm. The boundary layer predicted by the IDS was measured to be 2.83 mm. The authors are still evaluating the reasons for the resulting discrepancy.

5.3. Investigating the strength of the hypersonic leading-edge interaction

An important phenomenon that is observed within the flow field in the vicinity of the leading edge is a rise in pressure. This high pressure emanates at the edge and it leads to an induced pressure gradient along the rest of plate. Therefore, the assumption of zero pressure gradient condition through the boundary conditions is debatable [10]. Today’s experiments have shown that under hypersonic flow conditions, the leading-edge experiences a bow shock, a significant pressure rise, very large skin friction, and very large heat flux. Further, these experimental observations have shown that adverse conditions are confined only to the leading-edge. To fully explore the leading edge behavior, researchers [19] introduced the shock interaction parameter, $\chi$, which is defined as:

$$\chi ={\mathit{Ma}}_{\infty }^3\left(\frac{C_w}{{\text{Re}}_{\infty }}\right)$$E27

This parameter, $\chi$, is a function of the freestream Mach number, Reynolds number and the specific heats ratio, and it serves to quantify the strength of the flow field interaction at the leading edge. Two sets of the shock interaction parameter, $\chi$, ranges are of interest to this study; weak interactions as described by $\chi <<1$ and strong interactions as described by $\chi >>1$. Figure 13 shows the experimental data from [20] indicating the strength of the shock interactions, $\chi$ at Mach numbers ranging from 5 to 10 to 20, as indicated by the green, blue and cyan colors. It is important to mention that the horizontal axis is described by the inverse of $\mathit{A\chi }$, where A is an aerodynamic parameter defined by

$$A=\frac{1}{2}\left(\gamma -1\right)0.664\left(1+2.6T_w/T_0\right)$$E28

Figure 13.

Induced pressure near the leading edge.

This transformation is done in efforts to allow for small values of $\mathit{A\chi }$ to map with the strong interaction regions and for large values of $\mathit{A\chi }$ to map with the weak interaction regions.

Using the data provided by [20], the weak interaction curve, highlighted in black, is recovered. Finally, the IDS solution is also plotted in Figure 13 and it is depicted by the red curve. As observed in Figure 13, the IDS solution matches the predictions governed by the interaction parameter, $\mathit{A\chi }$, and as noted, lies closest to the Mach 5 interaction curves.

The IDS hypersonic flat-plate solution confirms that the scheme is capable of accurately resolving the complex flow physics within the vicinity of the hypersonic leading edge, and it correctly qualifies the inviscid-viscous interactions associated with this region. However, it is important to note that the IDS scheme is currently being upgraded with OpenMP and MPI capabilities in efforts to handle three-dimensional flows. It is the author’s opinion that the mechanism driving turbulence dynamics is three dimensional in nature, and therefore cannot be captured by two-dimensional flow fields. As a result, in 2D flow fields, vortex stretching and other fluid mechanisms important to the development of turbulence flows are absent. Once these capabilities are validated, the hypersonic flat-plate problem will be revisited.

Acknowledgments

The authors thank Maida Bermudez Bosch M.A., CCC-SLP for her help on editing and proofreading the written article as well as Professor Kenneth M Flurchick from the Department of physics and Computational Science and Engineering program at North Carolina A&T State University (NCAT) for his guidance and advice on developing the parallel version of IDS. The computational resources used in this research were made available in part by a grant from the National Science Foundation Grant No.#1429464, entitled MRI: Acquisition of CRAY XC-30 HPC Cluster to Support NCAT Research Computing, Education and Outreach. In Addition, the authors want to acknowledge the financial support from AFRL under the Grant/Cooperative Agreement Award #FA8651-17-2-0001.