Mathematical Modeling of Aerodynamic Heating and Pressure Distribution on a 5-Inch Hemispherical Concave Nose in Supersonic Flow

Abstract

During the design of an aircraft, a significant parameter that is taken into consideration is aerodynamic heating. Aerodynamically induced heating affects both the structure of the aircraft and its vulnerability to heat-seeking missiles in modern warfare. As a result, the ability to calculate efficiently the heating produced as well as the pressure distribution in such flows is crucial. Therefore, in this present study, the PHOENICS CFD code as modified by DRA Farnborough in order to calculate heat transfer and pressure measurement on a 5-inch hemispherical concave nose at a Mach number of 2.0 is evaluated and customized in order to produce faster and more accurate results. Apart from numerical alterations, different turbulence models are being examined as well as different discretization schemes. Numerical solutions show improvement up to 6% in comparison with the original model, both in terms of convergence rate and in terms of agreement with the available experimental data. With the new modeling suggested in the present work, the significance of both the discretization scheme and the choice of the turbulence modeling is demonstrated for the flows under consideration. The use of a high-order discretization scheme is suggested for more acute areas of the body modeled in order to improve results further.

1. Introduction

Aerodynamic heating has been a major study issue for a long time now and many investigations have been conducted by researchers in order to achieve high-quality simulation results, both on the theoretical and the numerical level. Experimental studies of the above phenomenon are rare, as conducting such research is prohibitively expensive.

Aerodynamic heating is a phenomenon worth studying as its significance lies in two major scientific fields. At first, aerodynamic heating is the most significant parameter of the emitted infrared radiation from aircraft, a fact that is of direct interest in military studies. According to statistics, the majority of aerial vehicles downed in modern warfare was due to weapons based on detecting Infrared signature level of those vehicles. Their percentage is as high as 89% of total aircraft shot downs between 1967 and 1993, according to Thompson J. et al. [1]. Apart from the emitted infrared radiation, the phenomenon of aerodynamic heating has to be taken into consideration when designing aerospace vehicles that fly at speeds comparable to that of sound. Ultrasonic aircraft and aircraft that re-enter the atmosphere can develop temperatures high enough that have to be carefully considered when designing the aircraft, in order to avoid unpleasant situations.

Aerodynamic heating is a phenomenon that takes place mostly due to heat from the warm part of the boundary layer to the body’s surface. However, when approaching the speed of sound, aerodynamic heating phenomena are enhanced by the radiation emitted from molecules that are being decomposed and ionized due to the temperature rising behind the shockwave. At speeds greater than 5000 m/sec, the emitted radiation becomes significant and plays a measurable role in the phenomenon. Finally, another way that aerodynamic heating takes place concerns aircraft that fly in the higher atmospheric layers, where the distance between the air molecules is comparable to the aerial vehicle dimensions. However, the present study is focused only on the convection phenomena that lead to aerodynamic heating, as the speeds examined do not exceed Mach numbers of 2.

As mentioned above, the significance of aerodynamic heating is a major issue in both the preliminary design of aerial vehicles and the calculation of the infrared radiation emitted from the vehicle in flight. Apart from experimental methods, which may offer reliable results but require enormous budgets, numerous simulations have been carried out via means of computational mechanics. Computational mechanics simulations approach the problem mainly through inviscid-viscous methods thus promising reduced computational time. AEROHEAT and INCHES are two of the simpler methods to approach the phenomenon of aerodynamic heating [2]. Both methods use the axisymmetric analog concept that allows axisymmetric boundary layer techniques to be applied to three-dimensional flows, provided that the surface streamlines are known. AEROHEAT calculates approximate surface streamlines based solely on the body geometry. INCHES uses an approximate expression for the scale factor in the windward and leeward planes, which describes the spreading of surface streamlines. An empirical relation then generates circumferential heating rates.

DRA Farnborough developed the SAPPHIRE software that is capable of calculating the surface temperatures and the characteristics of the flow field for a given geometry and defined flight characteristics. The above software is a series of codes that use PHOENICS [3] as the main CFD solver. The main advantage of CFD techniques is that aerodynamic heating can be related to other phenomena, as well, such as conjugating heat transfer and providing more detailed results. In order for SAPPHIRE to provide reliable results on the thermal trace, aerodynamic heating of the body has to be calculated with an error that does not exceed 5% [4]. For this purpose, a PHOENICS-based code was developed that has the possibility to validate results with available experimental data. The purpose of the present work is to make the model more time efficient and more accurate, producing results closer to the available experimental data.

According to finite volume or codes with FEM, attempts have been made, with acceptable results, by means of conjugating CFD and FEM. Murakami et al. [5] used CFD coding for the Navier-Stokes equations to calculate the outer flow, while the finite element method was used to perform the thermal analysis via equations describing unsteady heat conduction for a 3D field. The results were impressive, as they were very close to the experimental data measured in a wind tunnel experiment. Lu Jianwei and Wang Qiang [6] in order to calculate the emitted IR radiation suggested a division of the computational grid into two regions: the fluid and the solid. The fluid region is modeled by the Navier-Stokes equations, while turbulence is modeled by the RNG k-ε model. In the solid region, the energy transfer equation is used. In a similar study for calculating the emitted IR radiation, Marlene Johansson and Mats Dalenbring [7] developed the SIGGE software. SIGGE uses inputs of temperature, pressure, and species distribution that have been roughly estimated beforehand, together with the grid of the geometry being studied. The Navier-Stokes equations are solved by the VOLSOL solver, while turbulence is described by the k-ε model.

Apart from CFD methods, attempts have been made to simulate phenomena that involve aerodynamic heating via other methods as well. Rodrigo C. Palharini and Wilson F. N. Santos [8] used the Direct Simulation Monte Carlo (DSMC) method to simulate a two-dimensional steady supersonic flow in an orthogonal cavity. Main interests of this study were the characteristics of the aerodynamic surface, such as heat transfer coefficients, etc. Another method introduced by Georgia Tech Research Institute [9] named GTSIG (Georgia Tech Signature Inviscid Model) involves the calculation of temperature distribution by creating thermal meshes and partial differential equations. Finally, Christopher J. Riley et al., in their study [2], developed a 3D technique without viscous calculations that has the ability to calculate the degree of surface heating. Surface streamlines are calculated from the inviscid solution and the axisymmetric analog. Then, they are used in combination with an equation for convective heat transfer, to calculate the total surface heat transfer. The results of this method are satisfying as they are close to experimental data and the results of both the Navier-Stokes and the VSL equation solutions.

