Discrete Boltzmann Modeling of Compressible Flows

Abstract

Mathematically, the typical difference of discrete Boltzmann model (DBM) from the traditional hydrodynamic one is that the Navier-Stokes (NS) equations are replaced by a discrete Boltzmann equation. But physically, this replacement has a significant gain: a DBM is roughly equivalent to a hydrodynamic model supplemented by a coarse-grained model of the thermodynamic non-equilibrium (TNE) effects, where the hydrodynamic model can be and can also be beyond the NS. Via the DBM, it is convenient to perform simulations on systems with flexible Knudsen number. The observations on TNE are being obtaining more applications with time.

1. Introduction

Generally, compressible flow is frequently referred also to as gas dynamics which is the branch of fluid mechanics that deals with flows having significant changes in fluid density. That is because gases, mostly, display such behavior. In fact, all materials are compressible. Besides liquids, gases and plasmas, to some extent, the plastic solids under strong shock can also be modeled as compressible flows. Because, in the last case, that the strength of material is negligibly smaller than that of the shock. Flows with a Mach number less than 0.3 are usually treated as being incompressible for that the variation of density due to velocity is less than 5% in that case. The study on compressible flow is relevant to various fields, such as high-speed aircraft, rocket motors, jet engines, high-speed entry into a planetary atmosphere, gas pipelines, etc. where shock wave and/or detonation play a significant role. In this chapter, the following several kinds of flows including high Mach number flows with combustion, multiphase flows with phase separation and complex flows with hydrodynamic instability¹ are taken as examples, and all flows are treated uniformly as being compressible. It is hopeful that the methods and ideas developed in this chapter may be adapted, more or less, to other kinds of complex flows.

Some common and typical features of these flows are as below: each of them possesses multi-scale structures and/or kinetic modes. There exist plenty of interfaces inside the system. The interfaces include material interfaces and mechanical interfaces. Each of them experiences complicated competitions between behaviors in various spatial-temporal scales. The forcing and responsive processes inside the system are very complicated. Such a system generally shows pronounced non-equilibrium behaviors.

1.1. Traditional models

The traditional macroscopic models of compressible flows are generally based on Navier-Stokes (NS) equations, even Euler equations. The model of Euler equations assumes that the system is always at its local thermodynamic equilibrium (LTE). The NS model considers the thermodynamic non-equilibrium (TNE) via the viscous stress and heat flux. The viscous stress and heat flux compose a set of convenient and effective description of the TNE. But the description is also quite dense and coarse-grained. Many specific TNE behaviors are invisible under NS description, even though they are helpful for understanding the specific TNE status. Besides, since it includes only the first-order term of Knudsen number in the Chapman-Enskog expansion [1], it is reasonable only when the Knudsen number is very small. It cannot be used to access deeper TNE behaviors. To access the complicated non-equilibrium behaviors, one possible solution is to use the molecular dynamics (MD) [2] or direct simulation Monte Carlo (DSMC) [3]. The MD simulation can help to understand some fundamental mechanisms from the atomic level. But the spatial and temporal scales it can access are too small to be comparable with experiments. The DSMC simulation has a similar constraint of computational cost.

The NS model is not enough to capture so complex non-equilibrium behaviors while the MD and DSMC cannot access spatial-temporal scales that are large enough. Under such conditions, a kinetic approach based on the non-equilibrium statistical mechanics (NESM) [1, 4, 5] is preferable.

1.2. Non-equilibrium statistical mechanics

Statistical mechanics is a branch of theoretical physics. It was developed to study the average behavior of a mechanical system with uncertain state by using probability theory. Microscopic mechanical laws do not contain concepts such as temperature, heat, or entropy; however, statistical mechanics shows how these concepts arise from the natural uncertainty about the state of a system when that system is prepared in practice. Statistical mechanics provides exact methods to connect thermodynamic quantities to microscopic behavior.

The NESM is based on mechanics and some necessary assumptions. The concept of macroscopic observation and assumption of coarse-grained density are cornerstones. The Liouville equation [4, 5] is the most fundamental governing equation when without quantum fluctuations. It describes the N-particle system using a partial-differential equation for the probability density function, f = f(ξ₁, ξ₂,  ⋯ , ξ_N, t), in a 6N-dimensional space,

$\frac{\partial f}{\partial t}+\sum _{i=1}^N\left(\frac{{\mathit{p}}_i}{m}\cdot \frac{\partial }{\partial {\mathit{q}}_i}-\frac{\partial }{\partial {\mathit{q}}_i}{\Phi }^N\cdot \frac{\partial }{\partial {\mathit{p}}_i}\right)f=0,$ E1

