Numerical Study on Remote Infrared Signal Fluctuation Properties of Gas Free Jet Based on Improved SGS-TRI Model

1. Introduction

The phenomenon of turbulence-radiation interaction (TRI) is pronounced in high-temperature turbulent gas jet. Turbulence induces fluctuations in temperature, composition and pressure of the non-uniform gas jet, which subsequently causes the fluctuations of local radiation absorption coefficient and the Planck function for radiation emission. These result in variations in internal radiation transmission and can even affect the fluctuations of radiation heat transfer, thereby influencing the turbulent flow of gas jet [1]. Notably, the gas radiation absorption coefficient and the Planck function for radiation emission are nonlinear functions of temperature. This nonlinearity leads to a significant disparity between the time-averaged characteristics of transient radiation of turbulent gas and the radiation characteristics of its time-averaged flow field. In many scientific contexts, particularly in engineering applications, the former is often more relevant. For instance, a collaborative study by Modest and Haworth examined the TRI characteristics of methane-air axisymmetric diffusion flames [1]. Their findings indicated a 150 K difference in the maximum flame temperature when TRI was considered compared to when it was not, as well as a sixfold variation in nitric oxide production. Furthermore, as current optical characteristic recognition algorithms evolve towards three-dimensional space-time domain, the fluctuating characteristics of infrared radiation from gas jets, driven by turbulence [2], are anticipated to emerge as significant target features for hydro-carbon fuel aircraft.

Most research of TRI has been conducted in the field of combustion. As early as the 1980s, the probability density function (PDF) method based on Reynolds-averaged Navier-Stokes equations (RANS) emerged [3,4]. The Snegirev method [5] emerged at the beginning of this century. Its fundamental concept involves expressing the fluctuation characteristics of the emission term in the radiative transfer equation (RTE) as a function of the local temperature fluctuation variance. Subsequently, this method allows for the determination of temperature fluctuation within the RANS flow field by establishing a specialized transport equation for variance distribution. The aforementioned TRI characteristic calculation method based on RANS is still under development [6,7]. However, due to the inability to model the fluctuation of absorption term of the radiative transfer equation, this method can only be used in the scenario meeting the optically thin fluctuation approximation (OTFA) [7], which restricts its applicability in many engineering contexts.

This article focuses on the gas-free jet, where the turbulent fluctuation scale at a distance from the nozzle outlet is comparable to the nozzle diameter, which does not conform to the OTFA assumption. The practical calculation of its TRI characteristics can only be addressed through the radiative transfer equation of the transient turbulent flow field, derived from large eddy simulations (LES). Blunck et al. conducted a numerical study on the TRI characteristics of the jet in the 4.3micron band for an exhaust system equipped with a pre-swirl device [8]. They discovered that turbulent fluctuation minimally affects the radiation emitted from the core area of the jet; however, at distances from the nozzle exit or in proximity to the jet mixing layer, turbulent fluctuation significantly impacts the time-averaged radiation. The fluctuating radiation characteristics of the gas jet exhibit many similarities to those of flames and when the Reynolds number is sufficiently high, these characteristics are not sensitive to variations in the Reynolds number. Additionally, computational study of Damien et al on the TRI characteristics of radiative heat transfer in a rectangular combustion chamber with a triangular flame stabilizer [9] revealed that the effects of temperature autocorrelation fluctuation, temperature and absorption coefficient correlation fluctuation and absorption coefficient autocorrelation fluctuation on flame emission radiation fluctuation are all substantial and cannot be overlooked. Notably, none of the aforementioned LES-based calculations account for the influence of sub-grid scale turbulence fluctuations on TRI characteristics. The Coelho team employed direct numerical simulation (DNS), ray tracing and narrow-band k-distribution models to evaluate the SGS-TRI characteristics of isotropic turbulence and free jets [10,11]. Their findings indicate that the impact of SGS-TRI on the emission term of the radiative transfer equation is significantly greater than that on the absorption term. Furthermore, they observed that as optical thickness increases, temperature fluctuations intensify, leading to a higher proportion of SGS-TRI in the overall radiation; specifically, when the optical thickness is 1, the proportion of SGS-TRI ranges from 13% to 25%. Additionally, the influence of SGS-TRI on the radiative transfer of free jets is primarily concentrated at the edge of jet, where the overall radiant brightness is minimal. Notably, the variation in radiant brightness calculations due to the consideration of sub-grid scale fluctuations in radiative transfer can exceed one-third. In the context of flame radiation heat transfer, the filtered density function (FDF) method [12], which is based on the probability density function (PDF) method, is pre-dominantly used to assess SGS-TRI characteristics. The results indicate that SGS-TRI is significantly influenced by factors related to sub-grid scale chemical reactions and temperature fluctuations. In flames with large optical thickness, the contribution of SGS-TRI to radiant heat is considerably greater than that observed in gas jets, it can even surpass the contribution from turbulent fluctuations at a resolvable scale [13,14,15].

As the significance of SGS-TRI became increasingly apparent, Poitou et al. proposed a modeling method in 2007 [16]. This method extends the previously mentioned Snegirev method into large eddy simulations by using the same expression of emission term fluctuation in the RTE against gas temperature fluctuation, combined with a sub-grid scalar variance model [17]. The Coelho team assessed the accuracy of this modeling method in gas radiative heat transfer calculation for a simple axisymmetric flame, using the FDF model as benchmark [18]. The results indicate that the existing sub-grid Snegirev expression significantly underestimates the emission term fluctuation in flame RTE. Additionally, the fixed coefficient sub-grid scalar variance model considerably underestimates the sub-grid temperature fluctuations. Consequently, the combined effect of these two factors leads to the predictions of the SGS-TRI characteristics of diffusion flames being considerably lower than the actual measurements.