2. Mathematical modeling

2.1 Governing equations

In order to describe the phenomena under investigation, use is made of the conservation equations describing a three-dimensional, single phase, turbulent, steady state flow. The above equations are formulated in such a manner as to incorporate the physical modeling of the case study under the following assumptions:

• One phase, steady-state flow

• Compressible fluid

• Newtonian fluid

• Absence of chemical reaction

It is also clear that because of the general geometries considered, the equations must be modified in a way to be suitable for use in a body-fitted coordinate (BFC) system.

2.1.1 Modification of conservation equations from Cartesian coordinate system to BFC

Assuming that $\overrightarrow{\mathit{X}}=\left({\mathit{X}}_1,{\mathit{X}}_2,{\mathit{X}}_3\right)$ the Cartesian coordinate system, where X_j are the Cartesian Coordinates and $\overrightarrow{{\mathit{i}}_{\mathit{j}}}$ the direction unit vectors, the velocity $\overrightarrow{\mathit{V}}$ can be expressed as:

$$\overrightarrow{\mathit{V}}=\sum _{\mathit{j}=1}^3{\mathit{u}}_{\mathit{j}}\overrightarrow{{\mathit{i}}_{\mathit{j}}}={\mathit{u}}_1\overrightarrow{\mathit{i}}\mathit{+}{\mathit{u}}_2\overrightarrow{\mathit{j}}\mathit{+}{\mathit{u}}_3\overrightarrow{\mathit{k}\mathit{.}}$$

Furthermore,

$$\overrightarrow{\mathrm{\nabla }}=\sum _{\mathit{i}=1}^3\overrightarrow{{\mathit{i}}_{\mathit{j}}}\frac{\mathit{\partial }}{{\mathit{\partial x}}_{\mathit{j}}}$$ E1

Assume a non-orthogonal coordinate system, defined as ${\mathit{\xi }}_{\mathit{\eta }}={\mathit{\xi }}_{\mathit{\eta }}\left({\mathit{X}}_1,{\mathit{X}}_2,{\mathit{X}}_3\right)$ , η = 1,2,3, where ξ_η are the curvilinear coordinates at any point and X₁, X₂, X₃ its Cartesian coordinates. Assume $\overrightarrow{{\mathit{g}}_{\mathit{\eta }}}$ to be the non-orthogonal unit vectors along the directions of the curvilinear coordinates; then the components of the velocity vector tangent to the unit vectors $\overrightarrow{{\mathit{g}}_{\mathit{\eta }}}$ can be defined as:

$${\mathit{V}}_{\mathit{\eta }}=\overrightarrow{{\mathit{g}}_{\mathit{\eta }}}\overrightarrow{\mathit{V}}=\sum _{\mathit{j}=1}^3{\mathit{g}}_{\mathit{\eta },\mathit{j}}{\mathit{U}}_{\mathit{j}},\mathrm{\eta }=\mathrm{1,2,3}$$ E2

where U_j are the Cartesian components and the unit tangent vectors alongside the three directions are:

$$\overrightarrow{{\mathit{g}}_{\mathit{\eta }}}={\overrightarrow{\mathit{X}}}_{,\xi \eta }/\left|{\overrightarrow{\mathit{X}}}_{\mathit{\xi \eta }}\right|=\sum _{\mathit{j}=1}^3{\overrightarrow{\mathit{g}}}_{\mathit{\eta }\mathrm{,}\mathit{j}}{\overrightarrow{\mathit{i}}}_{\mathit{j}},\mathrm{\eta }=\mathrm{1,2,3}$$ E3

where $\overrightarrow{\mathit{X}}$ is the position vector of a point and the index ξ_η denotes the partial derivative to the curvilinear coordinate ξ_η.

2.1.2 Conservation of momentum

In order to obtain the resolute [10] of the momentum equation in each of the coordinate directions, we form the dot product of tangent unit vectors with the momentum equation (Eq. (1)). The momentum equation, considering the assumptions and the physics of this simulation, is derived as follows:

$$\mathit{\rho }\left[\left(\overrightarrow{\mathrm{\nabla }}\overrightarrow{\mathit{V}}\right)\overrightarrow{\mathit{V}}\right]-\overrightarrow{\mathrm{\nabla }}·\stackrel{\sim }{\mathit{\tau }}=-\overrightarrow{\mathrm{\nabla }}\mathit{p}+\overrightarrow{\mathit{F}},$$ E4

where

