2.1. The gas kinetic model

While the dynamics of monatomic gases is described by the Boltzmann equation, that of the polyatomic gas is described by the Wang-Chang & Uhlenbeck (Reference Wang-Chang and Uhlenbeck1951) equation. To reduce the computational cost, simplified kinetic equations are usually adopted (Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015; Li et al. Reference Li, Zeng, Su and Wu2021). Here, two reduced VDFs $f_0(t, \boldsymbol {x}, \boldsymbol {\xi })$ and $f_1(t, \boldsymbol {x}, \boldsymbol {\xi })$ are used to describe the translational and rotational states of molecular gas, where $t$ is the time, $\boldsymbol {x}=(x_1,x_2,x_3)$ is the spatial coordinate and $\boldsymbol {\xi }=(\xi _1,\xi _2,\xi _3)$ is the molecular velocity. The rotational degree of freedom is a constant $d_r$, and it is assumed that the vibrational and electronic energy levels are not excited, although the model can be easily extended to include these non-equilibrium and even multiphysics effects such as radiation (Li et al. Reference Li, Zeng, Huang and Wu2023). Macroscopic quantities, such as the density $\rho$, flow velocity $\boldsymbol {u}$, deviatoric stress $\boldsymbol {\sigma }$, translational and rotational temperature $T_t$ and $T_r$ and translational and rotational heat flux $\boldsymbol {q}_t$ and $\boldsymbol {q}_r$, are obtained by taking the moments of VDFs:

(2.1a)$$\begin{gather} \left(\rho,\rho\boldsymbol{u},\boldsymbol{\sigma},\frac{3}{2}\rho RT_t,\boldsymbol{q}_{t}\right) =\int\left(1,\boldsymbol{\xi},\boldsymbol{c}\boldsymbol{c}-\frac{c^2}{3}\boldsymbol{\mathsf{I}},\frac{c^2}{2},\frac{c^2}{2}\boldsymbol{c} \right) f_0 \,\mathrm{d}\boldsymbol{\xi}, \end{gather}$$

(2.1b)$$\begin{gather}\left(\frac{d_r}{2}\rho RT_r,\boldsymbol{q}_{r}\right) =\int\left(1,\boldsymbol{c}\right)f_1\,\mathrm{d}\boldsymbol{\xi}, \end{gather}$$

where $\boldsymbol {c}=\boldsymbol {\xi }-\boldsymbol {u}$ is the peculiar (thermal) velocity and $\boldsymbol{\mathsf{I}}$ is the identity matrix. The total temperature $T$ is defined as the equilibrium temperature between the translational and internal modes $T=(3T_t+d_rT_r)/(3+d_r)$. The pressure related to the translational motion is $p_t=\rho R T_t$, whereas the total pressure is $p=\rho RT$.

The evolution of VDFs is governed by the following Boltzmann kinetic equations:

(2.2a)$$\begin{gather} \frac{\partial f_0}{\partial t}+\boldsymbol{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla} f_0 = \frac{g_{0t}-f_0}{\tau}+\frac{g_{0r}-g_{0t}}{Z_r\tau}, \end{gather}$$

(2.2b)$$\begin{gather}\frac{\partial f_1}{\partial t}+\boldsymbol{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla} f_1 = \frac{g_{1t}-f_1}{\tau}+\frac{g_{1r}-g_{1t}}{Z_r\tau}, \end{gather}$$

where $\tau =\mu /p_t$ and ${Z}_{r}\tau$ are the elastic and inelastic collision characteristic time, respectively, with ${Z}_{r}$ being the rotational collision number and $\mu$ being the shear viscosity. The elastic collision conserves the kinetic energy, whereas the inelastic collision exchanges the translational and rotational energies. In order the make the kinetic model mimic the behaviours (such as recovering the shear viscosity, bulk viscosity, translational and internal heat conductivities) of the Wang-Chang & Uhlenbeck (Reference Wang-Chang and Uhlenbeck1951) equation as closely as possible, the reference VDFs to which the gas system will relax are designed as follows:

(2.3a)$$\begin{gather} g_{0t}= \rho\left(\frac{1}{2{\rm \pi} RT_t}\right)^{3/2}\exp\left(-\frac{c^2}{2RT_t}\right)\left[1+\frac{2\boldsymbol{q}_{t}\boldsymbol{\cdot}\boldsymbol{c}}{15RT_tp_t} \left(\frac{c^2}{2RT_t}-\frac{5}{2}\right)\right], \end{gather}$$

(2.3b)$$\begin{gather}g_{0r}= \rho\left(\frac{1}{2{\rm \pi} RT}\right)^{3/2}\exp\left(-\frac{c^2}{2RT}\right)\left[1+\frac{2\boldsymbol{q}_{0}\boldsymbol{\cdot}\boldsymbol{c}}{15RTp} \left(\frac{c^2}{2RT}-\frac{5}{2}\right)\right], \end{gather}$$

(2.3c)$$\begin{gather}g_{1t}=\frac{d_r}{2}RT_rg_{0t} + \left(\frac{1}{2{\rm \pi} RT_t}\right)^{3/2}\frac{\boldsymbol{q}_{r}\boldsymbol{\cdot}\boldsymbol{c}}{RT_t}\exp\left(-\frac{c^2}{2RT_t}\right), \end{gather}$$

(2.3d)$$\begin{gather}g_{1r}=\frac{d_r}{2}RTg_{0r} + \left(\frac{1}{2{\rm \pi} RT}\right)^{3/2}\frac{\boldsymbol{q}_{1}\boldsymbol{\cdot}\boldsymbol{c}}{RT}\exp\left(-\frac{c^2}{2RT}\right), \end{gather}$$

with $\boldsymbol {q}_{0},\boldsymbol {q}_{1}$ being linear combinations of translational and internal heat fluxes (Li et al. Reference Li, Zeng, Su and Wu2021),

(2.4)\begin{equation} \begin{bmatrix} \boldsymbol{q}_{0} \\ \boldsymbol{q}_{1} \end{bmatrix}=\begin{bmatrix} (2-3A_{tt})Z_r+1 & -3A_{tr}Z_r \\ -A_{rt}Z_r & -A_{rr}Z_r+1 \\ \end{bmatrix} \begin{bmatrix} \boldsymbol{q}_{t} \\ \boldsymbol{q}_{r} \end{bmatrix}, \end{equation}

where $\boldsymbol {A}=[A_{tt},A_{tr},A_{rt},A_{rr}]$ is determined by the relaxation rates of heat flux. For nitrogen, we choose $Z_r=3.5$ and $\boldsymbol {A}=[0.786,-0.201,-0.059,0.842]$.

The gas kinetic model describes the gas–gas interaction; in order to fully determine the gas dynamics in wall-bounded problems, the gas–surface interaction should be specified. In this work, we adopt the Maxwell gas–surface boundary condition, where the gas molecules hitting the solid wall will be reflected diffusely.

