2.1. Direct simulation Monte Carlo

Generally, the Knudsen number ($Kn$), which is defined as the ratio of the molecular mean free path to the characteristic length scale of system, is used to classify various flow regimes, with the continuum regime being in the range of $Kn \le 0.01$, the slip regime being $0.01 < Kn \le 0.1$, the transition regime being $0.1 < Kn \le 10$ and the free molecular flow being $Kn > 10$. While the computational fluid dynamics (CFD) methods based on the NSF equations have been successfully applied to the continuum regime and can also be used for the slip regime with appropriate slip boundary conditions, they cannot accurately simulate transitional and free molecular gas flows. On the contrary, the Boltzmann equation can describe gas flows in all flow regimes, which can be written as

(2.1)\begin{equation} \frac{{\partial f}}{{\partial t}} + {\boldsymbol{c}}\frac{{\partial f}}{{\partial {\boldsymbol{x}}}} + {\boldsymbol{G}}\frac{{\partial f}}{{\partial {\boldsymbol{c}}}} = \boldsymbol{S}. \end{equation}

Here, $f$ is the molecular velocity distribution function, which represents the density of molecules in phase-space (i.e. the space spanned by position and velocity). Here $\boldsymbol {c}$ is molecular velocity, $\boldsymbol {G}$ is external force and $S$ is a complex integral term, which represents the change of distribution function due to the interaction between molecules. Owing to the presence of $S$, the Boltzmann equation is difficult to solve numerically.

The DSMC is a stochastic molecule-based method that approximates the molecular velocity distribution function in the Boltzmann equation with simulation molecules (Oran, Oh & Cybyk Reference Oran, Oh and Cybyk1998). It has been demonstrated that DSMC converges to the solution of the Boltzmann equation in the limit of a large number of simulation molecules (Wagner Reference Wagner1992). In DSMC, each simulation molecule represents $F$ real molecules and $F$ is the so-called simulation ratio. The DSMC tracks the simulation molecules as they move, collide with other molecules and reflect from boundaries. The macroscopic physical quantities are obtained via sampling corresponding molecular information and making an average at the sampling cells. Specifically, density ($\rho$) and macroscopic velocity ($v_i$) are computed with

(2.2a,b)\begin{equation} \rho =\frac{mFN_{p}}{V_{cell} },\quad \mathrm{and}\quad v_{i} =\frac{1}{N_p}\sum_{cell}c_i, \end{equation}

where $m$ is molecular mass, $N_p$ is the number of simulation molecules in the sampling cell, and $V_{cell}$ is the volume of the sampling cell. The velocity of each molecule can be regarded as a sum of two parts, i.e. the macroscopic velocity of the cell to which it belongs and the molecular thermal velocity (which is defined as $C_i=c_i-v_i$). Based on the molecular thermal velocity, temperature ($T$) is computed with

(2.3)\begin{equation} T=\frac{1}{3k_{B}N_{p}}\sum_{cell}m\boldsymbol{C}^{2}, \end{equation}

where $k_B$ is the Boltzmann constant, and $\boldsymbol{C}^2=C_1^2+C_2^2+C_3^2$. Pressure tensor ($p_{ij}$) is computed with

(2.4)\begin{equation} p_{ij}=\frac{mF}{V_{cell}}\left ( {\sum_{cell}c_ic_j}-N_pv_iv_j \right ), \end{equation}

and is generally split into its trace and trace-free parts, i.e. pressure ($p$) and stress tensor ($\tau _{ij}$), $p_{ij}=p\delta _{ij}+\tau _{ij}$. Here, $\delta _{ij}$ is Kronecker delta (Heinbockel Reference Heinbockel2001). Therefore, pressure ($p$) and stress tensor ($\tau _{ij}$) are computed with

(2.5a,b)\begin{equation} p=\tfrac{1}{3} (p_{11}+p_{22}+p_{33})\quad \mathrm{and}\quad \tau_{ij}=p_{ij}-p\delta_{ij}. \end{equation}

The heat flux ($q_i$) is computed with

(2.6)\begin{equation} q_i=\frac{F}{2V_{cell}}\sum_{cell}m\boldsymbol{C}^{2}C_i . \end{equation}

The local viscosity coefficient is computed with

(2.7)\begin{equation} \mu =\mu _0\left ( \frac{T }{T _0} \right)^\omega, \end{equation}

where $\mu _0 = 2.117\times 10^{-5}\ \mathrm {kg}\ {\rm m}^{-1}\ {\rm s}^{-1}$ is the reference viscosity at the reference temperature $T_0 = 273\ \mathrm {K}$. Here $\omega$ is the viscosity exponent and is determined by the molecular model employed. For the hard sphere (HS) model (Bird Reference Bird1994), the Maxwell molecule model (Bird Reference Bird1994) and the variable hard sphere (VHS) model (Bird Reference Bird1981), $\omega$ is equal to 0.5, 1.0 and 0.81, respectively. The local heat conductivity coefficient is computed according to the local viscosity coefficient using

(2.8)\begin{equation} \kappa =\frac{5}{2}\frac{\mu}{Pr}, \end{equation}

where $Pr$ is the Prandtl number and is equal to $2/3$ for a monatomic gas.

For a specific application, DSMC first initializes the simulation molecules according to the initial distributions of density, temperature and macroscopic velocity. Then, the molecular motions and intermolecular collisions are sequentially conducted in each computational time step. The molecular motions are implemented in a deterministic way. Each molecule moves ballistically from its original position to a new position, and the displacement is equal to the product of its velocity times the time step. If the trajectory crosses any boundaries, an appropriate gas–wall interaction model would be employed to determine the reflected velocity. Common gas–wall interaction models include specular, diffuse and Maxwell reflection models. In our simulations, the diffuse model is employed. The intermolecular collisions are implemented in a probabilistic way. Among the several algorithms for modelling intermolecular collisions, the no-time-counter technique (Bird Reference Bird1994) is the most widely used. In the no-time-counter, any two molecules in the same computational cell are selected as a collision pair with a probability that is proportional to the relative speed between them. Then, the post-collision velocities of molecules are determined depending on the molecular model employed. In this work we employ the HS model for the first two benchmark cases (diffusion flow and Taylor–Green vortex), the Maxwell molecule model for the shock wave case and the VHS model for the Poiseuille flow case. It is noteworthy that because the molecular motions and intermolecular collisions are conducted in a decoupled manner, the time step in DSMC needs to be smaller than the molecular mean collision time, and the cell size for the selection of collision pairs needs to be smaller than the molecular mean free path.

