Computation of Real-Fluid Thermophysical Properties Using a Neural Network Approach Implemented in OpenFOAM

1. Introduction

This introduction provides an overview of the historical and evolving role of ML in CFD, highlighting its transformative potential across numerous applications. As we delve deeper into the intersection of ML and CFD, we will explore specific methodologies, case studies, and emerging trends that underscore the symbiotic relationship between these two fields, ultimately reshaping our understanding of fluid dynamics and propelling innovation in diverse industries.

Admittedly, much progress has been made the last years to promote understanding of combustion, from computational sides, analyzing hydrogen high-pressure injection and mixing processes [1] including assessment of a model which accounts for real fluid thermodynamics and transport properties, making use of the Peng-Robinson equation of state (PR-EoS) and the Chung transport model [2]. Reitz and co-workers focused on a thermodynamic analysis of the mixture states [3]. The challenges primarily come from real-gas effects, which have significant effects on these processes, introducing non-linearities into the thermodynamic system and resulting in non-physical pressure oscillations that can substantially affect accuracy. Fluid properties in states near a pure component’s critical point are the most difficult to obtain from both experiments and from models such as EoS [4]. The principal experimental difficulty is that the density of a near-critical fluid is so extremely sensitive to variations in pressure and temperature, that maintaining homogeneous and stable conditions takes extreme care [5]. Even gravity influences the measurements. The principal model difficulty is that near-critical property variations do not follow the same mathematics as at conditions well-removed from the critical point [6].

Advanced propulsion and transportation systems, from rocket to diesel engines, often involve high-pressure injection of substances. In liquid-fueled rockets, propellants are introduced in a supercritical or transcritical state, exceeding fluid critical pressure and experiencing temperatures near or below the critical temperature. This deviates from traditional subcritical injection, causing significant changes in fluid properties due to temperature fluctuations. The assumption of negligible molecule volume breaks down, leading to deviations from ideal gas behavior.

Numerous theoretical and experimental studies, such as Puissant et al.'s investigations on shear coaxial injection [7] and Mayer's understanding of LOX/GH2 rocket engine combustion processes [8], have been conducted on high-pressure flows. These studies aim to gain insight into chemical and physical processes, providing a database for improved modeling and validation of simulation codes in combustor design and optimization [9]. Analyzing hydrogen high-pressure injection in the near-nozzle region, investigating the formation process and the structure of the Mach disks and the transition to turbulent jets, is done by Rahantamialisoa et al. [10], and Oschwald et al. [11] present experiments on inert binary injection and mixing processes.

While these works contribute significantly, their constraints in fully capturing information through experiments are acknowledged. Consequently, there is an increasing emphasis on modeling studies focused on transcritical and supercritical flows.

In the near future, the necessity will arise for designing adaptable energy production and propulsion systems that can efficiently utilize a wide range of fuels, including synthetic, bio-derived, fossil fuels, or various combinations thereof. Modeling these fuels, especially when their properties are not initially known, will demand flexible methods that can describe essential operational traits with minimal input data. To achieve this objective, following Koukouvinis et al. [12], we propose the application of machine learning (ML) to provide quantitative predictions under unknown conditions by leveraging existing research, be it from experiments or simulations.

Using tabulation, ML and artificial Neural Networks (NN) models for tackling the complex issue of transcritical sprays, which are relevant to modern compression-ignition engines can assist the speeding up the calculations. Direct simulations can incur significant CPU costs, whereas tabulation may be memory-intensive and challenging to expand as the number of chemical species grows [13,14]. In contrast, NN have emerged as a remarkably efficient technique for classification and response prediction. In high pressure, high temperature conditions, the NN methods are used to predict the thermodynamic properties [12] and the tabulation method offers high-precision calculation results in a wide temperature and pressure range [13,14]. ML techniques are rapidly advancing in this era of big data, and there is high potential in exploring the combination between ML and combustion research and achieving remarkable results. Much of this interest is attributed to the remarkable performance of deep learning architectures, which hierarchically extract informative features from data [15]. In the field of fluid mechanics, NNs are currently being investigated as a complementary tool to Computational Fluid Dynamics (CFD) for accelerating design processes [12]. A recent thorough review discusses a variety of applications in fluid mechanics, involving flow feature extraction, turbulence modelling, optimization and flow control [15]. In addition, an implementation of the NN models in OpenFOAM using PyFOAM tool, which is an in situ data analysis with Python and OpenFOAM is done by R. Maulik et al. [16,17,18].