2.2. The NS equations and the RANS model

Exploiting the relation between macroscopic quantities and mesoscopic VDFs in (2.1), the macroscopic moment equations can be derived from the kinetic equation (2.2) as follows:

(2.5a)$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u}) = 0, \end{gather}$$

(2.5b)$$\begin{gather}\frac{\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \boldsymbol{u}) =-\boldsymbol{\nabla} \boldsymbol{\cdot} (\rho R T_t \boldsymbol{\mathsf{I}} + \boldsymbol{\sigma}), \end{gather}$$

(2.5c)$$\begin{gather}\frac{\partial \rho e}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho e \boldsymbol{u}) =-\boldsymbol{\nabla} \boldsymbol{\cdot} (\rho R T_t \boldsymbol{u} + \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{u} + {\boldsymbol{q}}_{t} + {\boldsymbol{q}}_{r}), \end{gather}$$

(2.5d)$$\begin{gather}\frac{\partial \rho e_r}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho e_r \boldsymbol{u}) =-\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{q}}_{r} + \frac{d_r \rho R}{2}\frac{T-T_r}{{Z}_{r}\tau}, \end{gather}$$

where $e=(3RT_{t}+{\Vert \boldsymbol {u} \Vert }^2)/2+d_r R T_r/2$, $e_r=d_r R T_r/2$. The first three equations in (2.5) describe the conservation of mass, momentum and total energy, respectively, whereas the fourth equation accounts for the exchange between translational and rotational energies.

Equation (2.5) is not closed, since the stress and heat fluxes have not been expressed in terms of the density, velocity, temperature and their spatial gradients. The Chapman–Enskog expansion of the kinetic equation to the first order of $Kn$ yields the following NS constitutive relations (Chapman & Cowling Reference Chapman and Cowling1990):

(2.6a)$$\begin{gather} \boldsymbol{\sigma}^{{NS}} =-\mu (\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T}-\tfrac{2}{3}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\boldsymbol{\mathsf{I}}), \end{gather}$$

(2.6b)$$\begin{gather}\boldsymbol{q}_{t}^{{NS}} =-\kappa_t\boldsymbol{\nabla} T_t,\quad \boldsymbol{q}_{r}^{{NS}} =-\kappa_r\boldsymbol{\nabla} T_r, \end{gather}$$

where the (laminar or physical) shear viscosity is obtained as $\mu =p_t\tau$, whereas the translational and rotational thermal conductivities $\kappa _t$ and $\kappa _r$ are (Li et al. Reference Li, Zeng, Su and Wu2021)

(2.7) \begin{equation} \left[\begin{array}{@{}c@{}} \kappa_t \\ \kappa_r \end{array}\right]= \frac{\mu}{2} \left[\begin{array}{@{}cc@{}} A_{tt} & A_{tr} \\ A_{rt} & A_{rr} \end{array}\right]^{-1}\left[\begin{array}{@{}c@{}} 5 \\ d_r \end{array}\right]= \mu\left[\begin{array}{@{}c@{}} 3.546 \\ 1.435 \end{array}\right]. \end{equation}

Theoretically, turbulence can be investigated by the DNS of (2.5) with the NS constitutive relations (2.6). However, its prohibitive numerical cost prompts the use of turbulence models. Though the RANS models are derived on coarse grids under multiple assumptions, they are widely accepted in engineering due to the low computational cost and good accuracy on capturing the flow pattern. Here, the popular $k$–$\omega$ shear stress transport (SST) model, which combines the advantages of both the low- and high-Reynolds-number $k$–$\omega$ models and embeds Bradshaw's assumption for the boundary layer, is considered (Menter Reference Menter1994). It is good at predicting flows with adverse pressure gradients and separation.

In the SST model, the moment equations (2.5) are used, but the constitutive relations not only contains the physical viscosity (heat conductivity), but also the turbulent ones. That is, $\mu$ and $\kappa$ (irrespective of translational or rotational) in (2.6) are modified as

(2.8a)$$\begin{gather} \mu = \mu_{lam} + \mu_{turb}, \end{gather}$$

(2.8b)$$\begin{gather}\kappa = \kappa_{lam} + \kappa_{turb} = \frac{\mu_{lam} c_p}{Pr_{lam}} + \frac{\mu_{turb} c_p}{Pr_{turb}} = \left(\mu_{lam} + \frac{\mu_{turb}Pr_{lam}}{{Pr}_{turb}}\right) \frac{c_p}{{Pr}_{lam}}, \end{gather}$$

where the subscript $turb$ denotes the turbulent parts, and the laminar to turbulent Prandtl number ratio is set to 0.8, with the value of $c_p/{Pr}_{lam}$ for $\kappa _t$ and $\kappa _r$ obtained with (2.7) accordingly.

The dimensionless turbulent viscosity in (2.8) is given by (here $a_R$ is a constant from the Bradshaw assumption, $F_\mu$ is a blending function that takes the value of one for boundary layer flows and zero for free shear flows, $F_\phi$ is the blending function which transforms the model and its coefficients from $k$–$\omega$ model in near-wall region $(F_\phi =1)$, to the $k$–$\epsilon$ model in the main flow area $(F_\phi =0)$, the model coefficient $\phi =(\sigma _k, \sigma _{\omega }, \beta, \beta ^*, \kappa _{\phi }, \gamma )$ is blended as $\phi =F_\phi \phi _1 + (1-F_\phi ) \phi _2$; these model functions and coefficients are given in Menter (Reference Menter1994))

(2.9)\begin{equation} \mu_{turb} ={a_R \rho k}/{\max\left(a_R \omega, \frac{\Vert\varOmega\Vert F_{\mu}}{{Re}_{ref}}\right)}, \end{equation}

where $k$ represents the turbulent kinetic energy, $\omega =\rho k/\mu _t$ is the turbulent dissipation frequency, with $\Vert \varOmega \Vert$ being the vorticity magnitude. They can be calculated by the following equations (to be consistent with the normalisation in GSIS (Zhang et al. Reference Zhang, Zeng, Yuan, Liu, Li and Wu2024), reference quantities are defined based on free flow properties denoted by subscript 0, i.e. density $\rho _{ref}=\rho _0$, temperature $T_{ref}=T_0$, velocity $v_{ref}=\sqrt {RT_0}$, pressure $p_{ref}=\rho _0 R T_0$, heat flux $q_{ref}=p_{ref}v_{ref}$, viscosity $\mu _{ref}=\rho _{ref}v_{ref}L$, turbulent kinetic energy $k_{ref}=v^2_{ref}$ and turbulent dissipation frequency $\omega _{ref}=\rho _{ref}v^2_{ref}/\mu _{ref}$):

