The rocket exhaust plume is a free turbulent flow with high Reynolds numbers, which involves multi-component non-equilibrium chemical reaction wave interference, shear, and turbulence effects, accompanied by vortex coupling effects of different scales. In order to simulate turbulent eddies and capture finer turbulence pulsation information, the Large Eddy Simulation (LES) method [16] is used. The LES method decomposes the instantaneous pulsating motion in turbulence into large-scale vortices and small-scale vortices through filtering functions. Large scale eddies that have a significant impact on the average flow field can be directly calculated using the Navier Stokes (N-S) equation, and small-scale eddies are simulated using the Sub Grid Scale (SGS) model [17]. The instantaneous pulsating is expressed in Equation 1: q = q ¯ + q ′ (1) where q ¯ represents the large-scale vortices and q ′ denotes the small-scale vortices. q ¯ can be obtained from the decomposition using the filtering function as Equation 2: q ¯ x i , t = ∫ G x i , x i ′ q x i ′ , t d x i ′ (2) where x i is the component in the coordinate direction, and G x i is the filtering function. The compressible N-S equation filtered by a filter can be expressed as Equation 3 [18]: ∂ ρ ∂ t + ∂ ρ u i ∂ x i = 0 ∂ u i ∂ t + ∂ u i u j ∂ x i = − 1 ρ ∂ P ∂ x i + 1 R e ∂ 2 u j ∂ x i ∂ x j + ∂ τ i j ∂ x j (3) where u is the filtered velocity component, P is pressure, ρ is the density, R e is the Reynolds number, and τ i j is a subgrid stress tensor that represents the influence of small-scale eddies on large-scale eddies.

The Smagorinsky Lilly model is adopted [18], which takes the form of Equations 4, 5: τ i j − 1 3 τ k k σ i j = 2 v t S ¯ i j (4) S ¯ i j = 1 2 ∂ u ¯ i ∂ x j + ∂ u ¯ j ∂ x i (5) where v t is the turbulent viscosity at the grid scale, expressed as Equations 6, 7: v t = C 2 Δ 2 S ¯ = C 2 Δ 2 2 S ¯ i j S ¯ i j (6) Δ = Δ x Δ y Δ z 1 / 3 (7) where Δ is the grid size after grid filtering, and Δ x 、 Δ y and Δ z are the grid sizes.

The finite volume method (FVM) is used to discretize the governing equations. The diffusive term is processed using a second-order central difference scheme, and the convective term is discretized using a flux limited second-order upwind scheme. The chemical reaction adopts a 9-species 10-reaction finite rate kinetics with the H₂/CO system. The detailed chemical reaction kinetics are borrowed from Ref. [19].

The main radiating species in the gas-phase products generated by propellant combustion have H₂O, CO₂, CO, OH, and HCl. The vibrational and rotational transition emission of high-temperature gas-phase components is the main source of plume infrared radiation. Without considering particle scattering and emission, physical parameters of infrared spectroscopy can be calculated using the Statistics Narrow Band (SNB) method [20]. The spectral transmittance is expressed as Equation 8: τ ¯ = exp − 2 a ¯ 1 + κ ¯ u ˙ a ¯ 1 / 2 − 1 (8) where κ ¯ characterizes the ability to emit/absorb photons, u is the length of the pressure path, and a ¯ = π γ ¯ / d ¯ denotes the fine-structure parameter related to Doppler broadening and collisional broadening. The Doppler broadening parameter γ D , s (expressed as Equation 9) for gas component s at wavenumber η is provided based on the Single Line Group (SLG) [21] model. γ D , s = 5.94 × 10 − 6 η T 273.0 M s (9) where M represents the molecular weight, and T represents the temperature.

The NASA-SP-3080 database [22] provided spectral band calculation parameters for six gas molecules, including H₂O, CO₂, CO, HCl, OH, and NO, in the temperature range of 300–3,000 K. Based on the given spectral parameters, the SNB method can be used to obtain the spectral parameters for the infrared radiation calculation.

Without considering particle scattering effects, the absorption and scattering terms of the radiative transfer equation (RTE) can be simplified. By introducing the optical thickness τ λ = κ s λ L , the RTE is expressed as Equation 10: d I λ τ λ d τ λ = − I λ τ λ + I b λ τ λ (10) where λ is the wavelength. I λ and I b , λ represent spectral radiation intensity and blackbody spectral radiation intensity, respectively. The line of sight (LOS) method is used to solve the RTE.

Gaussian white noise is added to the plume flow field data, and the infrared thermal imaging data is generated through infrared radiation calculation and the SPOD denoising method. The infrared radiation data after adding noise is represented as Equation 11: Q ∼ = Q + Q Noise (11) where Q ∼ represents the radiation thermal imaging data with noise, which is arranged in matrix format from thermal imaging data at different times, Q represents the thermal imaging data without noise, and Q Noise represents the residual caused by noise.