This work discusses the capabilities of a new solver, called realFluidReactingNNFoam using NN algorithm for replacing complex real-fluid thermophysical properties evaluations, using the approach of coupling OpenFOAM and Python-developed models.

2. Materials and Methods

2.1. CFD Model

The governing equations for a fully conservative, homogeneous, multicomponent, and compressible non-reacting two-phase flow are the conservation equations for mass, momentum, energy, and species which are written below [2,19,20]:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{U}\right)=0$$

(1)

$$\frac{\partial \rho \mathbf{U}}{\partial t}+\nabla ·\left(\rho \mathbf{U}\mathbf{U}\right)=\nabla ·\left(-p\mathbf{I}+\mathit{\tau }\right)$$

(2)

$$\frac{\partial \left(\rho h_T\right)}{\partial t}+\nabla ·\left(\rho h_T\mathbf{U}\right)=\frac{\partial p}{\partial t}+\nabla ·\left(\mathbf{U}·\mathit{\tau }\right)-\nabla ·\mathbf{q}$$

(3)

$$\frac{\partial \left(\rho Y_k\right)}{\partial t}+\nabla ·\left(\rho Y_k\mathbf{U}\right)=-\nabla ·{\mathbf{J}}_k$$

(4)

In the above equations, $\rho$ is the mixture density, $\mathbf{U}$ is the mixture velocity vector, $p$ is the pressure shared by all species, $\mathit{\tau }$ is the viscous stress tensor of the mixture, $Y_k$ is the mass fraction of species k, and ${\mathbf{J}}_k$ and $\mathbf{q}$ refers to diffusion flux of species $k$ and heat flux, respectively. Also, $h_T=h+1/2{\mathbf{U}}^2$ is the total enthalpy of the mixture.

Additionally, for the evaluation of thermodynamic fluid properties under supercritical pressure, the PR-EoS is employed:

$$p\left(v,T\right)=\frac{RT}{v-b_m}-\frac{a_m}{v^2+2b_{m}v-b_m^2}$$

(5)

where $v$ is molar volume, T is the temperature, $b_m$ is the effective molecular volume, $a_m$ is the attractive force between molecules (note that the subscript m refers to the mixture) and $R$ is the gas constant. In such approach, local mechanical and thermal equilibrium between the two phases is always enforced. As for transport properties, Chung correlations are used for viscosity and thermal conductivity [21]. Non-linear mixing rules are used to calculate any mixture coefficients appearing inside the PR and Chung models. More details about all the aspects of the RFM thermophysical model can be found in reference [2].

It is important to emphasize that this modeling approach does not account for phase splitting when retrieving the temperature from the mixture enthalpy. Instead, it calculates the so called frozen adiabatic mixing temperature, without considering phase stability and splitting. In other words, phase volume fractions are not considered for the calculation of local temperature, which is an approximation in terms of model accuracy but it does not have any impact for the objective of the present work.

It is also worth stating explicitly that this paper does not focus on the physical accuracy of the numerical results; its sole objective is to evaluate the proposed NN model in contrast to the direct solution based on the cubic EoS for multicomponent systems.

2.2. Artificial Neural Networks (NN) and Multilayer Perceptrons (MLP)

The multilayer perceptron is a widely used NN model for function approximation consisting of multiple layers with each layer containing several neurons [22]. The multilayer perceptron models nonlinear relationships between input and output vectors, utilizing simple nonlinear functions in a feed-forward neural network [23]. The backpropagation algorithm is commonly used for training, adjusting weights based on the local gradient of the error surface. In summary, the network is initialized with weights, processes input vectors to generate output, calculates error signals, and iteratively adjusts weights to minimize errors until an acceptable level is reached [24].