(2.10)\begin{gather} \frac{\partial{\rho k}}{\partial{t}}+\frac{\partial{\rho u_j k}}{\partial{x_j}}=\frac{1}{{Re}_{ref}} \frac{\partial}{\partial{x_j}} \left[ (\mu_{lam}+\sigma_k \mu_{turb})\frac{\partial{k}}{\partial{x_j}} \right] +\text{Prod} -\underbrace{{Re}_{ref}\beta^* \rho \omega k}_{{\rm Diss}_k}, \end{gather}

(2.11)\begin{gather}\begin{aligned}\frac{\partial{\rho \omega}}{\partial{t}}+\frac{\partial{\rho u_j \omega}}{\partial{x_j}}&=\frac{1}{{Re}_{ref}} \frac{\partial}{\partial{x_j}} \left[(\mu_{lam}+\sigma_\omega \mu_{turb})\frac{\partial{\omega}}{\partial{x_j}} \right]+\gamma \frac{\rho}{\mu_{turb}}\tau_{turb,ij}\frac{\partial{u_i}}{\partial{x_j}}\nonumber\\&\quad -{Re}_{ref}\beta \rho \omega^2 +\frac{2 (1-F_\phi)}{{Re}_{ref}} \frac{\rho \sigma_{\omega 2}}{\omega} \frac{\partial{k}}{\partial{x_j}}\frac{\partial{\omega}}{\partial{x_j}}. \end{aligned}\end{gather}

Note that $\boldsymbol {\tau }_{turb} = \mu _{turb}(\boldsymbol {\nabla }{\boldsymbol {u}}+ \boldsymbol {\nabla }\boldsymbol {u}^{T}-\frac {2}{3}\boldsymbol {\nabla }\boldsymbol{\cdot}\boldsymbol {u}\boldsymbol{\mathsf{I}}) - \frac {2}{3} \rho k \boldsymbol{\mathsf{I}}$ is the Favre-averaged Reynolds stress computed under the Boussinesq eddy-viscosity hypothesis. To eliminate certain erroneous spikes of $\mu _{turb}$ rooted in two-equation turbulence models (Menter Reference Menter1993), the production term in (2.10) is $\text {Prod}= \min (\tau _{turb,ij}({\partial {u_i}}/{\partial {x_j}}), 20 {\rm Diss}_k)$. At solid surface, $k=0$ and $\omega =({60\mu _{lam}})/({{Re}_{ref}{\beta }_{1}\rho D_{1}^{2}})$ are specified, with $D_{1}$ being the distance of the first cell centre to the wall.

The SST model, originally developed for incompressible flows, does not account for density fluctuations. However, in hypersonic or high heat transfer scenarios, three terms in the density-weighted Favre-averaged equation for turbulent kinetic energy are not modelled: pressure-dilatation, pressure work, and curl-free dilatation dissipation. These terms are crucial as they contribute to the reduction of turbulent kinetic energy in flows with $Ma>5$. In this study, the pressure-dilatation and dilatation dissipation terms are taken into account as per Sarkar (Reference Sarkar1992), while the pressure work is discarded due to the lack of satisfactory correction model and proof of correctness for coupling differently established correction models of different terms. The final form of the implemented SST model reads:

(2.12)\begin{align} \frac{\partial{\rho k}}{\partial{t}}+\frac{\partial{\rho u_j k}}{\partial{x_j}}&=\frac{1}{{Re}_{ref}} \frac{\partial}{\partial{x_j}} \left[ (\mu_{lam}+\sigma_k \mu_{turb})\frac{\partial{k}}{\partial{x_j}} \right] +\text{Prod}-\underbrace{{Re}_{ref}\beta^* \rho \omega k}_{{\rm Diss}_k}\nonumber\\ &\quad -\underbrace{{Re}_{ref}\beta^* \xi^* M^2_t \rho \omega k}_{\textit{Dilatation dissipation}} +\underbrace{\alpha_2 \tau_{turb,ij}\frac{\partial{u_i}}{\partial{x_j}} M_t + {Re}_{ref} \alpha_3 \beta^* \rho \omega k M^2_t}_{\textit{Pressure dilatation}}, \end{align}

(2.13)\begin{align} \frac{\partial{\rho \omega}}{\partial{t}}+\frac{\partial{\rho u_j \omega}}{\partial{x_j}}&=\frac{1}{{Re}_{ref}} \frac{\partial}{\partial{x_j}} \left[ (\mu_{lam}+\sigma_\omega \mu_{turb})\frac{\partial{\omega}}{\partial{x_j}} \right] +\gamma \frac{\rho}{\mu_{turb}} \tau_{turb,ij}\frac{\partial{u_i}}{\partial{x_j}}\nonumber\\ &\quad -{Re}_{ref}\beta \rho \omega^2 +\frac{2 (1-F_\phi)}{{Re}_{ref}} \frac{\rho \sigma_{\omega 2}}{\omega} \frac{\partial{k}}{\partial{x_j}}\frac{\partial{\omega}}{\partial{x_j}} +\underbrace{{Re}_{ref}\beta^* \xi^* M^2_t \rho \omega^2}_{\textit{Dilatation dissipation}}, \end{align}

where $M_t=\sqrt {2k}/\sqrt {\gamma R T}$ is the turbulent Mach number. Such corrections work fine for jets and free shear flow at low Reynolds number. Other model coefficients ($\alpha _2, \alpha _3, \xi ^*$) are summarised by Wilcox (Reference Wilcox2006).

2.3. The GSIS-SST model for coexisting turbulent and rarefied flows

To investigate the rarefied gas dynamics, the kinetic equation (2.2) should be solved numerically. However, since the VDFs are defined in six-dimensional phase space, the computational cost is much higher than the NS equations. Fortunately, since in most engineering problems the rarefied gas flows are laminar, only the steady states are considered. Therefore, the implicit iteration methods are normally used. Here, we first briefly introduce the recently developed GSIS for efficiently simulating the rarefied laminar flows, and then propose the GSIS-SST model for coexisting turbulent and rarefied flows.

The GSIS is a multiscale method that alternatively solves the mesoscopic kinetic equation and macroscopic synthetic equation (Su et al. Reference Su, Zhu, Wang, Zhang and Wu2020a; Zhang et al. Reference Zhang, Zeng, Yuan, Liu, Li and Wu2024). The synthetic equation is given by (2.5), but now the stress and heat flux are constructed as

(2.14a)$$\begin{gather} \boldsymbol{\sigma}=-\mu \left(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T}-\frac{2}{3}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\boldsymbol{\mathsf{I}}\right) + \text{HoT}_{\boldsymbol{\sigma}}, \end{gather}$$

(2.14b)$$\begin{gather}\boldsymbol{q}_{t} =-\kappa_t\boldsymbol{\nabla} T_t + \text{HoT}_{\boldsymbol{q}_{t}},\quad \boldsymbol{q}_{r}=-\kappa_r\boldsymbol{\nabla} T_r + \text{HoT}_{\boldsymbol{q}_{r}}, \end{gather}$$