$$\overrightarrow{\mathrm{\nabla }}·\stackrel{\sim }{\mathit{\tau }}=\left(\overrightarrow{\mathrm{\nabla }}\mathit{\mu }\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{\mathit{V}}+\sum _{\mathit{\lambda }=1}^3\left[\sum _{\mathit{j}=1}^3\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\left(\mathit{\mu }\frac{\mathit{\partial }{\mathit{u}}_{\mathit{j}}}{\mathit{\partial }{\mathit{x}}_{\mathit{\lambda }}}\right)\right]\overrightarrow{{\mathit{i}}_{\mathit{\lambda }}}$$ E5

In order to obtain the resolute of the momentum equation in each of the coordinate directions, the dot product of the tangent unit vector and the above momentum equation is formed:

$$\begin{array}{c}\overrightarrow{g_n}\rho \left[\left(\overrightarrow{\nabla }\overrightarrow{V}\right)\overrightarrow{V}\right]-\overrightarrow{g_n}\overrightarrow{\nabla }·\tilde{\tau }=-\overrightarrow{g_n}\overrightarrow{\nabla }p\leftrightarrow \overrightarrow{g_n}\left(\rho \overrightarrow{V}·\overrightarrow{\nabla }\right)\overrightarrow{V}-\overrightarrow{g_n}·\left(\overrightarrow{\nabla }·\mu \overrightarrow{\nabla }\right)\overrightarrow{V}-{\overrightarrow{g}}_n·S_{\mu }\\ =-\overrightarrow{g_n}\overrightarrow{\nabla }p\end{array}$$ E6

where ${\mathit{S}}_{\mathit{\mu }}={\sum }_{\mathit{\lambda }=1}^3\left[{\sum }_{\mathit{j}=1}^3\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\left(\mathit{\mu }\frac{\mathit{\partial }{\mathit{u}}_{\mathit{j}}}{\mathit{\partial }{\mathit{x}}_{\mathit{\lambda }}}\right)\right]\overrightarrow{{\mathit{i}}_{\mathit{\lambda }}}$

From vector calculus, using the identity $\mathit{A}\mathrm{x}\left(\mathit{B}\mathrm{x}\mathit{C}\right)=\mathit{B}\left(\mathit{A}\mathit{·}\mathit{C}\right)\mathit{-}\mathit{C}\left(\mathit{A}\mathit{·}\mathit{B}\right)$ , we have:

$$\overrightarrow{\mathrm{\nabla }}\left({\overrightarrow{\mathit{g}}}_{\mathit{n}}·\mathit{\rho }\overrightarrow{\mathit{V}}\right)=\left(\mathit{\rho }\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}\right){\overrightarrow{\mathit{g}}}_{\mathit{n}}+\left({\overrightarrow{\mathit{g}}}_{\mathit{n}}·\overrightarrow{\mathrm{\nabla }}\right)\mathit{\rho }\overrightarrow{\mathit{V}}+\mathit{\rho }\overrightarrow{\mathit{V}}\mathrm{x}\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}{\overrightarrow{\mathit{g}}}_{\mathit{n}}\right)+{\overrightarrow{\mathit{g}}}_{\mathit{n}}\mathrm{x}\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}\mathit{\rho }\overrightarrow{\mathit{V}}\right)$$ E7

Also,

$$\overrightarrow{\mathit{V}}·{\overrightarrow{\mathit{g}}}_{\mathit{n}}\mathrm{x}\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}\overrightarrow{\mathit{V}}\right)=-{\overrightarrow{\mathit{g}}}_{\mathit{n}}·\left(\overrightarrow{\mathit{V}}\mathrm{x}\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}\overrightarrow{\mathit{V}}\right)\right)\overrightarrow{\mathit{V}}$$ E8

$$\left(\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{\mathit{V}}=\overrightarrow{\mathrm{\nabla }}\left(\frac{{\left|{\mathit{V}}_{\mathit{n}}\right|}^2}{2}\right)-\overrightarrow{\mathit{V}}\mathrm{x}\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}\overrightarrow{\mathit{V}}\right)$$ E9

$$\overrightarrow{\mathit{V}}·\overrightarrow{\mathit{V}}\mathrm{x}\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}{\overrightarrow{\mathit{g}}}_{\mathit{n}}\right)=\left(\overrightarrow{\mathit{V}}\mathrm{x}\overrightarrow{\mathit{V}}\right)·\left(\overrightarrow{\mathrm{\nabla }}\mathrm{x}{\overrightarrow{\mathit{g}}}_{\mathit{n}}\right)=0,\mathit{\text{since}}\overrightarrow{\mathit{V}}\mathrm{x}\overrightarrow{\mathit{V}}=0$$ E10

Thus,

$${\overrightarrow{\mathit{g}}}_{\mathit{n}}·\overrightarrow{\mathrm{\nabla }}\left(\frac{{\left|{\mathit{V}}_{\mathit{n}}\right|}^2}{2}\right)=\overrightarrow{\mathit{V}}·\left({\overrightarrow{\mathit{g}}}_{\mathit{n}}·\overrightarrow{\mathrm{\nabla }}\right)·\overrightarrow{\mathit{V}}$$

Forming the dot product of the above equation and (Eq. (7)), we have:

$$\left(\mathit{\rho }\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}\right){\overrightarrow{\mathit{V}}}_{\mathit{n}}=\overrightarrow{\mathit{V}}·\left(\mathit{\rho }\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}\right){\overrightarrow{\mathit{g}}}_{\mathit{n}}+{\overrightarrow{\mathit{g}}}_{\mathit{n}}·\left(\mathit{\rho }\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{\mathit{V}}$$ E11

Also, for a scalar variable φ_n:

$$\overrightarrow{\mathrm{\nabla }}·\left(\mathit{\rho }\overrightarrow{\mathit{V}}{\mathit{\phi }}_{\mathit{n}}\right)=\mathit{\rho }\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}{\mathit{\phi }}_{\mathit{n}}$$ E12

and