In the present case, the neural network used consists of one input layer including 3 properties (composition, temperature, and pressure), two hidden layers with 100 number of neurons, and one output layer including one parameter, as shown exemplary in Figure 1.a for predicting density. Figure 1.b shows a typical structure of a neuron.

The mathematical representation of a single hidden-layer unit, often referred to as a single layer unit, is straightforward. It involves a linear combination of inputs processed through a nonlinear 'activation' function, typically a simple elementary mathematical function. In general, we will denote such units as follows [25]:

$$f_j^{\left(1\right)}\left(\mathit{x}\right)=a\left(v_j\left(\mathit{x}\right)\right)$$

(6)

$$v_j\left(\mathit{x}\right)={\mathit{w}}_j^{\left(1\right)T}\mathit{x}=w_{0,j}^{\left(1\right)}+\sum _{n=1}^N{w}_{n,j}^{\left(1\right)}{x}_n,j=1,\dots ,U_1$$

(7)

in which $f_j^{\left(1\right)}$ is the first layer output, with $j$ the index of each unit in the layer, $v_j$ is the linear combination, ${\mathit{w}}_j^{\left(1\right)}=\left[w_{0,j}^{\left(1\right)},w_{1,j}^{\left(1\right)},w_{2,j}^{\left(1\right)},\dots ,w_{n,j}^{\left(1\right)},\dots ,w_{N,j}^{\left(1\right)}\right]$ is the weight vector where $w_{0,j}^{\left(1\right)}$ is the bias and $w_{n,j}^{\left(1\right)}$ are the weights from each single input, and $\mathit{x}=\left[1,x_1,x_2,\dots ,x_n,\dots ,x_N\right]$ is the input layer. Then, $a\left(v\right)$ is the nonlinear activation function applied to the value $v$, which is Rectified Linear Unit (ReLU) activation function in our case. The following equation shows the ReLU activation function:

$$a\left(v\right)=ReLU\left(v\right)=max(0,v)$$

(8)

As expressed by the formula, ReLU acts as a feature detector. To make the concept of extending from a single layer to a multilayer neural network clearer, let us break down the sequence into two key operations: combining inputs linearly and passing the result through a nonlinear activation. This approach is essential for creating single-layer perceptron units [25]. The following equations show the output value for the second layer and the L^th layer respectively. $f^{\left(2\right)}$ and $f^{\left(L\right)}$ show the outputs of the second and the L^th layers.

$$f_j^{\left(2\right)}\left(\mathit{x}\right)=a\left(w_{0,j}^{\left(2\right)}+\sum _{i=1}^{U_1}{w}_{i,j}^{\left(2\right)}{f}_i^{\left(1\right)}\left(\mathit{x}\right)\right)$$

(9)

$$f_j^{\left(L\right)}\left(\mathit{x}\right)=a\left(w_{0,j}^{\left(L\right)}+\sum _{i=1}^{U_{L-1}}{w}_{i,j}^{\left(L\right)}{f}_i^{\left(L-1\right)}\left(\mathit{x}\right)\right)$$

(10)

2.3. Training Database

As training tool, Keras and Tensorflow libraries are used which are open-source libraries that provide a Python interface for NN [26]. To effectively train the network, it requires data featuring properties relevant to the main problem, such as consistent thermophysical properties.