where the high-order terms (HoTs) describing the rarefaction effects are constructed as

(2.15a)$$\begin{gather} \text{HoT}_{\boldsymbol{\sigma}} = \int \left(\boldsymbol{c}\boldsymbol{c}-\frac{c^2}{3}\boldsymbol{\mathsf{I}}\right)f_0 \,\mathrm{d}\boldsymbol{\xi} -\boldsymbol{\sigma}^{{NS}}, \end{gather}$$

(2.15b)$$\begin{gather}\text{HoT}_{\boldsymbol{q}_{t}} = \int \boldsymbol{c} \frac{c^2}{2}f_0 \,\mathrm{d}\boldsymbol{\xi} -\boldsymbol{q}_{t}^{{NS}},\quad \text{HoT}_{\boldsymbol{q}_{r}} = \int \boldsymbol{c}f_1\, \mathrm{d}\boldsymbol{\xi} -\boldsymbol{q}_{r}^{{NS}}, \end{gather}$$

with $\boldsymbol {\sigma }^{{NS}}$, $\boldsymbol {q}_{t}^{{NS}}$ and $\boldsymbol {q}_{r}^{{NS}}$ calculated based on (2.6) using the macroscopic properties from the moments of the VDFs as well.

As sketched in figure 1, the kinetic solver provides HoTs (it must be emphasised that, as shown in figure 1, (2.14) and (2.15) are solved at different iteration steps, thus the NS constitutive relations cannot cancel each other out until the converged solution is found) to the synthetic equation (2.5), whereas the solution of synthetic equation, when solved to the steady state, steers the evolution of VDFs towards the steady state. As has been demonstrated in various numerical simulations, this kind of treatment facilities the fast convergence in the whole range of gas rarefaction (Su et al. Reference Su, Zhu, Wang, Zhang and Wu2020a; Liu et al. Reference Liu, Zhang, Zeng and Wu2024a; Zhang et al. Reference Zhang, Zeng, Yuan, Liu, Li and Wu2024). Furthermore, according to the Chapman–Enskog expansion, in the continuum limit, the HoTs are proportional to $Kn^2$ (Su, Zhu & Wu Reference Su, Zhu and Wu2020b), such that the constitutive relation (2.14) asymptotically preserves the NS limit when $Kn$ is low.

Figure 1.

Schematic of the GSIS and GSIS-SST. The kinetic equation is solved by the discrete velocity method iteratively (here $k$ is the iteration step). Between each iteration, the HoTs extracted from the kinetic solver are fed into the macroscopic synthetic equation, which is solved to the steady state for $m\approx 200$ steps. Then, the updated macroscopic quantities are fed into the kinetic solver through the modification of VDFs. Detailed implementation of GSIS are given by Liu et al. (Reference Liu, Zhang, Zeng and Wu2024a). Different to the pure GSIS, in GSIS-SST the turbulent viscosity and heat conductivities are added in the macroscopic synthetic equation.

To effectively simulate turbulence, even within the scope of the NS equations, it is essential to employ turbulence models. Given the significantly higher computational cost associated with the kinetic solvers, the incorporation of turbulence models is in high demand. Thanks to the explicit inclusion of the NS constitutive relations in (2.14), the RANS model can be integrated seamlessly into the GSIS framework. That is, once the turbulence viscosity from the SST model is obtained, it is superimposed into the physical viscosity in the NS equations. Consequently, the stress and heat fluxes are formulated as follows:

(2.16a)$$\begin{gather} \boldsymbol{\sigma} =-(\mu_{lam} + \mu_{turb})(\boldsymbol{\nabla}{\boldsymbol{u}}+\boldsymbol{\nabla}\boldsymbol{u}^{T} -\tfrac{2}{3}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\boldsymbol{\mathsf{I}}) + \text{HoT}_{\boldsymbol{\sigma}}, \end{gather}$$

(2.16b)$$\begin{gather}\boldsymbol{q}_{t} =-(\kappa_{t,lam}+\kappa_{t,turb})\boldsymbol{\nabla}T_t + \text{HoT}_{\boldsymbol{q}_{t}},\quad \boldsymbol{q}_{r} =-(\kappa_{r,lam}+\kappa_{r,turb})\boldsymbol{\nabla}{T_r} + \text{HoT}_{\boldsymbol{q}_{r}}. \end{gather}$$

The general algorithm of the GSIS-SST is sketched in figure 1.

3.1. High Reynolds number

According to the Chapman–Enskog expansion (Su et al. Reference Su, Zhu and Wu2020b; Su, Zhang & Wu Reference Su, Zhang and Wu2021), when $Kn$ is low, the HoTs are proportional to $Kn^2$, such that the GSIS-SST with the constitutive relation (2.16) asymptotically preserves the pure SST model. This is validated in the boundary layer problem in turbulent flows, which is a standard problem included in the NASA turbulence modelling resource for model development and benchmarking.

As shown in figure 2, we consider the hypersonic flow of $Ma=9.22$ passing through a wedge. Taking 1 m as a reference length and the free flow temperature of 64.5 K as a reference temperature, the Reynolds number is $4.7\times 10^{7}$ and the Knudsen number is $2.91\times 10^{-7}$. The solid wall temperature is 295 K. A total of 31 284 nonuniform spatial grids are used, where the height of the first layer of cells is 1 ${\rm \mu}$m. The two-dimensional velocity space is discretised uniformly into $60\times 40$ quadrilaterals in the range of $\arrowvert \xi _1\arrowvert \times \arrowvert \xi _2 \arrowvert = 30\times 20$. In the SST model, initial conditions on free flow boundary are specified using a turbulence intensity $I_t=0.3\,\%$ and a turbulent to laminar viscosity ratio $\mu _{r}=\mu _{turb}/\mu _{lam}=15$. Then $k$ and $\omega$ are derived as $k=1.5(I_t \Vert \boldsymbol {u} \Vert )^2$ and $\omega =\rho k /\mu _{turb}$.

Figure 2.

(a) Geometry and the computational grid in the turbulent boundary layer problem (unit: metres). Red line on the left-hand side stands for in-flow condition, blue lines on the top, bottom and right-hand side are out-flow condition, whereas green lines are stationary isothermal wall. (b) Velocity $u_x$ and Mach number. Experimental data are sampled at $x_1=0.076$ m on the upper flat surface (Marvin, Brown & Gnoffo Reference Marvin, Brown and Gnoffo2013). Reference solutions are obtained from conventional CFD solver. Here $Y_n$ is the vertical distance normalised by the boundary layer thickness.

Figure 2 shows that the boundary layer velocity and Mach number obtained from GSIS-SST fit well with the experimental data. Furthermore, the pure SST model is solved by the conventional CFD solver, which also agrees with the GSIS-SST result. Thus, the asymptotic behaviour of the GSIS-SST in turbulence modelling is confirmed.