Current research on TRI primarily examines its influence on narrow-band infrared signals used in non-contact measurement for flow field and on full-spectrum radiative heat transfer, which have different characteristics [19]. However, no published research has addressed the TRI and SGS-TRI characteristics of gas jet as well as their modeling methods for the atmospheric infrared window wide-bands, which are pertinent to the field of remote sensing (e.g., 3-5 micron and 8-14 micron band). This paper employs a Multi-Scale Multi-Group Wide-Band k-distribution model [20,21] (Section 2.1) to calculate the TRI characteristics for remote infrared signals of axisymmetric gas jets in the 3-5 micron wave band. It also demonstrates that the model can effectively perform radiative transfer calculations in scenarios where the correlation-k characteristic of the gas absorption spectrum is notably disrupted due to substantial differences in temperature and composition between the combustion gas and the ambient air as well as strong thermodynamic state parameter inhomogeneities in their mixing layer (Section 3.3).To address the low prediction accuracy of existing SGS-TRI models, this study introduces a term describing correlation between combustion gas species concentration and temperature fluctuations into the sub-grid Snegirev relationship (Section 2.2). Furthermore, a new ensemble average calculation method, which is more applicable in complex engineering contexts, is developed to resolve the coefficients of dynamic Smagorinsky model and temperature variance model (Section 2.3). The validity of these algorithms is confirmed through DNS calculations of isotropic turbulence (Section 3.1) and comparison with measurement data of high-temperature air jet flow field (Section 3.2). Based on these, the TRI and SGS-TRI characteristics of 3-5micron band infrared signal from an axisymmetric combustion gas jet at different detection distances are analyzed (Section 4.1), and the characteristics of the MSMGWB model are employed to elucidate the formation mechanism of the observed phenomena (Section 4.2).

2. Radiation Transfer Equation Based on k-Distribution Model for Solvable Scale Turbulent Flow Field

2.1. MSMGWB Model

For the 3-5 micron band infrared radiation emitted by combustion gases without particles and attenuated by ambient air before reaching optical detectors, the RTE (without considering gains from ambient air) can be written as [20,21]

$$\left\{\begin{array}{l}{\displaystyle \frac{d{I}_{\eta }}{ds}}=k_{\eta }(\underset{\_}{\varphi })I_{b\eta }(T)-\left(k_{\eta }(\underset{\_}{\varphi })+k_{s\eta }(\underset{\_}{\varphi })\right)I_{\eta }\\ I_{3-5\mu m}={\displaystyle {\int }_{2000c{m}^{-1}}^{3333c{m}^{-1}}{I}_{\eta }d\eta }\end{array}\right.$$

(1)

where, $\eta$ is the wave number, $k_{\eta }$ and $k_{s\eta }$ are the spectral extinction coefficient and the spectral scattering coefficient at wave number $\eta$, respectively. $\underset{\_}{\varphi }=(T,P,x)$ is the thermodynamic state parameters of the mixed medium that affect$k_{\eta }$ and$k_{s\eta }$ (comprising gas species mole fraction $x$ , pressure $p$ and temperature $T$), $I_{\eta }$ the spectral radiance, $I_{b\eta }(T)$ the blackbody’s spectral radiance at $T$ and $s$ the location along the radiation transfer path. In this work, $k_{\eta }={\displaystyle \sum _{n=1}^4{k}_{n\eta }}$, $k_{s\eta }=k_{4s\eta }$where $k_{n\eta }$ is the spectral absorption coefficient of n-th species, and subscript $n=1~4$ denote water vapor, carbon dioxide, carbon monoxide and aerosol particles, respectively. The values of ${\left.k_{n\eta }\right|}_{n=1~3}$ were obtained from the LBL approach based on HITEMP database (2010 version) [22], and the values of $k_{4\eta }$ and $k_{4s\eta }$ were calculated using Mie scattering theory [23,24]. All the parameters will be based on the local $\underset{\_}{\varphi }$ unless otherwise indicated.

In hot combustion gas remote detection scenario, water vapor and carbon dioxide are treated as participation in both emission and absorption process, while carbon monoxide and aerosols only contribute to the attenuation part, because they exist only in the atmosphere rather than in turbo gas. Besides, the emission and scattering gain from atmosphere are not considered here since they are commonly treated as background radiation and calculated separately. Therefore, based on its linearity, Equation (1) can be rewritten as [19,21,25],

$$\left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{I}_{n\eta }}{\mathrm{d}s}}=k_{n\eta }{I}_{\mathrm{b}\eta }(T)-k_{\eta }{I}_{n\eta };\\ I_{\eta }={\displaystyle \sum _{n=1}^2{I}_{n\eta }};\\ k_{3\eta }=k_{4\eta }=k_{4\mathrm{s}\eta }=0\begin{array}{cc}& in combustion gas;\end{array}\\ I_{\mathrm{b}\eta }=0\begin{array}{cc}\begin{array}{cc}& \begin{array}{cc}& \end{array}\end{array}& in atmosphere\end{array}.\end{array}\right.$$

(2)

where $I_{n\eta }$ results from the emission of H₂O or CO₂, but is subject to the absorption of both in combustion gases and the absorption and scattering of all four species in atmosphere.

The MSMGWB model employs a method of grouping spectral absorption coefficients to address the significant differences in temperature and composition between combustion gas and ambient air. These differences result in a lack of correlated-k characteristics in absorption spectra, which significantly reduces the calculation accuracy of conventional k-distribution models. The principle is to allocate the spectral subintervals, which are dominated by 'cold absorption lines' and 'hot absorption lines', into different subset, ensuring spectral absorption coefficients within each subset have similar thermodynamic state parameter dependences, i.e., exhibit correlated-k characteristics under the concerned gas thermodynamic state parameter combinations, as illustrated in [20].

By performing independent k-distribution transformations for each spectral subset, the final radiative transfer equation of the MSMGWB model is [20],