where m is the particle mass, ζ_i = (q_i, p_i), q_i and p_i are the coordinate and momentum of the ith particle, Φ^N(q₁,  … , q_N) is the interaction potential of all the N particles. By integration over part of the variables, the Liouville equation can be transformed into a chain of equations, which is referred to as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [4], where the first equation connects the evolution of one-particle probability density function with the two-particle probability density function, the second equation connects the two-particle probability density function with the three-particle probability density function, and generally the Sth equation connects the S-particle probability density function, f_s,with the (S + 1)-particle probability density function, f_s + 1. Specifically,

$\frac{\partial f_S}{\partial t}+\sum _{i=1}^S\left(\frac{{\mathit{p}}_i}{m}\cdot \frac{\partial }{\partial {\mathit{q}}_i}-\frac{\partial }{\partial {\mathit{q}}_i}\Phi \left({\mathit{q}}_1,\dots ,{\mathit{q}}_s\right)\cdot \frac{\partial }{\partial {\mathit{p}}_i}\right)f_S=\frac{N-S}{V}\int \sum _{i=1}^S\left({\nabla }_{{\mathbf{r}}_i}{\Phi }_{i,S+1}\right)\cdot \left({\nabla }_{{\mathit{p}}_i}{f}_{S+1}\right)d^6{\zeta }_{S+1},$ E2

where V is the volume of system, $\Phi \left({\mathit{q}}_1,\dots ,{\mathit{q}}_s\right)=\sum _{i<j\le s}{\Phi }_{i,j}\left({\mathit{q}}_i,{\mathit{q}}_j\right)$ is the total two-particle interaction potential, Φ_{i , S + 1} is the interaction potential between the particle i and particle S + 1.

Here

$V^n={\int }_{\Gamma }{f}_n\left({\zeta }_1,{\zeta }_2,\cdots ,{\zeta }_n,t\right)d^6{\zeta }_1\cdots d^6{\zeta }_n,$ E3

and

$f_S\left({\zeta }_1,{\zeta }_2,\cdots ,{\zeta }_S,t\right)=V^{\left(S-N\right)}\int f_N\left({\zeta }_1,{\zeta }_2,\cdots ,{\zeta }_N,t\right)d^6{\zeta }_{S+1}\cdots d^6{\zeta }_N.$ E4

When S = N, the above equation recovers to the Liouville equation. For a macroscopic system, the particle number, N, is in the order of Avogadro constant, 10²³. Generally, there is no way to solve Eq. (1) or Eq. (2).

We need to simplify the model via considering some simpler cases. If considering only the case where correlations among three and more particles are negligible, the two-particle interaction is relevant to their distance, and the two-particle distribution function can be written as the product of two single particle distribution functions, specifically,

${\Phi }_{1,2}\left({\mathit{q}}_1,{\mathit{q}}_2\right)={\Phi }_{1,2}\left(\left|{\mathit{q}}_1-{\mathit{q}}_2\right|\right),$ E5

$f_2\left({\mathit{q}}_1,{\mathit{p}}_1;{\mathit{q}}_2,{\mathit{p}}_2;t\right)=f_1\left({\mathit{q}}_1,{\mathit{p}}_1;t\right)f_1\left({\mathit{q}}_2,{\mathit{p}}_2;t\right),$ E6

we obtain the Boltzmann equation [1, 4],

$\frac{\partial f}{\partial t}+\mathit{v}\cdot \frac{\partial f}{\partial \mathit{r}}+\mathit{a}\cdot \frac{\partial f}{\partial \mathit{v}}=Q\left(f,f\right),$ E7

with

$Q\left(f,f\right)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{0}{\overset{4\pi }{\int }}\left(f^{\ast }{f}_1^{\ast }-f{f}_1\right)\mathrm{g}\sigma \mathit{d\Omega }d{\mathit{v}}_1.$ E8

Here is the single particle distribution function. Compared with the MD and Liouville equation, the Boltzmann equation is a much coarse-grained model. The purposes of establishing Boltzmann equation are to define and calculate the entropy in non-equilibrium state, to derive, prove, even modify the fundamental hydrodynamic differential equation [5].

1.3. Non-equilibrium statistical mechanics and macroscopic description

All hydrodynamic quantities are some kinds of kinetic moments of the distribution function and can be expressed as:

$\mathbf{W}=\int f\left[\begin{array}{l}1\hfill \\ \mathit{v}\hfill \\ \mathit{v}\cdot \mathit{v}/2\hfill \end{array}\right]d\mathit{v},$ E9

where the particle mass has been assumed to be 1, W is the vector of macroscopic quantity composed of density, momentum and energy, W = [ρ, ρu, ρE]^T. Besides, the total stress σ, viscous stress Π and heat flux q are related to f via