3.2. High Knudsen number

At high $Kn$, the turbulent viscosity $\mu _{turb}$ will be negligible compared with the physical viscosity $\mu _{lam}$, so that the GSIS-SST becomes the pure GSIS which exactly solves the Boltzmann kinetic equation. To prove this, we simulate the hypersonic flow of $Ma=25$ passing over a two-dimensional flat plate with blunt leading edge, by both the GSIS and GSIS-SST. The computational geometry is shown in figure 3(a). The thickness of the plate is 0.03 m, and it is chosen as the reference length. The physical domain ($-0.9\ {\rm m}\le x_1 \le 0.09\ {\rm m}$, $|x_2| \le 0.9$ m) is discretised into 136 240 cells, with 381 nodes along and 360 normal to the model surface. The determination of velocity space follows the $5\sigma$ criteria of the standard normal distribution, where $\sigma =\sqrt {RT_0}$ is calculated based on the stagnation temperature of free flow. The final range and grid spacing is adjusted based on a grid convergence study in velocity space. The refined velocity space is finally defined in the range of $-90 \le \xi _1 \le 90$ and $-70 \le \xi _2 \le 70$, and is cut into $360\times 280$ uniform cells.

Figure 3.

Geometry and simulation results from the GSIS and GSIS-SST in the hypersonic flow of $Ma=25$ passing over a blunt leading edge, where the Knudsen number is $Kn=1$: (a) geometry and grid (unit: metres); (b) total temperature; (c) velocity and viscosity at $x_1 = 0.014$ m; and (d) surface heat flux.

The nitrogen gas flow of Mach number 25 is coming from the left boundary, with a temperature equal to the plate temperature. Turbulent conditions on free flow boundary stay the same as the flat plate case. The Knudsen number is 1 under the free stream condition.

Figure 3(b) shows the temperature contours, where the shock width is a few mean free paths ($\lambda =0.03$) and is comparable to the plate thickness. The GSIS and GSIS-SST solutions overlap perfectly, enlightening that the turbulence model loses its effect in highly rarefied flows. The velocity $u_x$ at $x_1=0.014$ m is plotted in figure 3(c), which exhibits a large slip at the solid surface and changes rather smoothly over quite a long distance, indicating the absence of turbulent boundary layer and the dominance of the rarefaction effects. This subfigure also shows that the turbulent viscosity $\mu _{turb}$ is nearly zero, which is much lower than the physical viscosity $\mu _{lam}$. This is because, in such highly rarefied flow environment and small velocity gradients, the density of near-wall flow and the production of turbulent kinetic energy are both small, leading to insignificant level of $\mu _{turb}$ as per (2.9). Therefore, the GSIS-SST degenerates to the pure GSIS for highly rarefied flows. This is further supported by the results in figure 3(d), where the surface heat fluxes predicted by GSIS and GSIS-SST overlap. Meanwhile, the pure NS and NS-SST results overlap, again proving the absence of turbulent effects in such low Reynolds number (i.e. ${Re}\approx 25$). in this way, the difference between the NS and GSIS results shows that the rarefaction effects are important when $Kn=1$, where the NS constitutive relations lose validity.

3.3. Crossover between turbulent and rarefied flows

Having analysed the asymptotic behaviour of GSIS-SST well into the turbulent and rarefied regimes, it is interesting to test its performance at intermediate Knudsen numbers. Since the pure GSIS solves the Boltzmann kinetic equation in the steady state, it is expected that it is inaccurate above certain Reynolds number. By comparing the solutions of GSIS with GSIS-SST, this critical $Re$ can be roughly identified.

Thus, three additional Knudsen numbers ($Kn=0.01$, 0.001 and 0.0001) are simulated in the hypersonic flow passing over a blunt leading edge. The relative error in surface heat fluxes are compared between the GSIS and GSIS-SST in figure 4(a). It is seen that the two methods start to depart from each other sensibly at $Kn\approx 0.001$, where the largest relative error reaches around 7 % (slightly higher than the engineering accepted range of 5 %). At $Kn=0.0001$, the relative error reaches 30 %, indicating the turbulent effects are significant. This example shows that the influence of the SST model on surface quantities is effectively off in the GSIS-SST when $Re \lessapprox 25\,000$.

Figure 4.

(a) Relative difference in the surface heat fluxes between the GSIS and GSIS-SST solutions, i.e. $(q_{GSIS\text -SST}/q_{GSIS}-1)$, when $Kn=0.01$ (solid line), 0.001 (dash-dotted line) and 0.0001 (dashed line). (b,c) Normal heat flux $q_2$ obtained from GSIS-SST at $x_1=0.014$ m when $Kn=0.01$ and $Kn=0.0001$, respectively. Contributions from the laminar (solid lines) and turbulent (dashed lines) parts, as well as the HoTs (circles), are compared. (Since there is a cold wall, when Knudsen number drops and the temperature jump near the wall diminishes, $\mu _{lam}$ experiences a sharp decrease. Therefore, a peak of $q_{2,lam}$ is produced even though the temperature gradient increases near the wall.)

We further compare in figure 4(b,c) the three components of heat flux obtained from GSIS-SST, i.e. $q_{2,lam}$ from the NS constitutive relations, $q_{2,HoT}$ from the high-order terms for rarefaction effects and $q_{2,turb}$ from the turbulence modelling. At $Kn=0.01$, the laminar part constitutes most of the heat flux, whereas the turbulent part is zero. A small contribution from the HoT, which represents the rarefaction effects, is seen in the vicinity of the solid wall. When $Kn$ is reduced to 0.0001, the laminar part of heat flux still dominates. The turbulent part of heat flux arises from zero at the wall as per the boundary condition of the SST model, quickly reaches a peak value of about one-third of the translational heat flux, and then decay to zero away from the wall. Since the HoT vanishes (which is consistent with the Chapman–Enskog analysis), the rarefaction effects are absent, and the GSIS-SST is effectively the SST model.

In the intermediate region, e.g. $Kn=0.001$, figure 5(a–c) shows that, along the blunt head, the laminar part dominates the surface heat flux $q_{1}$, the turbulent part also plays a noticeable role, whereas the HoT has negligible contribution around the blunt head, as the high-density stagnation point is not that rarefied. Since the GSIS does not take the turbulent effects into account, it predicts a surface heat flux 7 % smaller than that of the GSIS-SST at the stagnation point, see figure 4(a). Along the plate, the laminar and turbulent parts are absent due to the absence of horizontal temperature gradient, however, the HoT exists due to the rarefaction effects. This is because this region is downstream the flow expansion, it is much more rarefied than the stagnation region. As for $q_{2}$, figure 5(d–f) shows that the laminar and turbulent parts still make up most of the heat flux along the blunt head. Along the plate, the laminar part also contributes greatly to $q_{2}$, whereas the turbulent part is still absent, due to the reduction in density and turbulence production in this expansion region. The HoT part makes a constant presence on the whole model surface.