The DSMC has been successfully applied to simulate gas flows in the whole flow regimes, and the results have been well validated (Stefanov, Roussinov & Cercignani Reference Stefanov, Roussinov and Cercignani2002; Sun & Boyd Reference Sun and Boyd2002; Zhang, Fan & Fei Reference Zhang, Fan and Fei2010; Gallis et al. Reference Gallis, Bitter, Koehler, Torczynski, Plimpton and Papadakis2017). The advantage of the DSMC method is that there is no need to assume any macroscopic governing equations a priori. The macroscopic quantities, such as density and velocity, are obtained by sampling molecular information and making an average at the sampling cells. In this work, the aim of using data generated by DSMC is to prove the ability of DHC-GEP to discover the governing equations that could be unknown.

2.2. Dimensional homogeneity constrained gene expression programming

2.2.1. General framework

The DHC-GEP is an enhanced GEP method, which compiles mathematical expressions into chromosomes and iteratively discovers the equation that is suitable for describing the training data by mimicking the natural law of evolution. The flowchart of DHC-GEP is shown in figure 2(a). It starts with creating $N_i$ random individuals of initial population. Each individual has two forms: genotype (chromosome (CS)) and phenotype (expression tree (ET)). Phenotype is the expression product of genotype. Each CS is composed of one or more genes. A specific schematic diagram of a gene is shown in figure 2(b). There are two parts in a gene, i.e. head and tail. The head consists of the symbols from the function set or terminal set, while the tail consists of only symbols from the terminal set. The function set and terminal set are both predefined according to the specific problem. For the problems in this work, the function set includes basic mathematical operators ($+,\ -,\ \times,\ \div$). If necessary for other problems, the function set can also include nonlinear functions such as sin, cos and even user-defined functions. The terminal set includes the symbols of variables and physical constants. Taking the first symbol of the gene as the root node, we can translate the gene into the ET (shown in figure 2b) through level-order traversal according to the argument requirement of each function. Note that the final four symbols are not expressed, and the region before them is called ORF. The length of ORF is variable, which is determined by the argument requirement of the function nodes in the head. On the contrary, the length of gene is fixed, preventing the individuals from bloating. Assuming the length of the head is $n$, then the length of the tail must be $n \times ( {h - 1} ) + 1$, where $h$ is the maximum number of arguments in the functions. Four genes with different ORFs are shown in figure 3. Although the four genes have the same length, the complexities of their expression products are different. Therefore, ORF is the decisive factor that ensures the diversity of expression products and the strong expressivity of DHC-GEP. This is also the reason that the length of the head determines only the upper limit of complexity, not the lower limit. When discovering unknown governing equations, a relatively long head would be a good choice. Because with ORF, a long gene can represent either a simple equation or a complex equation.

Figure 2. Main characteristics of DHC-GEP. (a) Flowchart of DHC-GEP. (b) Schematic diagram of a gene and its corresponding expression tree and mathematical expression. (c) Strategy of dimensional verification.

Figure 3. Four genes with different ORFs.

Then, via dimensional verification (which will be introduced in § 2.2.2), all the individuals are classified into valid individuals and invalid individuals according to whether they satisfy the dimensional homogeneity. The valid individuals would be translated into mathematical expressions (MEs) and evaluated for losses with input data according to the loss function, which is defined using the mean relative error (MRE) as

(2.9)\begin{equation} {{L_{MRE}}} = \frac{1}{N}\sum_{i = 1}^N {\left| {\frac{{{{\hat Y}_i} - {Y_i}}}{{{Y_i}}}} \right|} . \end{equation}

Here, the variable with a superscript $\wedge$ is the predicted variable, and $N$ is the total number of data points. The invalid individuals would be directly assigned a significant loss and eliminated in the subsequent iterations. The best individual (with the lowest loss) is replicated to the next generation directly. Other individuals are selected as superior individuals with probabilities that are inversely proportional to their losses. Based on the superior individuals, the population of the next generation is reproduced via genetic operators. Common genetic operators include replication, mutation, transposition and recombination. The schematic diagrams of them are shown in figure 4; the details are introduced as follows.

1. (i) Replication is responsible for copying the selected superior individuals to the next generation. After replication, the population maintains the same size as the initial population.

2. (ii) Mutation is a change in some symbols of the genes. To maintain the structural organization of the genes, the symbols in the head can be changed into either functions or terminals, while the symbols in the tail can only be changed into terminals. In addition, note that mutations in genes do not necessarily lead to changes in expression products, unless the mutation occurs in ORF.

3. (iii) Transposition is the replacement of a gene segment with another position. Specifically, it includes insertion sequence (IS) transposition, root insertion sequence (RIS) transposition and gene transposition. For IS transposition, a short segment with a function or terminal in the first position is transposed to the head of the genes. The length of the transposed segment and the insertion site are randomly selected. The sequence before the insertion site remains unchanged. The sequence after the insertion site is shifted backwards as a whole, and the last symbols (as many as the length of the transposed segment) of the head are deleted. The RIS transposition is similar to IS transposition, except that the first position of the transposed segment must be a function, and the insertion site must be the root. Gene transposition occurs in the CSs that consist of multiple genes. During gene transposition, one gene switches position with another gene.

4. (iv) Recombination is an exchange of gene segments between two selected parent CSs, including one-point recombination, two-point recombination and gene recombination. For one-point recombination, one bond is randomly selected as the crossover point. Then, both selected CSs are cut at this point, and exchange the segments after the crossover point. For two-point recombination, two crossover points are selected, and the exchanged segments are between these two points. For gene recombination, the exchanged segment is an entire gene.

Figure 4. Schematic diagrams of genetic operators.

Note that the genetic operators, except replication, are not necessarily invoked in every CS, but are invoked with certain probabilities. In this work, the probabilities of the genetic operators being invoked refer to Ferreira (Reference Ferreira2006), listed in table 1, which were concluded from various examples via the trial-and-error approach.

Table 1. Probabilities for genetic operators being invoked.

The above processes are iteratively conducted until a satisfactory individual is obtained. Specifically, in benchmark studies, we manually terminate the evolution when the target equation is discovered. In cases of discovering the unknown constitutive relations, the evolution is terminated automatically if the optimal equation remains unchanged for 1000 consecutive generations.

Other advanced knowledge is provided in the appendices, including the Dc domain (Appendix B) for generating numerical constants and the linear scaling (Appendix C) for discovering numerical coefficients.

2.2.2. Dimensional verification

The dimensional verification is the additional process to implement the constraint of dimensional homogeneity, and is introduced as follows (a simplified version is shown in figure 2c). The general form of a governing equation is assumed to be

(2.10)\begin{equation} \boldsymbol{Y} = f(\boldsymbol{X},\boldsymbol{\vartheta}). \end{equation}

