Data-driven nonlinear constitutive relations for rarefied flow computations

To overcome the defects of traditional rarefied numerical methods such as the Direct Simulation Monte Carlo (DSMC) method and unified Boltzmann equation schemes and extend the covering range of macroscopic equations in high Knudsen number flows, data-driven nonlinear constitutive relations (DNCR) are proposed first through the machine learning method. Based on the training data from both Navier-Stokes (NS) solver and unified gas kinetic scheme (UGKS) solver, the map between responses of stress tensors and heat flux and feature vectors is established after the training phase. Through the obtained off-line training model, new test cases excluded from training data set could be predicated rapidly and accurately by solving conventional equations with modified stress tensor and heat flux. Finally, conventional one-dimensional shock wave cases and two-dimensional hypersonic flows around a blunt circular cylinder are presented to assess the capability of the developed method through various comparisons between DNCR, NS, UGKS, DSMC and experimental results. The improvement of the predictive capability of the coarse-graining model could make the DNCR method to be an effective tool in the rarefied gas community, especially for hypersonic engineering applications.


Introduction
Due to non-equilibrium effects inside the shock wave with large macroscopic gradients on a scale of several mean free paths, traditional hydrodynamic methods such as Navier-Stokes equations are no longer valid. To resolve the steady shock wave structure of monatomic and diatomic gases physically, particle methods such as Direct Simulation Monte Carlo (DSMC) [1] and numerical schemes for Boltzmann equations have been extensively employed. However, the computational cost for both microscopic method and mesoscopic method is large in comparison with hydrodynamic methods, especially in the near-continuum flow regime with multiscale characteristics. The same defects are also encountered in multi-dimensional rarefied flows, especially hypersonic flows around a space vehicle. To overcome current defects, the most popular numerical methods which have been employed in shock wave structure or rarefied hypersonic flow computations are reviewed subsequently. DSMC is a stochastic particle method to simulate rarefied non-equilibrium flows by tracking and counting representative molecules per cell. Not only for one-dimensional shock wave structure calculations, it has been widely applied in engineering applications after proposed because it could be easily implemented in the hypersonic reentry flow with thermal and chemical non-equilibrium effects and yields accurate results when the Knudsen number is large. However, the computational efficiency decays significantly in near-continuum flow or high-density regions in a multiscale flow field due to the surging number of particles per cell and limitations of time step and mesh size. In the mesoscopic method, the Boltzmann equation describes the evolution of gas velocity distribution function in theory. However, the complex collisional term in the right-hand side prevents the full Boltzmann equation from being solved widely. Linearized or simplified Boltzmann model equations such as the Bhatnagar-Gross-Krook (BGK) [2], Shakhov-BGK [3] and ellipsoidal statistical BGK (ES-BGK) [4] model equations were proposed to reduce the numerical complexity. In order to get the solutions of these approximations, discrete velocity and ordinate methods (DVM) such as gaskinetic unified algorithms (GKUA) [5] are widely employed and a new framework named unified gas kinetic scheme (UGKS) [6] proposed by Xu and collaborators a few decades ago is regarded as a multiscale scheme for both hydrodynamic and rarefied regimes. In comparison with other kinetic solvers, UGKS could reveal accurate solutions in continuum flow regime with a large cell size and time step, which overcomes the cell size barrier requiring a mesh size on the order of the mean free path for DSMC and other kinetic schemes. However, to represent the velocity distribution function in velocity space, all these discrete velocity methods need huge memory consumption, which is too expensive and sometimes inaccessible in hypersonic flows. On the other hand, hydrodynamic methods also attract a lot of attention by extending its covering range in various approaches. Burnett equations [7] derived from the second order Chapman-Enskog expansion of VDF are classical high-order moment equations in literature. But the original formulation was proved unstable due to violation of the second law of thermodynamics, although some remedies have been reported in recent years [8]. Besides Burnett equations, Grad equations [9] and regularized moment equations, e.g. R-13 [10] and R-26 [11] equations, and other moment equations form the community of extended hydrodynamic equations (EHE). Eu [12] proposed generalized hydrodynamic equations (GHE) from a different viewpoint in 1992 and developed nonlinear coupled constitutive relations (NCCR) with Myong [13,14] based on GHE subsequently. As a simplified approximation of GHE, NCCR includes the evolution equations of conservative variables and nonlinear coupled closure of stress tensor and heat flux, which is strictly consistent with the second law of thermodynamics [15]. Except DSMC, Boltzmann equation and moment equation solvers, many coupling models have also been proposed recently such as traditional CFD-DSMC hybrid method [16], unified gaskinetic wave-particle method (UGKWP) [17], general synthetic iterative scheme (GSIS) [18], unified stochastic particle method based on the BGK model (USP-BGK) [19] and simplified unified wave-particle method (SUWP) [20]. These methods try to couple hydrodynamic and microscopic descriptions of the rarefied flow to construct an efficient and accurate hybrid scheme.
The study of shock waves by the above methods has been widely carried out to compare with experimental results in recent decades. Both monatomic and diatomic gases are investigated in a wide range of upstream Mach numbers. Thermal non-equilibrium in shock waves with both translational and rotational temperatures is also investigated. In general, kinetic methods yield much more accurate macroscopic parameter profiles than hydrodynamic solvers especially at high Mach numbers. However, the computational cost of solving various Boltzmann equations is much larger than extended hydrodynamic equations due to the discretization in velocity space separately. Unfortunately, this deficiency will become unaffordable in three dimensional engineering applications and that could be the main reason why mesoscopic methods are still not successfully applied in predictions of aerodynamic and thermal loads on hypersonic vehicles.
While traditional fluid dynamics mainly focused on developing new physical models to increase their accuracy, another effective way arising in recent years to improve the predictive capability is data-driven modeling method. Especially in turbulence modeling, many research studies have been carried out by utilizing high-fidelity direct numerical simulation (DNS) or experimental data to construct a predictive data-driven Reynolds-averaged Navier-Stokes (RANS) model. Duraisamy et al. [21] summarized and reviewed these methods in detail and Xiao [22,23] proposed a physical-informed machine learning approach (PIML) for reconstructing Reynolds stress modeling discrepancies based on DNS data. Ling et al. [24] proposed a tensor basis neural network (TBNN) by combining traditional neural network with Pope's Reynolds averaged turbulence model [25]. Using TBNN, they successfully predicted the separation phenomenon of angular vortex and wavy wall flow in pipe flow. Not only in turbulence modelling, Kutz [26] reviewed deep neural networks (DNNs) in fluid dynamics and predicted that it was only a matter of time before deep learning made its mark in the general area of high-dimensional, complex dynamical systems. In fact, data methods have been utilized widely in computational fluid dynamics successfully such as proper orthogonal decomposition (POD) or dynamic mode decomposition (DMD), which reduces computational costs significantly. Sekar [27] proposed a fast prediction method for airfoil flow field based on deep convolutional neural networks (CNN) and multi-layer perceptron (MLP). Bar-Sinai [28] proposed an optimal approximation method of partial differential equations based on the actual solutions of known flow governing equations, using neural networks to estimate the spatial derivatives. Mao [29] successfully modeled the high-speed air flow based on Euler equation by using the physical-informed neural networks (PINN) [30], and solved the one-dimensional Euler equation by using smooth solution and two-dimensional oblique shock wave solution with contact discontinuity. Raissi [31] tried to combine the PINN method with the convection diffusion equation of particles to propose the HFM method, which describes the technical characteristics and application prospects of the PINN method. All in all, using limited high-fidelity data to capture features that dominate physics and then improve the predictive capability of coarse-graining models will provide an effective tool in the fluid community. In rarefied gas dynamics regime, Zhang [32] utilized machine learning methods to derive various macroscopic governing equations for fluid dynamics and all training data comes from DSMC simulations. It is proved that the data-driven method to discover governing equations for a non-equilibrium flow is promising.
In this paper, data-driven nonlinear constitutive relations (DNCR) will be developed by machine learning method to provide a predictive approach for rarefied nonequilibrium flows. Both NS equations solver and UGKS solver provide training data during training phase to map feature vectors and responses of stress tensor, heat flux and flow parameter derivatives between NS and UGKS data. Based on this off-line model, test cases excluded from training data set could be predicated rapidly and accurately by solving conventional equations with modified stress tensor and heat flux. The whole paper will be divided into four sections. "Data-driven nonlinear constitutive relations" will provide a detailed construction process of data-driven non-linear relations and the solving methodology. Numerical cases of shock wave structure and hypersonic flow passing through a blunt circle cylinder are presented in "Numerical simulations" to validate and assess the capability of the current DNCR method. And finally, conclusions and future work are summarized in the last section.
2 Data-driven nonlinear constitutive relations

