Multiscale modelling and simulation of coexisting turbulent and rarefied gas flows

Abstract

Simulating complex gas flows from turbulent to rarefied regimes is a long-standing challenge, since turbulence and rarefied flow represent contrasting extremes of computational aerodynamics. We propose a multiscale method to bridge this gap. Our method builds upon the general synthetic iterative scheme for the mesoscopic Boltzmann equation, and integrates the $k$–$\omega$ model in the macroscopic synthetic equation to address turbulent effects. Asymptotic analysis and numerical simulations show that the macroscopic–mesoscopic coupling adaptively selects the turbulence model and the laminar Boltzmann equation. The multiscale method is then applied to opposing jet problems in hypersonic flight surrounding by rarefied gas flows, showing that the turbulence could cause significant effects on the surface heat flux, which cannot be captured by the turbulent model nor the laminar Boltzmann solution alone. This study provides a viable framework for advancing understanding of the interaction between turbulent and rarefied gas flows.

1. Introduction

Numerical simulation of aerodynamics has two distinct branches. When the Reynolds number ($Re={\rho u L}/{\mu }$) is large, the Navier–Stokes (NS) equations are used to study turbulence. When the Knudsen number ($Kn$, the ratio of the molecular mean free path $\lambda$ to the characteristic flow length $L$) is large, the Boltzmann equation is employed to study rarefied flows. The Reynolds number and the Knudsen number are connected through the Mach number ($Ma={u}/{\sqrt {\gamma RT}}$) as

(1.1)\begin{equation} Re =\sqrt{\frac{{\rm \pi}\gamma}{2}} \frac{Ma}{Kn}, \end{equation}

where $u$ is the flow velocity, $\gamma$, $R$, $\mu$, $p$ and $\rho$ are the specific heat ratio, gas constant, shear viscosity, pressure and density of the gas, respectively. Therefore, it is generally believed that the turbulent and rarefied flows do not coexist (Tennekes & Lumley 1972).

From the physical perspective, the Boltzmann equation is elemental, as the NS equations can be derived from it when $Kn$ is low (Chapman & Cowling 1990). When $Kn$ increases, the gas flow gradually falls into the transition and free-molecular flow regimes, where the constitutive relations in the NS equations become inapplicable, necessitating the use of the Boltzmann equation. However, from the numerical perspective, since the velocity distribution function (VDF) is defined in the six-dimensional phase space, the Boltzmann equation, which is solved by the stochastic direct simulation Monte Carlo (DSMC) method (Bird 1994) or deterministic discrete velocity method (Aristov 2001), is typically applied in laminar rarefied flows.

There are many length scales in a turbulent flow. If $L$ is chosen to be the smallest Kolmogorov length scale, the local Knudsen number can be large, and the NS equations may fail to describe the multiscale turbulence. Tennekes & Lumley (1972) estimated the turbulent Mach number and Reynolds number based on fluctuation velocity and Kolmogorov length scale, and suggested that except under extreme conditions such as large-scale hydrogen nebula, the smallest eddies in turbulence still run in a continuum. Due to the growth in computational capabilities, many researchers have been trying to find the rarefaction effects in turbulent flows. For instance, Stefanov & Cercignani (1993) used the DSMC method to simulate the two-dimensional Taylor–Couette flow of a rarefied gas, and found that ‘the gas clearly exhibits an instability for a Taylor number of the order of the critical value for fluids governed by the NS equations; beyond that value, a new instability and a transition to chaotic dynamics is witnessed’. Komatsu et al. (2014) conducted molecular dynamics simulations of freely decaying turbulence, but ‘the energy spectrum is observed to scale reasonably well according to Kolmogorov scaling, even though the Kolmogorov length is of the order of the molecular scale’. Gallis et al. (2017) simulated the Taylor–Green vortex flow by the DSMC method, and found it ‘reproduces the Kolmogorov law and agrees well with the turbulent kinetic energy and energy dissipation rate obtained from direct numerical simulation (DNS) of the NS equations’. Li (2018) simulated the decaying isotropic turbulence by solving the Boltzmann kinetic equation deterministically, and found that the kinetic equation predicts the energy-decaying exponent about 10 % larger than the NS equations. However, the difference in statistical quantities is only a few per cent. McMullen et al. (2022) discovered that the sole divergence in the energy spectrum from the DSMC and NS simulation results occurs at large wave numbers, where thermal fluctuations induce the energy spectrum to scale proportionally to the square of wave number. (This is also observed in the molecular dynamics simulation by Komatsu et al. (2014), and can be predicted by the fluctuating NS equations (Landau & Lifshitz 1959). Therefore, these results cannot be attributed to the rarefaction effects in turbulence.) They also reported that ‘thermal fluctuations have little impact on the large-scale evolution of the flow’.