Here, $Y$ is the target variable, and $f(\cdot )$ is the mathematical expression coded by a CS. Here $X$ and $\vartheta$ represent variables and physical constants, respectively. Dimensional homogeneity means that $f(\cdot )$ should have the same dimension as $Y$, and every part in $f(\cdot )$ should satisfy the constraint that the parameters for operators $+$ or $-$ must have the same dimensions.

Generally, dimensional verification includes calculating the dimension of $f(\cdot )$ and comparing it with the dimension of $Y$. For calculating dimensions, there are five principles that should be noted.

1. (i) Computers cannot directly deal with symbolic operations, but only numerical operations.

2. (ii) In SI units, there are seven base dimensions: length ($L$); mass ($M$); time ($T$); electric current ($I$); thermodynamic temperature ($\varTheta$); amount of substance ($N$); luminous intensity ($J$). The dimensions of any other physical quantities can be derived by powers, products or quotients of these base dimensions.

3. (iii) Base dimensions are independent of each other. Any dimension cannot be derived by other base dimensions.

4. (iv) Physical quantities with different dimensions cannot be added or subtracted. Adding or subtracting the physical quantities with the same dimension does not change the dimension.

5. (v) When physical quantities are multiplied or divided, the corresponding dimensions are multiplied or divided equally.

Considering the above principles, we propose to assign number tags to physical quantities and replace the dimensional calculations with numerical calculations. For the seven base dimensions, the tags are 2, 3, 5, 7, 11, 13 and 17, respectively. These numbers are prime numbers, ensuring that the base dimensions are independent of each other. The tags for other physical quantities are derived according to their dimensions. For example, as the dimension of velocity is length divided by time ($L/T$), its tag is defined as $2/5$. It is noteworthy that the tags are always in the form of fractions rather than floats to avoid truncation errors. The tags for numerical constants and dimensionless quantities are set to $1$. In addition, according to the last two principles, we modify the operation rules of mathematical operators. The pseudocodes for addition ($+$) and multiplication ($\times$) are provided in Algorithm 1 and Algorithm 2, respectively. Subtraction ($-$) and division ($\div$) are defined in similar ways. With the number tags and modified mathematical operators, we can compute the dimension of each node in the ET from the bottom up. If the parameters for one function node ($+$ or $-$) do not have the same dimension, the individual is directly identified as invalid one. Otherwise, we can calculate the specific dimension of the individual, i.e. the number tag for the root node. Comparing it with the tag for the target variable, we can conclude whether the individual is dimensionally homogeneous.

Algorithm 1 Modified addition function

Algorithm 2 Modified multiplication function

Note that the datasets and source codes used in this work are available on GitHub at https://github.com/Wenjun-Ma/DHC-GEP.

3.1. Diffusion equation

We employ the DSMC method to simulate a two-dimensional diffusion flow in a square domain, as shown in figure 5(a). The side length of the square is $L = 100\lambda$, and $\lambda$ is the mean free path of argon gas molecules under standard conditions (pressure $p = 1.01 \times {10^5}\ {\rm {Pa}}$ and temperature $T = 273\ {\rm K}$), i.e. $\lambda =6.25\times 10^{-8}\ \mathrm {m}$. As a result, the Knudsen number ($Kn = \lambda / L$) is 0.01. All the four sides are set to periodic boundaries. The simulation includes only argon gas, which is divided into two species, denoted as species-A and species-B, respectively. At the initial instant, the density distributions of these two species are ${\rho _A} = 0.5{\rho _0}( {1 - \cos ( {2{\rm \pi} y/L} )} )$ and ${\rho _B} = 0.5{\rho _0}( {1 + \cos ( {2{\rm \pi} y/L} )} )$, respectively. Here, ${\rho _0}=1.78\ {\rm kg}\,{{\rm m}^{-3}}$ is the density of argon gas under standard conditions, and $y$ is the spatial coordinate in the vertical direction. The initial macroscopic velocity is uniformly zero. During simulation, the macroscopic velocity and the total density of the two species remain uniform and unchanged, while the respective densities of these two species vary with space and time.

Figure 5. (a) Schematic diagram of the computational domain for the diffusion flow. (b) Distributions of density at the instants of $t= 0,\ 200\tau,\ 400\tau,\ 600\tau, 800\tau,\ 1000\tau$.

The computational domain is divided into $256 \times 256$ cells, which means that the cell size (${\rm \Delta} x,{\rm \Delta} y$) is approximately $0.4\lambda$. In each cell, 4000 simulation molecules are randomly distributed at the initial instant, with one simulation molecule representing $4\times 10^{6}$ real molecules. The initial thermal velocity for each simulation molecule is sampled randomly from a Maxwell distribution at 273 K. The computational time step (${\rm \Delta} t$) is $0.1\tau$, where $\tau$ is the mean collision time of argon gas molecules under standard conditions with the definition of $\tau =\lambda /\bar {c}$. Here, $\bar {c}$ is the molecular mean velocity, i.e. $\bar {c} =\sqrt {8k_{B}T/{\rm \pi} m}$, where $k_B$ and $m$ are the Boltzmann constant and molecular mass, respectively.

Note that the density varies only along the vertical direction. We average the data along the horizontal direction to reduce the statistical errors. During simulation, we sample the macroscopic densities of species-A at the instants of $t = 0, 10\tau, 20\tau,\ldots, 1000\tau$. The temporal evolution of density distributions is shown in figure 5(b). The final dataset is composed of 25 856 ($256 \times 101$) data points that are distributed in 256 cells along the vertical direction at 101 instants.

The theoretical governing equation of this flow is the diffusion equation, i.e.

(3.1)\begin{equation} \frac{{\partial \rho }}{{\partial t}} = D\frac{{{\partial ^2}\rho }}{{\partial {y^2}}}. \end{equation}

Here, $D$ is the diffusion coefficient. According to the Chapman–Enskog theory (Chapman & Cowling Reference Chapman and Cowling1990), the diffusion coefficient for an HS gas at the first-order approximation is

(3.2)\begin{equation} D=\frac{3{\rm \pi} }{16} \frac{\lambda ^{2} }{\tau }. \end{equation}

Therefore, for the argon gas in our simulation, $D$ equals $1.40 \times {10^{-5}}\ {{\rm m}^2}\ {{\rm s}^{ - 1}}$.