Data-driven modeling approach
Data-driven approach is a suitable tool for building high-dimensional and complex regression functions between input features and output responses. Constitutive relations of macroscopic equations in rarefied regime are just impossible to specify analytically in math despite the fact that traditional moment equations try to extend their covering range by various expansion or simplification of velocity distribution function (VDF). Therefore, the intention to obtain non-equilibrium constitutive relations through a data-driven modeling approach is the main task for this paper.
In order to construct regression functions, we need to obtain training data first before adopting a machine learning method. Two numerical methods including NS equations solver and UGKS solver will be employed in current research to account for coarsegraining model and high-fidelity model respectively. The covering range or size of training cases is also important for the capability of the regression functions. Both local continuum and rarefied regions should be included inside the training cases to make it comprehensive. Another crucial issue is the choice of flow features as regression inputs. Due to non-equilibrium characters of target flows, we will choose local Knudsen number Kn GLL (ρ, T) and gradients of pressure ∇p, density ∇ρ and other favorable flow-field variants. Following the basic principles of a data-driven method, the input parameters should be obtained from the coarse-graining model directly, such as NS solver in current framework, because the high-fidelity data from UGKS solver is unavailable in the test cases. Secondly, the physical coordinates x should not be chosen as features in the DNCR method to overcome limitation that the test case should have the same geometry and grid as the training case. And only local flow quantities are considered as well-chosen input parameters. All these adopted features in current research strictly satisfy the above principles.