$$\begin{array}{c}\left\{\begin{array}{c}{\displaystyle \frac{d{I}_{\mathrm{j}}}{ds}}={\dot{k}}_{\mathrm{j}}^{*}(\varphi )A_{\mathrm{j}}(\varphi )I_{b,2000c{m}^{-1}-3333c{m}^{-1}}(T)-{\dot{k}}_{\mathrm{j}}(\varphi )I_{\mathrm{j}}\\ I_{3-5\mu m}={\displaystyle {\sum }_{\mathrm{j}=1}^{64}{I}_{\mathrm{j}}}\end{array}\right.\end{array}$$

(3)

where, ${\dot{k}}_{\mathrm{j}}^{*}=\rho x{k}_{\mathrm{j}}^{*}$ , ${\dot{k}}_{\mathrm{j}}=\rho k_{\mathrm{j}}$, $A_{\mathrm{j}}$ and $I_{b,2000c{m}^{-1}-3333c{m}^{-1}}$ represent the emission coefficient, absorption coefficient, stretching coefficient and blackbody radiance, respectively, corresponding to the jth radiative transfer equation of MSMGWB model.

The reliability of MSMGWB model in remote sensing scenario is preliminary validated through two types of 1D cases, as sketched in Figure 1, using calculation results of line-by-line (LBL) approach as benchmark. The results shown in Figure 2 and Figure 3. indicate that MSMGWB model can almost perfectly deal with radiation transfer in gases with temperature, participating component mole ratio, as well as pressure inhomogeneities and v aerosols with various sizes and components. Where radiation characteristics of various aerosols were calculated by Mie theory [26,27]. More calculation results of the MSMGWB model in the 1D cases can be found in reference [21].

2.2. Dynamic Smagorinsky Model Based on Relaxation Factor

The governing equations of the solvable scale flow field can be derived by filtering the NS equation, so the governing equations of large eddy simulation can be obtained as follows,

$$\begin{array}{c}{\displaystyle \frac{\partial \overline{\rho }}{\partial t}}+{\displaystyle \frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial x_i}}=0\end{array}$$

(4)

$$\begin{array}{c}{\displaystyle \frac{\partial \left(\overline{\rho }{\tilde{u}}_i\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\overline{\rho }{\tilde{u}}_i{\tilde{u}}_j\right)}{\partial x_j}}=-{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+{\displaystyle \frac{\partial {\tilde{\sigma }}_{ij}}{\partial x_j}}+{\displaystyle \frac{\partial {\tau }_{ij}}{\partial x_j}}\end{array}$$

(5)

$$\begin{array}{c}{\displaystyle \frac{\partial \left(\overline{\rho }\tilde{E}\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\overline{\rho }\tilde{H}{\tilde{u}}_j\right)}{\partial x_j}}={\displaystyle \frac{\partial \left({\tilde{u}}_i{\tilde{\sigma }}_{ij}\right)}{\partial x_j}}+{\displaystyle \frac{\partial \left(\frac{\mu c_p}{\mathrm{Pr}}\frac{\partial \tilde{T}}{\partial x_j}\right)}{\partial x_j}}+{\displaystyle \frac{\partial \left(C_p{Q}_j\right)}{\partial x_j}}+{\displaystyle \frac{\partial {\Pi }_j}{\partial x_j}}+{\displaystyle \frac{\partial D_j^H}{\partial x_j}}\end{array}$$

(6)

$$\begin{array}{c}{\displaystyle \frac{\partial \left(\overline{\rho }{\tilde{y}}_{O_2/H_{2}O/C{O}_2/CO}\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\overline{\rho }{\tilde{u}}_j{\tilde{y}}_{O_2/H_{2}O/C{O}_2/CO}\right)}{\partial x_j}}={\displaystyle \frac{\partial {\Gamma }_{O_2/H_{2}O/C{O}_2/CO,j}}{\partial x_j}}\end{array}$$

(7)

where, $\rho$,$u_i$,$p$,$y$,$E$,$H$and$\tau$are represent density, velocity, pressure, mass fraction of each component, total energy per unit mass of fluid, total enthalpy per unit mass of fluid, and subgrid-scale stress, respectively.

Where,

$$\begin{array}{c}{\tilde{\sigma }}_{ij}=2\mu {\tilde{S}}_{ij}\end{array}$$

(8)

$$\begin{array}{c}{\tilde{S}}_{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial {\tilde{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\tilde{u}}_j}{\partial x_i}}\right)-{\displaystyle \frac{1}{3}}{\displaystyle \frac{\partial {\tilde{u}}_k}{\partial x_k}}{\delta }_{ij}\end{array}$$

(9)

$$\begin{array}{c}{\tau }_{ij}=\overline{\rho }\left({\tilde{u}}_i{\tilde{u}}_j-\widetilde{u_i{u}_j}\right)\end{array}$$

(10)

$$\begin{array}{c}\overline{\rho }\tilde{E}=\overline{\rho }\tilde{H}-\overline{p}=\overline{\rho }\tilde{h}+{\displaystyle \frac{\overline{\rho }{\tilde{u}}_i{\tilde{u}}_i}{2}}-{\displaystyle \frac{{\tau }_{kk}}{2}}\end{array}-\overline{p}$$

(11)