$\mathbf{\sigma }=\int f\left(\mathit{v}-\mathit{u}\right)\left(\mathit{v}-\mathit{u}\right)d\mathit{v}=\mathit{\Pi }+p\mathit{I},$ E10

$\mathit{\Pi }=\int \left(f-f^{\mathit{eq}}\right)\left(\mathit{v}-\mathit{u}\right)\left(\mathit{v}-\mathit{u}\right)d\mathit{v},$ E11

and

$\mathit{q}=\frac{1}{2}\int \left(f-f^{\mathit{eq}}\right)\left(\mathit{v}-\mathit{u}\right)\left(\mathit{v}-\mathit{u}\right)\cdot \left(\mathit{v}-\mathit{u}\right)d\mathit{v},$ E12

respectively. Compared with Boltzmann equation, in Navier-Stokes equations,

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathit{u}\right)=0,$ E13

$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla \cdot \left(p\mathit{I}+\rho \mathit{uu}\right)=-\nabla \cdot \mathit{\Pi },$ E14

$\frac{\partial }{\partial t}\left(\mathit{\rho E}\right)+\nabla \cdot \left[\left(\mathit{\rho E}+p\right)\mathit{u}\right)]=-\nabla \cdot \left[\mathit{q}+\mathit{\Pi }\cdot \mathit{u}\right],$ E15

where Π = 2/3μ(∇ ⋅ u)I − μ(∇u)^T − μ∇u, q =  − κ∇T, T is the temperature, μ and κ are the coefficients of viscosity and heat conductivity, respectively. More details of molecular motions are neglected and the distribution function f disappears as well. What remained are the conserved kinetic moments, density ρ, momentum (density) ρu and energy(density) ρE, and a few non-conserved quantities, momentum fluxes ρuu and Π, energy fluxes ρEu and q. The relation between internal energy ρT and pressure p is given by the equation of state. In Euler equations, the viscous stress Π and heat flux q are further neglected.

From molecular dynamics to Boltzmann equation, to Navier-Stokes, and further to Euler equations, in each step, the description becomes coarse-grained, and the contained physical information becomes less [6]. The switching of model in each step corresponds to the state of system under consideration gets closer to its thermodynamic equilibrium, the behavior is simpler, consequently the system can be described by fewer physical variables. For a given non-equilibrium system, the switching of model in each step corresponds to the spatial-temporal scale that we use to observe the system becomes larger, consequently more smaller structures and quicker kinetic modes are invisible. What obtained are the remaining larger structures and slower kinetic modes. Based on Boltzmann equation, the most relevant TNE effects accompanying the hydrodynamic behaviors can be studied, in addition to the general hydrodynamic behaviors described by the hydrodynamic model.

2. Discrete Boltzmann theory

2.1. Discrete Boltzmann modeling

From Boltzmann equation to DBM, two steps of coarse-grained physical modelings are needed. The principle for coarse-grained modeling is that the physical quantities used to measure the system must keep the same values after simplification.

Step 1: Linearization of the collision term

Even though, compared with MD or Liouville equation, Boltzmann equation is a much coarse-grained model, its collision term is still too complicated to be solved in most practical cases. The simplest way to simplify is to introduce a local equilibrium velocity distribution function, f^eq, and write the collision term into the following linear form,

$\frac{\partial f}{\partial t}+\mathit{v}\cdot \frac{\partial f}{\partial \mathit{r}}+\mathit{a}\cdot \frac{\partial f}{\partial \mathit{v}}=-\frac{1}{\tau }\left(f-f^{\mathit{eq}}\right).$ E16

The physical meaning of the linearized collision model is thus: collisions of molecules result in that f approaches to thermodynamic equilibrium f^eq and the relaxation time is controlled by the parameter τ. If we do not aim to measure the system using the specific values of f, but using only some of its kinetic moments, then only if these kinetic moments keep invariable after the simplification, that will be OK. The linearization of the collision term requires

$\int -\frac{1}{\tau }\left(f-f^{\mathit{eq}}\right)\mathit{\Psi d}\mathit{v}=\int Q\left(f,f\right)\mathit{\Psi d}\mathit{v},$ E17

where Ψ = [1, v, vv, vvv, ⋯]^T contributes the kinetic moments used to measure the system. The specific form of f^eq depends on the terms that Ψ takes. According to the specific form of Ψ or f^eq, the linearized collision model may be referred to as the Bhatnagar-Gross-Krook (BGK) model [7, 8, 9, 10], ellipsoidal statistical BGK model [11], Shakhov model (for monatomic gas) [12], Rykov model (for diatomic gas) [13], Liu model [14], etc. When only the mass, momentum and energy conservation laws are kept, f^eq takes the simplest form, the Maxwellian. This is the so-called BGK model. Due to its simplicity, the BGK model is most extensively used.