Evolution equations of DNCR method
The basic macroscopic conservation laws of generalized hydrodynamics can be summarized as where Π, Q are stress tensor and heat flux, u are macroscopic velocity vectors. To close above conservation laws, constitutive relations that describe how the stress tensor, heat flux and mass diffusion respond to the gradients of velocity, temperature and mass fraction are needed. In the continuum region, the Navier-Stokes-Fourier constitutive relations provide a typical first order mathematical model, where the stress tensor and heat flux are expressed in terms of the first-order macroscopic velocity and temperature. However, NS equations are only valid in the continuum regime due to its intrinsic continuous hypothesis.
Beyond the constitutive relations, the conservation laws are regarded as a valid model in all scales such as continuum, transitional and rarefied flows. If we could find a set of proper constitutive relations in a desired flow, even discrete relations, the macroscopic quantities could be predicted physically by solving the closed mathematical model. However, the aim of finding analytical constitutive relations to describe the nonequilibrium characteristic becomes quite complicated for a complex system, where too many parameters need to be fitted, and the physical meaning of these terms becomes obscure. That's a key point for us to figure out how to obtain these exact constitutive discrete relations as Eq. (2) through the machine learning method.
Once the discrete constitutive relations in every grid cell center are obtained, the conservation equations could be solved to convergence by explicit or implicit timemarching method. Take a two-dimensional flow as an example, the stress tensor and heat flux substituted into the macroscopic conservation laws are obtained as Eq. (3), where the local flow parameters with a superscript tip come from the trained off-line model. With the constraint of above Eq. (3), the stress tensor and heat flux will automatically regress to the values of the off-line model as the gradients of flow parameters converge to a steady state. It is worth noting that the output of the trained off-line model not only includes stress tensor and heat flux, but also contains some gradients of key flow parameters. During the numerical simulation process, an identical multiblock structure grid is employed currently for both numerical solvers in finite volume method to guarantee the data structure of the flow field consistent.