The matrix Q ∼ is partitioned into smaller matrices through block matrix processing, which is referred to as blocks. After this partitioning, the matrix Q ∼ is divided into N b blocks, with each block containing N _f snapshots. A discrete Fourier transform is then applied to each block as Equation 12: Q ^ n = F F T Q n = q ^ 1 n , q ^ 2 n , ⋯ , q ^ N f n (12) where q ^ k n is the Fourier component at the k-th discrete frequency f _k in the n-th block. The data matrix Q ^ f k is composed of Fourier components with a frequency of f _k in each block matrix [23], expressed as Equation 13: Q ^ f k = q ^ 1 1 , q ^ 2 2 , ⋯ , q ^ k N b ∈ R N × N b (13) where N is the number of grid points. The SPOD mode and its corresponding eigenvalues can be obtained by constructing a cross spectral density matrix M f k = Q ^ f k Q ^ f k * and performing eigenvalue decomposition on it. The resulting modes are then sorted according to the magnitude of their eigenvalues.

At each frequency, the Fourier transform matrix Q ^ m ω k is reconstructed as Equation 14: Q ∼ f k = Φ f k n b l k Λ f k 1 / 2 Θ f k * (14)

The Fourier snapshot matrix of the l-th block is represented as Equation 15 [24]: Q ^ m ω l = ∑ i = 1 L a i k l ϕ i k = 1 , ∑ i = 1 L a i k l ϕ i k = 2 , ⋯ , ∑ i = 1 L a i k l ϕ i k = N f (15)

Data denoising can be achieved by rearranging the matrix Q ^ f k to form a reconstructed Fourier matrix Q ^ n and then applying the inverse Fourier transform, expressed as Equation 16: q k n = 1 w j F − 1 q ^ k n (16)

By the inverse Fourier transform, the denoised flow field Q _Denoise can be obtained. The detailed derivation process can be found in Ref. [12]. The denoising of infrared thermal images is achieved using the SPOD inversion algorithm. Figure 1 illustrates the schematic diagram of infrared thermal image denoising using the SPOD inversion.

The noise applied to the unsteady flow field can be amplified after infrared radiation calculation. Figure 2 show the infrared radiation intensity map of the plume in the 2.7 and 4.3 μm wavelength bands, respectively. The illustrations from top to bottom show the snapshots of original transient infrared thermal images the thermal image with noise and the difference distribution between the two. As shown in the figure, a series of Mach cells are present in the intrinsic core of the plume, showing a trend of higher infrared radiation intensity as they approach the nozzle exit. Due to the influence of shear effects, small-scale vortex structures are rolled up in the shear layer near the third cell from the nozzle (x ≈ 20 cm). These small-scale eddies cause unstable fluctuations in infrared thermal images. The noise causes the details of the thermal image to become blurred, and the shock wave structure in the intrinsic core area also exhibits irregular changes. Figures 2C1, C2 show the residual distribution of infrared radiation intensity with and without noise influence. It can be observed that the residuals are mainly distributed in the areas where the plume passes after applying Gaussian white noise, and the residuals in the core area are relatively large, while the residuals in the turbulent mixing area are relatively small.

The relative error is used to evaluate the difference between the noisy and real data, and its formula is expressed as Equation 17: E r r o r = norm Q ∼ − Q norm Q (17) where norm (▪) represents the norm of a matrix or vector.

Figure 3A shows the error distribution between the thermal images with and without noise in the 2.7 and 4.3 μm bands. From Figure 3A, it can be seen that the noise increases the error level of infrared radiation intensity to varying degrees across different bands. The average errors of infrared radiation intensity in the 2.7 μm and 4.3 μm bands are 21.1% and 17.6%, respectively. It means that the infrared radiation intensity in the 2.7 μm band is more affected by noise compared with the 4.3 μm band.

The suppression of infrared noise is achieved through the SPOD frequency domain reconstruction. Gaussian white noise has different effects on the structure of the flow field at different scales. Therefore, the threshold truncation method is used to filter out the small-scale coherent structures in the plume that are greatly affected by noise, in order to preserve the large-scale flow structure. By using the SPOD inversion, it is possible to achieve radiation denoising while preserving the original data of the thermal imaging field to the greatest extent possible. Figure 3B shows the error distributions of different bands before and after infrared thermal imaging denoising. From Figure 3B, it can be seen that the error of the denoised infrared thermal image has a trend of first decreasing and then increasing with the increase of the truncation threshold. When the truncation threshold is low, the dominant noise structure in high frequency is not filtered out, resulting in poor denoising quality. An excessively large truncation threshold can result in the filtering out of low-frequency large-scale flow structures with high energy consumption in infrared thermal images, and the noise reduction quality is also poor. Numerical analysis reveals that truncation occurs around λ = 0.001, and the flow field denoising reaches its minimum value. Infrared thermal imaging errors in the 2.7 and 4.3 μm bands have been reduced from 21.1% and 17.6% to 9.4% and 9.2%, respectively.

Figures 4A, B show the 4.3 μm and 2.7 μm infrared thermal images obtained by using the threshold truncation and the inversion SPOD algorithm, respectively. The upper illustration shows a thermal image snapshot with noise. It can be seen that the application of noise causes blurring of thermal image. The lower illustration shows the distribution of the thermal image contours without noise. It can be observed that the key features of the thermal image after denoising are more prominent, indicating that the proposed infrared thermal image denoising algorithm is reasonable.