To date, no evidence of coexisting turbulence and rarefied flows with engineering significance have been identified. There might be two reasons for the absence of such cases. First, the DNS based on the Boltzmann equation is computationally intensive, which limits the range of flow scenarios and parameters explored. Second, in previous studies the rarefaction effects are searched on the Kolmogorov scale. Even if such rarefaction effects do exist, their effect on the overall macroscopic flow phenomena is limited since they occur in the smallest scale. Therefore, in this paper, we are not going to search for the rarefaction effects inside the bulk turbulent flow, but instead study a turbulent flow inside a bulk rarefied flow. In such a scenario, the intense turbulent flow is confined to a relatively small region, whereas the weaker rarefied flows extend over a much larger area. This spatial distribution ensures that their respective powers or influences are on par with each other. Consequently, their interaction could potentially result in substantial deviations from the expected behaviour of isolated turbulent and rarefied flows in simulations. Such scenarios emerge with the rapid advancement of trans-atmospheric hypersonic vehicle technology, where new challenges in fluid dynamics have arisen, particularly in active flow control of aerodynamic force and heat. In higher-altitude and higher-Mach-number rarefied flow, employing a high-density jet may be the only feasible approach for thermal management and flight control. However, there is a current lack of understanding of the interactions between high-Mach-number rarefied flows and high-density turbulence.

This paper is dedicated to developing a multiscale method to simulate the gas flow spanning from the turbulent to rarefied conditions, and to identifying engineering scenarios where turbulence and rarefaction coexist and interact within a single flow configuration.

2. A multiscale model from the turbulent to rarefied flows

Although the Boltzmann equation encapsulates the gas dynamics across turbulent to rarefied regimes, the computational demands of DNS are overwhelming. In regions of high $Re$, the spatial grid must be finely resolved, whereas in regions of high Knudsen number $Kn$, the discrete velocity grid requires dense sampling. Given that the DNS of NS equations is already a time-intensive process for simulating turbulence, the prospect of performing the DNS of the Boltzmann equation in scenarios where turbulence and rarefaction coexist is beyond the reach of current computational capability. To address this challenge, we introduce a multiscale modelling approach that offers a viable alternative. Our idea is to couple the Reynolds-averaged Navier–Stokes (RANS) turbulence model into the Boltzmann equation, to ease the requirement of fine spatial grids in turbulent regime. The coupling should process the asymptotic-preserving property, i.e. it automatically recovers the RANS model when $Re$ is high, and the Boltzmann equation when $Kn$ is large. Moreover, from the numerical perspective, since the RANS produces the time-averaged steady state, the Boltzmann equation will also be solved by the general synthetic iteration scheme (GSIS) tailored to efficient steady-state simulations (Su et al. 2020a; Zhang et al. 2024).

In the following, we first introduce the Boltzmann kinetic model which describe the gas dynamics from the continuum to rarefied flow regimes. Second, we introduce the two-temperature NS equations derived from the kinetic equation when the Knudsen number is small, and describe the RANS model for efficient simulation of turbulence on a coarse spatial grid. Third, we couple the GSIS and RANS to efficiently find the solution of the kinetic equation from the turbulent to rarefied regimes.

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 (1951) equation. To reduce the computational cost, simplified kinetic equations are usually adopted (Wu et al. 2015; Li et al. 2021). 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. 2023). 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 (1951) 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. 2021),

(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 1990):

(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. 2021)

(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 1994). 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 (1994))

(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. 2024), 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 1993), 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 (1992), 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 (2006).

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. 2020a; Zhang et al. 2024). 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. 2020a; Liu et al. 2024a; Zhang et al. 2024). Furthermore, according to the Chapman–Enskog expansion, in the continuum limit, the HoTs are proportional to $Kn^2$ (Su, Zhu & Wu 2020b), 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. (2024a). 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. Validation of the GSIS-SST model

In this section, we validate the GSIS-SST model deep into the turbulent and rarefied flow regimes, respectively. In addition, we show that the GSIS-SST provides reasonable modelling of the crossover between turbulent and rarefied conditions.

3.1. High Reynolds number

According to the Chapman–Enskog expansion (Su et al. 2020b; Su, Zhang & Wu 2021), 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 2013). 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. Opposing jet in a rarefied environment