Solution methodology
2) Establish the feature vector field q(x), e.g., local Knudsen numbers, pressure gradient, density gradient and flow parameters, based on the NS-predicted flow fields in Step 3) Construct regression functions f : q ↦ ΔΦ through a machine learning method based on the training data. The extremely randomized trees model [33] is employed in current research to train the regression relations.
The choice of the feature vectors is a key step for constructing regression functions and will affect the final precision of prediction greatly. The input flow features and output responses in the regression and descriptions are presented in Tables 1 and 2 as below.
In the above training procedure, the training set has to include different flow types to guarantee the applicable capability of the model. However, in the current framework of DNCR, only a few different typical inflow conditions are needed due to the multi-scale feature of transitional flow. Every grid cell will provide a basic training set and the local non-equilibrium degree will be totally different inside a multi-scale flow field. In other words, it is even possible to predict a flow field utilizing the DNCR method based on a training set that   That means if the training flows consist of various elementary cases such as compression and expansion flows in transition regime, complex engineering flow which has a subset or all of them could be predicted well ultimately. This challenging task is the final goal of this method and is outside the scope of this study. We mainly focus on the demonstration of feasibility first and will achieve that goal gradually. 3) Solve conservation equations with baseline stress tensor, heat flux and flow parameter gradients in Step 2 until converging to a steady state and finish the whole procedure.
Finally, the test flow field with high-fidelity will be obtained without carrying out complex and expensive simulations through UGKS or DSMC methods. The whole flow diagram is shown in Fig. 1 with both training phase and predicting phase.
The above predicting phase is separate from the training phase because the maps between input features and target responses have already been established before prediction. Only NS feature vectors are needed in order to obtain the final discrete stress tensor and heat flux for test flow. Unlike linear NSF relations or nonlinear Burnett constitutive relations, there is no explicit formula in DNCR to close the conservation equations. And actually, the universality and simplicity of explicit constitutive relations to describe non-equilibrium phenomenon is also questionable. However, the constitutive Combined with the final baseline stress tensor, heat flux and flow parameter gradients, solving convention equations is much more efficient than UGKS or DSMC solvers because the total computational time is on the same order of magnitude as NS equations solver. The goal of combining the accuracy of microscale description and the efficiency of macroscale equations could be achieved, especially for hypersonic rarefied engineering flows.

One-dimensional shock structure
In this section, the monatomic gas shock structure with non-equilibrium effects limited to the translational energy mode is presented. Shock structure computation mainly focuses on the accuracy of the non-equilibrium model without considering complex wall boundary conditions first. The viscosity of Argon is calculated by using the power law: where s = 0.72 and other physical properties are given in Table 3.
The training set in one-dimensional shock structure is constructed by different shock wave cases with Mach numbers ranging from 1.2 to 10 which includes Ma = 1.2, Ma =  To compare the discrepancy of NS, UGKS and DNCR results, the distributions of temperature and density inside the shock wave are plotted for Mach 1.6 in Fig. 2, where the x coordinates are nondimensionalized by the upstream molecular mean free path. Due to the low Mach number condition, the gradient of flow parameter inside the shock is not too steep, which results that the local distribution function is not far from equilibrium state. In spite of that, the DNCR results are still in good agreement with that of UGKS, especially high-order moments such as stress tensor and heat flux as shown in Fig. 3. It is also demonstrated that the local stress tensor and heat flux in DNCR solver could converge to the values of UGKS as Eq. (3) when the steady state is reached.
When the Mach number grows to 8, the gap between NS and UGKS results becomes enormous as shown in Figs. 4 and 5. The non-equilibrium characteristic inside the strong shock wave makes Navier-Stokes-Fourier constitutive relations no longer validate. To overcome this defect, various numerical algorithms from the Boltzmann equation or statistical method were proposed in literature to describe the flow structure inside the shock. UGKS is one of these successful implementations and its accuracy has been demonstrated in one-dimensional shock wave computations. As shown in Figs. 4 and 5, the DNCR method yields comparable results with that of UGKS. The temperature rising point predicted by DNCR and UGKS is about 4 times earlier than that of NS, which indicates a thicker and weaker shock compression process in reality. To assess the applicability of DNCR for different Mach number conditions, Fig. 6 presents the comparison of inverse density thickness of shock wave with experimental data using different numerical methods. Both UGKS and DNCR results are in good agreement with most experimental data with Mach number ranging from 1.2 to 10.