It should be note that the Single-Relaxation-Time model works when all the kinetic modes approaching to thermodynamic equilibrium share more or less the same relaxation time. For more complicated cases where the relaxation times of different kinetic modes approaching to thermodynamic equilibrium are significantly different, the multiple-relaxation-time (MRT) collision model is needed [15].

Step 2. Discretization of the particle velocity space

We first consider the case without the force term. The discrete Boltzmann equation reads,

$\frac{\partial f_i}{\partial t}+v_{\mathit{i\alpha }}\frac{\partial f_i}{\partial r_{\alpha }}=-\frac{1}{\tau }\left(f_i-f_i^{\mathit{eq}}\right),$ E18

where i is the index of discrete velocity. Since the common schemes for discretizing the space and time do not work for discretizing the particle velocity space. To find an effective means to discretize the particle velocity space, we go back to consider what we really need. Here, we do not aim to describe the system using specific values of the discrete distribution function f_i, but its kinetic moments. So, only if these kinetic moments, originally in integral form, can be equally calculated in summation form, that will be acceptable. Specifically,

$\int f{\mathit{\Psi }}^{\text{'}}\left(\mathit{v}\right)d\mathit{v}=\sum _i{f}_i{\mathit{\Psi }}^{\text{'}}\left({\mathit{v}}_i\right),$ E19

where the left hand side gives the kinetic moments of f needed to describe the system and v_i at the right hand side is the discrete particle velocity. According to the Chapman-Enskog analysis, a kinetic moment of f can be finally calculated via an appropriate kinetic moment of f^eq. Therefore the requirement of Eq. (18) can be further written as

$\int f^{\mathit{eq}}{\mathit{\Psi }}^{\text{'}\text{'}}\left(\mathit{v}\right)d\mathit{v}=\sum _i{f}_i^{\mathit{eq}}{\mathit{\Psi }}^{\text{'}\text{'}}\left({\mathit{v}}_i\right),$ E20

where the left hand side gives the kinetic moments of f^eq being involved in the process of constructing DBM. The requirement of Eq. (20) can be rewritten as a matrix equation,

${\widehat{\mathit{f}}}^{\mathit{eq}}=\mathit{M}\cdot {\mathit{f}}^{\mathit{eq}},$ E21

where ${\widehat{\mathit{f}}}^{\mathit{eq}}\equiv \left[{{\widehat{f}}_k}^{\mathit{eq}}\right]$ , f^eq ≡ [f_i^eq], M ≡ [m_ki], ${{\widehat{f}}_k}^{\mathit{eq}}$ is the kth kinetic moment of f^eq. The elements, m_ki, are determined by the discrete velocities. That is to say, the discrete velocities should be chosen in such a way that the requirement of Eq. (20) is satisfied. Under such a constraint, the discretization of particle velocity space is flexible. In fact, M can also be non-squared rectangular matrix. The choice of discretization scheme depends on the compromise between the numerical cost, stability and physical gain (such as physical symmetry).

Via the same idea, it is straight forward to formulate MRT-DBM. To ensure the relaxation times have clear physical meanings, in the MRT-DBM, the collision term is first calculated in the kinetic moment space, and then transformed back to the discrete velocity space. To ensure DBM describe reasonable flow behaviors, a correction term is needed [15].

According to the Chapman-Enskog analysis, to access system which is deeper into thermodynamic non-equilibrium, higher order terms in Knudsen number should be considered. As a result, the requirement of f in Eq. (19) will lead to more kinetic moment relations of f^eq in Eq. (20). Consequently, more discrete velocities are needed. For the case with inter-particle force, one should generally first finish the derivative calculation of f with respect to v, and then perform the discretization of particle velocity space.

2.2. Non-equilibrium: definition and measuring

Once the concept of equilibrium is defined, the concept of non-equilibrium is clear. In classical mechanics, a particle is in mechanical equilibrium if the net force on that particle is zero. By extension, a physical system made up of many parts is in mechanical equilibrium if the net force on each of its individual parts is zero. In addition to defining mechanical equilibrium in terms of force, there are many alternative definitions for mechanical equilibrium which are all mathematically equivalent. In terms of momentum, a system is in equilibrium if the momentum of its parts is all constant. In terms of velocity, the system is in equilibrium if velocity is constant. In a rotational mechanical equilibrium, the angular momentum of the object is conserved and the net torque is zero. More generally in conservative systems, equilibrium is established at a point in configuration space where the gradient of the potential energy with respect to the generalized coordinates is zero. Similarly, a fluid system is in fluid mechanical equilibrium or hydrodynamic equilibrium if the net force on each of its ‘fluid particles’ (small fluid elements) is zero and without temperature gradient around the ‘fluid particle’.