Before using the dataset generated by DSMC, we test DHC-GEP and Original-GEP on a clean dataset. The only difference between the DSMC dataset and the clean dataset lies in the target variable (${{\partial \rho }}/{{\partial t}}$). In the clean dataset, ${{\partial \rho }}/{{\partial t}}$ is directly computed based on $D$ and ${{{\partial ^2}\rho }}/{{\partial {y^2}}}$, according to (3.1). We define the function set as $\{ { + , - , \times, \div } \}$, and the terminal set as $\{ {\mu,D, {{\partial \rho }}/{{\partial y}}, {{{\partial ^2}\rho }}/{{\partial {y^2}}}, {{{\partial ^3}\rho }}/{{\partial {y^3}}}} \}$. Here, $\mu$ is the viscosity coefficient, and is introduced as a distraction term. According to the Chapman–Enskog theory (Chapman & Cowling Reference Chapman and Cowling1990), $\mu$ equals $({5{\rm \pi} }/{32}) ({\lambda ^{2} }/{\tau })\rho = 2.08 \times {10^{ - 5}}\ {\rm kg}\ {{\rm m}^{-1}}\ {{\rm s}^{ - 1}}$. All the hyperparameters are set to be the same for both GEP methods across all cases (see Appendix E). The derived equations are shown in table 2. It can be observed that the derived equation of DHC-GEP is consistent with the theoretical equation, while the derived equation of Original-GEP is partially correct. Specifically, Original-GEP identifies the correct function form and coefficient ($1.40 \times {10^{ - 5}}$). In the view of data fitting, this is a correct result. However, Original-GEP fails to recognize that the coefficient ($1.40 \times {10^{ - 5}}$) preceding the diffusion term represents the diffusion coefficient. Consequently, while the resulting equation may be numerically accurate, it lacks complete physical significance. For another diffusion flow with a different diffusion coefficient, the derived equation of Original-GEP is no longer valid, while the derived equation of DHC-GEP is generally applicable.

Table 2. Derived equations of DHC-GEP.

Subsequently, we test DHC-GEP and Original-GEP on the DSMC dataset. All settings are consistent with those in the test on the clean dataset. The derived equations are shown in table 2. The coefficient of the equation discovered by DHC-GEP has a minor deviation from that of the theoretical equation, and the loss is not zero. This is because DSMC is a stochastic molecule-based method, the data of which are inherently noisy. It is impossible for DSMC to simulate a flow with the diffusion coefficient being exactly $1.40 \times {10^{ - 5}}\,{{\rm m}^2}\,{{\rm s}^{ - 1}}$. Small fluctuations around the theoretical value are acceptable. In addition, calculating derivatives also introduces errors. Therefore, we believe that DHC-GEP has discovered the correct equation, while the derived equation of Original-GEP is clearly wrong for not satisfying the fundamental dimensional homogeneity. Besides, we generate another dataset of diffusion flow with a diffusion coefficient being $3.0 \times {10^{ - 5}}\,{{\rm m}^2}\,{{\rm s}^{ - 1}}$ and calculate the losses of the two derived equations based on the new dataset. The losses for the derived equations of Original-GEP and DHC-GEP are 0.55 and 0.06, respectively. Therefore, the equation obtained by Original-GEP is clearly overfitting, while the equation obtained by DHC-GEP is generally applicable in different diffusion flows.

To test the sensitivity of DHC-GEP to hyperparameters, we conduct a parametric study on the length of the head, the number of genes in a CS and the number of individuals in a population. Specifically, the parameter spaces for these three hyperparameters are set to $\{ 10, 11, 12, 13, 14, 15 \}$, $\{1, 2 \}$ and $\{800, 1000, 1200, 1400, 1600, 1660 \}$, respectively, resulting in a total of 72 ($6\times 2\times 6$) distinct parameter combinations. Other hyperparameters are kept unchanged, including the probabilities of invoking the genetic operators and the maximum number of evolution generations. In all 72 sets of parameter combinations, DHC-GEP consistently discovers the correct equation. In contrast, Original-GEP discovers 71 completely different results, all of which are overfitting and fail to satisfy the dimensional homogeneity constraint despite their lower losses than the correct equations.

Considering that DSMC is a simulation method at mesoscale, which coarse grains the molecular description to the hydrodynamic regime (Hadjiconstantinou Reference Hadjiconstantinou2000), the influence of the coarse graining on the derived equations of DHC-GEP is discussed. For the same computational domain, we employ seven sets of sampling cells with different resolutions ($256\times 256$, $128\times 128$, $64\times 64$, $32\times 32$, $16\times 16$, $8\times 8$, $4\times 4$) to acquire macroscopic physical quantities. The derived equations based on these datasets are summarized in table 3. It is encouraging to note that DHC-GEP is still capable of discovering the correct equation when the resolution is $16\times 16$ and the size of the sampling cells is up to $6.3\lambda$. From the results, it can be concluded that DHC-GEP is robust to coarse graining, as long as the sampled data can accurately capture the gradients of macroscopic quantities.

Table 3. Derived equations of DHC-GEP based on the datasets with different resolutions.

In DHC-GEP, the overfitting results are automatically filtered out due to not satisfying the dimensional homogeneity. On the contrary, Original-GEP always favours the equations with smaller loss, so it tends to converge to overfitting results if data are noisy. Our numerical experiences show that Original-GEP may discover the correct equations only when tuning the hyperparameters repeatedly and terminating the evolution at a proper generation (when overfitting equations have not appeared).

We compare the computational cost needed per 1000 generations of evolution. Based on the 3.5 GHz Intel Xeon E5-1620 processor, the average CPU runtime of Original-GEP (443 s) is almost twice that of DHC-GEP (216 s). The main reason is that DHC-GEP can identify some individuals as invalid through dimensional verification and then skip the process of evaluating losses for these individuals. It can still save computational time despite the extra expense of dimensional verification. Additionally, as the complexity of the problem increases, the number of invalid individuals also increases, so the advantage of the computational efficiency of DHC-GEP becomes more significant.

Furthermore, we conduct a much more challenging test in which we do not include the transport coefficients in the terminal set. Specifically, the terminal set is $\{ {\tau,\lambda, {{\partial \rho }}/{{\partial y}}, {{{\partial ^2}\rho }}/{{\partial {y^2}}}, {{{\partial ^3}\rho }}/{{\partial {y^3}}}} \}$, which involves only the fundamental physical property parameters of gas (molecular mean collision time ($\tau$) and molecular mean free path ($\lambda$)), but excludes the diffusion coefficient ($D$). Other hyperparameters are the same as those in the original test. The derived equation is

(3.3)\begin{equation} \frac{{\partial \rho }}{{\partial t}} = 0.579\frac{\lambda ^{2} }{\tau } \frac{{{\partial ^2}\rho }}{{\partial {y^2}}}. \end{equation}