Precision in tuning or optimizing parameters poses a challenge. Hyperparameters, especially the learning rate (α), play a crucial role, impacting computational cost, execution time, and network structure. The learning rate is vital in minimizing the loss function during backpropagation, where it iteratively updates biases and weights in hidden layers. This process, involving multiplying the learning rate by the output error of neurons, continues until an acceptable error is achieved for the specified data in each iteration, influencing prediction accuracy and generalization capability. As it is shown in Table 1, the considered inputs are composition ($Y_{H2}$), temperature ($T$), and pressure ($p$) by order, for predicting density ($\rho$), viscosity ($\mu$), thermal diffusivity ($\alpha$), enthalpy ($h$) and the compressibility ($\psi$). For predicting temperature ($T$) as output, composition ($Y_{H2}$), pressure ($p$) and enthalpy ($h$) are considered as the inputs. The full definitions and units of outputs are shown in Table 2.

In this specific test case, the properties exhibit highly non-linear distributions. Consequently, various data transformations are experimented to enhance network performance. Figure 2 illustrates that the data distributions for enthalpy and temperature are notably smoother, indicating fewer challenges in training these specific properties.

2.3. Model Training

In the following, all the methods tested to improve the efficiency of the network are elaborated. For all the properties, the min-max scalar is applied for inputs and outputs, to facilitate training. As mentioned before, ReLU is considered as the activation function. The optimizer algorithm is Adam, and the loss metrics is the Mean Squared Error (MSE). The final values of other hyperparameters like hidden layer size, epochs, and batch size is obtained by experience for each output parameter. Table 3 reports the number of epochs, and the loss metrics for each output parameter. Learning rate is 0.0015, the batch size is 256, and the hidden layer size is considered 100×100 for all of the outputs. The influences of choosing different hyperparameters to train the network can be seen in Figure 3 and Figure 4. As it is shown in Figure 3, by increasing the number of epochs the relative error decreases for the density, but for a greater number of epochs there will be no change, and also it is probable to end up in overfitting. Also, the effects of using different learning rates is shown in Figure 4. A good way for choosing the hyperparameters is to use the previous tested ones, or by experience.

The training algorithm aims to minimize the sum of the square errors [22]. This implies that the training minimizes absolute errors across all MLP outputs, regardless of relative errors. However, upon closer examination, as an example, within the compressibility output value range (1.12e-6 - 1.178e-3 kg/(m³ Pa)), cf. full-range chart in Figure 5, it becomes evident that relative errors are notably high, reaching up to 100 % in certain points. It is crucial to note that these suboptimal outcomes do not necessarily indicate inadequate MLP training. Instead, they underscore the importance of selecting an appropriate indicator of the MLP performance. The substantial relative errors for smaller target outputs could significantly impact the application of MLPs in real simulation [22].

For instance, in our model, temperature serves both as an input and as an output. A substantial relative error in temperature propagates errors in calculating other variables, leading to rapid divergence after a few time steps due to the cumulative effect. To mitigate this, it becomes imperative to diminish the prediction errors of NN for outputs with small magnitudes. One viable approach is to partition the data into ranges and train the network separately for each output within those ranges.

As illustrated in Figure 5, for the compressibility values within the range of 1.12×10^-6 -1.178×10^-3 kg/(m³ Pa), certain points, despite multiple attempts to optimize performance, exhibit relative errors reaching up to 100 %, as already mentioned, on the small value side. However, the objective is to achieve the least relative error. To address this, the data is split into two ranges: the first encompasses 90 % of the data, and the second covers the remaining 10 %. A log transformation is employed before scaling the data within 0-1, in the first sub-range (1.12×10^-6 - 6×10^-5 kg/(m³ Pa)). By normalizing the input and output parameters according to their respective characteristics, we can optimize the learning process, resulting in more efficient predictions and improved overall performance [27]. Also, the logarithmic distribution accommodates the specific characteristics and value ranges associated with output parameters, contributing to enhanced model performance and accuracy [27].