$$\left\{\begin{array}{c}\overline{p}=\overline{\rho }R\tilde{T}\\ \tilde{h}=C_p\tilde{T}={\displaystyle \frac{\gamma }{\gamma -1}}R\tilde{T}\end{array}\right.$$

(12)

$$\begin{array}{c}{Q}_j=\overline{\rho }\left({\tilde{u}}_j\tilde{T}-\widetilde{u_{j}T}\right)\end{array}$$

(13)

$$\begin{array}{c}{\Pi }_j={\displaystyle \frac{\overline{\rho }\left({\tilde{u}}_j\widetilde{u_k{u}_k}-\widetilde{u_j{u}_k{u}_k}\right)}{2}}\end{array}$$

(14)

$$\begin{array}{c}{D}_j^H=-({\Gamma }_{j,O_2}+{\Gamma }_{j,H_{2}O}+{\Gamma }_{j,C{O}_2}+{\Gamma }_{j,CO}){\tilde{h}}_{N_2}+{\Gamma }_{j,O_2}{\tilde{h}}_{O_2}+{\Gamma }_{j,H_{2}O}{\tilde{h}}_{H_{2}O}+{\Gamma }_{j,C{O}_2}{\tilde{h}}_{C{O}_2}+{\Gamma }_{j,CO}{\tilde{h}}_{CO}\end{array}$$

(15)

$$\begin{array}{c}{\Gamma }_j\begin{array}{c}=\overline{\rho }\end{array}D{\displaystyle \frac{\partial \tilde{y}}{\partial x_j}}+{\Theta }_j=\overline{\rho }\end{array}D{\displaystyle \frac{\partial \tilde{y}}{\partial x_j}}+\overline{\rho }\left({\tilde{u}}_j\tilde{y}-\widetilde{u_{j}y}\right)$$

(16)

where, $\gamma$,$R$,$C_p$and$D$are represent specific heat ratio, universal gas const, specific heat at constant pressure, and diffusion coefficient, respectively.

In equations (4)-(16), the superscript "-" denotes the volume-averaged quantity within a cube with characteristic side length$\overline{\Delta }$ , termed the "filtering length". The superscript "~" signifies the Favre average. In equations (4)-(7), there exist unclosed terms, i.e. the sub-grid terms, that require modeling for closure. This paper employs the Smagorinsky sub-grid model:

$$\begin{array}{c}{\tau }_{ij}=\overline{\rho }\left({\tilde{u}}_i{\tilde{u}}_j-\widetilde{u_i{u}_j}\right)=2\overline{\rho }{C}_s^2{\overline{\Delta }}^2\sqrt{2{\tilde{S}}_{ij}{\tilde{S}}_{ij}}{\tilde{S}}_{ij}\end{array}$$

(17)

$$\begin{array}{c}{Q}_j=\overline{\rho }\left({\tilde{u}}_j\tilde{T}-\widetilde{u_{j}T}\right)=\end{array}\overline{\rho }{C}_q{\overline{\Delta }}^2\sqrt{2{\tilde{S}}_{ij}{\tilde{S}}_{ij}}{\displaystyle \frac{\partial T}{\partial x_j}}$$

(18)

$$\begin{array}{c}{\Pi }_j={\displaystyle \frac{\overline{\rho }\left({\tilde{u}}_j\widetilde{u_k{u}_k}-\widetilde{u_j{u}_k{u}_k}\right)}{2}}={\tilde{u}}_k\end{array}{\tau }_{jk}$$

(19)

$${\Theta }_j\begin{array}{c}=\end{array}\overline{\rho }\left({\tilde{u}}_j\tilde{y}-\widetilde{u_{j}y}\right)=\overline{\rho }{C}_d{\overline{\Delta }}^2\sqrt{2{\tilde{S}}_{ij}{\tilde{S}}_{ij}}{\displaystyle \frac{\partial \tilde{y}}{\partial x_j}}$$

(20)

When applying dynamic Smagorinsky model of Lilly [28], the calculation method for $C_s$ is as follows,

$$\begin{array}{c}{C}_s^2={\displaystyle \frac{〈M_{ij}{L}_{ij}〉}{〈M_{ij}{M}_{ij}〉}}\end{array}$$

(21)

$$\begin{array}{c}{L}_{ij}=\widehat{\rho }{\widehat{u}}_i{\widehat{u}}_j-\widehat{\overline{\rho }{\tilde{u}}_i{\tilde{u}}_j}\end{array}$$

(22)

$$\begin{array}{c}{M}_{ij}=2\widehat{\rho }{\widehat{\Delta }}^2\sqrt{2{\stackrel{⌣}{S}}_{ij}{\stackrel{⌣}{S}}_{ij}}{\stackrel{⌣}{S}}_{ij}-\widehat{2\overline{\rho }{\overline{\Delta }}^2\sqrt{2{\tilde{S}}_{ij}{\tilde{S}}_{ij}}{\tilde{S}}_{ij}}\end{array}$$

(23)

Here, $\widehat{\Delta }$ denotes a filtering scale larger than $\overline{\Delta }$. The superscript '^' signifies the volume-averaged value of this parameter within a cube with characteristic side length $\widehat{\Delta }$, while the superscript $\cup$ indicates the Favre average at this scale.

In equation (21),

denotes the ensemble average. For two-dimensional flows, the ensemble average of equation (21) can be spatially averaged in statistically homogeneous directions. However, in complex three-dimensional flows, there are no statistically homogeneous directions, preventing the use of spatial averaging to calculate dynamic eddy viscosity coefficients. The direct temporal averaging method proposed by Meveneau is computationally expensive due to excessive memory requirements; whereas the fluid particle averaging trajectory method necessitates solving additional transport equations, which is complicated and has low computational stability [29]. To enhance the applicability of both the dynamic Smagorinsky model and temperature variance model in Large Eddy Simulation (LES) calculations for complex three-dimensional flow fields, this paper establishes an ensemble averaging calculation method based on relaxation factors ($\alpha$) as follows (using an arbitrary variable $f$ as an example):

$$〈f〉={\displaystyle \frac{{\displaystyle \sum _{n=-\infty }^0{\alpha }^{-n}{f}_n}}{{\displaystyle \sum _{n=-\infty }^0{\alpha }^{-n}}}}={\displaystyle \frac{f_0+\alpha {\displaystyle \sum _{n=-\infty }^0{\alpha }^{-n}{f}_{n-1}}}{{\displaystyle \sum _{n=-\infty }^0{\alpha }^{-n}}}}$$

(24)

here,$0<\alpha <1$, $n$ is serial number of calculation time step, and $n=0$ denotes current time step. This algorithm eliminates the need to solve additional transport equations and only requires the storage of values ${\displaystyle \sum _{n=-\infty }^0{\alpha }^{-n}{f}_n}$ on each computational grid, thus resulting in minimal memory requirements.

2.3. Improved SGS-TRI Model

Applying a box filter to equation (3) in the same manner, we obtain the following equation,

$$\begin{array}{c}{\displaystyle \frac{d{\overline{I}}_{\mathrm{j}}}{ds}}=\overline{\rho x{k}_{\mathrm{j}}^{*}(\varphi )A_{\mathrm{j}}(\varphi )I_{b,{\eta }_L-{\eta }_U}(T)}-\overline{\rho k_{\mathrm{j}}(\varphi )I_{\mathrm{j}}}\end{array}$$

(25)

Extensive research indicates [30,31] that when the filter scale $\overline{\Delta }$ is optical thin, the SGS-TRI characteristics of its absorption term are not pronounced. As a result, equation (25) can be reformulated as follows,

$$\begin{array}{c}{\displaystyle \frac{d{\overline{I}}_{\mathrm{j}}}{ds}}=\overline{\rho x{k}_{\mathrm{j}}^{*}(\varphi )A_{\mathrm{j}}(\varphi )I_{b,{\eta }_L-{\eta }_U}(T)}-\overline{\rho }{\tilde{k}}_{\mathrm{j}}(\varphi ){\overline{I}}_{\mathrm{j}}\end{array}$$

(26)

For the first term on the right-hand side of the equation (26), based on the method proposed by Snegrive [5], expanding it in the Taylor series and neglecting higher-order terms yields,

$$k_{\mathrm{j}}^{*}(\varphi )A_{\mathrm{j}}(\varphi )I_{b,{\eta }_L-{\eta }_U}(T)=k_{\mathrm{j}}^{*}({\varphi }_0)A_{\mathrm{j}}({\varphi }_0)I_{b,{\eta }_L-{\eta }_U}(T_0)+{\displaystyle \frac{\partial k_{\mathrm{j}}^{*}{A}_{\mathrm{j}}{I}_{b,{\eta }_L-{\eta }_U}}{\partial T}}(T-T_0)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2{k}_{\mathrm{j}}^{*}{A}_{\mathrm{j}}{I}_{b,{\eta }_L-{\eta }_U}}{\partial T^2}}{(T-T_0)}^2$$