According to (3.2), the coefficient ($0.579({\lambda ^{2} }/{\tau })$) is essentially the diffusion coefficient. This result demonstrates that DHC-GEP has the capability to derive accurate equations with minor prior knowledge of the essential transport coefficients, and can deduce the essential transport coefficients from fundamental physical property parameters using dimensional homogeneity as a constraint.

3.2. Vorticity transport equation

We employ DSMC to simulate the temporal evolution of a Taylor–Green vortex, as shown in figure 6(a). The simulation model is approximately the same as that of diffusion flow, including simulation conditions, geometry and boundary conditions. The major difference is that the initial distribution of macroscopic velocity is a Taylor–Green vortex,

(3.4)\begin{equation} \begin{cases} {\begin{array}{*{20}{@{}l}} {u = {v_0}\cos (x)\sin (y)}, \\ {v ={-} {v_0}\sin (x)\cos (y)}, \end{array}}\end{cases} \end{equation}

where $u$ and $v$ are the velocities in the horizontal ($x$) and vertical ($y$) directions, respectively, and ${v_0}$ is the initial amplitude of velocity and equals $30\,{\rm m}\,{{\rm s}^{ - 1}}$ in our simulation.

Figure 6. Temporal evolution of vorticity. (a) Taylor–Green vortex at the instants of $t = 0,100\tau, 200\tau$. Here, $\tau$ and $\lambda$ are the mean collision time and mean free path of argon gas molecules at standard condition, respectively. (b) Viscous flow around a cylinder at the instants of $t = 0\,{\rm s},\ 3\,{\rm s},\ 6\,{\rm s}$.

The computational cell and time step are set the same as those in the diffusion flow. Additionally, we introduce sampling cells. In this case, each sampling cell consists of 16 computational cells, i.e. the whole computational domain has $64 \times 64$ sampling cells. The macroscopic physical quantities obtained with sampling cells are cleaner than those obtained with computational cells, due to more molecules being sampled in each cell. In addition, to further reduce the statistical errors, 10 independent simulations are conducted with different random number sequences to make an ensemble average. During simulation, we sample the macroscopic velocities ($u,v$) and vorticities (${\omega _z}$) at the instants of $t = 0, 10\tau, 20\tau,\ldots, 310\tau$. The final dataset is composed of 131 072 ($64 \times 64 \times 32$) data points distributed in $64 \times 64$ sampling cells at 32 instants.

The theoretical governing equation is the vorticity transport equation,

(3.5)\begin{equation} \frac{{\partial {\omega _z}}}{{\partial t}} = \upsilon \left({\frac{{{\partial ^2}{\omega _z}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\omega _z}}}{{\partial {y^2}}}} \right) - \left( {u\frac{{\partial {\omega _z}}}{{\partial x}} + v \frac{{\partial {\omega _z}}}{{\partial y}}} \right). \end{equation}

Here, ${\omega _z}$ is the vorticity in the $z$ direction. $\upsilon$ is the kinematic viscosity and approximately equals $1.17 \times {10^{ - 5}}\,{{\rm m}^2}\,{{\rm s}^{ - 1}}$. Compared with the diffusion equation (3.1), this equation is more complex, involving multiple variables and nonlinear terms.

In this case, we define the target variable as ${{\partial {\omega _z}}}/{{\partial t}}$, the function set as $\{ { + , - , \times, \div } \}$, and the terminal set as $ \{ \upsilon,u,v,\omega _z, {{\partial u}}/{{\partial x}}, {{{\partial ^2}u}}/{{\partial {x^2}}}, {{\partial u}}/{{\partial y}}, {{{\partial ^2}u}}/{{\partial {y^2}}}, {{\partial v}}/{{\partial x}}, {{{\partial ^2}v}}/{{\partial {x^2}}}, {{\partial v}}/{{\partial y}}, {{{\partial ^2}v}}/{{\partial {y^2}}}, {{\partial \omega _z }}/{{\partial x}}, {{{\partial ^2}\omega _z }}/{{\partial {x^2}}}, {{\partial \omega _z }}/{{\partial y}}, {{{\partial ^2}\omega _z }}/{{\partial {y^2}}}\}$.

We also test the performances of DHC-GEP on the clean dataset and the DSMC dataset sequentially. The results are provided in table 2. An interesting observation is that the derived equations of DHC-GEP miss two convective terms $( {u ({{\partial {\omega _z}}}/{{\partial x}}), v({{\partial {\omega _z}}}/{{\partial y}})})$. This is caused by the unique feature of the Taylor–Green vortex. The analytical solution of the Taylor–Green vortex is

(3.6)\begin{equation} \begin{cases} \begin{array}{*{35}{@{}l}} u=v_0\, \cos(x)\sin(y) {\exp ({ - 2\upsilon t})}, \\ v={-}v_0\, \sin(x)\cos(y) {\exp ({ - 2\upsilon t})}, \\ {{\omega _z}} ={-}2v_0\, \cos(x)\cos(y) {\exp ({ - 2\upsilon t})}. \end{array} \end{cases} \end{equation}

Substituting (3.6) into the convective terms $( {u ({{\partial {\omega _z}}}/{{\partial x}}), v({{\partial {\omega _z}}}/{{\partial y}})})$, it is clear that the sum of the convective terms is automatically zero. Therefore, the equations without convective terms are also correct for the Taylor–Green vortex. It can be considered as a specialized variant of the complete vorticity transport equation in the Taylor–Green vortex.

To discover the complete vorticity transport equation and validate the performance of DHC-GEP on discovering nonlinear terms, we further consider a viscous flow around a cylinder at Reynolds number being 100, as shown in figure 6(b). The dataset is the open access dataset provided by Rudy et al. (Reference Rudy, Brunton, Proctor and Kutz2017), and the theoretical governing equation is

(3.7)\begin{equation} \frac{{\partial {\omega _z}}}{{\partial t}} = 0.01 \,{\cdot}\, \left( {\frac{{{\partial ^2}{\omega _z}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\omega _z}}}{{\partial {y^2}}}} \right) - \left( {u\frac{{\partial {\omega _z}}}{{\partial x}} + v \frac{{\partial {\omega _z}}}{{\partial y}}} \right). \end{equation}