In addition, from the learning curves which are shown in Figure 6, the higher efficiency of this split approach can be seen. Learning curves serve as a frequently employed diagnostic tool for algorithms that iteratively learn from a training dataset. Throughout the training process, the model's performance is assessed on both the training dataset and a separate validation dataset. Plots illustrating the performance measurements are created to display these learning curves. Analyzing the shape and trends within a learning curve can provide insights into the behavior of a machine learning model and potentially offer suggestions for configuration adjustments to enhance learning and performance. The ultimate objective is to achieve a well-fitting model, indicated by a training and validation loss that decreases and stabilizes with minimal disparity between their final loss values. Typically, the model's loss will be lower on the training dataset than on the validation dataset. Consequently, it is expected to observe a gap between the learning curves of training and validation losses like the curve related to the first range of compressibility which is shown in Figure 6, but not too much loss like the curve related to the whole range of compressibility. So, for compressibility after data splitting the curves related to two different ranges show the training development. Learning curves can also be utilized to diagnose characteristics of a dataset, particularly in terms of its representativeness. An unrepresentative dataset implies that it may not accurately capture the statistical characteristics compared to another dataset from the same domain, such as a comparison between a training and a validation dataset. This situation often arises when the dataset's sample size is insufficient relative to another dataset. An example of this behavior is the curve of compressibility in the second range (6×10^-5 - 1.17×10^-3 kg/(m³ Pa)) in Figure 6, in which the oscillations of the loss show that the training data set is unrepresentative which means that it does not provide enough information for effective learning in comparison to the validation dataset used for evaluation [28,29].

However, this approach poses concerns as there are no predefined criteria for guessing a-priori the range of interest. While the analytical thermophysical library of OpenFOAM could be utilized to get the output value and make the choice, this is a heavy calculation, therefore it cannot be considered as a general method. Another possibility is using information provided by the NN models from a previous timestep, but finding a consistent criterion for subsequent iterations and time steps proves challenging. Another approach would be calling several NN per cell to begin with as it is mentioned in reference [22], in which the training data was clustered in to 400 subdomains, and each subdomain was fitted by an individual MLP. Data clustering using the K-mean method, which is a non-supervised clustering algorithm is done by Xi Chen et al. [29], .Consequently, this method is set aside, and an alternative approach involving the use of a single network for each property is adopted for the present CFD test case. In essence, for the accuracy and stability of the overall solution algorithm it has been experimented that it is easier to have just one network for each property, but with careful selection of the range needed.

As depicted in the Figure 7, the relative error consistently meets less than 1 % for each specified range. This successful outcome is observed across density, compressibility, diffusivity, and viscosity. Here the logarithm transformation is applied, together with an adapted data range close to the needs of the CFD test case (detailed in Section 3). Temperature and enthalpy result in low values of the relative errors without logarithm transformation. An essential lesson learned from experience is the careful selection of input and output ranges, particularly those proximate to the test case where the modules will be implemented. In these modules, various data manipulations such as min-max scaling and log transformation are performed. The interplay of these data treatments hinges on the chosen min-max values. Therefore, the judicious selection of a suitable range for training proves pivotal to the overall performance of the network, in the following implementation step within the OpenFOAM framework.

2.4. Implementation of the NN Models in OpenFOAM

In this study, Machine Learning is used to estimate changes in thermodynamic properties of fluids during the numerical simulation of injection in modern fuel systems and to create surrogate models for two phase flow predictions. To accomplish this, an NN is established for each fluid property and trained using data derived from dedicated codes, such as 0D thermodynamic simulation models. In the realFluidReactingNNFoam solver, evaluations via Python-trained NN models replace direct evaluations via the original thermophysical libraries, for each thermophysical and transport properties. To achieve the implementation of the NN Models in OpenFOAM, the “Neural Network Inference in C made Easy” (NNICE) code is employed for inferences [31] NNICE is a simple C++ library designed to provide neural network inference capabilities without the complexities associated with compiling traditional ML library C++ APIs. Initially developed for seamless integration into 3D CFD codes, this library is utilized in scenarios where calls to neural networks are made in a constrained parallelization environment [31].

