Detailed simulation of LOX/GCH4 flame-vortex interaction in supercritical Taylor-Green flows with machine learning

For two-dimensional TGV with a small Reynolds number, an analytical solution can be obtained. In this configuration, liquid oxygen at $150$ K and $10$ MPa is selected as the fluid with the kinematic viscosity of $\nu=7.7548\times 10^{-8}$ m²/s. The analytical solution for the velocity is given by

$$u(x,y,t)=u_{0}\times sin(\frac{x}{L_{0}})\times cos(\frac{y}{L_{0}})\times e^{\frac{-2\nu t}{L_{0}^{2}}},$$

(2)

$$v(x,y,t)=-u_{0}\times cos(\frac{x}{L_{0}})\times sin(\frac{y}{L_{0}})\times e^{\frac{-2\nu t}{L_{0}^{2}}}.$$

(3)

Figure 2 shows the comparison between the analytical solution and simulation results at $t=10\tau_{ref}$, using different grid resolutions and spatial schemes. The results suggest that the linear scheme (2nd order central differencing) provides a satisfactory accuracy as in the 4th order cubic scheme. Also, resolutions of $N=128$ and $256$ give no obvious difference, indicating a good numerical grid convergence.

This step aims to further investigate the fluid processes in TGV with the real gas of oxygen at $150$ K and $10$ MPa. In Fig 3, the temporal evolution of the normalised volume-averaged kinetic energy and enstrophy is compared to the non-dimensional results obtained with a pseudo-spectral method in [27, 14]. The figure suggests that the progression of non-dimensional enstrophy under extreme pressure and thermodynamic conditions in cold flow scenarios exhibits a similar pattern to that in subcritical cases. A notable difference is that the enstrophy increases faster in the supercritical case, which is mainly due to the larger density gradients resulting in stronger vortex formation as observed in previous mixing layer studies [18].

Initial fields of temperature and $Q$-criterion for the three-dimensional multi-species reacting case are shown in Fig 4. This setup closely follows that of [14] to allow for a reasonable comparison with the ideal gas benchmark.

The differences lie in the thermodynamic states, for which a supercritical pressure of $10$ MPa is specified uniformly across the domain. The initial temperature and species distributions are similar to those in [14], but for a CH₄/LO_x mixture with unburnt oxygen at $150$ K and methane at $300$ K. The high temperature region is interpolated from the equilibrium state calculated from a two-step five-species  [28] mechanism. To compare the differences of the flame-vortex interactions using the ideal gas model and PR model, the isovolume of temperature and logarithm of $Q$ fields are depicted from $t=1$ to $4\tau_{ref}$ for both cases in Fig 5. The $Q$-criterion is defined as $Q=\frac{1}{2}~{}(\|\boldsymbol{\Omega}\|^{2}-\|\boldsymbol{S}\|^{2})$, where $\boldsymbol{\Omega}$ is the vorticity tensor and $\boldsymbol{S}$ is the strain rate tensor. The positive values of $Q$ are indicative of areas in the flow field where the vorticity dominates and negative values of $Q$ are indicative of strain rate or viscous stress dominated areas. Thus, logarithm of $Q$ is used in Fig 5 to visualise vortical structures. It is seen that similar vortical structures can be observed in both cases, and the initial lumpy vortices are rolled up into vortex tubes surrounded by the spreading flame surface from $t=1$ to $2\tau_{ref}$. Parts of the tubes are then stretched and twisted into thinner ones in following time instants. Considering the flame-vortex interaction, the ideal gas case has the identical flame structures and evolving stages as compared to the hydrogen TGV-flame at subcritical conditions simulated using DNS with a similar Reynolds number [17]. At the same reference times, although the temperature under extreme pressure is much higher, the deformations of the flame and the stretching of the vortices are similar. By contrast, the interactions between the flame and vortex in real gas case start to behave differently from $t=2\tau_{ref}$. Compared to the ideal gas case, more pronounced deformations of the flame surface occur due to stronger vortices and higher strain rate. Moreover, breaking and quenching of the flame are observed especially in the $y$-direction.

To further understand the differences appearing from $t=2\tau_{ref}$, the contour plots for the $x$-$y$ cross-section plane is presented in Fig. 6. Several evident differences can be observed in this figure. Firstly, the real gas case exhibits significant density gradient near the flame due to the substantial difference in density between the LOX and GCH4 mixtures. Moreover, there are two regions with notably high $U_{y}$ around $y$ central line (marked by black dashed lines), swirling and stretching the flame away from the centre of the domain. As seen in Fig. 6b, flame quenching occurs in this region due to the high strain rates, which is significantly different to the ideal gas case in Fig. 6a. To eliminate the possible influence of chemical mechanism, the simulation performed using a detailed 17-species/44-step mechanism [29, 30] is presented in Fig. 6c (see detailed quantitative comparison in the Supplementary Material). Despite the considerably low peak temperature, the detailed chemistry case exhibit similar flow and flame patterns as the 2-step mechanism case. In particular, the flame quenching induced by the high strain near the central line is also observed, suggesting that this behaviour is mainly driven by the real gas thermodynamic effects.