(27)

By setting $T_0=\tilde{T}$ and ${\varphi }_0=\tilde{\varphi }$ in equations (27), we can obtain the filtered radiative transfer equation:

$$\begin{array}{c}{\displaystyle \frac{d{\overline{I}}_{\mathrm{j}}}{ds}}=\overline{\rho }\tilde{x}{k}_{\mathrm{j}}^{*}(\tilde{\varphi })A_{\mathrm{j}}(\tilde{\varphi })I_{b,{\eta }_L-{\eta }_U}(\tilde{T})+\left({\displaystyle \frac{\partial x}{\partial T}}{\displaystyle \frac{\partial k_{\mathrm{j}}^{*}{A}_{\mathrm{j}}{I}_{b,{\eta }_L-{\eta }_U}}{\partial T}}+{\displaystyle \frac{1}{2}}\tilde{x}{\displaystyle \frac{{\partial }^2{k}_{\mathrm{j}}^{*}{A}_{\mathrm{j}}{I}_{b,{\eta }_L-{\eta }_U}}{\partial T^2}}\right)\left(\overline{\rho T^2}-\overline{\rho }\tilde{T}\tilde{T}\right)-\overline{\rho }{\tilde{k}}_{\mathrm{j}}{\overline{I}}_{\mathrm{j}}\end{array}$$

(28)

where $\partial x/\partial T$and$\overline{\rho T^2}-\overline{\rho }\tilde{T}\tilde{T}$are terms to be modeled. The latter can refer to the subgrid-scale scalar variance model by Pierce and Moin [17]:

$$\overline{\rho T^2}-\overline{\rho }\tilde{T}\tilde{T}=C_T\overline{\rho }{\overline{\Delta }}^2{‖\nabla \tilde{T}‖}^2$$

(29)

the model coefficient $C_T$ can also be determined based on Lilly's dynamic model [28].

$$C_T=={\displaystyle \frac{〈G_A{G}_B〉}{〈G_B{G}_B〉}}$$

(30)

$$G_A=\widehat{\overline{\rho }\tilde{T}\tilde{T}}-\widehat{\rho }\stackrel{⌣}{T}\stackrel{⌣}{T}$$

(31)

$$G_B=\widehat{\rho }{\widehat{\Delta }}^2{‖\nabla \stackrel{⌣}{T}‖}^2-\widehat{\overline{\rho }{\overline{\Delta }}^2{‖\nabla \tilde{T}‖}^2}$$

(32)

The ensemble average in equation (30) is also calculated using the method proposed in section 2.1, namely equation (24). It is noteworthy that the impact of the term $\partial x/\partial T$ in equation (28) is neglected in the existing subgrid-scale Snegirev method [5]. However, turbulence theory and extensive DNS simulation results have confirmed a strong correlation between fluctuations in component concentration ($x$) and local temperature fluctuations within chemically frozen turbulent mixing layers [32] (see section 3.1.4). Furthermore, since ${\dot{k}}_{\mathrm{j}}^{*}=\rho x{k}_{\mathrm{j}}^{*}$, in most cases, the existing subgrid-scale Snegirev method underestimates the influence of subgrid-scale turbulent fluctuations on the emission terms of the RTE. The term $\partial x/\partial T$ in equation (28) can be estimated based on local components and temperatures at different times by fitting their slope using least squares. Following a principle similar to that used for equation (24), wherein higher weights are assigned to current time and lower weights to more distant time, we obtain,