$$\left(\overrightarrow{\mathrm{\nabla }}·\mathit{\mu }\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{\mathit{V}}·{\overrightarrow{\mathit{g}}}_{\mathit{n}}={\overrightarrow{\mathit{g}}}_{\mathit{n}}·\left(\overrightarrow{\mathrm{\nabla }}·\mathit{\mu }\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{\mathit{V}}+\overrightarrow{\mathit{V}}·\left(\overrightarrow{\mathrm{\nabla }}·\mathit{\mu }\overrightarrow{\mathrm{\nabla }}\right){\overrightarrow{\mathit{g}}}_{\mathit{n}}$$ E13

Using Eqs. (11)–(13), we can write (Eq. (6)), as follows:

$$\overrightarrow{\mathrm{\nabla }}·\left(\mathit{\rho }\overrightarrow{\mathit{V}}{\mathit{V}}_{\mathit{n}}\right)-\left(\overrightarrow{\mathrm{\nabla }}·\mathit{\mu }\overrightarrow{\mathrm{\nabla }}\right){\mathit{V}}_{\mathit{n}}-\overrightarrow{{\mathit{g}}_{\mathit{n}}}·\overrightarrow{{\mathit{S}}_{\mathit{\mu }}}=-{\mathit{g}}_{\mathit{n}}\overrightarrow{\mathrm{\nabla }}\mathit{p}+\overrightarrow{\mathit{V}}\left(\mathit{\rho }\overrightarrow{\mathit{V}}·\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{{\mathit{g}}_{\mathit{n}}}-\overrightarrow{\mathit{V}}\left(\overrightarrow{\mathrm{\nabla }·}\mathit{\mu }\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{{\mathit{g}}_{\mathit{n}}}$$ E14

and

$$\left(\overrightarrow{\mathit{V}}\mathrm{·}\overrightarrow{\mathrm{\nabla }}\right){\mathit{V}}_{\mathit{n}}=\overrightarrow{\mathit{V}}\left(\overrightarrow{\mathit{V}}\mathrm{·}\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{{\mathit{g}}_{\mathit{n}}}\mathrm{+}\overrightarrow{{\mathit{g}}_{\mathit{n}}}\left(\overrightarrow{\mathit{V}}\mathrm{·}\overrightarrow{\mathrm{\nabla }}\right)\overrightarrow{\mathit{V}}$$

From Eqs. (1) and (3), the partial derivative in the direction ξ_n is:

$$\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{\xi }}_{\mathit{n}}}=\overrightarrow{{\mathit{g}}_{\mathit{n}}}·\overrightarrow{\mathrm{\nabla }}=\sum _{\mathit{j}=1}^3{\mathit{g}}_{\mathit{n},\mathit{j}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}$$ E15

With inverse

$$\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{n}}}=\sum _{\mathit{j}=1}^3{\mathit{G}}_{\mathit{n},\mathit{j}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{\xi }}_{\mathit{j}}}$$ E16

where G_n,j the elements of the inverse matrix of the elements g_n,j.

2.1.3 Conservation of mass

Since the case study is for steady flow, the conservation of mass equation is hence simplified in the following form for the Cartesian coordinate system:

$$\mathit{div}\left(\mathit{\rho V}\right)=\frac{\mathit{\partial }\left(\mathit{\rho }\overrightarrow{\mathit{V}}\right)}{\mathit{\partial x}}+\frac{\mathit{\partial }\left(\mathit{\rho }\overrightarrow{\mathit{V}}\right)}{\mathit{\partial y}}+\frac{\mathit{\partial }\left(\mathit{\rho }\overrightarrow{\mathit{V}}\right)}{\mathit{\partial z}}=0$$

Using body-fitted coordinates, the continuity equation must be solved in terms of the velocity components along the coordinate directions. These components, u_n, would be $\mathit{V}={\sum }_{\mathit{n}=1}^3{\overrightarrow{\mathit{g}}}_{\mathit{n}}{\mathit{u}}_{\mathit{n}}$ and are calculated from the velocity resolutes by solving the simultaneous equations

$${\mathit{V}}_{\mathit{n}}=\sum _{\mathit{j}=1}^3{\overrightarrow{\mathit{g}}}_{\mathit{n}}·\left({\overrightarrow{\mathit{g}}}_{\mathit{j}}{\mathit{u}}_{\mathit{j}}\right),\mathit{n}=1,2,3.$$ E17

Velocity components in terms of the resolutes are

$${\mathit{u}}_{\mathit{n}}=\sum _{\mathit{j}=1}^3{\mathit{G}}_{\mathit{n},\mathit{j}}{\mathit{V}}_{\mathit{j}},\mathrm{n}=1,2,3$$ E18

where

$$\raisebox{1ex}{${\mathit{G}}_{\mathit{n},\mathit{j}}$}\!\left/ \!\raisebox{-1ex}{$\mathit{\sigma }$}\right.=\left[\begin{array}{ccc}\left(1-{\mathit{c}}^2\right)& \left(\mathit{bc}-\mathit{a}\right)& \left(\mathit{ac}-\mathit{b}\right)\\ \left(\mathit{bc}-\mathit{a}\right)& \left(1-{\mathit{b}}^2\right)& \left(\mathit{ab}-\mathit{c}\right)\\ \left(\mathit{ac}-\mathit{b}\right)& \left(\mathit{ab}-\mathit{c}\right)& \left(1-{\mathit{a}}^2\right)\end{array}\right]$$ E19

where

$$\mathit{a}={\mathit{g}}_1·{\mathit{g}}_2,\mathrm{b}={\mathit{g}}_1·{\mathit{g}}_3,\mathrm{c}={\mathit{g}}_2·{\mathit{g}}_3$$ E20

and

$$\mathit{\sigma }=\left[1-\left({\mathit{a}}^2+{\mathit{b}}^2+{\mathit{c}}^2-2\mathit{abc}\right)\right]$$ E21

Using Eqs. (1) and (18), the following expression for the velocity is derived:

$$\overrightarrow{\mathit{V}}=\sum _{\mathit{n}=1}^3{\overrightarrow{\mathit{i}}}_{\mathit{n}}\sum _{\mathit{j}=1}^3{\mathit{G}}_{\mathit{n},\mathit{j}}{\mathit{V}}_{\mathit{j}}=\sum _{\mathit{j}=1}^3{\overrightarrow{\mathit{Q}}}_{\mathit{j}}{\overrightarrow{\mathit{V}}}_{\mathit{j}}$$ E22

where

$$\overrightarrow{{\mathit{Q}}_{\mathit{j}}}=\sum _{\mathit{j}=1}^3{\mathit{G}}_{\mathit{n}\mathrm{,}\mathit{j}}{\overrightarrow{\mathit{i}}}_{\mathit{n}}$$ E23

the vector used to obtain the resolutes (projections) of the velocity vector $\overrightarrow{\mathit{V}}$ .

From Eqs. (1), (16), and (20):

$$\overrightarrow{\mathrm{\nabla }}=\sum _{\mathit{j}=1}^3{\overrightarrow{\mathit{Q}}}_{\mathit{j}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{\xi }}_{\mathit{j}}}$$ E24