Figure 5.

(a–c) Contours of horizontal heat flux $q_{1}$: (a) $q_{1,lam}$, (b) $q_{1,turb}$ and (c) $q_{1,HoT}$. (d–f) Contours of vertical heat flux $q_{2}$: (d) $q_{2,lam}$, (e) $q_{2,turb}$ and (f) $q_{2,HoT}$. The contour levels are limited in the range of $-$40 to 40 for a better view. The Knudsen number is $Kn=0.001$.

In summary, the asymptotic analysis and numerical simulations show that the GSIS-SST can automatically select the SST model and the laminar Boltzmann equation based on the flow conditions.

4.1. Numerical configurations

The physical model is shown in figure 3(a), except that a 0.004-m-wide slot is now placed at the head of the blunt leading edge, from which pressurised nitrogen is injected towards the incoming free-stream flow of Mach number 25. The jet is characterised by the pressure ratio $P_{ratio}$ between the jet flow and the post-shock value of free flow (a constant pressure ratio could ensure a quite similarly located jet flow field, even if all free flow conditions change significantly (Tian, Duan & Chen Reference Tian, Duan and Chen2023)):

(4.1)\begin{equation} P_{ratio} = \frac{p_{j,0}}{p_{f,s,0}}=2.5, \end{equation}

where $p_{j,0}$ is the stagnation pressure of the jet and $p_{f,s,0}$ is the stagnation pressure of the incoming shock in the downstream (calculated with the normal–shock relation based on free-stream quantities). The jet has a static temperature of 250 K and an exit Mach number of 1. As a consequence, the jet density is about 800–1000 times higher than the incoming free-stream density. The turbulent boundary condition for free flow is the same as that used in § 3.1. A characteristic length scale $l_t=0.118$ (non-dimensionalised) is provided to the jet, and $I_t=3\,\%$ for the jet. Then $k$ and $\omega$ are derived as $k=1.5(I_t \Vert \boldsymbol {u} \Vert )^2$ and $\omega =\sqrt {k}/l_t$.

We consider the hypersonic flow at altitudes of 80 and 62 km, where the Knudsen numbers are 0.125 and 0.01, respectively, when the reference length is chosen as the diameter of the model leading edge. The velocity grid is the same as that used in the validation section. Since in the jet problem complex flow structures exist ahead of the model, much higher grid density is required. To benefit from symmetry and reduce the total cell count, the spatial domain for the jet study is half ($x_2\geq 0$) that in figure 3(a). A study on spatial grid independence is carried out using the coarse (comprising 68 210 quadrilateral cells), moderate (with 119 682 cells) and fine (containing 162 797 cells) grids. The primary distinction among these grids is the number of nodes in the wall-normal direction and along the flat section of the model. The results for surface quantities are found to be consistent across all grids. Considering both computational cost and accuracy, the moderate grid resolution of 119 682 cells (with 20 nodes for jet boundary, 200 for wall and 550 normal to wall) was selected for further analysis.

4.2. Jet flow field under a rarefied environment

We use the GSIS and GSIS-SST to investigate the role of a turbulent jet in a rarefied environment. The NS and NS-SST equations are not used since they do not describe the rarefaction effects, e.g. see figure 3(d).

Figure 6 shows the temperature contours when $Kn=0.01$ and 0.125. The under-expanded jet is injected into the flow field from the nozzle at $x_1=-0.0150$ m and is confined by the barrel shock. It expands and accelerates vigorously until it encounters the free flow, resulting in a pronounced Mach disk. The Mach disk acts to compress and decelerate the now over-expanded jet flow. Beyond the Mach disk, the jet flow achieves a local equilibrium with the post-shock free flow, subsequently flowing laterally and in reverse, reattaching to the model surface. A recirculation zone is present between the barrel shock, the reversed jet flow and the model surface. This zone arises from two separation phenomena: the first caused by the jet's expansion at the nozzle and the second induced by the reversed jet flow's reattachment to the model surface. The reattachment stimulates the development of a recompression shock, and the change in surface curvature from the blunt to the flat section leads to flow expansion near $x_1=0$ m. The jet displaces the high-temperature post-shock region, creating a low-temperature sheath that exists between the hot outer shock and the model. This cool sheath is roughly delineated by the streamlines emanating from the jet exit. Consequently, this cool sheath is also known as the jet-controlled area.

Figure 6.

Temperature contours from GSIS (upper half) and GSIS-SST (lower half) and velocity streamlines in GSIS-SST: (a) $Kn=0.01$ and (b) $Kn=0.125$.

When $Kn$ increases, the outer shock gets thicker due to the increased rarefaction effects. In the jet-controlled area, flow structures are sharp and clear. The position of Mach disk moves slightly closer to the jet exit, while the shape of Mach disk becomes flatter. When $Kn=0.01$, the post-shock stagnation pressure of the free flow is actually lower than lateral positions, since the oblique outer shock is weaker laterally, thus a curved Mach disk is formed. However, when $Kn=0.125$, the rarefaction effect weakens the outer shock and reduces such lateral difference, thus producing a flatter Mach disk.

Figure 7 reveals another change of flow pattern under the rarefied environment: the production of an additional inflection point (around $x_1=-0.027$ m) along the centreline. The free flow coming from the left experiences a drop in $u_x$, whereas the under-expanded jet flow coming from $x_1=-0.015$ m experiences an increase of magnitude in $u_x$. Normally, this change of velocity is monotonic. However, figure 7(b) shows $u_x$ follows a decrease–increase–decrease process. Such effect becomes stronger when $Kn$ increases. We believe this counterintuitive change of velocity profile is the outcome of a microscopic mixing process, in which molecules from both directions flow into each other for a certain range without sufficient collision when $Kn$ is large. In contrast, in a continuum regime such layer is too thin to be noticed.

Figure 7.

(a) Contour of horizontal velocity with streamlines from GSIS at $Kn=0.01$ (with $x_2\leq 0$ mirrored from upper half). (b) Horizontal velocity along the centreline at $Kn=0.01$ and 0.125, obtained from the GSIS and GSIS-SST. Inset: enlarged view showing the diffusion layer.

On the influence of turbulence, GSIS and GSIS-SST contours and plots in figures 6 and 7 look alike in general shape under the same $Kn$ condition. Since the introduced turbulence functions as an extra viscosity, structures with steep changes are generally smoothed in GSIS-SST. Another alternation is the shift of Mach disk, which stands slightly closer to the model in GSIS-SST. Due to the enhanced mixing of the jet and free flow, the thickness of the low-temperature layer normal to the wall shrinks; see the blue regions around the model surface in figure 6. These changes finally build up the difference on surface quantities distributions, as analysed in the following.