To investigate how these differences are generated, vorticity and density fields at $t=1\tau_{ref}$ are compared between the real gas and ideal gas cases in Fig 7. As mentioned in Section 2, the unburned oxygen lies in the ‘liquid-like’ region while the high-temperature mixture are more ‘gas-like’ during the TGV evolution. Similar to that observed in [31], the significant density gradients lead to a dramatic change in the momentum of the two regions, and thus the dense oxygen acts as a ‘solid obstacle’ against the rotating internal gaseous part of the flow. Comparing the two cases in this figure, it can be seen that due to the constraining effect of the ’solid obstacle’, when the vortex rolls up the flame surface to the centre, the flame appears in a narrower region near the centre of the $x$-$y$ plane with higher vorticity. The stronger vorticity and strain rate in the real gas case leads to further stretching and quenching of the flame as shown earlier in Fig 6.

To shed further light into the flame-vortex interaction, Fig. 8 presents the Probability Density Function (PDF) of curvature and tangential stain rate for the two cases at $t=4\tau_{ref}$. The positive tangential strain rate of real gas with higher PDF is seen in Fig. 8a for the real gas case indicating a stronger flame stretching, which is consistent with the curvature PDF in Fig. 8b. The high peak of curvature distribution around zero in the real gas case also implies that the flame is globally less deformed but with substantial local wrinkles (a second peak exists) compared to the more distributed shape in the ideal gas case.

The current investigation adopts a multilayer perceptron (MLP) framework, incorporating $3$ hidden layers with the Gaussian Error Linear Unit (GeLU) activation function. Inputs and outputs are tailored for the CFD solver used, where the real fluid properties as model outputs are updated after solving the energy and species equations, while inputs contain enthalpy, pressure, and Bilger’s mixture fraction from the previous iteration. Seven dedicated networks are designed for inferring temperature, density, viscosity, compressibility factor ($\psi=\frac{\partial\rho}{\partial p}$ at constant enthalpy), thermal diffusivity, mixture-averaged diffusion coefficients of each species, and nondimensional standard state enthalpies. Should the variable exhibits inter-species variations, the network is configured with multiple outputs; otherwise, a single-output neural network suffices. Among these variables, the accuracy of temperature, density, and compressibility factor significantly impacts the convergence and stability of CFD simulations.

The training dataset is obtained from the initial 10 time steps of a Taylor-Green Vortex (TGV) simulation. To capture a wide range of operational states during TGV evolution, additional stochastic datasets are generated from the collected results by

$$t=c+a_{1}\frac{max(h)-min(h)}{a_{2}},$$

(4)

where $t$ is the training data point, $c$ is the data sampled from CFD, $h$ is the physical variable, and $a_{1}$, $a_{2}$ are two random numbers that can be chosen to the fit specific case. This process requires constraints to prevent non-physical outcomes, ensuring mixture fractions to range from 0 to 1, and densities remain positive. Recognising the multi-scale nature of the broad density spectrum in supercritical conditions, pre-processing techniques are used to optimise the neural network performance. This includes applying the Box-Cox Transformation (BCT) to both density and mixture fraction, highlighting their significant but subtle features, particularly the clustering towards smaller magnitude values.

Subsequent to sample data generation, the models are trained using the machine learning frameworks PyTorch and Scikit-learn.The L1 norm is designated as the loss function to minimise the absolute error during training. The training of each model is limited to $1,000$ epochs, beyond which no significant improvement in loss reduction is observed. The wall-clock time required for training each network approximates one hour. This duration is kept short to ensure that the expeditiousness imparted by the deployment of the DNN is not nullified by the training process itself.

The fidelity of network predictions is evaluated by a priori assessment and comprehensive comparisons with flow field results obtained from CFD simulations. Figure 9 depicts a scatter plot comparing neural network (NN) predictions with labelled data accompanied by a colour-coded gradient to illustrate the distribution of local relative errors. Each plot includes an R2-score, a widely recognised measure of regression model performance, which affirms the strong alignment between NN outputs and the intended targets, suggesting an overall good model accuracy for the thermophysical quantities considered. Subsequent application of the NN framework to a three-dimensional reacting TGV simulation (discussed in Section 3.3) yields the density $\rho$ and temperature fields presented in Fig. 10. Both the field and line graphs show satisfying agreement between the DNN and Cantera real fluid package. Employing the relative Mean Squared Error (MSE) defined in [20] as a metric for validation, the relative errors in density $\rho$, pressure $p$, and temperature $T$ are confined to $2.1\%$, $1.2\%$, and $9\%$, respectively. This suggests a reasonable predictive accuracy of the NN model for the thee-dimensional supercritical reactive TGV simulation considered in this study.

Finally, computational acceleration provided by the machine learning approach is of primary interest. Figure 11 compares computation time required to advance one time step of the supercritical reactive TGV simulation using Cantera and DNN for thermophysical property calculations with 16 CPU cores. The DNN inference was performed using one additional GPU card. A speed-up factor of 13 is achieved for the thermodynamics when the DNN model is used. A further acceleration strategy is applied for the other time-consuming part $-$ chemical reaction source term integration [22], resulting in an overall 12 times speed-up for the entire CFD simulation.