Real-fluid thermodynamics and transport are integrated into a solver which utilizes the PIMPLE algorithm. This algorithm is a pressure-based segregated approach that combines elements of both SIMPLE and PISO methods [2]. To tackle severe pressure and velocity oscillations a modified reactingFoam solver in OpenFoam is used, called realFluidReactingFoam, which includes a modification of the pressure equation to include linear formulation of the density and changes in the frequency of updates of enthalpy and species [19]. This solver with direct evaluation of real-fluid properties will be considered as the reference.

In the following, a description of the overall CFD solution algorithm is provided, focusing on the points where properties need to be evaluated, and therefore when also NN inference replaces direct evaluation.

• In the initialization stage, $\left(Y_{H2},T,p\right)$ fields are read from dictionaries and values are used to infer each property of the flow through the Python-trained NN models via NNICE.
• Then, as shown in the flow-chart in Figure 8, at the beginning of a time-step, the continuity equation is first solved, and the PIMPLE iteration begins with the momentum predictor step.
• With entering the PISO loop the specie and the energy equations are solved. Then the thermodynamic properties, inferred through the Python-trained NN models with NNICE, replace the original thermo.correct() OpenFOAM function call. More in detail, temperature $T$ is first retrieved after the enthalpy equation is solved using $\left(Y_{H2},p,h\right)$ information; then all other properties (density $\rho$, viscosity $\mu$, thermal diffusivity $\alpha$, and compressibility $\psi$) are updated using most recent values of $\left(Y_{H2},T,p\right)$.
• Other inferences are then needed inside each PISO loop. Density is first updated with the corresponding NN model before the pressure equation is solved (details of the pressure equation can be found in Figure 8.b). Afterwards, density is explicitly updated solving the continuity equation, and, after checking the continuity error, the velocity field is updated.
• The last step with the PISO loop concerns a new recalculation of the density field via its NN model.
• Outside the PISO loop turbulence equations are solved, but no additional inferences are needed there. In addition, in this paper we do not include a turbulence model and assume a laminar flow, to show the application of the machine learning method in OpenFOAM.

2.5. CFD Case Setup

The assessment of the numerical framework is carried out on a cryogenic two-species mixing layer under trans -critical state, which incorporates liquid oxygen (LOX) and gaseous hydrogen (GH2). The case study has already been considered by many researchers, and serves a well-known benchmark test [2,9,32,33,34,35,36,37,38].

Specifically, we analyze a two-dimensional mixing layer with an injector lip that divides streams of liquid oxygen (LOX) and gaseous hydrogen (GH2), both at supercritical pressure. As shown in Figure 9, this setup serves as a general representation of a cryogenic coaxial injector, where a central dense oxygen jet interfaces with high-speed coaxial hydrogen flow, creating conditions to effective control mixing and flame stability in actual rocket engines [37].

The simulation setup involves an injector lip with a height (h) of 0.5 mm, segregating high-speed gaseous hydrogen (GH2) in the surroundings from dense liquid oxygen (LOX) in the center. A two-dimensional layout is considered. The area of interest of the computational domain spans 15h along the axial direction and 10h in the transversal direction. An additional 15h streamwise region serves as sponge layer, designed to mitigate pressure oscillations arising from the outlet. For this study, a relatively coarse mesh is employed with a uniform grid spacing of ∆x/h = ∆y/h = 20. The specified resolution is applied up to 15h in the x-direction and extends to 2.5h on both sides of the lip center in the y-direction. The remaining downstream section (from 15h to 30h) of the domain undergoes stretching in the axial direction with a global expansion ratio of 10 [2]. The boundary conditions and computational domain are based on the description provided in [37], but unlike the original benchmark description, no turbulence modelling is employed. As a matter of fact, the present calculations are not meant for assessing model accuracy, but for proving the effectiveness of the Neural Network approach for fast evaluation of real-fluid properties, therefore direct numerical simulation, not including turbulence model, makes the test case even stronger and more probing.

At the injector lip, we implement an adiabatic no-slip wall condition, while the top and bottom boundaries are treated as adiabatic slip walls. Furthermore, we apply a Dirichlet pressure boundary condition at the outlet [2].