For ideal gas system, in thermodynamic equilibrium there are no net macroscopic flows of matter or energy either within a system or between systems. In non-equilibrium systems, by contrast, there are net flows of matter or energy. Global thermodynamic equilibrium means that the relevant intensive parameters are homogeneous throughout the whole system, while local thermodynamic equilibrium means that those intensive parameters are varying in space and time, but are varying so slowly that, for any point, one can assume thermodynamic equilibrium in some neighborhood about that point. Rarefied gases at ordinary temperatures behave very nearly like ideal gas and the Maxwell speed distribution is an excellent approximation for such gases. Thus, it forms the basis of the kinetic theory of gases.

It is evident that Euler equations are used to investigate fluid flows which are at local thermodynamic equilibrium but mechanical non-equilibrium, while NS, Burnett and Super-Burnett equations are used to investigate fluid flows which are at mechanical and thermodynamic non-equilibrium. Only when all kinds of kinetic moments of (f − f^eq) are zero, the system is at thermodynamic equilibrium. hydrodynamic non-equilibrium (HNE) drives the evolution of f and results in viscous stress and heat flux.

In many practical cases, it is neither necessary to know all the details of f nor necessary to know all the kinetic moments of (f − f^eq).Therefore, one can (i) care only the main feature of f, specifically, neglect the high order terms of Knudsen number, (ii) care only the kinetic moments of (f − f^eq) which are most relevant to the macroscopic behaviors under consideration, specifically, those involved in constructing the discrete Boltzmann model.

A centrally important motivation of DBM is to check, measure and analyze the non-equilibrium state and effects [6, 16, 17, 18]. The DBM presents two sets of measures for the TNE. One set is dynamically from the difference of the kinetic moments of f and f^eq,

${\mathit{\Delta }}_m={\mathit{{\rm M}}}_m\left(f\right)-{\mathit{M}}_m\left(f^{\mathit{eq}}\right),$ E22

where M_m(f) is the mth rank moment of velocity tensor of f, ${\Delta }_m^{\ast }$ can be obtained when M_m(f) is the mth rank central moment. The other set includes the viscous stress and heat flux. The former describes the specific TNE status, while the latter describes the influence of those TNE to the macroscopic control equations. The former is local, while the latter is non-local. The former is finer, while the latter is coarser. At each step of DBM simulation, both the f and f^eq are calculated. Therefore, the TNE effects are naturally included in each step. The captured TNE effects are just those being most relevant to the macroscopic flow behaviors under consideration.

Entropy production is a highly concerned quantity in both physics and engineering studies. From the physics side, it is helpful for understanding the complex non-equilibrium behaviors. From the engineering side, a process with lower entropy production may have higher energy transformation efficiency. Following the way of defining entropy equilibrium equation in the non-equilibrium thermodynamics, a new entropy equilibrium equation can be obtained as follows [19],

$\frac{\partial s}{\partial t}=-\nabla ·{\mathit{J}}_s+\sigma ,$ E23

where s is the entropy density,

${\mathit{J}}_s=s\mathit{u}+\frac{1}{T}{\mathit{\Delta }}_{3,1}^{\ast }$ E24

and

$\sigma ={\mathit{\Delta }}_{3,1}^{\ast }·\mathit{\nabla }\frac{1}{T}-\frac{1}{T}{\mathit{\Delta }}_2^{\ast }:\mathit{\nabla }\mathit{u}+\rho \frac{Q}{T}F\left(\lambda \right)$ E25

are the entropy flux and entropy production rate, respectively. $F\left(\lambda \right)=\dot{\lambda }$ is a rate function describing the process of chemical reaction occurring in the system, where λ is defined as the ratio of the product density to the overall density. From Eq. (24) one can find that the entropy production results from three kinds of resources, the non-organized energy flux (NOEF), non-organized momentum flux (NOMF), and chemical reaction. From the relation, one can study the various mechanisms resulting in increase of entropy and their relative importance [19].

The TNE behaviors are very complex and difficult to quantitatively investigate. Finding a convenient and efficient method to characterize the TNE status and effects is the corner stone, DBM presents such a potential approach [16, 17, 18, 19, 20, 21, 22, 23, 24, 25].

2.3. DBM versus CFD