For the sake of simplicity, we regard the data as dimensional data with kinematic viscosity $\upsilon = 0.01\,{{\rm m}^2}\,{{\rm s}^{ - 1}}$. We randomly sampled 200 spatiotemporal data points from the Taylor–Green vortex dataset and another 200 spatiotemporal data points from the flow around cylinder dataset, forming a hybrid dataset consisting of 400 spatiotemporal data points. The hyperparameters of DHC-GEP are consistent with those in the case of the Taylor–Green vortex. The results are provided in table 2. The DHC-GEP discovers the correct and complete vorticity transport equation. This demonstrates that in some special cases, where the overfitting equations may satisfy the dimensional homogeneity and fit well with the training data, employing multiple datasets can help mitigate overfitting. Overfitting equations are generally dataset specific. While an overfitting equation may adequately describe a particular dataset, it cannot be expected to generalize well to multiple datasets. Besides, note that the training dataset consisting of 400 spatiotemporal data points is relatively small. Generally, if the size of dataset is small, the information carried by data is sparse, leading to the derived equations overfitting to specific phenomena. However, in DHC-GEP, the overfitting equations are automatically filtered out by the constraint of dimensional homogeneity.

4.1. One-dimensional shock wave

In the research community of non-equilibrium flows, the one-dimensional shock wave is a benchmark to validate solvers and formulations in numerical computations. We employ DSMC to simulate the one-dimensional shock wave of argon gas, the general structure of which is shown in figure 7.

Figure 7. General structure of one-dimensional shock. Each variable is normalized with its maximum value in the computational domain.

The one-dimensional shock essentially connects two equilibrium states. In this work, for upstream (i.e. the free stream), the temperature (${T_1}$) is $300\,{\rm K}$, the density (${\rho _1}$) is $1.067 \times {10^{ - 4}}\,{\rm kg}\,{{\rm m}^{ - 3}}$ and the mean free path (${\lambda _1}$) is $1.114 \times {10^{ - 3}}\,{\rm m}$. For downstream, the state variables are computed with Rankine–Hugoniot relations

(4.2) \begin{equation} \begin{cases} \begin{array}{*{35}{l}} {T_2} = \dfrac{{({5Ma_\infty ^2 - 1})({Ma_\infty ^2 + 3})}}{{16Ma_\infty ^2}}{T_1}, \\ {\rho _2} = \dfrac{{4Ma_\infty ^2}}{{Ma_\infty ^2 + 3}}{\rho _1}, \\ {v_2} = \dfrac{{Ma_\infty ^2 + 3}}{{4Ma_\infty ^2}}{v_1}. \end{array} \end{cases} \end{equation}

Here, $M{a_\infty }$ and ${v_1}$ are the Mach number and velocity of the free stream, respectively. Note that for this one-dimensional problem, the velocity direction of each position in the flow field is along the $x$ direction, and all the flow variables vary along only the $x$ direction. Therefore, we set one cell along the $y$ and $z$ directions, but 120 cells along the $x$ direction. A schematic diagram is shown in figure 8. The lengths of the computational domain and computational cell are $30{\lambda _1}$ and $0.25{\lambda _1}$, respectively. The time step is $0.1\tau$.

Figure 8. Schematic diagram of the computational domain for one-dimensional shock.

During simulation, we sample the macroscopic velocities ($u$), densities ($\rho$), temperatures ($T$) and heat flux in the $x$ direction (${q_x}$) at the instants of $t = 10\,000\tau, 10\,010\tau,\ 10\,020\tau,\ldots,50\,000\tau$. Then, we obtain the final data by averaging all the sample data.

In this case, the target variable is defined as the heat flux in the $x$ direction (${q_x}$). The function set is defined as $\{ { + , - , \times, \div } \}$ as usual, while the terminal set is meticulously constructed to embed physics knowledge as follows.

1. (i) Constitutive relations describe the local transport mechanisms of momentum and energy, and hence we regard local non-equilibrium parameters as key factors. Specifically, the GLL Knudsen number defined as $K{n_Q} = ({{{\lambda _l}}}/{Q})({\partial Q}/{\partial x})$ is selected. The local non-equilibrium characteristics intensify as the increase of the absolute value of $K{n_Q}$. Here, ${\lambda _l}$ is the local mean free path, and $Q$ represents state variables, including temperature ($\theta$) and density ($\rho$).

2. (ii) The transports of momentum and energy are driven by the gradients of state variables. As a result, the gradient terms are selected, including ${{\partial u}}/{{\partial x}}, {{\partial \rho }}/{{\partial x}}, {{\partial \theta }}/{{\partial x}}, {{{\partial ^2}u}}/{{\partial {x^2}}}, {{{\partial ^2}\rho }}/{{\partial {x^2}}}, {{{\partial ^2}\theta }}/{{\partial {x^2}}}, {{{\partial ^3}u}}/{{\partial {x^3}}}, {{{\partial ^3}\rho }}{{\partial {x^3}}}$ and ${{{\partial ^3}\theta }}/{{\partial {x^3}}}$.

3. (iii) The state variables themselves are important factors. Besides, it is noteworthy that the constitutive relations should satisfy Galilean invariance (Han et al. Reference Han, Ma, Ma and E2019; Li, Dong & Wang Reference Li, Dong and Wang2021), which means that the constitutive relations cannot contain velocity ($u$) explicitly outside the partial differential operators. The proof is provided in Appendix F. Therefore, the state variables ($\theta$ and $\rho$) excluding velocity are selected. This is also the reason that we do not select the GLL Knudsen number of velocity.

4. (iv) The parameters representing the physical properties of gas are selected, including local viscosity ($\mu$), local heat conductivity ($\kappa$), heat capacity ratio ($\gamma = 5/3$ for argon gas) and viscosity exponent ($\omega = 1$ for Maxwell molecules).

Finally, the terminal set is $\{K{n_\theta },K{n_\rho }, {{\partial u}}/{{\partial x}}, {{\partial \rho }}/{{\partial x}}, {{\partial \theta }}/{{\partial x}}, {{{\partial ^2}u}}/{{\partial {x^2}}}, {{{\partial ^2}\rho }}/{{\partial {x^2}}}, {{{\partial ^2}\theta }}/{{\partial {x^2}}}, {{{\partial ^3}u}}/{{\partial {x^3}}}, {{{\partial ^3}\rho }}/{{\partial {x^3}}}, {{{\partial ^3}\theta }}/{{\partial {x^3}}},\rho,\theta,\mu,\kappa, \gamma,\omega \}$.

Moreover, constitutive relations ought to satisfy the second law of thermodynamics. Theoretically, whether a constitutive relation satisfies the second law of thermodynamics can be determined through the Clausius–Duhem inequality (Comeaux Reference Comeaux1995)