As discussed in reference [2] regarding the time discretization, a first-order Euler scheme is applied, while second-order accurate Gauss limited-linear schemes are employed for both advection and diffusion.

3. Results and Discussion

This paper extensively explains the development of a model that seamlessly combines real-fluid model (RFM) with NN. In terms of surrogate models for approximating thermodynamic functions, NNs demonstrated good performance in predicting properties as a function of pressure, temperature, and mass fraction of a two-component mixture. When used in simulations, since the NN regression requires minimal time during evaluation, the computational cost is small. On the other hand, the only alternatives are either evaluating the complex thermodynamic models through the RFM, which is rather demanding on the fly, or performing interpolations from tables, the latter being the only practical way in complex cases. However, tables have a large storage footprint and become very cumbersome as the interpolation dimensionality increases (interpolating data among pressure, temperature and the mass fractions of multiple components) [12,27]. Therefore, NN seems to be a very attractive option, as the trained network has the size of KBs, which is over a thousand or million times smaller than a table. Further applications where NN can be used, is the development of surrogate models that can be trained using existing, validated data. In the current demonstration case, networks with two hidden layers and 100 neurons each can describe adequately well the properties as a function of pressure, temperature and mass fraction for a LO2/GH2 transcritical mixing layer. It is important to mention that the training for a single vector output network takes about 30-60 minutes (depending on the number of epochs). Its evaluation happens instantaneously much faster than any simulation could produce.

Another noteworthy aspect is the remarkable versatility of ML. The methods discussed here, even when employed in tandem with simulations, exhibit flexibility in terms of input types. ML can effectively handle diverse forms of data presented in suitable formats. It is crucial to emphasize that this work does not assert to have exhaustively explored the full potential of ML. Admittedly, there are numerous avenues for integrating ML, given its expansive range of techniques adaptable to a multitude of problems. It is only very recently that such ML techniques have begun to be explored in the context of multiphase flows [10].

The open-source code OpenFOAM-v2112 is employed to build distinct solvers, whose results and performance are then presented and discussed. Specifically, such CFD model variants are termed as follows.

• The original CFD model is referred to as realFluidReactingFoam solver. As already explained, it is a real-fluid pressure-based multi-species solver using PR-EoS and Chung models with non-linear mixing rules coded in OpenFOAM [2,10,19,38].
• The second CFD model is referred to as realFluidReactingNNFoam_1 which adopts a neural network approach (here referred to as NN_1 variant) wherein Python-trained NN models and NNICE inference are used to replace original calculations through the RFM.
• The third CFD solver is referred to as realFluidReactingNNFoam_2 which is similar to the second model, but incorporates log transformations for density, compressibility, viscosity, and thermal diffusivity to achieve target accuracy criteria (here referred to as NN_2 variant).
• The last model is referred to as realFluidReactingNNFoam. It is like the previous one, but using the adapted data range, which means that the data range for training is chosen to be as close as possible to the case study needs, to improve networks accuracy (here referred to as NN variant).

We start by showing in Figure 10 a comparison of the results obtained at 0.1 ms by the codes above listed, namely, the NN variants vs. the original code (realFluidReactingNNFoam_1, realFluidReactingNNFoam_2, realFluidReactingNNFoam, and the reference realFluidReactingFoam). Contours are shown for H2 mass fraction, temperature, velocity, and density, at 0.1 ms. Overall, all fields appear reasonable, and the three ML variants have similar features compared to the original code. Of course, no exact correspondence can be expected, as each simulation produces its own instantaneous results. However, by further scrutinizing these results, it is shown that the first neural network variant (NN_1) starts generating noisy fields for temperature, velocity and density, which are attributed to insufficient ML model accuracy. The second variant (NN_2), with transformed input variables, improves the velocity field quality but does not improve other fields. Conversely, the final neural network variant (NN) preserves clean and smooth fields for all quantities.