4.3. Interaction between turbulent and rarefied flows

4.3.1. Velocity and viscosity profiles

The boundary layer is examined to study the interaction between the jet, the rarefaction effect and the turbulence in the near-wall region. Furthermore, the turbulent viscosity $\mu _{turb}$ is calculated to understand how the turbulence shapes the flow field.

When $Kn=0.01$, figure 8(a) shows that the velocity profile from GSIS is well layered. From top to bottom, the free flow region, outer shock, free-flow-dominated post-shock layer, shear/mixing layer between the free and jet flow, jet-flow-dominated layer and near-wall boundary layer could be identified. The typical boundary layer still exists, where the flow velocity increases from 0 across the boundary layer to 2 and remains at the jet flow velocity until $x_2\approx 0.019$ m. Then, the jet flow starts to mix with free flow and $u_x$ keeps increasing. After a turning point at $x_2\approx 0.035$ m, the free-flow-dominated layer is reached and the rate of increase of $u_x$ slows down. Through the outer shock with large gradients, $u_x$ finally reaches the undisturbed free flow velocity of 29.6 (Mach number 25). The GSIS-SST result shows generally the same pattern as GSIS, but the transitions between layers are smoothed. In addition, $u_x$ from GSIS-SST increases faster inside the boundary layer than the GSIS.

Figure 8.

Profiles of horizontal velocity and viscosity at $x_1=0.006$ m: (a) $Kn=0.01$, (b) $Kn=0.01$, (c) $Kn=0.125$ and (d) $Kn=0.125$. Inset: velocity and viscosity in the vicinity of a solid wall.

It is the turbulence viscosity that shapes this difference. The viscosity profile from GSIS-SST at $Kn=0.01$ is depicted in figure 8(b), which shows that $\mu _{turb}$ increases from 0 at the wall, rises up and crosses the physical viscosity $\mu _{lam}$ curve around $x_2=0.0155$ m; within this range $\mu _{turb}$ is produced due to large velocity gradient of $u_x$ in $x_2$ direction. A small peak of $\mu _{turb}$ is found at $x_2\approx 0.016$ m, this portion of $\mu _{turb}$ is produced from the vortex in the recirculation region and is transported downstream. The turbulence viscosity dominates over the physical viscosity in the range of $0.0155< x_2<0.0308$ m (turbulent jet flow and shear/mixing layer) as the jet flow is a source of turbulence. It also forms a local peak at around $x_2\approx 0.05$ m (right behind the outer shock) due to large velocity gradients across shock. The high turbulence viscosity in these regions boosts the momentum and energy exchange. High-speed and high-temperature post-shock free flow introduce more kinetic and internal energies through intense turbulent mixing into the lower layers (shear/mixing layer and subsequently the jet flow). Therefore, the abrupt transitions around $x_2\approx 0.0225$ m and $x_2\approx 0.035$ m in GSIS's horizontal velocity profile is smoothed in GSIS-SST, see figure 8(a).

When $Kn$ goes up to 0.125, velocity profiles from GSIS and GSIS-SST become smoother, so that the transition of different velocity layers in figure 8(a) is now blurred in figure 8(c). This is the result of rarefaction effects: since the free flow is rather rarefied and the density of the jet flow is also decreased (as the pressure ratio $P_{ratio}$ in (4.1) is fixed), the diffusion between different layers is intensified, disturbance travels further into each other, leading to higher temperature and momentum around the model. In GSIS-SST, adding together the rarefaction effects and turbulence, the velocity profile is no doubt smoother. The viscosity profiles figure 8(d) show similar pattern as those of $Kn=0.01$, yet the turbulence-viscosity-dominated area in GSIS-SST is now restricted to jet-flow-dominated area, i.e. when $x_2<0.0243$ m. The local peak of $\mu _{turb}$ at $x_2\approx 0.05$ in $Kn=0.01$ still exists at $x_2\approx 0.0675$ in $Kn=0.125$, but is now now sharper than that at lower $Kn$, and smaller than the corresponding laminar viscosity for a wider range as $x_2$ increases. This is the result of the higher level of rarefaction at $Kn=0.125$, since $\mu _{turb}$ is proportional to density.

Figure 9 further depicts the ratio between the turbulent and laminar viscosities in GSIS-SST. When $Kn=0.01$, $\mu _{turb}$ exhibits a pronounced concentration towards the Mach disk, the outer shock (around $x_1=-0.035$ m) and the recompression shock (around $x_1=-0.01$ m). These structures are associated with high-velocity gradients which promotes production of turbulent kinetic energy and, hence, $\mu _{turb}$. Unlike $\mu _{turb}$, the viscosity ratio $\mu _r$ has only high values in the cooler jet-controlled area. Since the jet temperature is much lower and, hence, lower $\mu _{lam}$, $\mu _{turb}$ could be dominant in these regions. When $Kn=0.125$, $\mu _{turb}$ still shows high concentration in the region between the outer shock and the Mach disk, although the recompression shock contributes now less to the formation of $\mu _{turb}$ compared with the $Kn=0.01$ case. The viscosity ratio $\mu _r$ at $Kn=0.125$ shows similar distribution as that at $Kn=0.01$, yet the high $\mu _r$ region shrinks at $Kn=0.125$. It is clear that the high $\mu _{turb}$ region (both at $Kn=0.01$ and $Kn=0.125$) behind the outer shock does not turn into actual effect, as $\mu _r$ there is quite low. This is because the high post-shock temperature largely increases $\mu _{lam}$, diluting the influence from $\mu _{turb}$.

Figure 9.

Contour of turbulent viscosity $\mu _{turb}$ (top half) and the turbulent to laminar viscosity ratio $\mu _r$ from GSIS-SST (bottom half): (a) $Kn=0.01$ and (b) $Kn=0.125$. The initial jet, the post-Mach disk region and the reattachment region are the main sources of turbulence effect.

4.3.2. Three contributions to shear stress

Figure 10 shows the three contributions of shear stress in GSIS-SST simulation. At $Kn=0.01$, $\sigma _{12,lam}$ starts from 0 in the free flow region and rises at $x_2\approx 0.07$ m, reaching rapidly to its peak at the core of outer shock ($x_2\approx 0.065$ m). Its magnitude then falls quickly to 0 in the free-flow-dominated post-shock layer. Then $\sigma _{12,lam}$ rises again when $x_2$ drops below 0.04 m to the shear layer. With a secondary peak in the centre of shear layer, $\sigma _{12,lam}$ drops again and stays near 0 in the jet flow layer. Finally, in the boundary layer, $\sigma _{12,lam}$ exhibits absolute dominance, especially in the viscous sublayer. Here $\sigma _{12,HoT}$, which measures the rarefaction effects, is significant in the shock layer. In addition, a slight positive peak around 0.0315 m in the shear layer could be identified. The profile of $\sigma _{12,turb}$ runs in a similar manner as $\sigma _{12,lam}$. Notably, $\sigma _{12,turb}$ dominates over $\sigma _{12,lam}$ below $x_2\approx 0.03$ m and above boundary layer, respecting the assumption of turbulence domination in jet-controlled area.