Finally, the continuity equation can be written as $\overrightarrow{\mathrm{\nabla }}\left(\mathit{\rho }\overrightarrow{\mathit{V}}\right)=0$ where

$$\overrightarrow{\mathit{V}}=\sum _{\mathit{n}=1}^3{\overrightarrow{\mathit{g}}}_{\mathit{n}}{\mathit{u}}_{\mathit{n}}$$ E25

2.2 Turbulence modeling

In order to evaluate the use of different turbulence models in the cases under consideration, a parametric study was accomplished by trying three “popular” turbulence models: (i) standard k-ε model, (ii) RNG k-ε model, and (iii) k-ω model. The k-ε model simulation was chosen in order to reproduce the original simulation accomplished by Steve Rooks [4]. In order to take it one step further and include the effect of the swirling of the flow, another simulation was performed using the RNG k-ω model. Last but not least, the k-ω turbulence model was used as well, as preliminary studies proved its suitability for wall-bounded flows.

2.2.1 The general ϕ-equation

By generalizing Eq. (14), the appropriate equation for any scalar variable can be provided [11]:

$$\overrightarrow{\mathrm{\nabla }}·\overrightarrow{\mathit{Q}}=\mathit{S}{\mathit{\phi }}_{\mathit{n}}$$ E26

where

$$\overrightarrow{\mathit{Q}}=\mathit{\rho }\overrightarrow{\mathit{V}}{\mathit{\phi }}_{\mathit{n}}-\overrightarrow{{\mathit{\Gamma }}_{{\mathit{\phi }}_{\mathit{n}}}}\overrightarrow{\mathrm{\nabla }}{\mathit{\phi }}_{\mathit{n}}$$ E27

and φ a general scalar variable like H (enthalpy), k, ε, etc. Regarding k (turbulence kinetic energy) and ε (turbulence dissipation rate):

$${\mathit{S}}_{\mathit{k}}=\mathit{G}-\mathit{\rho \epsilon }$$ E28

$${\mathit{S}}_{\mathit{\epsilon }}=\raisebox{1ex}{${\mathit{C}}_1\mathit{\epsilon }$}\!\left/ \!\raisebox{-1ex}{$\mathit{kG}$}\right.-\raisebox{1ex}{${\mathit{C}}_2\mathit{\rho }{\mathit{\epsilon }}^2$}\!\left/ \!\raisebox{-1ex}{$\mathit{\kappa }$}\right.$$ E29

$$\mathit{G}=\mathit{\mu }\left(\raisebox{1ex}{${\mathit{\partial u}}_{\mathit{j}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{\partial x}}_{\mathit{i}}$}\right.+\raisebox{1ex}{${\mathit{\partial u}}_{\mathit{i}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{\partial x}}_{\mathit{j}}$}\right.\right)\raisebox{1ex}{${\mathit{\partial u}}_{\mathit{j}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{\partial x}}_{\mathit{j}}$}\right.,\mathit{\text{in Cartesian}}-\mathit{\text{tensor form}}$$ E30

where u_i, u_j the velocity components in the direction of the Cartesian coordinates x_i, x_j with i = 1,2,3 and j = 1,2,3 [10]. It must also be noted that the effective exchange coefficients for k, ε are as follows:

$${\mathit{\mu }}_{\mathit{l}}+{\mathit{\mu }}_{\mathit{t}}=\mathit{\mu }$$ E31

$${\mathit{\mu }}_{\mathit{t}}=\mathit{\rho }{\mathit{C}}_{\mathit{\mu }}\frac{{\mathit{k}}^2}{\mathit{\epsilon }}$$ E32

$${\mathit{\Gamma }}_{\mathit{k}}=\raisebox{1ex}{$\left({\mathit{\mu }}_{\mathit{l}}+{\mathit{\mu }}_{\mathit{t}}\right)$}\!\left/ \!\raisebox{-1ex}{${\mathit{\sigma }}_{\mathit{k},\mathit{t}}$}\right.$$ E33

$${\mathit{\Gamma }}_{\mathit{\epsilon }}=\raisebox{1ex}{$\left({\mathit{\mu }}_{\mathit{l}}+{\mathit{\mu }}_{\mathit{t}}\right)$}\!\left/ \!\raisebox{-1ex}{${\mathit{\sigma }}_{\mathit{\epsilon },\mathit{t}}$}\right.$$ E34

where μ_l is the laminar viscosity; C₁, C₂, and C_μ are turbulence model constants; and σ_k,t and σ_ε,t are the turbulence Prandtl numbers for k, ε [10].

In total seven differential equations are solved, namely, the continuity equation, the three momentum equations, the enthalpy equation, and the two equations for the turbulence variables.

2.2.2 The near-wall region

Special attention was given to the treatment of the near-wall region, as it is of high importance for the type of flows studied and the 2-equation models used to simulate turbulence do not, in general, produce results of high accuracy in near-wall regions. The near-wall region was treated via a log-law wall function that also played a major role in deciding the computational grid as well as ensuring computational economy since there is no need to solve throughout the viscous sub-layer. The use of this method was introduced by Launder and Spalding [12], and the formula used is:

$$\frac{\mathit{ln}\left(\mathit{E}{\mathit{Y}}^{+}\right)}{\mathit{k}}=\frac{\mathit{U}}{{\mathit{U}}_{\mathit{\tau }}}$$ E35

where

U: the absolute velocity parallel to the wall in the first node of the grid

U_τ: the friction velocity

k: the Von Karman constant

E: the surface’s roughness parameter

Y⁺: the dimensionless distance from the wall

The above function is applied to nodes whose dimensionless distance from the wall is between 30 and 130. In this region of the boundary layer, the effects of turbulence and viscosity are equally important.

2.3 Boundary conditions

Because of the selection of body-fitted coordinates to design the grid of the problem, the velocity must be broken down to its components and transferred from the Cartesian coordinate system to the BFC used in the present study.