Hypersonic flow around a blunt circular cylinder
The hypersonic flow of Argon passing through a blunt cylinder at Ma = 10 is considered in this section. In this case, the radius of the cylinder is 0.5 m and the computational domain is divided with 100 × 120 (normal direction × circular direction) structured cells as Fig. 7. In UGKS computations, the velocity space is discretized with rectangle cells of 200 × 200. The detailed inflow conditions to be predicted are listed below.  out to obtain all these features and output responses in order to train the nonlinear constitutive relations model. Finally, the DNCR solver calculates the case with conditions as Eq. (5) and the comparisons are presented here to assess the improvements of the DNCR method.
Flow parameter contours such as pressure, temperature and velocity in x-direction are presented in Figs. 8, 9 and 10. It is obvious that DNCR yields a very similar flow field as UGKS does under these inflow conditions. All rarefied characteristics such as the thicker shock wave and the weaker compression behind the bow shock are captured accurately by the DNCR method. By checking the inflow Knudsen number with a characteristic length of 1.0 m, which is the diameter of the circle, the NS solution will be no longer validated in this case because Kn ∞ = 0.085. This could also be observed clearly in the figures below where the NS solutions are different from the other two approaches.
To investigate the flow field more carefully, the distribution lines of pressure, temperature and x-velocity along the stagnation line are plotted in Figs. 11, 12 and 13. These lines indicate that the results of DNCR are in good agreement with that of DSMC and the discrepancies between NS and DNCR exist from the bow shock to the solid surface. In order to quantify the improvement degree of the proposed DNCR method, a parameter ξ would be introduced as Eq. (6).
where y denotes a flow parameter and the superscript i means the index of the grid cell. As we can see from the definition, the formulation 1 -ξ indicates the ratio of the discrepancy between DNCR and UGKS to the discrepancy between NS and UGKS. Therefore, the number of ξ ranges from 0 to 1 which indicates that the DNCR results converge to NS or UGKS solutions respectively. In this case, the ξ of pressure, temperature and x-velocity are 0.4908, 0.7680 and 0.6614 respectively and the improvement of temperature is more profound than the other two parameters. Besides physical parameters in the flow field, the heat flux and shear stress distributions on the wall are also investigated in the test case to validate our approach. Due to the symmetry, only the shear stress coefficient and heat flux on the upper surface are plotted in Figs. 14 and 15.  The results of heat flux and shear stress on the wall yielded by the DNCR model are still in good agreement with that of UGKS, which indicates that the DNCR method has a promising capability in engineering application to predict aerodynamic heating and drag for hypersonic vehicles.
To check the computational efficiency of different numerical methods, the time and memory consumption is listed in Table 4 as below.
The computational time in Table 2 denotes the time to get the converged steady state solution. It is worth noting that the training time of DNCR is not included in the above computational time because the final aim of any machine learning algorithm is to establish an off-line trained model which could be utilized in different configurations or inflow conditions. Therefore, it is not convenient or necessary to account for the whole training time in a single computational case to be predicted. For the implicit UGKS computation, the computational time is about 17 h and for DNCR is around 13 min. The memory cost of UGKS is around 54 GB and for DNCR it is about 20 MB. Both simulation time and memory consumption of the DNCR solver are comparable to NS because the framework of both methods is very similar. To compare the computational efficiency of DNCR and UGKS, the time consumption of DNCR is about 78 times faster than the implicit UGKS method, and memory cost is 2760 times less than UGKS.

Conclusions
This paper describes the development of a data-driven non-linear constitutive relation model for multiscale flow simulation based on the extremely randomized trees model. Based on the training set constructed by the NS and UGKS data, an off-line regression model between input features and output responses is established first. In the prediction phase, NS simulation is carried out to compute the feature vectors, which could be utilized to yield baseline stress tensor, heat flux and flow parameter gradients by querying the trained regress functions. With this baseline flow field parameters, the macroscopic conservation laws are solved implicitly to obtain the final results of the test flow. 1-D shock wave structure and 2-D hypersonic flow around a blunt circle cylinder are investigated by the DNCR method to assess the capability, accuracy and efficiency. The simulation results of DNCR are in good agreement with that of UGKS in all test cases when the distribution of flow parameters and high-order moments are compared. And the computational efficiency of DNCR is much higher than UGKS but with less memory consumption. As a result, the developed DNCR method is proved to be a promising method and expected to pave a new way to predict complex three dimensional hypersonic rarefied gas flows in future.