Figure 10.

(a,b) Three contributions of shear stress $\sigma _{12}$ at $x_1=0.006$ m from GSIS-SST, namely, $\sigma _{12,lam}$ from molecular viscosity, $\sigma _{12,HoT}$ from rarefaction effects and $\sigma _{12,turb}$ from turbulence model: (a) $Kn=0.01$ and (b) $Kn=0.125$. The plot of $Kn=0.125$ is limited to $x_2<0.1$ m as the shear stress decays to zero above this point. (c–h) Contours of the three contributions of $\sigma _{12}$ at $Kn=0.01$ and 0.125, respectively: (c) $\sigma _{12,lam}$, (d) $\sigma _{12,turb}$, (e) $\sigma _{12,HoT}$, (f) $\sigma _{12,lam}$, (g) $\sigma _{12,turb}$ and (h) $\sigma _{12,HoT}$. The arrow marks the $\sigma _{12}$ sampling location $x_1=0.006$ m.

The contour of three contributions in the second row of figure 10 shows that the laminar part concentrates to locations with high velocity gradients, such as shock, shear layer and boundary layer. It is biased towards free flow in the shear layer due to higher temperature. The turbulent part looks similar to the laminar part but with higher magnitude, especially in the jet-controlled area. The HoT is negligible when compared with the laminar and turbulent parts, as it is proportional to $Kn^2$ in the bulk flow as per the Chapman–Enskog expansion.

When $Kn$ increases from 0.01 to 0.125, the pattern of rise and drop remains similar for $\sigma _{12,lam}$. However, the near-zero parts of $\sigma _{12,lam}$ curves (at $Kn=0.01$, $x_2\approx 0.05$) disappear now at $Kn=0.125$, due to the increased mean free path so that the two peaks in $\sigma _{12,lam}$ mix. The $\sigma _{12,turb}$-dominated area is reduced, representing the degraded influence from turbulence at high $Kn$, even in the less-rarefied jet layer. In contours at $Kn=0.125$ shown in the third row of figure 10, the laminar part is less concentrated compared with that at $Kn=0.01$, as the whole flow field is more rarefied. The largest difference is that the relative magnitude of HoT compared with the other two components increases sensibly, and it also takes effect more at the outer shock ($x_1=-0.035$ m).

4.3.3. Three contributions to heat flux

Figure 11 reveals the three contributions of heat flux in the jet-controlled area when $Kn=0.01$. The laminar part is significant in the shear layer, Mach disk and boundary layer. It also has noticeable magnitude around the recirculation region. In contrast to the laminar part, the high-magnitude turbulent part occupies a larger area. Therefore, significant change on heat flux is expected on the blunt head. The HoT is relatively lower than in the turbulent part, but outmatches the laminar part in certain regions, such as the barrel shock. It also produces a large negative contribution just before the jet flow enters the shear layer, where the laminar and turbulent parts are positive.

Figure 11.

Contours of the three contributions of heat fluxes at $Kn=0.01$: (a) $q_{1,lam}$, (b) $q_{1,turb}$, (c) $q_{1,HoT}$, (d) $q_{2,lam}$, (e) $q_{2,turb}$ and (f) $q_{2,HoT}$.

When $Kn=0.125$, figure 12 shows that the pattern of laminar heat flux is close to that at $Kn=0.01$, except being more dispersed. The turbulent part loses its influence in most of the area. It gathers now only around the jet exit, Mach disk and shear layer. The recirculation region is cleared of turbulent component. The reattachment point is under the influence of a turbulence contribution from the shear layer rather than that from recirculation zone at $Kn=0.01$. The HoT becomes more visible especially in the boundary layer.

Figure 12.

Same as figure 11, except that here $Kn=0.125$: (a) $q_{1,lam}$, (b) $q_{1,turb}$, (c) $q_{1,HoT}$, (d) $q_{2,lam}$, (e) $q_{2,turb}$ and (f) $q_{2,HoT}$.

4.4. Changes in surface quantities

Surface pressure and heat flux are of great engineering significance, and the differences between the GSIS and GSIS-SST are analysed here. Figure 13(a) shows that, when the Knudsen number is fixed, the surface pressure is small around the jet nozzle, grows through the barrel shock and drops in the recirculation region, increases again to the peak value at the reattachment point and decreases almost monotonically after the peak. When the Knudsen number increases from 0.01 to 0.125, the shape of the pressure profile barely changes, but its magnitude increases. Comparison between the GSIS and GSIS-SST results show that the turbulence model increases the surface pressure, e.g. surface pressure from GSIS-SST is higher than that of GSIS by 10 % and 7 % at $Kn=0.01$ and 0.125, respectively.

Figure 13.

(a) Surface pressure. Results from GSIS and GSIS-SST remain constant after $x_1=0.02$ m and, hence, are not shown. (b) Surface heat flux.

Figure 13(b) shows the distributions of surface heat flux $q$. When $Kn=0.01$, $q$ obtained from the GSIS is small in the recirculation region near the nozzle exit, and increases to peak value at the reattachment point (around $x_1=-0.008$ m). After the peak, $q$ drops fast along the surface of the blunt leading edge, and then stays nearly constant along the flat plate. When $Kn=0.125$, the variation of $q$ curve, obtained from the GSIS, is similar to that at $Kn=0.01$, except that when $x_1>0$ m the surface heat flux from GSIS continues to drop, first sensible and then stays around a constant (from $x_1\approx 0.05$ m) when approaching the end of the model.

When the turbulence is considered, large deviation in surface heat flux between the GSIS and GSIS-SST could be spotted on the blunt part. At $Kn=0.01$, the GSIS-SST curve shows >20 % higher reattachment peak heat flux than the GSIS around $x_1=-0.01$ m. On the flat part of the model with $x_1>0.02$ m, GSIS-SST predicts a increasingly higher heat flux of 20–70 % as $x_1$ increases. At $Kn=0.125$, the GSIS-SST presents a 12 % higher reattachment peak heat flux and it greatly deviates from the GSIS one when $x_1>0$ m. The $q$ curve of GSIS-SST climbs almost linearly with $x_1$ on the flat part of the model, and resulting in a maximum relative deviation as high as 370 %. This implies a large degradation of the thermal protection effect, e.g. figure 6 shows that the low-temperature jet flow layer is much more limited in GSIS-SST than GSIS. With more dispersed jet flow layer and hotter environment at $Kn=0.125$, the rest of the model (from $x_1=0$ m) in GSIS-SST is thus heated quickly.