At entrance, a uniform flow is applied of Mach 2 speed. The rate of flow is consistently calculated for each control volume. The kinetic energy of turbulence and its rate of dissipation are considered constant with values of 21.78 m²/s² (10% of the mean velocity squared) and 2.765*10⁵ m²/s³ (assuming turbulence viscosity 1000 times the laminar value), respectively ( Figure 1 ).

Figure 1.

The plane of inlet boundary conditions specification.

The outlet boundary conditions describe the fact that there is no alteration in the rate of mass flow, as there is neither mass creation nor destruction. As for the turbulence kinetic energy, turbulence dissipation rate, and the enthalpy, Neumann conditions prevail, suggesting there is no alteration of values regarding these properties along the boundary ( Figure 2 ). Since the boundary conditions must be transported in the BFC system, the Jacobian of the transformation used is:

Figure 2.

The plane of outlet boundary conditions specification.

$$J=\frac{\partial x}{\partial g_1}\frac{\partial y}{\partial g_2}-\frac{\partial x}{\partial g_2}\frac{\partial y}{\partial g_1}$$

and the general form of the functions:

$$\frac{\partial f}{\partial x}=\frac{1}{J}\left(\frac{\partial f}{\partial g_1}\frac{\partial y}{\partial g_2}-\frac{\partial f}{\partial g_2}\frac{\partial y}{\partial g_1}\right)$$

$$\frac{\partial f}{\partial y}=\frac{1}{J}\left(\frac{\partial f}{\partial g_2}\frac{\partial x}{\partial g_1}-\frac{\partial f}{\partial g_1}\frac{\partial x}{\partial g_2}\right)$$

Finally, the boundary condition at the northern boundary is prescribed external pressure. Like the boundary conditions at the outlet, the kinetic energy of turbulence, its dissipation rate, and the enthalpy are described by Neumann boundary conditions, implying absence of factors capable of altering their values along the northern boundary.

2.4 Computational grid and method used

2.4.1 Computational grid generation

The study for grid independence came down to the same results as the one proposed by Steve Rooks [4]. The computational grid used is one with dimensions 48*45*1 and has been proved to lead to grid-independent results [4]. Briefly, the grid divides the computational field in three regions along the x-axis and three along the y-axis, as one can see in the following figure ( Figure 3 ).

Figure 3.

The computational grid.

Along the y-axis the three sub-regions correspond to the area below the lip, the area of the lip, and finally the area above the lip. A most important role in the dimensions of the grid played the fact that a denser or a sparser grid would result in improper use of the wall functions, as it would affect the dimensionless distance Y+.

2.4.2 Computational methods used

In order to solve the equations, the finite-volume method is employed. As for regular grids, the finite-volume equations can be derived by integrating over a control volume in a system with body-fitted coordinates. The use is made of the Gauss theorem that transforms a volume integral to a surface one:

$${\int }_{\mathit{\Delta V}}\overrightarrow{\nabla }·\overrightarrow{Q}\mathit{dV}={\int }_S\overrightarrow{n}·\overrightarrow{Q}\mathit{dS}={\int }_S\overrightarrow{n}·\left(\rho \overrightarrow{V}{\phi }_{\eta }\right)\mathit{dS}-{\int }_S\overrightarrow{n}·\left(\overrightarrow{{\Gamma }_{{\phi }_{\eta }}}\overrightarrow{\nabla }{\phi }_{\eta }\right)\mathit{dS}$$ E36

where $\overrightarrow{n}$ is the outward vector drawn normal to any of the six surfaces defining ΔV and S denotes the enclosing the volume surface area [10]. The first term represents convection, where properties are regarded as uniform over the cell faces, except in respect of the velocities for which the face-center values are stored. The values prevailing at the cell face are determined by using upwind interpolation. As for the second term, the one representing diffusion the property gradients and the transport properties that they multiply are uniform over cell faces. In order to calculate the gradients, it is supposed that properties vary linearly and also the transport properties are arithmetic averages of those on either side of the cell faces. In source terms, the nodal values are supposed to prevail over the whole of the cell volume.

Any dependent variable’s partial differential equation is represented by a coupled set of algebraic equations of the following form:

$$A_P{\phi }_P=\sum _i{A}_i{\phi }_i+S_{\phi }$$ E37

where the A_i’s represent the influence of convection and diffusion, S_φ is the source term, and P refers to the control volume under consideration. Σ_i indicates the summation over the neighboring nodes. All the above equations are derived through work that cannot be presented here due to limited space [10].

2.4.3 Numerical simulation

The numerical simulations in this work are realized by the fluid dynamics package PHOENICS. The SIMPLEST algorithm is used in the iterative procedure, whereas the equations are solved via the finite-volume method. The momentum equations are solved using the initial estimate of the pressure field, using velocities that do satisfy momentum but not in the general continuity. For each cell, the equation of (inlet-outlet) is formed, and an equation of pressure correction is solved where the coefficients are A (d(vel))/dp and the sources are the errors in continuity. Afterwards, the pressure field is corrected alongside with the velocity field. The iterative procedure continues until the errors in the equations of momentum and continuity are acceptable.

2.5 Relaxation methods used

In order to improve convergence and eliminate great fluctuations in values of variables in consecutive iterations that could lead to the divergence of the solution, relaxation techniques are used. In the present study, two techniques are used: the linear relaxation and the false time-step relaxation. Linear relaxation is used for the variables of pressure, enthalpy, turbulence kinetic energy, and turbulence dissipation rate as it is more suitable for scalar variables in the following form: φ^new = φ^old + α(φ^new − φ^old). In the case of the vector variables of U1 and V1, false time-step relaxation is used as it is more suitable for the velocity variables in the following form ρV_p/δt(Φ* − Φ_ρ) = β(Φ* − Φ_ρ) where V_P is the cell volume, Φ* the previous iteration value of variable Φ, Φ_ρ the value of φ in the last iteration, and δt_false the false time step, generally of the order of the cell residence time.

2.6 Discretization schemes used

In the original coding, the Hybrid discretization scheme was used. However, in order to investigate further the most efficient way to simulate the phenomenon, runs were also carried out using the upwind discretization scheme that, according to literature, is more efficient when describing directional flows but introduces more “false-diffusion” errors.