The primary objective of this work is to demonstrate the impact of employing NN methods in enhancing real-fluid model simulations by substituting the NN algorithm for the direct RFM model. Considering the final and more accurate solver version, named realFluidReactingNNFoam which as explained uses a tailored dataset for training, results reported in Table 4 show that the code based on neural networks is about 3 times faster on the total execution time, and about 7.5 faster on a time-step basis. This performance is measured considering 0.1 ms of simulated physical time on one compute node with 20 cores, as the test case is not very large in terms of mesh size, to avoid including scalability aspects related to node communication. The tolerances for solving governing equations are kept the same, and in particular 1e-8 is used for the velocity field, which is quite strict. The NN version of the code requires some more iterations per time-step to reach convergence, so despite more iterations are needed in general each one is much faster. Further work can be done on further improving the accuracy of the NN models and optimize the overall code coupling, so that the number of iterations does not increase, to get an even better speed-up, but the current result appears already promising.

To assess the simulations, result quality, comparisons with those from the original approach proposed in [2,20] are presented. Figure 11 shows a scatter plot of mixture temperature vs. hydrogen mass fractions for an instantaneous field, at 1 ms, which based on the low speed stream velocity (30 m/s) corresponds already to 4 flow through times (FTT) over a distance of 15h (cf. Figure 9). The NN solver reproduces exactly the thermodynamic states of the reference model, without any appreciable discrepancy.

The corresponding instantaneous fields for H2 mass fraction, temperature, velocity, and density at 1 ms are shown in Figure 12, for the NN solver realFluidReactingNNFoam and the original solver realFluidReactingFoam. The current numerical solution effectively captures the three prominent vortical structures in the velocity fields, anticipated within x/h = 10, along with the steep density gradient. Kelvin-Helmholtz mechanisms initiate the formation of the initial eddies in the mixing layer, downstream of the lip, prominently on the hydrogen side. These structures, due to interfacial instabilities, contribute to the development of thicker vortices in the oxygen stream [2]. Results of the NN solver are in good agreement with those obtained by the direct RFM solver. These findings showcase the proficiency of the current numerical NN framework in managing mixing processes involving real-fluid thermodynamic and transport properties, under transcritical conditions under the assumption of frozen temperature [21].

Additionally, as a further assessment, a quantitative comparison of statistics collected over time for axial velocity, temperature, and density are presented in Figure 13. These transversal profiles of mean and fluctuating root-mean-square statistics (rms) are obtained by time-averaging instantaneous fields from 1 to 4 FTTs, and are taken at various axial locations, namely at x/h = 3, 5, and 7. Overall, there is very good agreement with the reference trends provided by the realFluidReactingFoam solver.

5. Conclusions

This paper presents a comprehensive exploration of the integration of NN models within the RFM approach for efficiently simulating complex thermodynamic functions in a CFD code, for non-evaporating binary mixture under transcritical conditions. The study demonstrates the efficiency of NNs as surrogate models, emphasizing their ability to approximate complex thermodynamic properties with very affordable training time compared to traditional dataset derivation and significant savings during CFD runtime. The paper highlights the advantages of using NN models in CFD simulations, where computational costs, though present during evaluation, are notably smaller than the original method using direct evaluation from the equation of state. The versatility of machine learning (ML) is underscored, particularly in handling various data types and distributions. The work acknowledges the vast potential of ML in the context of multiphase flows and suggests numerous avenues for further exploration. The application of the developed NN models in OpenFOAM to simulate a transcritical mixing case reveals the potential of neural networks to replace traditional RFM models, offering faster execution times and very good fit, as evidenced by contour plots and averaged transverse profiles. The comparison of realFluidReactingNNFoam results with the original realFluidreactingFoam indicates that, despite the latter requiring fewer iterations, the NN model outperforms in terms of execution time at the same level of solution accuracy.

In summary, the successful application of NN in conjunction with RFM in simulating real-fluid thermodynamic and transport properties signifies a promising direction for future research in the field of accelerating computational fluid dynamics models. The study opens avenues for exploring the full potential of ML techniques in addressing complex problems related to multiphase flows with more than just two species.