The traditional CFD needs to first know the exact form of the hydrodynamic equations, then design numerical algorithm according to the properties of those equations. The DBM is a coarse-grained model derived from the Boltzmann equation. In principle, it can be formulated and applied to simulate flows without knowing the exact form of the hydrodynamic equations, only if necessary kinetic moments of f^eq are followed. From the perspective of physical application, a DBM is roughly equivalent to a hydrodynamic model supplemented by a coarse-grained model of TNE. Via the DBM, it is easy to perform multi-scale simulations in a wide range of Knudsen number.

3. Applications

3.1. Combustion system

Combustion has long been playing a dominant role in the transportation and power generation. To improve combustion efficiency and decrease pollution, in recent years, some new combustion concepts have been proposed. For example, pulsed detonation engine, spinning detonation engine, microscale combustion, nanopropellent, partially premixed and stratified combustion, plasma-assistant combustion, cool flames, etc. All these new combustion concepts involve complicated non-equilibrium chemical and transport processes [26].

The chemical reaction process is very complex and may include varieties of reaction mechanism. So far, most of the chemical reaction kinetic models are phenomenological. As the first step, we consider only the simplified form of Lee-Tarver chemical reaction rate law [21]. Considering the thermal initiation, the reaction kinetic is described by

$\frac{\mathit{d\lambda }}{\mathit{dt}}=\left\{\begin{array}{l}a\left(1-\lambda \right)+b\left(1-\lambda \right)\lambda ,T\ge T_{\mathit{th}}\mathrm{and}0\le \lambda \le 1\\ 0,\text{else},\end{array}\right.$ E26

where a and b are constants, λ is the concentration of the product and works as the reaction process parameter, T_th is the temperature threshold for chemical reaction. Consider the case where the chemical reaction is slow enough compared with the kinetic process of approaching thermodynamic equilibrium, so we can treat f as f^∗eq during the reaction process. The evolution equation of single-relaxation-time DBM for combustion reads

$\frac{\partial f_i}{\partial t}+v_{\mathit{i\alpha }}\frac{\partial f_i}{\partial r_{\alpha }}=-\frac{1}{\tau }\left(f_i-f_i^{\ast \mathit{eq}}\right),$ E27

where $f_i^{\ast \mathit{eq}}=f_i^{\mathit{eq}}\left(\rho ,u,T+\left(\gamma -1\right)\mathit{\tau QF}\left(\lambda \right)\right)$ is the local equilibrium distribution function taking into account chemical reaction effect. Q is the amount of heat released by the chemical reactant per unit mass. If the relaxation times of different kinetic modes approaching to thermodynamic equilibrium are significantly different, one needs the MRT model [22]. Some observations brought by DBM are as below.

Non-equilibrium quantities defined in Eq. (22) are used to study a simple case of detonation [21]. For the case of CJ detonation shown in Figure 1(a), the corresponding non-equilibrium quantities, ${\mathit{\Delta }}_2^{\ast }$ and ${\mathit{\Delta }}_3^{\ast }$ , are shown in Figure 1(b) and (c), respectively. Interestingly, at the position of von Neumann pressure peak, the system is not far from but is nearly in its thermodynamic equilibrium. The internal energies in two degrees of freedom deviate from their mean value in opposite direction with the same amplitude. The internal energy in each degree of freedom deviates from the mean value in opposite direction before and after the von Neumann pressure peak. The amplitude in front of the von Neumann pressure peak is larger. The physical reasons are as follows. The mechanical non-equilibrium is the driving force for TNE. When a strong shock comes, the density, temperature and flow velocity increases abruptly so that the system does not have enough time to relax to its thermodynamic equilibrium. With decrease of the gradients of density, temperature and flow velocity, the system gets relatively more time for thermodynamic relaxation. At the position of von Neumann pressure peak, the system has been close to its thermodynamic equilibrium. After the von Neumann pressure peak, the density and flow velocity decrease quickly so that the system does not have enough time for thermodynamic relaxation again. The total deriving force for TNE in front of the von Neumann pressure peak is larger than behind.

Figure 1.

Profiles of physical quantities for CJ detonation. (a) The density ρ, pressure P, temperature T, x-component of velocity u, and the reaction process λ. (b) and (c) The non-equilibrium quantities ${\Delta }_2^{\ast }$ and ${\Delta }_3^{\ast }$ , respectively (see Ref. [21] for more details).

Figure 2 gives the pressure P and corresponding non-equilibrium quantity Δ_2 , xx at a certain time [22]. From Figure 2(a)–(c), the relaxation time, 1/R_i, decreases by 10 times in each case so that the detonation wave changes from unsteady to steady. In Figure 2(d)–(f), the shaded area, enclosed by the curve Δ_2 , xx and the x-axis, presents a rough description on the global TNE effect around the detonation wave. It is clear that the viscosity (and/or heat conductivity) decreases the local TNE while increases the global TNE around the detonation wave.

Figure 2.

Profiles of physical quantities for different relaxation times. (a)–(c) The pressure P. (d)–(f) The corresponding non-equilibrium quantities Δ_2 , xx (from Ref. [22] with permission).

Figure 3 shows some numerical results aiming to investigate the main mechanisms resulting in entropy increase and their relative importance in the combustion system [19]. It is clear that, in the checked cases, the most pronounced contribution to entropy increase is from the chemical reaction, Δs_CHEM, the lest contribution is from the non-organized energy flux, Δs_NOEF, the contribution of non-organized momentum flux, Δs_NOMF, is in between. With the increasing of Mach number, the entropy production caused by non-organized momentum flux becomes more remarkable.

Figure 3.

Mechanisms for global entropy production in four cases (from Ref. [19] with permission).

3.2. Multiphase flow with phase separation

Phase separation is an important branch in the field of multiphase flows. It is also a kind of non-equilibrium phase transition. The key step for modeling phase separation is to incorporate the non-ideal gas effects into the discrete Boltzmann equation. Enskog equation can be regarded as an extension of Boltzmann equation under the hard-ball molecule model [1]. Although the specific treatments may be different, the aims are the same. Those are all to replace the equation of state of ideal gas with a more practical one.

In 2007 Gonnella, Lamura, and Sofonea (GLS) [27] introduced an appropriate inter-particle interaction to the Watari’s model [28] to describe van der Waals fluids. The evolution equation of GLS model reads:

$\frac{\partial f_{\mathit{ki}}}{\partial t}+v_{\mathit{ki\alpha }}\frac{\partial f_{\mathit{ki}}}{\partial r_{\alpha }}=-\frac{1}{\tau }\left(f_{\mathit{ki}}-f_{\mathit{ki}}^{\mathit{eq}}\right)+I_{\mathit{ki}},$ E28

where the external force term I_ki takes the following form:

$I_{\mathit{ki}}=-\left[A+B_{\alpha }\left(v_{\mathit{ki\alpha }}-u_{\alpha }\right)+\left(C+C_q\right){\left(v_{\mathit{ki\alpha }}-u_{\alpha }\right)}^2\right]f_{\mathit{ki}}^{\mathit{eq}}$ E29

In a recent work [20], the GLS model was further developed to be a kinetic model which can be used to access both the hydrodynamic non-equilibrium and the thermodynamic non-equilibrium. To roughly and averagely estimate the derivation amplitude from the thermodynamic equilibrium, a TNE strength can be defined as

$D=\sqrt{{\mathit{\Delta }}_2^{\ast 2}+{\mathit{\Delta }}_3^{\ast 2}+{\mathit{\Delta }}_{3,1}^{\ast 2}+{\mathit{\Delta }}_{4,2}^{\ast 2}}$ E30

Figure 4 shows that the maximum value point can work as a physical criterion to discriminate the two stages, spinodal decomposition and domain growth, of phase separation. The TNE strength increases with time in the first stage while decreases with time in the second stage. More details are referred to the original publication [20].

Figure 4.

Evolutions of the boundary length L and the TNE strength D for the phase separation process (see Ref. [20] for more details).

3.3. Rayleigh-Taylor interfacial instability

Rayleigh-Taylor instability (RTI) occurs at the interface between two fluids with different densities. The compressible RTI system can be described by [23, 24]

$\frac{\partial f_i}{\partial t}+v_{\mathit{i\alpha }}\frac{\partial f_i}{\partial r_{\alpha }}-\frac{a_{\alpha }\left(v_{\mathit{i\alpha }}-u_{\alpha }\right)}{\mathit{RT}}{f}_i^{\mathit{eq}}=-\frac{1}{\tau }\left(f_i-f_i^{\mathit{eq}}\right)$ E31

where a_α is the acceleration due to external force.

Generally, the depth of the mixing layer is an important parameter to measure the evolution of RTI. For incompressible RTI, the measurement is readily performed by tracing the constant density. However, for the compressible case, how to measure the mixing layer remains a thorny problem. Here we present two independent interface-tracking methods. One is by tracking the mean temperature of the upper and bottom fluids while the other is by tracking the maximum values of TNE characteristic quantities, such as ${\mathit{\Delta }}_{3.1,y}^{\ast }$ . The latter method is based on the fact that ${\mathit{\Delta }}_{3.1,y}^{\ast }$ takes its maximum value at the position of the interface along the y direction of the spike and bubble. The perturbation amplitudes developing with time obtained by the two methods are shown in Figure 5. The good agreement shows that the two approaches validate each other and the local TNE, ${\mathit{\Delta }}_{3.1,y}^{\ast }$ , can be used to track interfaces in numerical experiments.

Figure 5.

The perturbation amplitudes developing with time obtained by two different tracking approaches (see Ref. [23] for more details).

3.4. Compressible flows under shock

Figure 6 shows the profiles of physical quantities in the process where the shock wave passes outwards from the heavy medium to the light one. Figure 6(a) and (b) are for the case without and with initial perturbations at the interface, respectively. From left to right, one can find three kind of interfaces, the rarefaction wave, material interface and shock wave, which are indicated by dashed lines [25].

Figure 6.

Profiles of physical quantities in the process of shock wave passing outwards from the heavy medium to the light one. (a) Without initial perturbations at the interface. (b) With initial sinusoidal perturbation at the interface (from Ref. [25] with permission).

According to the TNE information, the main feature of actual particle velocity distribution function can be qualitatively recovered. Figure 7 shows an example, where the interface is not perturbed initially. The details are referred to Ref. [25]. DBM simulations [25] show that the shear stress exists only for the oblique shocking. As a consistent correspondence, MD results [29] show that fluctuating shear stresses exist if observed in a scale with a few angstroms, while their mean value becomes negligibly small when being averaged over a scale with several hundred angstroms.

Figure 7.

Sketches of the actual distribution functions in velocity space (v_r, v_θ). Panels (a)–(c) show the recovered distribution functions at the rarefaction wave, the material interface, and the shock wave, respectively (from Ref. [25] with permission).

A comparison of the DBM results and those of the MD is shown in Figure 8. Figure 8(a) and (b) shows the DBM results for the cases where the interface is not and is perturbed initially, respectively. Figure 8(c) shows the shear stresses from MD simulation. Since the MD uses particle description, the existence of locally fluctuating shear stress corresponds to the observation in the case with perturbed interface; the observation that mean value of shear stresses becomes negligibly small in a much larger scale roughly correspond to the case with non-perturbed interface for the DBM simulations.

Figure 8.

Shear stress within shock wave. (a) DBM results for the case without initial interface perturbation. (b) DBM results for the case with initial interface perturbation. (c) MD results (figures (a) and (b) are from Ref. [25] with permission, figure (c) is from Ref. [29] with permission).

3.5. Shock wave in plasma

Figure 9 shows an example for that the TNE effects can be used to physically discriminate shock wave in plasma from those in common fluid. From the first two rows, the two TNE quantities, ${\mathit{\Delta }}_2^{\ast }$ and ${\mathit{\Delta }}_{4,2}^{\ast }$ , show quite similar behaviors around shock wave and/or detonation wave in common fluid, even though they two have different physical meanings. However, the two quantities show qualitative difference around shock wave in plasma [30].

Figure 9.

TNE characteristics of three types of shock waves. The three rows, from top to bottom, correspond to pure shock wave (see Ref. [25] for more details), shock wave with chemical reaction (detonation) (see Ref. [21] for more details), and shock wave in plasma (see Ref. [30] for more details). The profiles on the left column show the values of ${\Delta }_2^{\ast }$ around the wave fronts while the profiles on the right column show the values of ${\Delta }_{4,2}^{\ast }$ .

4. Summary

Understanding compressible flows need more time. DBM presents a convenient way to model and simulate systems with trans-scale Knudsen number. Mathematically, the only difference of discrete Boltzmann from the traditional hydrodynamic modeling is that the NS equations are replaced by a discrete Boltzmann equation. But physically, this replacement has a significant gain: a DBM is roughly equivalent to a hydrodynamic model supplemented by a coarse-grained model of the TNE, where the hydrodynamic model can be and can also beyond the NS. The TNE provided by DBM has been used to investigate non-equilibrium effect during detonation process, to discriminate different stages of phase separation, to recover actual particle velocity distribution function qualitatively, to track the interfaces of different fluid, and to discriminate shock wave in plasma from those in common fluid. More use of those TNE quantities are further being discovered with the deeper investigation of the compressible and complex flows.

Acknowledgments

We warmly thank Profs. Wei Kang, Zhihua Chen, Yanbiao Gan, Feng Chen, Chuandong Lin, Huilin Lai, Zhipeng Liu, Bo Yan, et al. for helpful discussions. We acknowledge support of the National Natural Science Foundation of China (under Grant Nos. 11475028 and 11772064), Science Challenge Project (under Grant No. JCKY2016212A501 and TZ2016002), and Science Foundation of Laboratory of Computational Physics.