2.7 Modifications in the original coding

In order to improve the simulation results both in terms of convergence speed and agreement with the experimental data, numerous modifications were made to the original coding. These modifications include the calculation of the diffusion coefficient by means of the harmonic mean and not the numerical mean, and the calculation of reference kinematic viscosity for laminar flow was accomplished via Sutherland’s law. Moreover, the initial values of the variables for enthalpy and velocity were not chosen randomly, as in the original coding but were calculated based on the simulation’s data so as to help the convergence procedure with more “realistic” initial variable values. At the northern boundary, a boundary condition that suggests constant mass flow was used so as to ignore the effects of diffusion. Finally, in opposition to the initial coding, where false time-step relaxation was used for the variables of enthalpy, turbulence kinetic energy, and turbulence dissipation rate, in the modified version, for these variables, linear under-relaxation was used.

2.8 Different case studies

In terms of this present research, many simulations were carried out, and results of the most relevant eight of them were brought into comparison in order to decide which method is more efficient in simulating the phenomenon studied by Markley et al. [13]. The first simulation carried out in terms of this research was the reproduction of the simulation carried out by Steve Rooks [4] in an up-to-date version of the PHOENICS. The other seven include the modified original version, described in the above paragraph, although no alteration in either the turbulence model or the discretization scheme was made. The other six were the combinations of the three turbulence models (k-ε, RNG k-ε, and k-ω) with the two discretization schemes (Hybrid and Upwind) studied in this research paper.

3. Results

In order to examine the validity of the results of the simulation, a comparison had to be made with the available experimental data. The available experimental data are for temperature, as that was the main point of interest at the experiment conducted by NACA [13], measured behind the lip, and high enough so as the flow would be unaffected by boundary layer phenomena. As a result, measurements were taken along the line Y = 18 and behind the region dictated by cell X = 20.

The first simulation carried out in terms of this research was the reproduction of the simulation carried out in [4] using our modified version of the PHOENICS code, so as comparison can be carried out. The modifications made to the original coding that led to improved results, both in terms of agreement with the experimental data and in terms of decreased convergence time, include the calculation of the diffusion coefficient by means of the harmonic mean rather than the numerical mean and the calculation of reference kinematic viscosity for laminar flow that was accomplished via Sutherland’s law. Moreover, the initial values of the variables for enthalpy and velocity were not chosen randomly, as in the original coding, but were calculated based on the simulation’s data so as to help the convergence procedure with more “realistic” initial variable values. At the northern boundary, a boundary condition that suggests constant mass flow was used so as to ignore the effects of diffusion. Finally, different to the initial coding, where false time-step relaxation was used for the variables of enthalpy, turbulence kinetic energy, and turbulence dissipation rate, in the modified version, for these variables, linear relaxation was used.

For this initial modified version, no alterations were made as far as the turbulence modeling and the discretization scheme are concerned. The initial modeling with k-ε turbulence modeling and the Hybrid discretization scheme were used. As a result of the above changes, the convergence of the coding was faster by 48.33% as convergence was reached at 5427 iterations instead of 10,503 for the same level of errors remaining, and there were even slight improvements in the agreement with the experimental results as well.

The following diagram depicts the distribution of temperature along the x-axis above the airfoil at the height of the line Y = 18, as mentioned above ( Figure 4 ).

Figure 4.

Comparison diagram for temperature distribution along the x-axis.

As seen in the above diagram, the two simulations differ only slightly, and it can be noted that the modified version offers a marginal 0.04% better accuracy along with the faster convergence. More precisely, at the point where the temperature drop is noticed, which is the lip region, the accuracy is improved by 0.015%, whereas after the threshold of 5.3 cm, which is behind the lip region, the accuracy is improved by 0.27%.

Likewise, the diagrams depicting the main velocity and density distribution along the x-axis slightly differ from the diagrams of the initial simulation since the mean fault is 0.5847 and 0.9%, respectively ( Figures 5 and 6 ).

Figure 5.

Comparison diagram for velocity distribution along the x-axis.

Figure 6.

Comparison diagram for density distribution along the x-axis.

Regarding the distributions of temperature, velocity and density along the Y-axis, defined by the line X = 20 ( Figures 7 – 9 ), no major difference is noticeable in their main part, as one can notice from Figures 7 – 9 . The mean errors are 0.7, 3.06, and 0.9%, respectively. Yet, reaching the northern boundary, the results from the modified simulation differ from the initial one, as they do not alter the value they have settled in following the exit of the boundary layer. In the initial simulation, an increase can be spotted for the variables of temperature and density and a decrease in the value of velocity. The increase in the temperature value along with the decrease in velocity can be explained from the fact that enthalpy is assumed constant. However, in the modified simulation, the steady mass flow boundary condition neglects diffusion phenomena, and therefore there is stability in this variable. This approaches the nature of the original experiment, as no factors able to alter the values of these variables exist near the northern boundary.

Figure 7.

Comparison diagram for temperature distribution along the y-axis.

Figure 8.

Velocity distribution along Y-axis.

Figure 9.

Density distribution along Y-axis.

Figures 10 – 13 present the contour plots for temperature, density, and velocity distributions as well as the vector diagram of the velocity vectors.

Figure 10.

Contour plot of the temperature distribution on the computational field for the modified simulation (k-ε hybrid).

Figure 11.

Contour plot of the velocity distribution on the computational field for the modified simulation (k-ε hybrid).

Figure 12.

Contour plot with the density distribution on the computational field for the modified simulation (k-ε hybrid).

Figure 13.

Velocity vectors zoomed in the area of interest depicting the recirculations.

In order to investigate further the problem and any possible aspects that may improve the quality of the results produced, a parametric study was conducted. The parameters studied were both the turbulence model and the discretization scheme. Simulations were conducted using all three different turbulence models: k-ε, RNG k-ε, and k-ω. Each of the above turbulence models was used in combination with either the Hybrid discretization scheme or the Upwind.