$${\displaystyle \frac{\partial x}{\partial T}}={\displaystyle \frac{{\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{T}_i{x}_i-{\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{T}_i{\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{x}_i}{{\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{T}_i^2-{\left({\displaystyle \sum _{i=-\infty }^0{\alpha }^{-i}}{T}_i\right)}^2}}$$

(33)

3. Verification of Model Calculation Accuracy

This section, the reliability of the relaxation factor-based methods to solve coefficients of dynamic Smagorinsky model and SGS temperature variance model was verified using DNS for compressible isotropic turbulence flow, which was also used to assess the relevance between species concentration fluctuations and local temperature fluctuations in turbulent mixing layers as well as rationality of incorporating term $\partial x/\partial T$ within the SGS-TRI model. Additionally, by comparing the results of hot air free jet LES calculations with experimental data, the accuracy of the dynamic Smagorinsky model based on relaxation factors was further validated.

3.1. DNS for Decaying Isotropic Turbulence

3.1.1. Verification of DNS Method

Direct numerical simulations of compressible isotropic turbulence serve as a prime example for testing the resolution of small-scale fluctuations. This paper selects a computational case [33] involving isotropic turbulence with an initial Taylor-Reynolds number ${\mathrm{Re}}_{\lambda }=72$ and an initial turbulent Mach number ${\mathrm{Ma}}_t=0.5$, which decays under the influence of viscosity. The computational domain is $\left[0,2\pi \right]\times \left[0,2\pi \right]\times \left[0,2\pi \right]$ , and the number of grids is 256³. Convective terms are computed using the TENO-5 scheme [34], while viscous terms are discretized using a sixth-order central difference scheme. Time advancement is achieved through a TVD-type three-step third-order Runge-Kutta method [35], with a time step $\Delta t$ of 0.001s and a calculation duration of 3.5 times the large eddy turnover time $\tau$. The initial field conditions assume that the energy spectrum satisfies equation (34) and that the divergence of both velocity and temperature fields is zero.

$$E(k)=A{k}^4{e}^{(-{\displaystyle \frac{2k^2}{k_0^2}})}$$

(34)

where, $k$ is the wave number, $k_0$ is the wave number at the peak of the wave, and A is a constant value used to determine the specified initial turbulent kinetic energy. In this paper, the value is taken as $(A,k_0)=(0.00013,8)$.The calculation methods for turbulence statistics such as the Taylor scale Reynolds number, turbulent kinetic energy, turbulent Mach number, velocity derivative skewness, and large eddy turnover time are as follows,

$${\mathrm{Re}}_{\lambda }={\displaystyle \frac{u^{\text{'}}\lambda 〈\rho 〉}{〈\mu 〉}}$$

(35)

$$K_t={\displaystyle \frac{1}{2}}〈\rho u^2+\rho v^2+\rho w^2〉$$

(36)

$$M{a}_t={\displaystyle \frac{{〈u^2+v^2+w^2〉}^{\frac{1}{2}}}{〈c〉}}$$

(37)

$$skewness={\displaystyle \frac{〈{\left(\frac{\partial u}{\partial x}\right)}^3〉}{{〈{\left(\frac{\partial u}{\partial x}\right)}^2〉}^{\frac{3}{2}}}}$$

(38)

$$\tau =\sqrt{{\displaystyle \frac{32}{A}}}{(2\pi )}^{{\displaystyle \frac{1}{4}}}{k_0}^{-\frac{7}{2}}$$

(39)

where,

$${\rho }^{\text{'}}=\rho -〈\rho 〉$$

(40)

$$u^{\text{'}}=\sqrt{{\displaystyle \frac{u^2+v^2+w^2}{3}}}$$

(41)

$${\lambda }^2={\displaystyle \frac{{u\text{'}}^2}{〈{\left(\frac{\partial u}{\partial x}\right)}^2〉}}$$

(42)

After the flow field development time exceeds one large eddy turnover time, it can be considered that it is no longer affected by the initial conditions. As shown in Figure 4, the time evolution of the Taylor-scale Reynolds number, turbulent kinetic energy, turbulent Mach number, and root mean square of the normalized density in this article are in good agreement with the data from Samtanery [33]. The velocity derivative skewness is a high-order statistic that is very sensitive to small-scale motion in the flow field. The TENO-5 scheme reacts well to the development trend of this quantity, and the peak and valley are calculated.

.

3.1.2. Coefficient Verification for Dynamic Smagorinsky Model Based on Relaxation Factor

Scale filtering ($n=\overline{\Delta }/\Delta =4,8,16,32$) is applied to the aforementioned DNS data. Here, $\Delta$ represents the grid scale of the DNS computation, while $\overline{\Delta }$ denotes the filtering scale. When $n$ takes values of 4, 8, 16, and 32, the subgrid-scale energy accounts for 5%, 15%, 40%, and 70% of the total turbulent kinetic energy, respectively. Subsequently, $\alpha =0.99$ is employed for calculating equation (21). As illustrated in Figure 5, for $t/\tau =6~10$ and $n=4\mathrm{or}8$, the dynamic Smagorinsky model coefficients solved by the relaxation method are very close to the theoretical value 0.0384 (the notation 4-8$\Delta$ indicating $\overline{\Delta }/\Delta =4$ and $\widehat{\Delta }/\Delta =8$).

3.1.3. Correction of Temperature variance model coefficient

Based on the aforementioned DNS data, Figure 6 presents the time evolution curves for the coefficients of the dynamic temperature variance model calculated at different filtering scales, as well as the correlation coefficients between the calculating values and true values of temperature variance.

It can be observed that the calculated sub-grid temperature variance model coefficients $C_{\mathrm{T}}$ stabilize at 0.1, 0.2, and 0.6 for the three different filtering scales. This indicates that the coefficients are significantly influenced by the filtering scale, making it unsuitable to use a fixed coefficient model for parameterization. Furthermore, the values of $C_T=C_s^2$reported in reference [28]$(C_s=0.08~0.18)$ are considerably lower across various filtering scales. The correlation between the true sub-grid temperature variance and the values predicted by the dynamic model is quite strong, reaching 0.9 and 0.7 at filtering scales for$n=4$ and $n=8$ (when $n=16$, the reconstruction rate of turbulent kinetic energy falls below the fundamental requirements for LES), thereby demonstrating the accuracy of this parameterization method.

Figure 7 and Figure 8 display the variations of the true values of dimensionless sub-grid temperature variance and the predicted values from the dynamic model over time and filtering scales. It is evident that, when $n=4$ and $n=8$, the modeled values closely align with the true values as time changes, with a notable discrepancy only at larger filtering scales ($n=16$). This demonstrates the accuracy of the parameterization method within the scales of interest for large eddy simulations.

(left) and the correlation coefficient $CR(C_T)$(right) at different filtering scales.

, the variation of the dimensionless sub-grid temperature variance (true values) and the calculated values from the dynamical temperature variance model across different filter scales.

3.1.4. Verification and Analysis of the Correlation Between Gas Components and Temperature

To examine the correlation between gas component fluctuations and temperature fluctuations in chemically frozen turbulent flow fields, thereby to verify the rationality of adding term $\partial x/\partial T$ in the Snegirev SGS-TRI model, the $C(T,x)$ were calculated based on the DNS data from compressible isotropic turbulence flow field, as shown be follow,

$$C(T,x)={\displaystyle \frac{〈(T-〈T〉)(x-〈x〉)〉}{{\left(〈{(T-〈T〉)}^2〉〈{(x-〈x〉)}^2〉\right)}^{\frac{1}{2}}}}$$

(43)

The time evolution of $C(T,x)$ at different filter scales are shown in Figure 9. It is evident that, at various filter scales, both the component fluctuations and temperature fluctuations exhibit a very high degree of correlation. When $n=4$, the correlation between the two is even higher than the direct numerical simulation (DNS) results. This can be attributed to the filtering out of certain small-scale fluctuations, which enhances the correlation. At the scales of interest in large eddy simulation (LES), namely $n=8$ and $n=16$, the correlation remains at a high level, around 0.99 and 0.97 respectively. This supports the rationale for including the $\partial x/\partial T$ term in the SGS-TRI model to establish Eq. (28).

3.2. LES for High Temperature Air Jet Flow

The axisymmetric convergent nozzle has an outlet diameter of 50.8 mm, with a nozzle pressure ratio (the ratio of jet stagnation pressure to ambient pressure) of 1.357 and a static temperature ratio of jet at the nozzle outlet to the ambient temperature of 1.764. The ambient temperature is 288 K and the pressure is 101,325 Pa. The total grid number used in the LES computation was about 10 million, achieving a reconstruction rate of turbulent kinetic energy about 93%. The resulting instantaneous Mach number and temperature contour at jet symmetry plane are presented in Figure 10. The iso-surface of vorticity colored with local temperature is shown in Figure 11. The vorticity is caculated by criterion, as Equation (44). Furthermore, the calculated jet velocity field is compared to the PIV data from Bridges and Wernet [36], while the temperature field is validated against the rotating Raman scattering results from Locke et al. [37] and the Rayleigh scattering results from Mielke et al. [38], as shown in Figure 12 and Figure 13. Here, $r$denotes the radial position, $D_j$ specifies the diameter of the jet outlet, $U_j$ indicates the velocity at the center of the jet outlet, $T_j$ represents the temperature at the center of the jet outlet, $T_{\infty }$ refers to the temperature of the ambient air, $\Delta T_j=T_j-T_{\infty }$ and $k$ signifies the turbulent kinetic energy of the jet.

$$Q={\displaystyle \frac{1}{2}}\left({‖{\displaystyle \frac{\Delta V-\Delta V^T}{2}}‖}^2-{‖{\displaystyle \frac{\Delta V+\Delta V^T}{2}}‖}^2\right)$$

(44)

The numerical results obtained from the large eddy simulation presented in this paper accurately capture the turbulent fluctuation characteristics of the hot jet, as well as the flow direction and spanwise decay characteristics of the jet in the surrounding air. Given that the size and Reynolds number of the jet in this case are relatively small, the results of the large eddy simulation exhibit considerable sensitivity to the sub-grid eddy viscosity model. The dynamic Smagorinsky model, which is based on a relaxation factor developed in this study, effectively addresses these requirements. The differences between the calculated results of turbulent kinetic energy distribution along the axis and the experimental data may stem from the difference between the Favre-averaged value of turbulent kinetic energy output by the sub-grid model and its ensemble average.

4. Analysis of TRI Characteristics of Axisymmetric High-Temperature Gas Jet

4.1. Verification of Infrared Radiation Characteristic Calculation for Combustion Gas Jet

The calculation of the infrared radiation characteristics of an axial combustion gas jet in the 3-5micron band presented in case1 of this section is still based on the convergent nozzle described in Section 3.2, the only difference is to replace the jet medium by combustion gas. This adjustment takes into account the influence of the effect of real gas, resulting in a gas jet length that is slightly longer than that of the air jet. Case 2 has higher gas temperature, which is regulated through the combustion of aviation kerosene at varying equivalence ratios, simultaneously adjusting its components. Additionally, the nozzle size is enlarged to maintain the Reynolds number consistent with that of case 1, as detailed in the table below.

Table 1. Parameters of LES cases.

	Case 1	Case 2
Diameter of nozzle outlet/mm	50.8	97.2
Total temperature/K	554.62	794.65
Total pressure/Pa	137498
Ma	0.67
Mass fraction of oxygen	20.85%	18.63%
Mass fraction of water vapor	0.79%	1.54%
Mass fraction of carbon dioxide	2.05%	3.98%
Mass fraction of carbon monoxide	0.000015%	0.000015%

Table 1. Parameters of LES cases.

	Case 1	Case 2
Diameter of nozzle outlet/mm	50.8	97.2
Total temperature/K	554.62	794.65
Total pressure/Pa	137498
Ma	0.67
Mass fraction of oxygen	20.85%	18.63%
Mass fraction of water vapor	0.79%	1.54%
Mass fraction of carbon dioxide	2.05%	3.98%
Mass fraction of carbon monoxide	0.000015%	0.000015%

The ambient air is set to the equatorial mode in the MODTRAN 5.0 software (where the mass fraction of water vapor and carbon dioxide are1.44%and 0.06%, respectively, meaning the proportion of radiation participating components is completely different from that of the combustion gas jet), and the aerosol is set to the sea surface mode with 23km visibility. The extinction efficiency of aerosols in 3-5 micron band is shown in Figure 14. The infrared detection direction is perpendicular to the jet axis, and the detection distances are set to 0.2, 3.0, and 45km, respectively. For each physical time step, flow field LES is performed (computation time is about 1.1 seconds on a single computing node with 256-cores/512-threads), followed by an infrared imaging calculation for the entire flow field (computation time is about 0.4 seconds). Figure 15 shows the instantaneous infrared imaging and infrared radiance calculation error distribution of case 2 at a certain moment under 45km detection distance where the radiance error distribution of MSMGWB model is calculated by contrasting with results of the line-by-line model [39] (requiring approximately 6 hours of computation time) ignoring the influence of SGS-TRI. The error is found to be within 8%. This demonstrates that the MSMGWB model is adept at long-range radiative transfer computations in scenarios characterized by substantial disruptions of correlated-k characteristics of absorption spectra, stemming from inhomogeneities of gas temperature and composition.

4.2. TRI Characteristic Analysis for High Temperature Combustion Gas Jet

Figure 16 and Figure 17 illustrate 3-5micron band infrared images for case 2 at various distances derived by ray-trace [21] based on the time-average flow field (i.e., ignoring TRI effect) and the instantaneous flow fields followed by time-averaging. The latter has longer core areas and exhibits higher maximum radiance compared to the former. The gas jet infrared signal detected by the actual infrared detector is clearly the latter.

Figure 18 illustrates the distribution of percentage increases in 3-5micron band infrared image, attributed to TRI for case 2 at 0.2 km detection distance. It is evident that flow field fluctuation significantly affects infrared radiation, especially in the violent turbulent regions of gas-ambient mixing layer. Conversely, in the core area characterized by high temperature and high radiance, the TRI impact is relatively small, especially near the center of nozzle exit, which aligns with the findings presented in the literature [8]. Radiation intensity of pixels in the infrared images corresponding to various positions of combustion gas jet at 0.2 km detection distance is shown in Figure 19. Figure 20 illustrates the time evolution curve of the integrated radiation intensity of the whole image. The normalized fluctuation of the integrated radiation intensity is defined as follows,

$$M\text{'}=(M-〈M〉)/〈M〉$$

(45)

The relative amplitudes and frequencies of radiation intensity of each pixel are very high, but the integral of them has much lower relative amplitude and frequency, causing by canceling one another’s fluctuation.

Figure 21 and Figure 22 illustrate the percentage influence of TRI/SGS-TRI on the integrated radiation intensity of the jet at various detection distances. Notably, the effect of component fluctuation in SGS-TRI is determined by excluding the term ${\displaystyle \frac{\partial x}{\partial T}}{\displaystyle \frac{\partial k_{\mathrm{j}}^{*}{A}_{\mathrm{j}}{I}_{b,{\eta }_L-{\eta }_U}}{\partial T}}$ in equation (28).

The two figures illustrate that, at various jet temperatures, the TRI accounts for more than 1/3 of the total infrared signal, with its influence further amplifying as the detection distance increases. Among the effects of TRI, the resolvable scale TRI plays a dominant role; however, the contribution of SGS-TRI is also significant, accounting for over 10% of the total infrared signal of the jet, which corresponds to approximately 50% of the contribution from the resolvable scale TRI. This finding is consistent with the proportions of SGS-TRI reported in the literature based on DNS studies of chemically frozen free jet [10,11]. Within the SGS-TRI effect, the contribution from component fluctuation is approximately one-third, underscoring its importance.

To analyze why the influence of TRI on the jet infrared signal increases with detection distance, Figure 23 presents the temperature dependences of emission coefficients and corresponding atmospheric transmittances for RTE’s with a significant contribution ratio (W) to $I_{3-5\mu m}$ at various detection distances, their cumulative contribution ratio accounts for approximately 90%. For clarity, the minimum transmittance is set at no less than $1\times {10}^{-20}$.

At shorter transmission distances, where the atmospheric attenuation is minimal (red line in the right part of Figure 23), the combustion gas radiation spectrum is mainly contributed by the cold absorption lines of water vapor and carbon dioxide, their corresponding $k_{\mathrm{j}}^{*}$ values vary slowly with temperature increasing. Conversely, at longer transmission distances, the emitted radiation from the cold absorption lines of combustion gas is almost entirely absorbed by the same cold lines of water vapor and carbon dioxide in low-temperature atmosphere (blue line in the right part of Figure 23). As a result, the remaining infrared signals predominantly originate from hot absorption lines in the combustion gas, which abilities of radiation absorption (red line in the left part of Figure 23) and emission are very weak at low temperature states, but exhibiting a rapid increase with temperature increasing (high $\partial k_{\mathrm{j}}^{*}/\partial T$ value). Consequently, temperature fluctuations in turbulence flow induce a more pronounced positive correlation between $k_{\mathrm{j}}^{*}$ and $I_{b,{\eta }_L-{\eta }_U}$for the spectrum consisted of the hot absorption lines compared to that of cold absorption lines. This ultimately results in an increased relative contribution ratio of TRI as detection distance increases.

5. Discussion and Conclusions

This paper employs LES and a MSMGWB k-distribution gas radiation model to investigate the effects of high-temperature combustion gas jet turbulence fluctuation on its remote infrared radiation characteristics within the 3-5micron band. To address the issue that existing algorithms significantly underpredict the SGS-TRI, we introduced a term related to the correlation between fluctuations of radiation-participating components and temperature into the SGS Snegirev relationship. Additionally, we developed an ensemble average calculation algorithm to solve dynamic sub-grid model coefficient which is easy to implement in complex engineering calculations. Building on this foundation, we established a dynamic temperature variance model, and the validities of the abovementioned algorithms was confirmed through isotropic turbulence DNS and comparison with experimental data for high-temperature air jet. The conclusions drawn include,

(1) The ensemble average calculation algorithm presented in this article to solve coefficients of dynamic sub-grid models is characterized by its simplicity, low memory consumption, and robust adaptability to complex three-dimensional flow problems. It is competent to free jet LES by combining with dynamic Smagorinsky model. Moreover, the coefficients of the SGS temperature variance model exhibit significant variation with filtering scale, meaning the use of fixed coefficients is inappropriate for its modeling. Additionally, the fixed coefficient values used in existing temperature variance models are considerably low at each filtering scale, which is the main cause of underprediction for SGS-TRI characteristics. This issue can be substantially mitigated by employing the dynamic temperature variance model in conjunction with the aforementioned ensemble average calculation algorithm.

(2) The influence of the TRI effect on the 3-5micron band infrared characteristics of high-temperature combustion gas jets is primarily concentrated in the combustion gas-air mixing layer, which accounts for more than 1/3 of the overall infrared radiation intensity. This proportion is expected to increase further with greater detection distances. The reason is that as the detection distance increases, infrared radiation from the "hot absorption lines" in the combustion gas jet gradually takes the dominant position. The emitted radiation from the "hot absorption lines" increases more quickly with gas temperature increasing, resulting in stronger radiance fluctuation induced by temperature fluctuations than that from the "cold absorption lines". SGS-TRI accounts for about 10% of the total infrared radiation signal of the gas jet, so it must be reasonably modeled in the calculation of the infrared characteristics of the gas jet in 3-5micron band. The proportion of temperature and gas component-related fluctuations in SGS-TRI at 3-5micron band can also exceed 1/3. Based on the analysis of DNS data, considering its impact can significantly improve the modeling accuracy for SGS-TRI of combustion gas jet flow.

(3) Calculation accuracies of conventional gas radiation models mainly depend on the correlated-k characteristics of gas media long radiation transmission paths. For radiation emitted by hot combustion gas jet and attenuated by atmosphere, the CK characteristics are disrupted by severe discrepancies in temperature and composition between the combustion gas and ambient air. The MSMGWB k-distribution model can effectively address the problem and due to its small computational load, it is suitable for calculating the impacts of turbulence fluctuations on the remote infrared radiation characteristics of combustion gas jet.