We now employ the proposed GSIS-SST to identify engineering scenarios where turbulence and rarefaction phenomena coexist. Since the flow quantities do not change significantly (except in the shock region) in § 3.3, the coexisting turbulent and rarefied flows only occurs around the blunt edge. This is probably also the reason why rarefaction effects are not observed with significant difference from the NS equations (Komatsu et al. 2014; Li 2018; McMullen et al. 2022). We believe a strong evidence of the coexisting turbulent and rarefied flows is most likely to arise during trans-atmospheric hypersonic flight, when the active jet technology is used to reduce the drag and surface heat flux (Hayashi & Aso 2003; Daso et al. 2009). In this case, the strong turbulent flow and the weak rarefied flow occupies a small and large area, respectively, so that their powers/influences are comparable, and their interaction might lead to significant deviations from the simulation for the turbulent flow and rarefied flow alone.

Numerical research on opposing jets is abundant, where turbulence models are always present, especially the $k$–$\omega$ SST model (Ou et al. 2018). On the other hand, kinetic solvers are used to study the rarefied jet flows, e.g. the opposing jet at an altitude of over 70 km is simulated by DSMC (Raeisi, Mohammadi-Amin & Zakeri 2019; Guo, Luo & Wu 2024), and the microscale jets into quiescent space and vacuum plume are solved by the hybrid NS-DSMC solver (Virgile, Albert & Julien 2022; Liu, Yue & Lin 2024b), where only laminar flows are considered. Here, we apply the GSIS-SST model to investigate the interaction of turbulent and rarefied flows for the first time.

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 2023)):

(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.

5. Conclusions and outlook

A multiscale method, which is built upon the Boltzmann kinetic equation to describe the rarefied gas dynamics and the $k$–$\omega$ SST model to capture the turbulence effect, has been established to simulate gas flows from highly rarefied regime down to the fully turbulent continuum regime. Asymptotic analysis and numerical simulations have shown that the GSIS-SST adaptively selects the SST model, the laminar NS equations and the laminar Boltzmann equation. It should be noted that, while GSIS exactly solves the Boltzmann kinetic equation for laminar rarefied gas flows, the SST is an approximate model for turbulence viscosity. Therefore, the overall accuracy is affected by the SST model. Fortunately, the turbulence models are well studied, where abundant models/parameters can be optimised for specific problems.

With the proposed GSIS-SST model, the interaction between turbulent and rarefied gas flows has been analysed. That is, for opposing jet problem under hypersonic rarefied flows, the turbulence effect predicted by the SST model promotes the momentum and energy exchange between the jet and the free flow, considerably enhancing the jet diffusion and reducing the heat reduction range of the jet. For example, the surface heat flux predicted by the GSIS-SST may be approximately four times higher than that of the GSIS. The predicted turbulence also increases reattachment surface pressure and heat flux significantly, and such differences increase when the $Kn$ number decreases. These findings underscore the necessity of considering the turbulence effect of the jet even under rarefied free flow conditions, which is critical for optimising thermal protection systems and aerodynamic control strategies in hypersonic vehicles.

The significance of this study is twofold. First, despite many researchers’ unsuccessful attempts to find the rarefaction effects in turbulence, our research has presented, for the first time, numerical evidence of coexisting turbulent and rarefied gas flows in a single flow configuration. These flows defy characterisation by conventional turbulence models and cannot be accurately represented by laminar Boltzmann solutions. Second, this study has provided a viable framework for advancing our understanding of the interaction between turbulent and rarefied gas flows, based on which more new flow phenomena can be explored.

While the current model provides valuable insights into the interaction between turbulent and rarefied gas flows, there are also opportunities for further refinement. Although we have assumed that the jet is fully turbulent before interacting with the laminar rarefied flow, we acknowledge that incorporating a transition model could improve the spatial accuracy of predicting where the GSIS and GSIS-SST solutions diverge. However, based on our analysis, we anticipate that although this addition might enhance spatial predictions, its effect on the magnitude of surface quantities, such as pressure or heat flux, would likely be modest, as the turbulent mixing primarily affects the higher-density jet flow, leading to its rapid dispersion in the GSIS-SST solution, while the influence of the low-density free flow on surface quantities remains limited. Nevertheless, the observed significant differences suggest that the current model adequately capture key flow behaviours, even in the absence of a transition model. Moving forward, we plan to integrate more advanced, unsteady, high-fidelity turbulence modelling approaches with the ability to predict transition, to further enhance the accuracy and robustness of the GSIS-SST in capturing complex flow dynamics. In addition to its application in hypersonic aerodynamics, the GSIS-SST may be applied in inertial confinement fusion to understand the interaction of turbulence in low-temperature high-density region and rarefied flows in high-temperature region with large mean free path of electron/ions (Rinderknecht et al. 2018).