In total, six different simulations were performed. The results were compared to the available experimental data in order to evaluate their accuracy. The schematic of the comparison can be seen in the diagram that follows ( Figure 14 ).

Figure 14.

Comparison diagram of parametric simulations and experimental data.

The simulation results can also be seen in numbers in the following table, where the mean error is calculated ( Table 1 ).

	k-ε Hybrid	k-ε Upwind	RNG k-ε Hybrid	RNG k-ε Upwind	k-ω Hybrid	k-ω Upwind
Overall mean error (%)	4.9529	3.5055	4.9729	3.5741	5.034	3.3418

Table 1.

Overall mean error calculation.

As one can see, the most successful simulation in terms of the overall mean error is the one engaging the k-ω turbulence model in combination with the Upwind discretization scheme. If we take a closer look, it is obvious that it is the only simulation that seems to approach the steep temperature descent in the 4–5 cm region. In that region, all simulations fail to keep the mean error below 5%; however, as one can see in the following table, the k-ω Upwind simulation achieves a mean error of 5.6274% ( Table 2 ).

	k-ε Hybrid	k-ε Upwind	RNG k-ε Hybrid	RNG k-ε Upwind	k-ω Hybrid	k-ω Upwind
Mean error in region 4–5 cm (%)	8.42	6.763	8.4331	6.8564	8.0097	5.6274

Table 2.

Mean error calculation in the 4–5 cm region.

The difficulty in simulating realistically the behind the lip region is possibly due to heat conduction phenomena. In the simulations performed, no solution is run for conjugate heat transfer, and given that the airfoil geometry is very thin, heat conduction taking place on the edge of the lip is speculated to be one among the causes that lead to higher calculated temperatures than the experimental data.

When examining the region behind the lip, where the flow presents a more stabilized form than the one when encountering the lip, the turbulence model that achieves better accuracy is the RNG k-ε model with the Upwind discretization scheme. This can be seen in both the diagram and the table that follow ( Figure 15 and Table 3 ).

Figure 15.

Comparison diagram between simulation performed and experimental data in the behind the lip region.

	k-ε Hybrid	k-ε Upwind	RNG k-ε Hybrid	RNG k-ε Upwind	k-ω Hybrid	k-ω Upwind
Mean error in behind the lip region (%)	6.5465	0.3348	6.6621	0.6050	8.0915	1.9930

Table 3.

Mean error calculation in behind the lip region.

It must be noted that in all simulations performed, the discretization scheme that produced the best results was the Upwind scheme. That would be justified by the fact that the flow is basically unidirectional and is characterized by high velocity. As a result, the diffusion phenomena are of very low intensity, with the convective phenomena being the main terms in the discretization equation. In both k-ω and k-ε simulations ran, the benefits from altering the initial code were still valid.

4. Recommendations

The calculation of aerodynamic heating phenomena includes many aspects that can interfere with both the speed of convergence but the accuracy of the results as well. Apart from the choice of proper initial value and other numeric factors, attention has to be paid to the flow field characteristics in order to decide the suitability of the turbulence model and the discretization scheme that will be used in order to achieve better results. Despite the fact that improvements for the 4–5 cm region were also suggested here, future research could focus on simulating the lip region, as this study’s main interest lied in the behind the lip region. In order to achieve that, high-order discretization schemes could be used, like the Van Leer scheme as well as high-order turbulence models, like the large eddy simulation (LES) and the direct numerical simulation (DNS) or any combinations.

5. Conclusion

Heat transfer and pressure measurement on a 5-inch hemispherical concave nose at a Mach number of 2.0 is investigated via the use of PHOENICS to evaluate aerodynamic heating phenomena. In continuation of [4], the code is modified in order to achieve better results both in terms of faster convergence but also in terms of more accurate results as well. Finite-volume method is used, while turbulence models and discretization schemes are studied as parameters in order to optimize results. Before altering the above mentioned parameters, attempts were made to improve the speed of convergence. Altering the values of input variables in order to be closer to the problem considered as well as choosing relaxation methods suitable for each variable type, scalar or vector, along with a few other changes described above, managed to improve the convergence speed by up 48.33%. Moreover, an improvement in boundary conditions of the northern boundary so as to ignore dissipation phenomena corrected the results of all variables at the northern region of the computational field, as their values were stabilized, as expected. In the parametric study, the main observation was that the discretization scheme plays a major role in the quality of the results, as the flow is unidirectional and is characterized by high speed of 2.0 Mach. As a result, dissipation phenomena are not comparable to convection, which plays the major role. The choice of the Upwind discretization scheme seems to be more accurate, as it considers the fact that the flow is unidirectional and provides results improved by more than 1.5% for given turbulence model on the whole computational field. That improvement can go up to 6% when discussing regions behind the lip, where the flow has been stabilized.

The parametric study regarding the turbulence model produced interesting results, as well. In order to evaluate solely the impact of the turbulence model, the discretization scheme used was the Hybrid, as in the original. As it appears, there is not only one fitting turbulence model for this simulation, as their suitability changes with regions. When studying the flow field in the lip region (4–5 cm), the k-ω turbulence model gave the best results, improved by 0.5% when compared to the original simulation performed with k-ε turbulence model. The improvement introduced by the use of the k-ω model in that area is due to the fact that this model can calculate better the effect of the created vortices. When both k-ω turbulence model and Upwind discretization scheme were implied, the total improvement in that area of interest went up to almost 3% (5.6274% discrepancy from the experimental data, compared to 8.42%). However, the study showed that behind the lip (x > 5.2), where flow has been stabilized, the use of the RNG k-ε turbulence model combined with the Upwind discretization scheme improved the mean error by more than 6% (0.3348% compared to 6.5465).

Acknowledgments

We would first like to thank our colleagues in DRA Farnborough, who performed the original step in this modeling effort. Their inspired research was the basis of this current research paper. Moreover, we would also like to thank Dr. Mike Malin, Technical Support Manager of CHAM Ltd., whose assistance was truly helpful to go over technicalities that appeared along the way of this present research.

PHOENICS is licensed by CHAM Ltd., London, UK. The present model may become available from the present authors.