(4.3)\begin{equation} \rho \frac{{{\rm D}{s_{eq}}}}{{{\rm D} t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \frac{\boldsymbol{q}}{\theta } ={-} \frac{1}{\theta }\boldsymbol{\tau} :\boldsymbol{\nabla} \boldsymbol{v} - \frac{1}{{{\theta ^2}}}\boldsymbol{q} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta \ge 0. \end{equation}

The two terms on the left-hand side are the local increase rate of entropy and the reversible outflow of entropy, respectively. The sum of the two terms on the right-hand side is the total entropy production. The second law of thermodynamics demands that the total entropy production must be non-negative. Note that the total entropy production is the sum of subentropy productions of all high-order macroscopic variables (the viscous stress ${\tau _{ij}}$ and heat flux ${q_i}$). If each subentropy production is non-negative, the total entropy production is naturally non-negative. This work deals with the discovery of the constitutive equation of a single high-order macroscopic variable, rather than all of them. Hence, the entropy production discussed in this work refers to the contribution of a single high-order macroscopic variable's constitutive relation. Specifically, in the one-dimensional shock wave case, the target variable is the heat flux in the $x$ direction (${q_x}$), and thus the entropy production considered here is

(4.4)\begin{equation} S_{p,q_{x} } ={-} \frac{1}{{{\theta ^2}}}q_{x}\frac{\partial \theta }{\partial x}. \end{equation}

To ensure that the resulting constitutive equation satisfies the second law of thermodynamics, we introduce this law as a constraint into the loss function (2.9). Specifically, we incorporate the constraint in a soft form by adding a loss term (${L}_{S_p}$) that represents the ratio of the number of the data points with negative entropy production to the total number of data points. The modified loss function is

(4.5)\begin{equation} \mathrm{Loss}={L}_{{MRE}}+\alpha {L}_{S_p}=\frac{1}{N} \sum_{i=1}^{N}\left | \frac{\hat{Y}_i-Y_i}{Y_i}\right |+\alpha \frac{\langle S_p\rangle _{-} }{N}. \end{equation}

Here, the variable with a superscript $\wedge$ is the predicted variable, and $N$ is the total number of data points. Here $\langle S_p \rangle _{-}$ is the operator that counts the number of data points with negative entropy production, and $\alpha$ is the weighting factor that controls the importance of the constraint.

We employ the DHC-GEP, combined with the constraint of the second law of thermodynamics, to discover the underlying constitutive relation based on the data of two cases with the free stream Mach number $M{a_\infty } = 3.0$ and $4.0$. We sequentially set $\alpha$ to 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0. Each $\alpha$ corresponds to an independent training run and a resulting equation with relatively low ${L}_{{MRE}}$. Among them, the equation with the smallest ${L}_{S_p}$ is ultimately selected as the optimal equation. For this case, the optimal equation is

(4.6)\begin{equation} {q_x} ={-} 0.468\left( {6K{n_\rho }+\frac{{K{n_\theta }}}{{K{n_\rho }}} + \frac{{K{n_\rho }}}{{K{n_\theta }}}} \right)\kappa \frac{{\partial \theta }}{{\partial x}}, \end{equation}

which is derived with $\alpha =0.3$. The entropy production is

(4.7)\begin{equation} S_{p,q_{x} } = \frac{0.468\kappa}{\theta ^{2} } \left( {6K{n_\rho } +\frac{{K{n_\theta }}}{{K{n_\rho }}} + \frac{{K{n_\rho }}}{{K{n_\theta }}}} \right) \left(\frac{{\partial \theta }}{{\partial x}} \right ) ^{2} . \end{equation}

According to figure 7, for the one-dimensional shock wave case, the gradients of density and temperature are non-negative. Additionally, according to the definition of the local Knudsen number ($K{n_Q} = ({{{\lambda _l}}}/{Q})({\partial Q}/{\partial x})$), both $K{n_\rho }$ and $K{n_\theta }$ are thus non-negative. Therefore, from a mathematical perspective, the entropy production (4.7) must be non-negative. As a comparison, whether the Burnett equation, augmented-Burnett equation and super-Burnett equation satisfy the second law of thermodynamics is obscure (Comeaux Reference Comeaux1995).

We compare the results predicted by (4.6) with those predicted by the NSF equation, the Burnett equation, the augmented-Burnett equation and the super-Burnett equation in figure 9. For a quantitative comparison, we define the relative error as

(4.8)\begin{equation} \mathrm{error} = \frac{\displaystyle{\sqrt{\sum_{i = 1}^N {{{({{{\hat q}_{x,i}} - {q_{x,i}}})}^2}} } }}{\displaystyle{\sqrt {\sum_{i = 1}^N {q_{x,i}^2} } }}. \end{equation}

The relative errors for each constitutive relation are summarized in table 4. It can be found that the derived constitutive relation of DHC-GEP exhibits higher accuracy than other equations over a wide range of $M{a_\infty }$, and is applicable to cases beyond the parameter space of the training data, i.e. $2.0 \le M{a_\infty } < 3.0$ and $4.0 < M{a_\infty } \le 5.5$.

Figure 9. Comparisons of the NSF equation, Burnett equation, augmented-Burnett equation (Au-Burnett), super-Burnett equation (Su-Burnett), derived equation of DHC-GEP and DSMC (exact solution) at $M{a_\infty } = 2.5,3.0,3.5,4.0,4.5$ and $5.0$ for the one-dimensional shock wave case.

Table 4. Relative errors of the derived equation of DHC-GEP, NSF equation, Burnett equation, Au-Burnett equation and Su-Burnett equation for the one-dimensional shock wave case.

Furthermore, to investigate the sensitivity to hyperparameters in non-equilibrium cases, we conduct a parametric study. Specifically, we focus on the sensitivity to the length of the head, the number of genes in a CS and the number of individuals in a population. The parameter spaces for these three hyperparameters are set to $\{ 10, 11, 12, 13, 14, 15 \}$, $\{1, 2 \}$ and $\{800, 1000, 1200, 1400, 1600, 1660 \}$, respectively, resulting in a total of 72 ($6\times 2\times 6$) distinct hyperparameter combinations. Other hyperparameters, including the probability of invoking the genetic operators, are kept unchanged. Among the total of 72 distinct hyperparameter combinations, DHC-GEP successfully derives the same equation as (4.6) based on 38 of these hyperparameter combinations. In contrast, Original-GEP discovers 70 completely different results, all of which fail to satisfy even the fundamental dimensional homogeneity constraint, let alone other physical constraints. It can be concluded that in non-equilibrium cases, while DHC-GEP is not entirely insensitive to model hyperparameters, it does exhibit significantly reduced sensitivity compared with Original-GEP.

4.2. Rarefied Poiseuille flow

Poiseuille flow is a flow confined between two infinite, parallel and relative static plates. The gas is driven by an external force in the $x$ direction. The external force is uniformly distributed along the $y$ direction. A schematic diagram of Poiseuille flow is shown in figure 10. The global Knudsen number ($Kn$) is defined as

(4.9)\begin{equation} Kn = \frac{\lambda }{H}. \end{equation}

Here, $H$ is the distance between two plates, and $\lambda$ is the mean free path of argon gas molecules at standard condition. For different $Kn$, $\lambda$ remains the same, while $H$ changes according to (4.9). The uniform external force is implemented through a uniform acceleration in molecular simulations. The accelerations for different $Kn$ are listed in table 5.

Figure 10. (a) Schematic diagram of Poiseuille flow. (b) Profiles of velocity, density and temperature along the direction normal to the plate. Each variable ($U$) is normalized with ${U_{Normal}} = ( {U - {U_{\min }}} )/( {{U_{\max }} - {U_{\min }}} )$.

Table 5. Accelerations for different $Kn$.

All the flow variables vary along only the $y$ direction. Therefore, we set one cell along the $x$ and $z$ directions, but 3000 cells along the $y$ direction. The time step is $0.2\tau$. During simulation, we sample the macroscopic velocities ($u$), densities ($\rho$), temperatures ($T$) and viscous shear stress (${\tau _{xy}}$) at the instants of $t = 100\,000\tau, 102\,000\tau, 104\,000\tau,\ldots, 1\,000\,000\tau$. Then, we obtain the final data by averaging all the sample data. In addition, to avoid the influence of boundary conditions, we subsample the data points in the region of $0.1H \le y \le 0.9H$ to form the training dataset.

In this case, the target variable is defined as the viscous shear stress (${\tau _{xy}}$), and the function set is defined as $\{ { + , - , \times, \div } \}$. The terminal set is almost the same as that in the case of the one-dimensional shock wave, i.e. $\{ K{n_\theta },K{n_\rho }, {{\partial u}}/{{\partial y}}, {{\partial \rho }}/{{\partial y}}, {{\partial \theta }}/{{\partial y}},\rho,\theta,\mu, \gamma,\omega \}$, except that the high-order gradients are removed. The motivation for removing high-order gradients is that the derived constitutive relation should be applicable in CFD. If containing high-order gradients, the constitutive relation would be unstable and require additional boundary conditions (Bobylev Reference Bobylev1982; Zhong et al. Reference Zhong, MacCormack and Chapman1993; Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003; Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2004; Singh, Jadhav & Agrawal Reference Singh, Jadhav and Agrawal2017), which are the common problems with the kind of Burnett equations. Besides, note that the GLL Knudsen number in this case is defined as $K{n_Q} = ({{{\lambda _l}}}/{Q})({\partial Q}/{\partial y})$.

Similar to the one-dimensional shock wave case, we employ the DHC-GEP, combined with the constraint of the second law of thermodynamics, to discover the underlying constitutive relation based on the data for the cases of $Kn = 0.01,\;0.3$ and $0.4$. The optimal equation is

(4.10)\begin{equation} {\tau _{xy}} ={-} 0.247\left( {\frac{5}{{Kn_\rho ^2 + \gamma }} - \frac{{K{n_\theta }}}{{K{n_\rho }}}} \right)\mu \frac{{\partial u}}{{\partial y}}, \end{equation}

which is derived with $\alpha =0.7$. The ${L}_{S_p}$ for (4.10) is 0, which quantitatively indicates that the entropy productions at all data points are non-negative.

We compare the results predicted by (4.10) with those predicted by the NSF equation, the Burnett equation, the augmented-Burnett equation and the super-Burnett equation in figure 11. The relative errors for each constitutive relation are summarized in table 6. It can be found that the derived constitutive relation of DHC-GEP is much more accurate than other equations in a wide range of $Kn$, and is applicable to cases beyond the parameter space of the training data, i.e. $0.4 < Kn \le 0.6$.

Figure 11. Comparisons of the NSF equation, Burnett equation, augmented-Burnett equation (Au-Burnett), super-Burnett equation (Su-Burnett), derived equation of DHC-GEP and DSMC (exact solution) at $Kn = 0.2, 0.3, 0.35, 0.4, 0.5$ and $0.6$ for the rarefied Poiseuille flow case. Note that the NSF equation is equal to the Burnett equation in this case.

Table 6. Relative errors of the derived equation of DHC-GEP, NSF equation, Burnett equation, Au-Burnett equation and Su-Burnett equation for the rarefied Poiseuille flow case. Note that the NSF equation is equal to the Burnett equation in this case.

Moreover, note that for continuum flows, the GLL Knudsen numbers ($K{n_\rho }$ and $K{n_\theta }$) approach zero, and (4.10) reduces to

(4.11)\begin{equation} {\tau _{xy}} ={-} 0.247\left( {\frac{5}{\gamma } - \frac{{K{n_\theta }}}{{K{n_\rho }}}} \right)\mu \frac{{\partial u}}{{\partial y}} ={-} 0.247\left( {3 - \frac{{K{n_\theta }}}{{K{n_\rho }}}} \right)\mu \frac{{\partial u}}{{\partial y}}. \end{equation}

Here, we used the value of $\gamma =5/3$ for monatomic gas in the above equation. For the Poiseuille flow in the continuum regime, the pressure along the $y$ direction is observed to remain constant. Taking the flow of $Kn = 0.01$ as an example, the pressure oscillates approximately $0.5\,\%$ of its absolute value. Considering the ideal gas equation of state $\rho = {p}/{\theta }$, we can obtain

(4.12)\begin{equation} \partial \rho ={-} \frac{p}{{{\theta ^2}}}\partial \theta ={-} \frac{\rho }{\theta }\partial \theta . \end{equation}

Then, combining (4.12) with $K{n_\rho } = ({{{\lambda _l}}}/{\rho })({{\partial \rho }}/{{\partial y}})$ and $K{n_\theta } = ({{{\lambda _l}}}/{\theta })({{\partial \theta }}/{{\partial y}})$, we can conclude that ${{K{n_\theta }}}/{{K{n_\rho }}} = {{\partial \theta }}/{{\partial \rho }} {\cdot } {\rho }/{\theta } = - 1$. Substituting this relation into (4.11) yields

(4.13)\begin{equation} {\tau _{xy}} ={-} 0.247({3 + 1})\mu \frac{{\partial u}}{{\partial y}} \approx{-} \mu \frac{{\partial u}}{{\partial y}}. \end{equation}

Therefore, it can be concluded that in the continuum regime, the derived constitutive relation of DHC-GEP can be reduced to NSF equations. This makes it applicable to real complex flow problems, where flows tend to be multiscale, i.e. consisting of both continuum flows and non-equilibrium flows.

Finally, we emphasize that as containing only the first-order gradient of velocity, the derived constitutive relation is stable and requires the same boundary conditions as NSF equations. It is convenient to embed it into the well-developed CFD frameworks with minor modification.