Advances in the kinetics of heat and mass transfer in near-continuous complex flows

Abstract

The study of macro continuous flow has a long history. Simultaneously, the exploration of heat and mass transfer in small systems with a particle number of several hundred or less has gained significant interest in the fields of statistical physics and nonlinear science. However, due to absence of suitable methods, the understanding of mesoscale behavior situated between the aforementioned two scenarios, which challenges the physical function of traditional continuous fluid theory and exceeds the simulation capability of microscopic molecular dynamics method, remains considerably deficient. This greatly restricts the evaluation of effects of mesoscale behavior and impedes the development of corresponding regulation techniques. To access the mesoscale behaviors, there are two ways: from large to small and from small to large. Given the necessity to interface with the prevailing macroscopic continuous modeling currently used in the mechanical engineering community, our study of mesoscale behavior begins from the side closer to the macroscopic continuum, that is from large to small. Focusing on some fundamental challenges encountered in modeling and analysis of near-continuous flows, we review the research progress of discrete Boltzmann method (DBM). The ideas and schemes of DBM in coarse-grained modeling and complex physical field analysis are introduced. The relationships, particularly the differences, between DBM and traditional fluid modeling as well as other kinetic methods are discussed. After verification and validation of the method, some applied researches including the development of various physical functions associated with discrete and non-equilibrium effects are illustrated. Future directions of DBM related studies are indicated.

Graphical abstract

Keywords

near-continuous flow / non-equilibrium / kinetics / discrete Boltzmann method / complex physical field analysis

1 Introduction

Fluid is a general term encompassing both gases and liquids. The two fluids we are most familiar with are the atmosphere and water: the atmosphere surrounds the entire Earth, and 70% of the Earth’s surface is covered by water. Since ancient times, water and air, on the one hand, have nurtured life. On the other hand, it also brings natural disasters such as typhoons, floods and droughts to human beings from time to time. Human daily life, production and other activities are inseparable from fluid. In the struggle with nature and through production practice, human beings gradually accumulate experience and enhance understanding of fluid behavior. The continuous summarization and enhancement of these understandings gradually formed the discipline of fluid mechanics. The air and water we perceive in our daily lives appear to be continuous, and this perception is based on a macroscopic understanding of continuous flow. However, with the emergence of molecular kinetic theory, people gradually realized that the sense of continuity is because the distances between molecules are exceedingly small compared to the scales we typically perceive in our daily lives.

In addition to the seemingly continuous flow of water and air, there exist various flows in our daily lives and industrial production that exhibit distinct discrete characteristics. Examples include dust explosions, particle flow, traffic flow, and the movement of crowds. For a uniform description, we treat people and cars as particles. These flows are perceived as discrete because the distance between particles is no longer negligible relative to the scale of our everyday perception. However, when observed on a significantly larger scale (for example, exceeding 1000 times the distance between particles), these flows can also be considered continuous. It can be seen that the fluid system itself is objective, but the description of continuous and discrete is dependent on the scale of the behavior we are concerned with. Hsue-shen Tsien, a famous Chinese scientist, first suggested that according to the Knudsen (Kn) number, the fluid behavior could be divided into continuous flow, slip flow, transition flow and free molecular flow [1]. In addition, it should be noted that unlike microscopic molecular collisions, collisions between particles composed of a large number of molecules are dissipative.

The development of scientific theories typically follows a progression from simplicity to complexity. The theory of fluid mechanics begins with a macroscopic continuous flow such as air and water in which collisions between microscopic particles do not dissipate. Considering the stage of the study and space limitation of current review, the continuous and discrete flows discussed in the following parts are “simple” cases in which there is no dissipation in the collisions of microscopic particles which constitute the fundamental units. Interestingly, and somewhat, the study of heat and mass transfer in these “simple” cases has been, and is still going through a long process.

To show that the study of these simple scenarios also contains a wealth of yet unexplored complex behaviors, we begin with an everyday phenomenon: the fluid (interface) instability. Fluid (interface) instability refers to the phenomenon that after the fluid interface is disturbed, the interface cannot automatically return to its original state under the action of the system’s own force, and the disturbance amplitude increases with time. Two common types of fluid instabilities are the Kelvin−Helmholtz instability (KHI) [2, 3] and Rayleigh−Taylor instability(RTI) [4, 5]. KHI occurs when there is a difference in tangential velocity on both sides of the interface, leading to phenomena like rolling sea waves and cirrus clouds. RT instability refers to the phenomenon of interface instability that can be caused by any disturbance when acceleration is directed from a medium with high density to one with low density. Due to the relativity of motion, another equivalent image of RT instability is the phenomenon of interface instability caused by the acceleration of a denser medium by a less dense medium. Common examples of RT instability, such as interface instability when milk is poured into a glass of water, volcanic eruptions, and mushroom clouds. Fig.1 shows the mechanical causes of the formation of KHI and RTI. In addition, there is a class of fluid (interface) instability known as Richtmyer−Meshkov instability (RMI) [6, 7]. It refers to the phenomenon that the interface disturbance cannot naturally recover by the system itself under the action of shock wave, and the disturbance amplitude increases gradually with time. In fact, RMI can be regarded as a special case of RTI when the acceleration takes the form of a pulse. At the same time, during the evolution of RTI and RMI, a tangential velocity difference on both sides of the interface will inevitably induce the KHI.

Fig.1 (a) Sketch of the formation mechanism of KHI. (b) Sketch of the formation mechanism of RTI.

Full size|PPT slide

The three kinds of fluid instabilities, RTI, RMI, and KHI, have been the focus of high energy density physics laboratories and inertial confinement fusion (ICF) laboratories around the world [8-13]. In fact, fluid instability phenomenon widely exists in nature and engineering technology fields, such as astrophysics, energy and power, aerospace, chemical engineering and materials science [14-24]. As a physical phenomenon that must occur when certain conditions are met, fluid instability is also a double-edged sword. On one hand, it can accelerate the mixing of substances, which is conducive to the mixing and combustion of liquid fuels in internal combustion engines, aeroengines and supersonic ramjets. It also plays a critical role in initiating explosives. On the other hand, it is also a key factor that seriously affects the success of ICF ignition and some equipment performance. It poses a potential threat to the safety of explosives. The guiding principle here is straightforward: if it is beneficial, strengthen it; if it is harmful, try to suppress it. The basis of effective utilization and suppression is to have a clear understanding of the characteristics, mechanisms and laws in the process of its occurrence and development. These systems often have the following characteristics: (i) Despite their macroscopic scale, they encompass a large number of intermediate-scale spatial structures and kinetic patterns. The existence and evolution of these structures and models greatly affect the physical properties and functions of the system. (ii) Such systems often have numerous internal interfaces, including material interface and mechanical interface (shock wave, rarefaction wave, detonation wave, etc.). (iii) The internal force and response processes within the system are exceedingly complex. Furthermore, the study of these systems is faced with the following problems: The large-scale slow behavior of these systems can be reasonably described by the Navier−Stokes (NS) equations. However, in the description of flow behavior in some low-pressure rarefaction regions, the description of internal structure of shock wave or detonation wave, and the description of non-equilibrium behavior caused by fast changing flow or reaction, NS equations are insufficient in physical function. At the same time, for the flow behavior of interest, microscopic Molecular Dynamics (MD) simulations are often powerless due to the limitations of applicable (spatiotemporal) scales. As shown in Fig.2, for the mesoscale dilemma where the macroscopic continuous NS models suffer from inadequate physical capabilities, while micro MD simulations are unable to reached the scale, the corresponding research is extremely weak due to the lack of suitable models and methods. Compared with the macroscopic continuous case, the typical characteristics of mesoscale behavior are as follows: the discrete effect is more significant, and the non-equilibrium effect is more significant. Currently, mesoscale kinetic modes related to discrete and non-equilibrium phenomena are far from being fully understood, which significantly hinders the assessment of various effects in mesoscale behavior and the development of corresponding control techniques.

Fig.2 The mesoscale dilemma between macroscopic continuous and microscopic particle descriptions.

Full size|PPT slide

Similar problems are faced in the following situations: In the field of aerospace, spacecraft traverse the rarefied atmosphere during its launch and reentry. With the increase of the height in the rarefied region, the fluid behavior predicted by NS deviates increasingly from experimental observations. Even at the same altitude, the reasonability of NS describing practical flow behaviors in different parts of the same spacecraft may also be different [1, 25, 26]. In the field of micro-porous flow and microfluidics, behaviors such as heat transfer, mass transfer, and phase changes within microchannels often show significant deviations from macroscopic NS descriptions [27-30]. Descriptions of complex flow processes, such as jet atomization, bubble collapse under shock, and droplet fragmentation under laser ablation, challenge the rationality of the NS model [31-36]. In combustion environments, the phase transition and mixing of liquid and gas [37, 38], as well as complex flow processes like droplet vaporization under laser ablation and laser ablation propulsion, pose challenges to the rationality and physical capabilities of NS modeling [9, 39, 40]. The internal non-equilibrium phenomena of plasma system are very rich and complex, and hydrodynamic description is an important means to study the behavior of plasma. However, the in-depth study of a series of strong non-equilibrium behaviors poses a challenge to the rationality and capabilities of NS modeling [40-45]. It has been pointed out that in ICF, the non-continuous and highly non-equilibrium kinetic behaviors may significantly affect the success of ignition [40, 44, 45].

The physical reason for the above dilemma is that the traditional NS model is based on the assumption of continuous media and near equilibrium approximation. However, near sharp interfaces such as shock waves and boundary layers, the average molecular spacing is no longer a negligible small quantity relative to the scale of the structure we are concerned about. In terms of describing fast changing modes, the thermodynamic relaxation time is no longer a negligible small quantity compared to the time scale of the flow (or reaction) we are concerned with. During the process of flow (or reaction), the system no longer has enough time to return to near thermodynamic equilibrium. In addition, traditional fluid modeling only focuses on/describes the spatiotemporal evolution of the conserved moments of the distribution function (density, momentum, and energy). Subsequently, we will further discuss that this approach is inadequate for describing micro−mesoscale structures and fast-changing modes, and its shortcomings become more pronounced as the degree of discreteness and non-equilibrium increases.

There are only two ways to advance to the mesoscale science: from macro to meso and from micro to meso. The former involves scaling down from larger to smaller, while the latter involves scaling up from smaller to larger. Given that the fluid models predominantly used in engineering applications such as aerospace, ICF and some other equipment are mainly macroscopic continuous NS models, the issue we are facing is: how to seamlessly connect mesoscale modeling with macroscopic continuous modeling? Therefore, in the discussions within this article, we adopt a scheme that gradually increases the degree of discreteness and non-equilibrium from macro to meso, with a primary focus on behaviors closer to the macroscopic side within the mesoscale range.

2 State of art for complex flow

2.1 Problems and challenges

The discontinuity (or discreteness) of a fluid system is closely related to the degree of Thermodynamic Non-Equilibrium (TNE). Knudsen (Kn) number, on the one hand, can be regarded as a rescaled average molecular distance, describing the discontinuity (or discreteness) of the system. On the other hand, it can also be regarded as a rescaled thermodynamic relaxation time to describe the degree of TNE of the system. The concepts of non-equilibrium flow and discontinuous flow (discrete flow) overlap from some sense in physical connotation. The study of non-equilibrium complex flows has made great progress, but still faces some fundamental scientific challenges, such as:

(i) Cross-scale modeling and simulation has been a research hotspot for about 15−20 years. During this period, it has experienced fluctuations in popularity. The main methodology is to modify the macroscopic fluid/solid equations based on (quasi) continuity and (near) equilibrium assumption [46]. That is, the physical quantities used to describe the system’s state and behavior are still based on those used in traditional fluid/solid mechanics equations. A very natural fundamental scientific question arises: As the degree of discreteness/non-equilibrium increases, is it really no problem to focus only on the few physical variables used in a continuous and near-equilibrium model (e.g., NS)?

(ii) It is a typical feature that complex structures and behaviors generally take on different characteristics when viewed from different angles. “How to extract more valuable information” combined with “how to analyze it” determines our research capabilities and depth. Among them, the technical key is as below: when dealing with systems that are becoming increasingly complex, how to achieve an intuitive geometric correspondence for describing complex system states and behaviors?

(iii) Turbulence research remains an enduring and important research focus due to its significance. Turbulent mixing is an important aspect of hydrodynamic instability research. It can be thought that as the vortex and other structures concerned become smaller and smaller, the average molecular distance is no longer a negligible small quantity relative to the scale of the structure or behavior concerned, that is, the discreteness becomes more pronounced, and discrete effects become more significant. However, the early concepts and theoretical framework of turbulence are established on the continuum theory [47]. So, the third fundamental scientific question is as follows: Is the turbulence concept derived from pure mathematical deductions based on a macroscopic continuous image entirely consistent with that in pursuit of physical origin?

(iv) At present, the studies on non-linear systems [48, 49] and heat and mass transfer in small systems have made rich progress [50-58]. Heat and mass transfer in small systems exhibit some behaviors that seem “anomalous” (different from the macroscopic continuous cases) and have attracted significant attention in the fields of statistical physics and nonlinear science [59]. However, it is worth noting that the small systems extensively studied in the literature of statistical physics and nonlinear science are often much smaller than the mesoscale non-equilibrium situations of interest in the mechanical engineering field. The situation we face is that engineering applications largely rely on macroscopic continuous modeling, but also encounter unsatisfactory or even unsolvable problems. The question we face is how to seamlessly connect mesoscale modeling with macroscopic continuous modeling?

(v) In the past years, the main idea of cross-scale modeling and simulation is based on reductionism. As a scientific method, reductionism is significant in uncovering the internal mechanisms of complex systems and phenomena, determining causality between the whole and its parts, and achieving precise modeling of complex systems. However, a typical characteristic of complex systems is that as complexity increases, the system may exhibit emergent novel behavioral features that cannot be derived from any constituent or coexistence rules from simpler or simplest cases. (Emergence, the idea that the whole is greater than the sum of its parts, is a central theme in complexity science.) So, the fifth fundamental scientific question that DBM needs to face is: Is the reductionist approach alone sufficient for cross-scale modeling and simulation?

(vi) The sixth basic science problem relates to the most basic parameter of non-equilibrium flow description/cognition, that is the non-equilibrium strength. Given the complexity of non-equilibrium behavior, any definition of non-equilibrium strength depends on the research perspective. In other words, there may be the following situations: the non-equilibrium strength is increasing from a certain perspective, while from another perspective, it appears to be decreasing. How to describe the non-equilibrium strength of complex flows accurately?

2.2 Ideas and schemes

The problems discussed above are all problems where the complexity of the problem challenges the physical function of the model. The problem of physical function cannot be solved by improving algorithm accuracy. Multi-scale/cross-scale modeling and simulation are efforts to address these and similar problems.

As shown in Fig.3, simulation study of complex flow includes three major steps, (i) physical modeling, (ii) discrete format selection/design, (iii) numerical experiment and complex physical field analysis. As users of discrete formats, the focus of the physical research groups generally lies in the two ends of (i) and (iii), leaving the design of the optimal discrete format in step (ii) to more specialized computational mathematicians. The core tasks of the physical research group are as follows: (i) Ensure the validity of the physical model (theoretical model) for the problem under investigation while maintaining simplicity, and (ii) try to extract more valuable physical information for the complex physical field from the massive data.

Fig.3 The three main steps of complex flow simulation research.

Full size|PPT slide

People’s understanding of fluids begins with continuous and (near) equilibrium flow. As the degree of discreteness and nonequilibrium increases, the complexity of system behavior sharply increases, which poses challenges for the description and control of corresponding flow behavior. Our understanding of fluid behavior has two progression modes: one is the macroscopic cognitive progression mode centered on conserved quantities, and the other is the microscopic cognitive progression mode centered on molecular dynamics. People’s understanding of fluid behavior begins with macroscopic continuous images. The order of scientific research is generally from conserved quantities, to spatiotemporally slow variables, followed by gradually increasing fast variables. Early research on conserved quantities and slow variables prepares cognitive and technical foundations for gradually increasing fast variable research in the later stages. The space slowly changes, corresponding to the large structure, and the relative average molecular spacing is small enough to be negligible. At this point, the assumption of continuity is reasonable. A slow change in time indicates that the system has enough time to recover to the thermodynamic equilibrium state, and the thermodynamic (near) equilibrium assumption is reasonable at this time. Traditional fluid modeling is based on and describes the three conservation laws of mass, momentum, and energy. However, as the rate of spatial change increases, that is, with a decrease in the focus structure, the relative average molecular spacing gradually becomes less negligible, and the discrete effect becomes more significant. Alternatively, as the rate of change over time increases, during the flow or reaction process, the system may not have enough time to return to thermodynamic equilibrium, resulting in a higher degree of thermodynamic non-equilibrium. At this point, the rationality of the two basic assumptions of traditional fluid modeling is increasingly challenged. How to model and analyze has become a problem to be solved. The path from large to small scale has encountered a theoretical rationality bottleneck.

Molecular dynamics is theoretically complete, but it requires extremely high computational resources, so the applicable spatiotemporal scale is very small. The behaviors we focus on often occur at spatiotemporal scales that are beyond the reach of molecular dynamics. So, on the path from small to large scale, we encountered a bottleneck in scale.

Statistical physics is a bridge that connects microscopic and macroscopic aspects. It utilizes a probabilistic approach based on the understanding of the microstructure of matter and interactions between microscopic particles to provide a microscopic explanation for the physical properties and macroscopic laws of objects composed of a large number of particles. It is also known as statistical mechanics. Statistical physics is both a theory of physics and a methodology. The microscopic particles described in its framework can be molecules in gas, atoms in crystals, photons in laser beams and other conventional microscopic particles, as well as stars in the Milky Way galaxy, cars on highways, sheep in flocks, people in social groups, and so on. The kinetic theory in non-equilibrium statistical physics is the foundation of macroscopic fluid mechanics theory. The kinetic theory based on the Boltzmann equation is a relatively complete part of non-equilibrium statistical physics [60]. In response to the “mesoscale” dilemma of insufficient physical functionality of macro models and limited applicability of micro models, as shown in Fig.2, developing kinetic models based on the Boltzmann equation has become the primary approach.

2.3 Statistical physics: Coarse-grained modeling and phase space description

The fundamental equation of non-equilibrium statistical physics is the Liouville equation:

(1)${\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }t}}+{\displaystyle \sum _{i=1}^N\left({\dot{\mathit{q}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{q}}_i}}+{\dot{\mathit{p}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{p}}_i}}\right)F=0.}$

It uses the $N$-body distribution function $F=$$F({\mathit{q}}_1,{\mathit{q}}_2,\cdots ,{\mathit{q}}_N,{\mathit{p}}_1,{\mathit{p}}_2,\cdots ,{\mathit{p}}_N,t)$ and its evolution equation to describe the state and behavior of the system, where $({\mathit{q}}_i,{\mathit{p}}_i)$ are the generalized coordinate and generalized momentum of the ith particle. The Bogoliubov−Born−Green−Kirkwood−Yvon (BBGKY) hierarchy chain is equivalent to Liouville equation. Equivalent to the description of molecular dynamics without any simplification, it can be regarded as the holographic description based on the molecular level, which is the starting point of theoretical thinking [60]. For a macroscopic system, $N$ is the order of magnitude of the Avogadro constant, i.e. ${10}^{23}$. In general, the Liouville equation (1) is fundamentally unsolvable and cannot be directly applied. So we have to simplify it. It is impossible and, in most cases, unnecessary to capture all the details of a system’s state and behavior simultaneously. We need to break down complex problems, according to the research needs, to grasp the main contradiction. The direction of simplification involves making the governing equation depend on a distribution function with fewer bodies. What about discarded information? According to the mean-field approach, a correction term is added to the simplified governing equation to compensate for the lost information, ensuring that the kinetic properties to be studied remain unchanged.

Specifically, since the $N$-body distribution function evolution equation (1) cannot be processed, we first thought of simplifying it to the ($N$−1)-body distribution function evolution equation:

(2)${\displaystyle \frac{\mathrm{\partial }{f}_{N-1}}{\mathrm{\partial }t}}+{\displaystyle \sum _{i=1}^{N-1}\left({\dot{\mathit{q}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{q}}_i}}+{\dot{\mathit{p}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{p}}_i}}\right)f_{N-1}=C_{N-1},}$

where $C_{N-1}$ is a correction term describing the contribution of the $N$-body distribution function, whose function is to recover some of the information lost due to the simplification of the model, so as to ensure that the kinetic properties we want to study will not change due to the simplification of the model. Since the evolution equation of the ($N$−1)-body distribution function is still unmanageable, we simplify it to the ($N$−2)-body distribution function evolution equation. Finally, the following 2-body distribution function evolution equation is obtained:

(3)${\displaystyle \frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }t}}+{\displaystyle \sum _{i=1}^2\left({\dot{\mathit{q}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{q}}_i}}+{\dot{\mathit{p}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{p}}_i}}\right)f_2=C_2,}$

where the modified term $C_2$ describes the contribution of the 3-body distribution function $F_3$. It is important to note the fact that the probability of 3 particles colliding together is much less than the probability of 2 particles colliding together at the same time. Thus, the effect of the 3-body interaction is much smaller than the effect of the 2-body interaction, we can approximate $F_3$ as zero. Equation (3) simplifies to

(4)${\displaystyle \frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }t}}+{\displaystyle \sum _{i=1}^2\left({\dot{\mathit{q}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{q}}_i}}+{\dot{\mathit{p}}}_i\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{p}}_i}}\right)f_2=0.}$

Let us make two further assumptions: (i) The interaction between particles depends only on the distance between them, not on the direction:

(5)$U_{12}\left({\mathit{q}}_1,{\mathit{q}}_2\right)=U_{12}\left(\left|{\mathit{q}}_1-{\mathit{q}}_2\right|\right).$

(ii) The 2-body distribution function can be written in the form of the product of two monomer distribution functions, i.e,

(6)$F_2({\mathit{q}}_1,{\mathit{p}}_1;{\mathit{q}}_2,{\mathit{p}}_2;t)=F_1({\mathit{q}}_1,{\mathit{p}}_1;t)F_1({\mathit{q}}_2,{\mathit{p}}_2;t).$

Then, we finally obtain the evolution equation of the monomer distribution function $F_1\equiv f$. That is the Boltzmann equation:

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

where

(8)$Q\left(f,f\right)={\displaystyle \underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \underset{0}{\overset{4\pi }{\int }}(f^{{}^{\prime }}{f}_1^{{}^{\prime }}-f{f}_1)g\sigma \mathrm{d}\mathrm{\Omega }\mathrm{d}{\mathit{u}}_1}}$

describes the rate of change of the distribution function caused by molecular collisions.

From Liouville equation to Boltzmann equation, the coarse-grained physical modeling has been experienced ($N$−1) times. It can be seen that compared with Liouville equation, Boltzmann equation is already a highly coarse-grained physical model. As the foundation of non-equilibrium statistical physics, the Liouville equation is in principle equivalent to molecular dynamics without any simplification, and is a “holographic” description of $N$-body systems, applicable to both solid and fluid systems. The Boltzmann equation, however, is limited to fluid systems.

In statistical physics, there are two commonly used phase spaces: $\mu$ space and $\mathrm{\Gamma }$ space. For a system composed of $N$ identical particles without internal degrees of freedom, if a 6-dimensional phase space is formed based on generalized coordinates and generalized momenta (or velocities) describing the microscopic states of particles, then a state of the system is a distribution of $N$ points in the space. Fig.4 shows the schematic of three common phase spaces, where position or velocity space shown in Fig.4(a) can be regarded as a sub-space of $\mu$ space shown in Fig.4(b), and $\mu$ space shown in Fig.4(b) can be regarded as a sub-space of $\mathrm{\Gamma }$ space shown in Fig.4(c). Fig.4(b) shows two states of the $N$ particle system. If the generalized coordinates and generalized momenta of these $N$ particles are used as the basis to open phase space, i.e., the $\mathrm{\Gamma }$ space, then a state of the system appears as a single point in the space. Fig.4(c) shows two states of the $N$ particle system. Obviously, the usual coordinate space and velocity space are subspaces of two kinds of phase spaces. In coordinate (velocity) space, it is impossible to know the velocity (coordinate) information of a particle by a distribution alone, as shown in Fig.4(a). In the coordinate space, the particle velocity information needs to be roughly obtained by comparing the position information of the two moments before and after. In velocity space, the position information of a particle needs the position information of adjacent two moments and the time step to be roughly known.

Fig.4 Schematic of phase space description method.

Full size|PPT slide

The process from Boltzmann equation to hydrodynamic equations is also a physical coarse-graining modeling process that grasps the main contradiction according to the research demand. In macroscopic fluid theory, fluid states are described using density $\rho$, temperature $T$, velocity $\mathit{u}$, and pressure $p$. The fundamental governing equations are the Navier−Stokes equations:

(9)$\{\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot (\rho \mathit{u})=0,\\ {\displaystyle \frac{\mathrm{\partial }(\rho \mathit{u})}{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot (\rho \mathit{u}\mathit{u}\text{+}p\mathit{I})=-\mathrm{\nabla }\cdot \mathbf{\Pi },\\ {\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot [(E+p)\mathit{u}]=-\mathrm{\nabla }\cdot (\mathit{q}+\mathbf{\Pi }\cdot \mathit{u}),\end{array}$

where $\mathit{I}$ is the unit tensor. In kinetic theory, the state of a system is described using the distribution function $f(\mathit{x},\mathit{v},t)$ and its kinetic moments. Among them, the three conserved moments of the distribution function $f$, (density $\rho$, momentum $\rho \mathit{u}$, and energy $E$) and the two non-conserved moments (viscous stress $\mathbf{\Pi }$ and heat flux $\mathit{q}$) are described into traditional fluid theory. The relationship between these moments and the distribution function $f$ is as follows:

(10)$\mathit{W}={\displaystyle \int f\mathbf{\Psi }\mathrm{d}\mathit{v},}$

where $\mathit{v}$ is the molecular velocity, and

(11)$\mathit{W}={\left[\rho ,\rho \mathit{u},E=(D+n)\rho T/2+\rho \mathit{u}\cdot \mathit{u}/2\right]}^{\mathrm{T}},$

(12)$\mathbf{\Psi }={\left[1,\mathit{v},\mathit{v}\cdot \mathit{v}/2\right]}^{\mathrm{T}},$

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

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

where $f^{eq}$ is the corresponding equilibrium distribution function. Here, $D$ is the spatial dimension and $n$ represents the extra degree of freedom. The higher order non-conserved moments of the distribution function $f$ are not included in the description of traditional fluid theory, which is naturally a double-edged sword: on the one hand, it brings the simplicity of traditional fluid theory, but at the same time, it also brings constraints to describe the more discrete and non-equilibrium cases. In a more reasonable kinetic description, the emergence of more high-order moments of $f$ is an inevitable requirement driven by the requirement/original intention of “the complexity of system behavior increases dramatically + the physical description/control ability does not decrease”. With the increase of the degree of discreteness and the degree of thermodynamic non-equilibrium, the complexity of the system behavior increases sharply, and more physical quantities are needed to describe the state and behavior. A lot of research can be done without increasing the number of physical quantities, but the obvious consequence is a sharp decline in the actual control of system behavior [61]. This point seems to have not received sufficient attention!

It should be noted here that, like NS, Burnett and other hydrodynamic equations, DBM describes a physical image, providing a series of physical constraints on the model for studying physical problems, naturally expressed by equations and relations, but without specific discrete formats.

To establish the connection between the macroscopic and microscopic scales, we must try to cross the threshold of mesoscale. The complexity of “mesoscale” behavior is far greater than expected by reductionist approach. When dealing with a complex system, the general cognitive process is to start with a general overview of its macroscopic (system) scale behavior, and then gradually explore the mechanism of its microscopic (unit) scale, and gradually establish the relationship between system scale behavior and unit scale behavior. However, directly establishing such a connection is often extremely challenging. The research shows that there may be a general dominant principle in the mesoscale between the unit scale and the system scale, that is, the coordination of different control mechanisms in competition. This is highly consistent with the idea that in addition to reductionism paradigm, emergentism paradigm is also needed in physics. Under the initiative of scientists like Academician Jinghai Li and many others, the field of meso science, which focuses on multi-scale and cross-scale research, is attracting more interest with time [62, 63].

3 Discrete Boltzmann method

The nature of an entity comprises both specific individual characteristic and commonality. As the individual characteristic gradually diminishes, what remains is more common. The further one moves away from the original specific system, the stronger the universality becomes. For a model to fully recover all the properties of the system, it must be as complex as the actual system. Perfect in theory, but too complicated to handle, it becomes impractical. Therefore, it is necessary to “choose something to do and ignore others” and simplify the model to the point where it can be handled. When it comes to mechanistic research, further simplification is required to the extent that the mechanism can be analyzed clearly, which is exactly what some mechanism studies do at the beginning. Research based on simple models provides cognitive basis and technical basis for subsequent research on more complex or realistic situations. Knowledge from practice ultimately returns to practice. The research of complex systems requires multi-level and multi-perspective methods and cognition.

Statistical physics is a bridge between micro and macro scales, and the kinetic theory is a bridge between micro and macro descriptions of fluid systems. Coarse-grained modeling is a basic means to grasp the main contradiction according to the research needs of non-equilibrium statistical physics. The related study on Boltzmann equation is a relatively mature part in the theory of kinetics. In the context of fluid dynamics, “cross-scale” typically refers to scaling across Knudsen numbers.

Complex flows are typically associated with unsteady behavior, but the study of steady-state behavior is also an important aspect. From a descriptive perspective, steady-state conditions can be achieved by setting the time derivative terms to zero in the governing equations. In contrast, the usual DBM is designed for unsteady flow [61], but steady-state DBM and its application are also important contents of DBM research.

For the sake of clarity, we refer to DBM with a nonzero time span as time-dependent DBM and DBM that only focuses on the behavior at a single moment as time-independent DBM. The simplest case of time-independent DBM is steady-state DBM [64]. In this review, we will first introduce time-dependent DBM designed for typical unsteady situations and then provide a brief overview of steady-state DBM. The physical behavior description functions of the two are complementary. But the design of the latter is significantly different from that of the former.

3.1 Brief review of non-steady DBM

The Discrete Boltzmann Method (DBM) is proposed to (partially) solve a series of fundamental problems, especially those presented in Section 2. It is a kind of (coarse-grained) physical model construction method and complex physical field analysis method developed based on the discrete Boltzmann equation,

(15)${\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }t}}+{\mathit{v}}_i\cdot {\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }\mathit{x}}}+(\text{force}\text{term}{)}_i=(\text{collision}\text{term}{)}_i,$

where the subscript “i” of the distribution function $f$ is the index of the discrete velocity, corresponding to the discrete velocity ${\mathit{v}}_i$. Obviously, the discrete velocity selection rule is the key technology in DBM modeling. Fig.5 shows the schematic of several commonly used discrete velocity sets in current literature.^†

¹⁾ Please note that this is by no means a standard or optimal choice that can be applied to every situation. The selection of the optimal discrete velocities depends not only on the discrete formats of time integral and spatial derivative, but also on the specific fluid behavior. This is still an open topic in computational mathematics and is beyond the scope of this article.

Fig.5 Schematic of several commonly used two-dimensional discrete velocity sets.

Full size|PPT slide

In terms of physical model construction, DBM includes two steps of coarse-grained physical modeling: modification and simplification of Boltzmann equation and discretization of particle velocity space. The basic principle of coarse-grained physical modeling is that the behaviors of the system to be studied cannot be changed due to the simplification of the model. In kinetic theory, in addition to the distribution function itself, system properties are described by the kinetic moments of the distribution function. So, both before and after simplifying the collision term and discretizing the particle velocity space, the control equation should provide the same moments of the distribution function that correspond to the desired kinetic properties. This is represented mathematically as

(16)${\displaystyle \int f\mathbf{\Phi }\left(\mathit{v}\right)\mathrm{d}\mathit{v}={\displaystyle \sum _i{f}_i\mathbf{\Phi }\left({\mathit{v}}_i\right).}}$

Due to the lack of analytical solutions for the kinetic moments of the distribution function $f$, except for the three conserved moments, it is not convenient to obtain the physical constraints required for discrete velocity selection using Eq. (16). So, we turned to examine all the kinetic moments of the equilibrium distribution function $f=f^{eq}$ involved in the calculation of the left side of Eq. (16). Furthermore, it is required that the results of these kinetic moments remain unchanged when converted to summation for calculation. Mathematically, this can be expressed as

(17)${\displaystyle \int f^{(0)}\mathbf{\Psi }\left(\mathit{v}\right)\mathrm{d}\mathit{v}={\displaystyle \sum _i{f}_i^{(0)}\mathbf{\Psi }\left({\mathit{v}}_i\right).}}$

Please note that moving from Eq. (16) to Eq. (17) provides sufficient conditions for the required physical constraints. Both $\mathbf{\Phi }$ and $\mathbf{\Psi }$ are column vectors, whose elements correspond to the kinetic moments that should keep values. The elements of the former $\mathbf{\Phi }$ are partial elements of the latter $\mathbf{\Psi }$. The latter $\mathbf{\Psi }$ contain elements with higher powers of particle velocity. As the degree of discreteness/non-equilibrium of the system increases, elements with higher particle velocity powers are gradually added to $\mathbf{\Psi }$.

There are a few points to highlight here. (i) Many behaviors of complex flows exceed the descriptive power of the original Boltzmann equation. The Boltzmann equation in DBM is actually a modified Boltzmann equation based on specific situations combined with mean field theory, and its applicability may, in certain aspects, exceed that of the original Boltzmann equation, as shown in Fig.6 [61, 65]. So, DBM is actually a coarse-grained modeling and analysis method that combines kinetic theory with mean field theory^†

²⁾ Physics study has two research paradigms, reductionism and emergentism. The latter is used to describe the part that cannot be included in the former. Statistical physics is a one of the disciplines that do not rely on reductionism. However, in the study of fluid physics, it is often necessary to use kinetic theory to investigate the mesoscopic kinetic behavior ignored by the traditional macroscopic modeling, so the role of kinetic theory here is similar to that of reductionism. DBM modeling is originated from a reasonable combination of kinetic theory and mean field theory, where the latter is used to describe the part that the former cannot describe and consequently its role here is similar to that of emergentism.

. (ii) When dealing with the external force term of the Boltzmann equation, in order to avoid the difficulty of making differentiation to the discrete velocity, the DBM modeling approach is to first approximate the external force term, use the property that the equilibrium distribution function is derivable to the particle velocity. Once the differentiation of particle velocities is accomplished, the expression is then rewritten in terms of discrete velocities. (iii) The model equations of DBM include the system behavior evolution equation(s) and discrete velocity constraint equations, where the evolution equations, in addition to discrete Boltzmann equation, depending on the specific situation, may also include phase field evolution equation, chemical reaction evolution equation, electromagnetic field evolution equation, etc. (iv) DBM only gives the physical constraints that need to be followed when choosing discrete velocities, but does not contain the specific discretization format. In order to reduce the influence of the dramatic change of flow field on the discrete velocity distribution function, the discrete velocity constraint equations can be written in the form of central moments and used to select the discrete velocities on the fluid element, i.e., the discrete peculiar velocities ${\mathit{c}}_i={\mathit{v}}_i-\mathit{u}$.

Fig.6 Schematic diagram of the application scopes of the Boltzmann equation, the original BGK model, and the BGK-like model in kinetic methods.

Full size|PPT slide

In terms of complex physical field analysis method, DBM is a specific application and further development of statistical physics coarse-grained description method, non-equilibrium behavior description method, and phase space description method within the framework of discrete Boltzmann equation. From a historical perspective, DBM has evolved from the physical modeling branch of the Lattice Boltzmann Method (LBM), with some abandoning and new addition. It preserves the use of discrete velocity, but is no longer limited to specific discrete format. It provides only the most necessary physical constraints for the discrete format. It is no longer based on the continuity assumption and near equilibrium approximation of traditional fluid modeling, and no longer uses the “lattice gas” image of standard LBM. It adds new detection, presentation, description, and analysis schemes for non-equilibrium states and resulting effects based on phase space, and introduces more information extraction techniques and complex physical field analysis techniques over time. Fig.7 depicts the expansion from the phase-space description of the non-conserved kinetic moments

Fig.7 Phase space description methods: from non-conserved moments of $f-f^{(0)}$ to any set of system characteristics.

Full size|PPT slide

(18)$\begin{array}{rl}{\mathbf{\Delta }}_n^{\ast }=& {\mathit{M}}_n^{\ast }\left(f-f^{(0)}\right)\\ =& \int \mathrm{d}\mathit{v}\left(f-f^{(0)}\right)\underset{n\text{-th order tensor}}{\underbrace{(\mathit{v}-\mathit{u})(\mathit{v}-\mathit{u})\cdots (\mathit{v}-\mathit{u})}}\\ =& \sum _i\left(f_i-f_i^{(0)}\right)\underset{n\text{-th order tensor}}{\underbrace{({\mathit{v}}_i-\mathit{u})({\mathit{v}}_i-\mathit{u})\cdots ({\mathit{v}}_i-\mathit{u})}}\end{array}$

to phase space description method based on any set of behavioral features,

(19)$\mathit{X}=\left\{X_1,X_2,X_3,\cdots \right\}.$

In Eq. (18), $\mathit{u}$ is the local flow velocity. In the phase space or its subspace, we can define the non-equilibrium strength (or the strength of corresponding behavioral features) of the corresponding perspective by the distance from the state point to the coordinate origin, such as $D$ in Fig.7. Using the concept of distance between two points to describe the differences between two non-equilibrium states (or corresponding behavioral features), such as $d$ in Fig.7. Therefore, a more complete name for DBM is the Discrete Boltzmann modeling and analysis method [61]. The use of non-conserved kinetic moments of $f-f^{(0)}$ to detect and describe the specific way in which a system deviates from its thermodynamic equilibrium state and the various effects caused by this are key techniques for DBM in complex physical field analysis.

The behavior of actual systems is often complex. Coarse-grained modeling is a process of losing information, but this losing information has a bottom line. The bottom line is that the properties under investigation cannot be changed due to the simplification of the model. DBM, based on the research needs of the problem, selects a perspective to study a set of kinetic properties of the system, and requires the kinetic moments describing these properties to be preserved during the model simplification process. As the degree of discreteness and non-equilibrium increases, the complexity of system behavior increases sharply. Then incorporating some higher-order kinetic moments into the descriptive perspective becomes an inevitable outcome to ensure that the control ability does not decrease. Therefore, using more physical quantities to describe the system state and behavior is a typical feature that distinguishes DBM from traditional fluid modeling and other current kinetic methods. From the perspective of Kinetic Macro Modeling (KMM),^†

³⁾ Corresponding to the physical function of DBM is the Extended Hydrodynamic Equations (EHE), which include not only the conservation moment evolution equations corresponding to the three conservation laws of mass, momentum, and energy, but also some of the most closely related non-conservation moment evolution equations. For convenience of description, we refer to the modeling method of deriving EHE based on kinetic equations as Kinetic Macro Modeling (KMM).

this is a requirement for obtaining more accurate constitutive relationships. From the perspective of the kinetic theory, this is a requirement for obtaining more accurate distribution function $f$ under given conditions. Different perspectives lead to the same goal.

The starting point of DBM lies within the mesoscale range, approaching the macroscopic side. Due to the stage nature of development, the primary situations currently considered by DBM still fall within the realm where the Chapman−Enskog (CE) multiscale analysis theory is valid. Therefore, the CE multi-scale analysis theory is a mathematical guarantee for the reasonable and effective DBM approach [61].

3.2 Rules for selecting discrete velocities

By simultaneously taking the density moment, momentum moment, and energy moment on both sides of the Boltzmann equation, a set of fluid dynamic equations that are consistent with the NS form can be obtained,

(20)$\{\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot (\rho \mathit{u})=0,\\ {\displaystyle \frac{\mathrm{\partial }(\rho \mathit{u})}{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot (\rho \mathit{u}\mathit{u}\text{+}p\mathit{I})=-\mathrm{\nabla }\cdot {\mathbf{\Delta }}_2^{\ast },\\ {\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot [(E+p)\mathit{u}]=-\mathrm{\nabla }\cdot ({\mathbf{\Delta }}_{3,1}^{\ast }+{\mathbf{\Delta }}_2^{\ast }\cdot \mathit{u}).\end{array}$

For convenience of description, here we refer to this as the generalized NS. The advantage of this generalized NS is that there is no approximation. Its viscous stress

(21)$\begin{array}{rl}{\mathbf{\Delta }}_2^{\ast }=& {\mathit{M}}_2^{\ast }\left(f-f^{(0)}\right)={\mathit{M}}_2^{\ast }\left(f^{(1)}+f^{(2)}+f^{(3)}+\cdot \cdot \cdot \right)\\ =& {\mathbf{\Delta }}_2^{\ast (1)}+{\mathbf{\Delta }}_2^{\ast (2)}+{\mathbf{\Delta }}_2^{\ast (3)}+\cdot \cdot \cdot \end{array}$

and heat flux

(22)$\begin{array}{rl}{\mathbf{\Delta }}_{3,1}^{\ast }=& {\mathit{M}}_{3,1}^{\ast }\left(f-f^{(0)}\right)={\mathit{M}}_{3,1}^{\ast }\left(f^{(1)}+f^{(2)}+f^{(3)}+\cdot \cdot \cdot \right)\\ =& {\mathbf{\Delta }}_{3,1}^{\ast (1)}+{\mathbf{\Delta }}_{3,1}^{\ast (2)}+{\mathbf{\Delta }}_{3,1}^{\ast (3)}+\cdot \cdot \cdot \end{array}$

include contributions from various orders of non-equilibrium distribution functions $f^{(j)}$, where $j=1,2,3,\cdot \cdot \cdot$, corresponding to the discrete/non-equilibrium order described by the power of Kn number. The subscript “3,1” signifies that a 3rd-order tensor has undergone 1 contraction operation to become a 1st-order tensor, which is a vector. The meanings of other subscripts below are similar, and so on. However, the limitation of this approach is that it cannot provide specific expressions for stress and heat flux directly. To derive analytical expressions for constitutive relationships, Chapman and Enskog developed the CE multi-scale analysis method later named after them. CE multiscale analysis is actually a generalized Taylor expansion and analysis method: the independent variable here is the Kn number. Not only the distribution function $f$, but also the temporal and spatial derivatives are Taylor expanded at the point, $\text{Kn}=0$ (i.e., continuity, thermodynamic equilibrium). More physical image explanations for CE multiscale analysis can be found in Ref. [61]. The CE multi-scale analysis informs us the following:

(i) The Boltzmann equation, under the quasi-continuous and near-equilibrium conditions, corresponds to the system of fluid dynamics equations that are the NS equations. The viscous stress and heat flux considering only the first order non-equilibrium $f^{(1)}$ are the NS stress $\mathbf{\Pi }$ and NS heat flux $\mathit{q}$ in Eqs. (13) and (14), respectively. As the degree of discretization/non-equilibrium increases, the contribution of second-order $f^{(2)}$ or even higher-order non-equilibrium $f^{(j)}(j>2)$ should be considered in viscous stress and heat flow.

(ii) Among all the kinetic moments of the distribution function $f$, only the three conserved moments (density, momentum, and energy) and two non-conserved moments (viscous stress and heat flux) enter the NS description. The remaining non-conserved moments have not entered the NS theory, which is a double-edged sword: on the one hand, they bring the simplicity of traditional fluid mechanics theory; however, on the other hand, it creates obstacles for NS to describe situations with higher degrees of discreteness/non-equilibrium, which is the physical reason why NS cannot describe well situations with higher degrees of discreteness/non-equilibrium.

(iii) The ability to recover corresponding levels of macroscopic hydrodynamic equations (such as NS and Burnett equations) is only part of the physical function of DBM; Corresponding to the physical functions of DBM is the Extended Hydrodynamic Equations (EHE), which include not only the three common hydrodynamic equations but also some of the most closely related non-conservative moment evolution equations. The extended part (i.e., the evolution equations for relevant non-conserved moments) becomes increasingly necessary as the degree of discreteness/non-equilibrium increases.

(iv) As the degree of discreteness/non-equilibrium characterized by Kn number increases, the complexity of DBM simulation increases more slowly compared to KMM simulation (i.e., deriving and solving EHE), thus allowing it to go further. If the order of the Kn number to be considered increases by 1, then the number of kinetic moments in DBM that need to be preserved increases by 2, while the complexity of deriving EHE from KMM sharply increases. Furthermore, in the case of multiple media, due to the existence of different media flow velocities and average flow velocity, and the definition of temperature depends on the referenced flow velocity, the form of EHE is not unique. These different forms of EHE correspond to descriptions from different perspectives of complex flow dynamics in multiple media. It can be seen that in the case of multiple media, the correspondence between DBM and physically equivalent EHE is one-to-several.

DBM simulation does not require knowing the specific form of extremely complex EHE with physical equivalence, but we can borrow KMM’s ideas and physical images to quickly screen out the kinetic moments involved in the more accurate stress and heat flux calculations (without strictly deriving equations, only checking the order of the kinetic moments). This part of the kinetic moments are converted from integral to summation for calculation, and the results need to remain unchanged. At this point, we have obtained the most necessary physical constraints to follow for the selection of discrete velocity: the linear equations of $f_i^{(0)}$ represented by Eq. (17).

3.3 Description of non-equilibrium behavior and effect

Non-equilibrium strength/degree/extent/intensity is one of the most fundamental parameters for describing and recognizing non-equilibrium flows. In addition to commonly used Kn numbers, spatial gradients of macroscopic physical quantities such as density, flow velocity, temperature, and pressure, DBM uses physical quantities based on non-conserved kinetic moments of $f-f^{(0)}$ to describe the way and magnitude of system deviation from equilibrium. Furthermore, DBM uses the (independent components of) non-conserved moments of $f-f^{(0)}$ [e.g., $\{{\mathbf{\Delta }}_n^{\ast },n=2,(3,1),3,(4,2),(5,3),\cdot \cdot \cdot \}$] as the bases to open phase space, providing an intuitive geometric correspondence for complex system states and behavioral characteristics, as shown in Fig.7. It is easy to understand that each non-conserved moment itself and any independent component of $f-f^{(0)}$ describe the way and magnitude of the system deviating from equilibrium from its own perspective, thus offering a unique view of non-equilibrium features. These non-equilibrium strengths from different perspectives are interrelated, complementary, and cannot be replaced by each other. Together, they form a more complete description. It should be emphasized that as a later and more fundamental description method, DBM naturally inherits all traditional methods of describing non-equilibrium behavior.

The definition of any non-equilibrium strength depends on the research perspective. The non-equilibrium strength may be increasing from a certain perspective; while from another perspective, it may be decreasing. This is one of the concrete manifestations of the complexity of non-equilibrium behavior. Studying non-equilibrium strength and effects from only one perspective is often one-sided or even incorrect. Being incorrect refers to the misconception that the conclusions obtained are independent of the research perspective and are universal. In view of this, DBM introduces a multi-perspective non-equilibrium description scheme based on non-conserved moments of $f-f^{(0)}$, gradients of macroscopic quantities, thermodynamic relaxation time $\tau$, Kn number, morphological description, etc., to provide a cross-positioning of non-equilibrium intensity of complex flows from multiple perspectives. For the convenience of description, the concept of non-equilibrium strength vector is further introduced, where each component of the vector represents a non-equilibrium intensity from a different perspective. For example, the composition of a non-equilibrium strength vector may be

(23)$\mathit{S}=\left\{D_2,D_3,|\mathbf{\nabla }\rho |,|\mathbf{\nabla }T|,\tau ,\text{Kn},\cdots \right\},$

where

(24)$D_2=\sqrt{{\mathrm{\Delta }}_{2,xx}^{\ast 2}+{\mathrm{\Delta }}_{2,xy}^{\ast 2}+{\mathrm{\Delta }}_{2,yy}^{\ast 2}}$

is the distance in the ${\mathbf{\Delta }}_2^{\ast }$ sub-space of non-conservative moment phase space, and is the non-equilibrium strength from the perspective of ${\mathbf{\Delta }}_2^{\ast }$. $T$ is the local temperature. In some cases, $D_2$ is defined as

(25)$D_2=\sqrt{{\mathrm{\Delta }}_{2,xx}^{\ast 2}+{\mathrm{\Delta }}_{2,xy}^{\ast 2}+{\mathrm{\Delta }}_{2,yx}^{\ast 2}+{\mathrm{\Delta }}_{2,yy}^{\ast 2}}.$

Similar but different options exist for defining non-equilibrium strengths such as $D_3$. These different choices correspond to different research perspectives.

Just as the description of phase space does not require its coordinates with the same unit, the various components of the non-equilibrium intensity vector also do not need to have the same unit. The non-equilibrium strength from different perspectives, as a set of behavioral characteristics, can also be described using the phase space method, providing an intuitive geometric correspondence for the non-equilibrium strength of complex flows from different perspectives. Of course, as a coarse-grained modeling and description method, the physical accuracy of DBM, including the precision of describing behavioral characteristics such as non-equilibrium intensities, should be adjusted based on specific research requirements [61].

3.4 Additional remarks on DBM

At a report on the 2nd U.S.-Japan Joint Seminar, held at Santa Barbara, California, from 7−10 August, 1996, the late eminent scholar, Mr. Chang-Lin Tien, et al. commented the LBM as follows: “Many physical phenomena and engineering problems may have their origins at molecular scales, although they need to interface with the macroscopic or ‘human scales’. The difficulty arises in bridging the results of these models across the span of length and time scales. The lattice Boltzmann method attempts to bridge this gap”. More details are referred to Ref. [66]. Many people’s understanding on LBM is directly or indirectly affected by this or similar comments.

In fact, starting from its predecessor, the lattice gas (cellular automata) method, there are two types of LBMs in the literature: coarse-grained physical model construction method [67-69] and numerical solving method for certain control equations [70-72]. The latter accounts for the vast majority of existing literature, to the extent that LBM has almost become the abbreviation or synonym for the latter. However, the LBM method mentioned by Mr. Chang-Lin Tien in the above review was clearly a cross scale physical model construction method, which falls into the former category. Coarse-grained physical modeling is a fundamental means for non-equilibrium statistical physics to address the main contradiction based on research needs. Cross scale, here referred to as cross Kn number, is a physical function inherited from the Boltzmann equation description. These two different types of LBM work on complementary dimensions, with different goals and construction rules, each of which is reasonable.

DBM is often referred to as the Discrete Boltzmann Method, Discrete Boltzmann Modeling, and Discrete Boltzmann Model without causing ambiguity. As mentioned earlier, the more complete term for the Discrete Boltzmann Method is Discrete Boltzmann modeling and analysis Method. When compared with other methods, DBM’s identity is “physical model construction method + complex physical field analysis method”, or rather a physical model construction method with built-in complex physical field analysis function. When compared with other physical models, DBM is a physical model with built-in complex physical field analysis capabilities. Because the commonly used physical models currently do not provide complex physical field analysis capabilities for simulated data, when referring to the term “discrete Boltzmann model or modeling”, people often do not realize that it also includes complex physical field analysis function. This is why our research group has to refer to it as the “physical model construction and complex physical field analysis method” in some situations.

The physical model construction part of DBM includes determining specific control equations and providing the most necessary physical constraints to be followed for selecting discrete velocities. Just like in NS simulation, in DBM simulation, the specific discretization formats for spatial derivatives, time integrals, and particle velocities should be chosen reasonably based on the specific operating conditions and flow patterns.

According to Eq. (18), from a dimensional perspective and physical image, the physical meaning of non-equilibrium feature quantities ${\mathbf{\Delta }}_{n+1}^{\ast }$ is the non-organized flux of ${\mathbf{\Delta }}_n^{\ast }$. It should be noted that in macroscopic description, a physical quantity (such as density $\rho$) multiplied by flow velocity $\mathit{u}$, its physical meaning is the flux of this physical quantity, which is an organized flux. The combination of organized flux and non-organized flux constitutes a more comprehensive description of complex flow kinetic behavior. Among all the non-organized fluxes composed of non-conserved central moments of $f-f^{(0)}$, only Non-Organized Momentum Flux (NOMF) ${\mathbf{\Delta }}_2^{\ast }$ and Non-Organized Energy Flux (NOEF) ${\mathbf{\Delta }}_{3,1}^{\ast }$ have counterparts in the NS description, corresponding to viscous stress and heat flux, respectively. The remaining non-conserved central moments are not directly included in NS, but as mentioned earlier, their physical images are clear. The absence of these physical quantities in the NS description is precisely the physical reason of NS’s inability to describe well situations with higher degrees of discreteness/non-equilibrium. As the degree of discreteness/non-equilibrium increases, these non-organized fluxes should receive more attention.

Since the traditional fluid mechanics theory is based on NS, we are most familiar with the NS equation. Comparing with NS is an important way to understand new models. Therefore, it is necessary to clarify the connection and difference between DBM and NS again: NS is just a physical model, and it does not inherently possess complex physical field analysis capabilities for data, it only considers first-order (in terms of the Kn number perspective) discrete/non-equilibrium effects, which is only applicable to quasi-continuous and near-equilibrium situations. In terms of numerical simulation research, NS is only responsible for pre-simulation and not for post-simulation. In contrast, DBM not only provides modeling capabilities equivalent to the corresponding level of EHE but also provides a set of complex physical field analysis methods for data. In terms of numerical simulation research, DBM is not only responsible for pre-simulation, but also for post-simulation.

DBM is beyond traditional fluid modeling in both its applicable range and physical capabilities, focusing more on the “mesoscale” dilemma where macroscopic continuous models lack physical functionality and microscopic molecular dynamics simulations are powerless due to limited scale applicability, as shown in Fig.8. In the schematic, RANS is the abbreviation for Reynolds Average Navier Stokes, which is a further coarse-grained physical modeling based on NS and describes large-scale behavior within the NS range. The “NS−RANS” in the schematic represents the portion that RANS fails to well describe within the scope of NS application, where “−” is a minus sign, indicating subtraction from the scope of application. The “DBM−NS” indicates the part within the scope of DBM that is not well described by NS and is the key focus of DBM.

Fig.8 Schematic for application ranges of several physical models.

Full size|PPT slide

The greater challenge in the study of mesoscale behavior lies in the modeling, simulation, and analysis of higher degrees of discreteness/non-equilibrium. The necessity of the most closely related non-conserved moments increases as the degree of discreteness/non-equilibrium increases. Although the goal is for situations with higher degrees of discreteness/non-equilibrium, almost everyone has to go through quasi continuous and near equilibrium situations because he or she requires cognitive and technical foundations. This is why progress in the field of mesoscale behavior research can be relatively slow, especially for those just entering the field.

Compared with macroscopic behavior, the typical characteristics of “mesoscale” behavior are significant discrete effects and significant thermodynamic non-equilibrium effects. The description of mesoscale behavior requires more physical quantities. The “mesoscale” characteristics of DBM include: (i) in terms of scale, between micro and macro, connecting micro and macro; (ii) in terms of functionality, surpassing NS in describing discrete/non-equilibrium effects, but weaker than MD. In DBM, the corresponding solutions to the six basic scientific problems proposed in Section 2 are as follows: For problem (i), add physical quantities based on the definition of non-conserved kinetic moments according to the degree of discreteness/non-equilibrium. In response to problem (ii), phase space description methods based on non-conserved moments of $f-f^{(0)}$, morphological and other behavioral feature quantities are introduced. In response to problem (iii), the model construction no longer relies on continuity and near equilibrium assumptions. To address problem (iv), develop a description method that starts from macro continuous near equilibrium and gradually increases the degree of discreteness and non-equilibrium. For problem (v), develop mesoscale modeling based on the combination of kinetic theory and mean field theory. For problem (vi), introduce the concept of non-equilibrium strength vector, where each component represents a non-equilibrium intensity, for a cross localization description. Tab.1 summarizes the developmental milestones in the evolution of DBM.

Tab.1 Milestones in the development of DBM.

Year	Schemes for detecting and describing discrete/TNE behavior and effects
Before 2012	There was no significant difference in physical function between the two types of LBMs.
2012	Proposed to use non-conserved moments of $f-f^{(0)}$ to detect and describe discrete/TNE states and effects, which is the starting point of current DBM method [73].
2015	Proposed to use the non-conserved moments of $f-f^{(0)}$ as bases to open phase space, and use the distance from a state point to the origin to define the TNE strength of one perspective. This is the starting point of the phase space description method in DBM [74].
2018	Proposed to use the distance between two points in the phase space to describe the difference between two discrete/TNE states, and use the mean distance between two points in a given time interval to describe the difference of two kinetic processes [75].
2021	Extended the phase space method to describe any set of system features [76].
2022	Proposed the concept of non-equilibrium strength vector, each of whose components is one TNE strength of a perspective, to multi-perspective cross-locate the non-equilibrium strength of complex flow [61, 77].

3.5 Modeling examples and flowchart

According to the fundamental theory of non-equilibrium statistical physics, there are infinite ways in which the distribution function $f$ deviates from the equilibrium distribution function $f^{(0)}$. There are infinitely many non-conserved moments of $f-f^{(0)}$, and there are infinite possibilities for which moments have larger amplitudes (or intensities) and which moments have smaller amplitudes. The approach of DBM involves breaking down complex problems, selecting a perspective, studying a set of kinetic properties of the system, and requiring the non-conserved moments that describe these properties to remain invariant during model simplification. Based on the independent components of the set of kinetic moments of $f-f^{(0)}$, construct a phase space and use it and its subspaces to describe the discrete and non-equilibrium behavior characteristics of the system. The research perspective and modeling accuracy adjust as the research progresses. Therefore, in principle, DBM modeling itself can be independent of CE multiscale analysis. As long as it is possible to obtain which kinetic moments of $f^{(0)}$ are necessary for grasping system behavior in a certain way, it is only required that the corresponding kinetic moments remain unchanged when converted into summation for calculation.

Given the staged nature of the research, we will only provide the discussion on DBM modeling based on CE multi-scale analysis theory. By utilizing CE multi-scale analysis and preserving the TNE effect of different orders, the Boltzmann equation can recover, which is actually reduced to, the corresponding order of fluid dynamic equations. In principle, it is not necessary, but most of the published DBM modeling is still within the framework of CE multi-scale analysis theory. Unlike the LBM numerical solving scheme for some given control equations, in the DBM modeling process, the role of CE multiscale analysis is not to help derive the corresponding EHE functions, but only to quickly identify and confirm the kinetic moments of $f^{(0)}$ that need to be preserved based on the tensor order of particle velocity $\mathit{v}$. Whether the specific form of the corresponding EHE is known or not does not affect DBM modeling and simulation. The CE multi-scale analysis theory serves as a reasonable theoretical basis for such kinetic models.

In fact, even in the case of a simple single-medium system, the theoretical derivation of 2nd-order or higher discrete and non-equilibrium effects in the KMM or EHE is already quite complex, let alone in the case of multi-medium scenarios. Even if EHE, which is functionally equivalent to DBM, is ultimately derived, the feasibility of numerical simulation remains a huge challenge. Because it involves stronger nonlinearity and higher-order spatiotemporal partial derivatives, and its number of terms increases sharply with the increase of discreteness and non-equilibrium degree, these pose substantial challenges to practical numerical simulation. Therefore, as the degree of discreteness or non-equilibrium increases, the KMM modeling and simulation approach quickly becomes infeasible.

Unlike KMM, DBM belongs to the kinetic direct modeling method. Except for explaining and verifying the functions of the model, functionally equivalent EHE is not required in actual DBM simulations. As the degree of discreteness or non-equilibrium increases, the complexity of DBM modeling and simulation approach also increases, but it is relatively slow, making it possible to go further.

Below, we illustrate the process of DBM modeling for an ideal monatomic gas in a single medium, showing the process of DBM modeling with different discrete or non-equilibrium degrees. The general form of CE multi-scale expansion is

(26)$\begin{array}{rl}& f=f^{(0)}+\text{Kn}{f}^{(1)}+{\text{Kn}}^2{f}^{(2)}+\cdot \cdot \cdot (\mathrm{a}),\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }t}=\text{Kn}\frac{\mathrm{\partial }}{\mathrm{\partial }{t}_1}+{\text{Kn}}^2\frac{\mathrm{\partial }}{\mathrm{\partial }{t}_2}+\cdot \cdot \cdot (\mathrm{b}),\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{r}}=\text{Kn}\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_1}+\cdot \cdot \cdot (\mathrm{c}).\end{array}$

To facilitate the interpretation of the compositions of distribution function, temporal rate of change, and spatial rate of change, we absorb the coefficients $K{n}^j(j=1,2,3,\dots )$ into the corresponding non-equilibrium distribution function $f^{(j)}$, temporal rate of change $\mathrm{\partial }/\mathrm{\partial }{t}_j$ and spatial rate of change $\mathrm{\partial }/\mathrm{\partial }{\mathit{r}}_j$. Thus, the CE multi-scale expansion becomes

(27)$\begin{array}{rl}& f=f^{(0)}+f^{(1)}+f^{(2)}+\cdot \cdot \cdot (\mathrm{a}),\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }{t}_1}+\frac{\mathrm{\partial }}{\mathrm{\partial }{t}_2}+\cdot \cdot \cdot (\mathrm{b}),\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{r}}=\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_1}+\cdot \cdot \cdot (\mathrm{c}).\end{array}$

In this way, the meaning of $f^{(j+1)}$ is the newly added $(j+1)$th order more detailed description on basis of the previous $j$th order level coarse-grained description, $f=f^{(0)}+$$f^{(1)}+f^{(2)}+\cdot \cdot \cdot +f^{(j)}$. For the description of temporal change rate, the meaning of $\mathrm{\partial }/\mathrm{\partial }{t}_{j+1}$ is, based on the previous coarse-grained description $\mathrm{\partial }t=\mathrm{\partial }{t}_1+$$\mathrm{\partial }{t}_2+\cdots +\mathrm{\partial }{t}_j$, as the temporal scale decreases from $t_j$ to $t_{j+1}$, the higher-frequency component that we observe. For the description of spatial change rate, the meaning of $\mathrm{\partial }/\mathrm{\partial }{\mathit{r}}_{j+1}$ is the $(j+1)$th order smaller structure that we newly observe on the basis of the previous $j$th order coarse-grained description, $\mathrm{\partial }\mathit{r}=\mathrm{\partial }{\mathit{r}}_1+\mathrm{\partial }{\mathit{r}}_2+\cdots +\mathrm{\partial }{\mathit{r}}_j$, when the spatial scale decreases from ${\mathit{r}}_j$ to ${\mathit{r}}_{(j+1)}$. It should be noted that in the existing CE theory in the literature, the spatial change rate only retains the first order. As mentioned earlier, the corresponding scenario of $f^{(0)}$ is continuous and equilibrium. So, what do $t_0$ and ${\mathit{r}}_0$ correspond to? Reference [61] points out that $t_0$ corresponds to the total time span of the system behavior that we aim to investigate. When we use the time scale $t_0$ to examine the rate of change of system behavior over time, we cannot see the change, that is $\mathrm{\partial }{t}_0=0$. This is the reason why this term does not appear in CE expansion. Similarly, ${\mathit{r}}_0$ corresponds to the total space span of the system behavior that we aim to investigate. When we use the space scale ${\mathit{r}}_0$ to examine the rate of change of system behavior over space, we cannot see the change, that is $\mathrm{\partial }{\mathit{r}}_0=0$. This is the reason why this term does not appear in CE expansion. In the following discussions on CE multiscale analysis, the spatial change rate $\mathrm{\partial }\mathit{r}$ is only discussed at one spatial scale ${\mathit{r}}_1$.

Substituting CE expansion (27) into the BGK-Boltzmann equation,

(28)${\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }t}}+\mathit{v}\cdot {\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\mathit{r}}}=-\frac{1}{\tau }\left(f-f^{(0)}\right),$

gives

(29)$\begin{array}{rl}& \left(\frac{\mathrm{\partial }}{\mathrm{\partial }{t}_1}+\frac{\mathrm{\partial }}{\mathrm{\partial }{t}_2}+\cdots \right)(f^{(0)}+f^{(1)}+f^{(2)}+\cdots )\\ & +\mathit{v}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_1}(f^{(0)}+f^{(1)}+f^{(2)}+\cdots )\\ & =-\frac{1}{\tau }(f^{(1)}+f^{(2)}+\cdots ),\end{array}$

where $\tau$ is the relaxation time. By equating the coefficients of the terms with the same order of Kn number on both sides of Eq. (29), we can obtain

(30)$\begin{array}{rl}{f}^{\text{(1)}}=& \tau \left[-{\displaystyle \frac{\mathrm{\partial }{f}^{(0)}}{\mathrm{\partial }{t}_1}}-{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_1}}(f^{(0)}\mathit{v})\right]\\ =& \tau \left[-{\displaystyle \frac{\mathrm{\partial }{f}^{(0)}}{\mathrm{\partial }{t}_1}}-{\mathbf{\nabla }}_1\cdot (f^{(0)}\mathit{v})\right],\end{array}$

(31)$f^{\text{(2)}}=-\tau \left[{\displaystyle \frac{\mathrm{\partial }{f}^{(0)}}{\mathrm{\partial }{t}_2}}+{\displaystyle \frac{\mathrm{\partial }{f}^{(1)}}{\mathrm{\partial }{t}_1}}+{\mathbf{\nabla }}_1\cdot (f^{(1)}\mathit{v})\right],$

(32)$f^{(3)}=\cdots .$

Assume we want to propose a DBM that only includes first-order discrete and non-equilibrium effects. It is easy to find

(33)$\begin{array}{rl}{\mathbf{\Delta }}_2^{\ast (1)}=& {\mathbf{\Delta }}_2^{(1)}={\mathit{M}}_2(f_{}^{(1)})\\ =& \tau {\mathit{M}}_2\left(-\frac{\mathrm{\partial }{f}_{}^{(0)}}{\mathrm{\partial }{t}_1}-{\mathbf{\nabla }}_1\cdot (f_{}^{(0)}{\mathit{v}}_{})\right)\\ \Rightarrow & {\mathit{M}}_2\left(f^{\left(0\right)}\right),{\mathit{M}}_3\left(f^{\left(0\right)}\right),\end{array}$

(34)$\begin{array}{rl}{\mathbf{\Delta }}_{3,1}^{(1)}=& {\mathit{M}}_{3,1}(f_{}^{(1)})\\ =& \tau {\mathit{M}}_{3,1}\left(-\frac{\mathrm{\partial }{f}_{}^{(0)}}{\mathrm{\partial }{t}_1}-{\mathbf{\nabla }}_1\cdot (f_{}^{(0)}{\mathit{v}}_{})\right)\\ \Rightarrow & {\mathit{M}}_{3,1}\left(f^{\left(0\right)}\right),{\mathit{M}}_{4,2}\left(f^{\left(0\right)}\right),\end{array}$

(35)$\begin{array}{rl}{\mathbf{\Delta }}_{3,1}^{\text{*}(1)}=& {\mathbf{\Delta }}_{3,1}^{(1)}-\mathit{u}\cdot {\mathbf{\Delta }}_2^{\ast (1)}\\ =& \cdots \\ \Rightarrow & {\mathit{M}}_{3,1}\left(f^{\left(0\right)}\right),{\mathit{M}}_{4,2}\left(f^{\left(0\right)}\right).\end{array}$

Among them, the meaning of the symbol “$\Rightarrow$” is “involving”. Therefore, the set of kinetic moments involved in the viscous stress and heat flux calculation formulas are the following seven kinetic moments, $\left\{{\mathit{M}}_0,{\mathit{M}}_1,{\mathit{M}}_{2,0},{\mathit{M}}_2,{\mathit{M}}_{3,1},{\mathit{M}}_3,{\mathit{M}}_{4,2}\right\}$. Consequently, DBM requires that these seven kinetic moments be strictly preserved during the process of converting from integrals to summations:

(36)${\mathit{M}}_0=\stackrel{N}{{\displaystyle \sum }_{i=1}}{f}_i^{eq}=\rho ,$

(37)${\mathit{M}}_1=\stackrel{N}{{\displaystyle \sum }_{i=1}}{f}_i^{eq}{\mathit{v}}_i=\rho \mathit{u},$

(38)${\mathit{M}}_{2,0}=\stackrel{N}{{\displaystyle \sum }_{i=1}}{\displaystyle \frac{1}{2}}{f}_i^{eq}(v_i^2+{\eta }_i^2)={\displaystyle \frac{1}{2}}\rho [(D+n)T+u^2],$

(39)${\mathit{M}}_2=\stackrel{N}{{\displaystyle \sum }_{i=1}}{f}_i^{eq}{\mathit{v}}_i{\mathit{v}}_i=\rho (T\mathit{I}+\mathit{u}\mathit{u}),$

(40)${\mathit{M}}_{3,1}=\stackrel{N}{{\displaystyle \sum }_{i=1}}{\displaystyle \frac{1}{2}}{f}_i^{eq}(v_i^2+{\eta }_i^2){\mathit{v}}_i={\displaystyle \frac{1}{2}}\rho \mathit{u}[(D+n+2)T+u^2],$

(41)${\mathit{M}}_3=\stackrel{N}{{\displaystyle \sum }_{i=1}}{f}_i^{eq}{\mathit{v}}_i{\mathit{v}}_i{\mathit{v}}_i=\rho (T\mathbf{\Theta }+\mathit{u}\mathit{u}\mathit{u}),$

(42)$\begin{array}{rl}{\mathit{M}}_{4,2}=& \stackrel{N}{\sum _{i=1}}{f}_i^{eq}\frac{v_i^2+{\eta }_i^2}{2}{\mathit{v}}_i{\mathit{v}}_i=\rho [\left(\frac{D+n+2}{2}T+\frac{u^2}{2}\right)T\mathit{I}\\ & +\left(\frac{D+n+4}{2}T+\frac{u^2}{2}\right)\mathit{u}\mathit{u}],\end{array}$

with $\mathbf{\Theta }=(u_{\alpha }{\delta }_{\beta \gamma }+u_{\beta }{\delta }_{\alpha \gamma }+u_{\gamma }{\delta }_{\alpha \beta }){\hat{\mathit{e}}}_{\alpha }{\hat{\mathit{e}}}_{\beta }{\hat{\mathit{e}}}_{\gamma }$, $({\hat{\mathit{e}}}_{\alpha },{\hat{\mathit{e}}}_{\beta },{\hat{\mathit{e}}}_{\gamma })$ denote unit vectors along the $\alpha$, $\beta$, and $\gamma$ axes of a fixed coordinate system. In Eqs. (36)−(42), the left side represents the analytical solution for the kinetic moment of the equilibrium distribution function $f^{eq}=f^{(0)}$. These Eqs. (36)−(42) can be written in matrix form as

(43)${\hat{f}}_k^{eq}=m_{ki}{f}_i^{eq},$

or

(44)${\hat{\mathit{f}}}^{eq}={\mathit{m}\mathbf{\cdot }\mathit{f}}^{eq}.$

Equations (36)−(42), or the matrix equation (43) or (44) represent the physical constraint that need to be satisfied for the discrete velocity selection in the DBM evolution equation.

If we want to achieve a DBM that includes second-order discrete/non-equilibrium effects, we only need to consider: which additional kinetic moments of $f^{(0)}$ should be preserved after adding the second-order non-equilibrium contributions to viscous stress and heat flux? Because

(45)$\begin{array}{rl}& {\mathbf{\Delta }}^{(2)}={\mathit{M}}_2(f_{}^{(2)})\Rightarrow {\mathit{M}}_2(f_{}^{(1)}),{\mathit{M}}_3(f_{}^{(1)}),\\ & \text{where}{\mathit{M}}_3(f_{}^{(1)})\Rightarrow {\mathit{M}}_4(f_{}^{(0)}),\end{array}$

(46)$\begin{array}{rl}& {\mathit{M}}_{3,1}(f_{}^{(2)})\Rightarrow {\mathit{M}}_{3,1}(f_{}^{(1)}),{\mathit{M}}_{4,2}(f_{}^{(1)}),\\ & \text{where}{\mathit{M}}_{4,2}(f_{}^{(1)})\Rightarrow {\mathit{M}}_{5,3}(f_{}^{(0)}).\end{array}$

It can be seen that there are two more kinetic moments, ${\mathbf{\Delta }}_4$ and ${\mathbf{\Delta }}_{5,3}$, that need to be conserved. It is clear that, in a DBM containing only up to the second-order discrete/non-equilibrium effects, the set of kinetic moments that need to be preserved are $\left\{{\mathit{M}}_0,{\mathit{M}}_1,{\mathit{M}}_{2,0},{\mathit{M}}_2,{\mathit{M}}_{3,1},{\mathit{M}}_3,{\mathit{M}}_{4,2},{\mathit{M}}_4,{\mathit{M}}_{5,3}\right\}$.

Fig.9 shows the kinetic moments that need to be preserved in the step-by-step DBM construction as the degree of discreteness/non-equilibrium represented by the Kn number increases and the corresponding hydrodynamic equations. It is easy to find that when considering up to the $j$-th order discreteness/non-equilibrium, the number of kinetic moments that need to be preserved is $3+2(j+1)$. For more details, please refer to references [65, 77-79]. It is emphasized again that the ability to recover the corresponding level of hydrodynamic equations is only a part of the physical functionality of DBM. Roughly equivalent to the physical functions of DBM is the extended hydrodynamic equation systems, which should include some of the most closely related non-conserved moment evolution equations in addition to the usual density, momentum, and energy evolution equations. Among them, the extended part is referred to the set of evolution equations of corresponding non-conserved moments, and its necessity increases with the increase of system discreteness/non-equilibrium degree. We provide a DBM simulation flowchart in Fig.10 corresponding to the work positioning shown in Fig.3.

Fig.9 Kinetic moments that should keep value in various levels of DBM.

Full size|PPT slide

Fig.10 Flowchart for DBM simulation.

Full size|PPT slide

3.6 Brief review of steady DBM

When conducting research on complex flow behaviors, the research perspective includes at least two dimensions of span, as shown in Fig.11: the span of characteristics that make up this set of behaviors and the time span of these behaviors. Among them, the span of kinetic behavior often correlates with the discrete/non-equilibrium degree of the system being described. In addition to the kinetic characteristics that evolve over time, people often need to conduct detailed research on some steady-state behavior of the system. The steady DBM is designed for this requirement.

Fig.11 Two dimensions of the perspective of complex flow behavior research.

Full size|PPT slide

Formally, in the discrete Boltzmann equation (15), setting the time derivative term to zero yields the control equation for steady-state DBM,

(47)${\mathit{v}}_i\cdot {\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }\mathit{x}}}+(\text{force}\text{term}{)}_i=(\text{collision}\text{term}{)}_i.$

However, the selection rules for discrete velocities here are completely different from the case of unsteady DBM, and the physical meaning of the discrete distribution function $f_i$ is the actual distribution function of corresponding discrete velocity ${\mathit{v}}_i$.

Considering that the distribution function $f$ has a small amplitude when $\left|\mathit{v}-\mathit{u}\right|$ is large, the simplest way in steady DBM is to uniformly discretize the particle velocity space in a finite region centered on the local flow velocity $\mathit{u}$. In the case of low particle number density, a more reasonable approach is to use non-uniform grids with local refinement around the local velocity field $\mathit{u}$ [80]. In Ref. [64], the discretization of particle velocity space introduces a Lagrangian grid that moves with local flow velocity $\mathit{u}$. In this way, as shown in Fig.12, in the three-dimensional case, the direct uniform discretization of particle velocity can be written as

Fig.12 Schematic of discretization of particle velocity space. The velocity discretization in $x$ direction is taken as example.

Full size|PPT slide

(48)$\begin{array}{rl}& v_{ix}=\left(i-\frac{N_{vx}+1}{2}\right)/\left(\frac{N_{vx}-1}{2}\right)v_{x,max},\\ & v_{jy}=\left(j-\frac{N_{vy}+1}{2}\right)/\left(\frac{N_{vy}-1}{2}\right)v_{y,max},\\ & v_{kz}=\left(k-\frac{N_{vz}+1}{2}\right)/\left(\frac{N_{vz}-1}{2}\right)v_{y,max},\end{array}$

where $(N_{vx},N_{vy},N_{yz})$ are the total numbers of grid in the particle velocity space at three degrees of freedom, $(v_{x,max},v_{y,max},v_{z,max})$ are the maximum values of particle velocity in the three dimensions. Because it is a uniform grid, the weight coefficients of each grid point are the same. In many cases, it is reasonable to take $v_{x,max}=v_{y,max}=v_{z,max}=v_{max}$. But in the case of local refinement, the weight coefficients of grid points are no longer the same. Taking direction of ${\mathit{v}}_y$ as an example, we provide an example of using a set of non-uniform grid $v_{jy}$ and weight coefficient $w_j$:

(49)$\begin{array}{rl}& v_{jy}={\left(j-\frac{N_{vy}+1}{2}\right)}^{\lambda }/{\left(\frac{N_{vy}-1}{2}\right)}^{\lambda }{v}_{y,max},\\ & w_j=\lambda {\left(j-\frac{N_{vy}+1}{2}\right)}^{\lambda -1}/{\left(\frac{N_{vy}-1}{2}\right)}^{\lambda }{v}_{y,max}.\end{array}$

The discretization of particle velocity space in steady DBM given here is only an example, The purpose is to highlight the differences from non-steady-state DBM. For more details, please refer to more professional literature.

A brief methodological review on the capabilities and characteristics of the two types of DBMs, time-dependent DBM and steady-state DBM, is provided as follows:

(i) Due to limitations in simulation capabilities, we often need to compromise and trade-offs between the two spans, time span and behavior span, of the descriptive perspective. The kinetic behavior span is closely related to the degree of discreteness/non-equilibrium. If we compress the time span of our attention, it is possible to exchange it for broadening the behavioral span of our attention. If the time span of the focus is compressed to a single point, the span of the behavior of the focus can be expanded to all the kinetic moments of the distribution function $f$, i.e., all the kinetic features of the complete and real $f$. This is the limit extension of the degree of discreteness/non-equilibrium that can be described under the effective within the framework of BGK-like model and Chapman−Enskog multi-scale expansion theory.

(ii) At this limit, the number of kinetic moments to be preserved becomes infinite. In the ideal case, $f_i$ becomes the real distribution function corresponding to ${\mathit{v}}_i$, and particle velocity ${\mathit{v}}_i$ becomes the real particle velocity. At this limit, the linear equations about $f^{(0)}$ that were originally designed to preserve only a finite number of kinetic moments, such as Eq. (43), become infinitely high-dimensional. Therefore, the original method of selecting discrete velocities is no longer applicable. When complexity and difficulty reach a certain level, they often turn around and become simpler. At this point, the discretization of particle velocity space can return to a direct discretization scheme similar to spatial and temporal discretization.

(iii) In time-dependent DBM, only a finite number of kinetic moments of $f$ are preserved, and what has a clear physical meaning is the kinetic moment of $f_i$, but $f_i$ itself does not have a clear physical meaning. The recovered distribution function is only approximately accurate, and is only accurate in terms of the features of interest. When it comes to the steady DBM, the focus is on the actual distribution function, which means that essentially all kinetic moments need to be preserved. It is not just the kinetic moments of that have clear physical interpretations, but also the distribution function itself, which represents the real distribution function at that particular particle velocity.

(iv) Both time-dependent DBM and time-independent DBM can be used together, achieving the ultimate or maximum descriptive power in two dimensions: the maximum in terms of time span (sacrificing behavior span or degree of discreteness/non-equilibrium, etc.) and the maximum in terms of degree of discreteness/non-equilibrium (sacrificing time span, etc.). The time-dependent DBM and steady-state DBM, similar to the time-dependent Schrödinger equation and the steady-state Schrödinger equation, also show significant differences in numerical simulation methods. Of course, it also contains the same or similar parts. It is noted that the Schrödinger equation has a clear physical image correspondence and is a physical model description equation. It does not contain numerical schemes for simulation.

(v) Compared to steady-state DBM, time-dependent DBM compresses the system’s behavior description ability from all (infinite) kinetic moments of $f$ to a finite set. The gain in the system’s behavior description ability is that the time span of the behavior is amplified infinitely from 0. However, there is a trade-off involved in this. The gain is “mesoscopic” and the cost is also “mesoscopic”. Both types of DBMs are typical representatives of “choosing something to do” or “doing something and abandoning the others”.

(vi) The transformation of coordinate systems from Eulerian to Lagrangian during the construction process of steady-state DBM in Ref. [64] has not yet involved discrete formats and belongs to the theoretical modeling part. Please note that MD simulation is Lagrangian, and the coordinate system follows the molecule. The reason why the coordinate system can run with molecule is because in MD, the position of molecule at the next moment is deterministic! In kinetic theory, the direction and speed $\mathit{v}$ of molecule at the next moment are uncertain and are probability events. So, in Lagrangian DBM, the coordinate system cannot follow the (uncertain and probabilistic) molecule, but can and should follow fluid element with a macroscopic flow the velocity $\mathit{u}$. This is the technical key in the construction of Lagrangian DBM.

4 Verification and validation

In this section we illustrate some verification and validation of the method. Fig.13-Fig.29 are for unsteady DBM, and Fig.30 is for steady DBM. It should be mentioned that the model verification and validation are based on known results, aims to verify and validate the chosen discrete format, program implementation, parameter settings, etc. This is a basic process but does not have any contribution to the new physics cognition.

Fig.13 Comparison of DBM simulation and DSMC simulation of a shock structure. Reproduced with permission from Ref. [77].

Full size|PPT slide

Fig.14 Schematic of velocity slip and temperature jump in Knudsen layer.

Full size|PPT slide

Fig.15 Comparison between DBM simulation results and analytical solutions: slip velocity. Reproduced with permission from Ref. [92].

Full size|PPT slide

Fig.16 Comparison between DBM simulation results and analytical solutions: velocity difference profile. Reproduced from Ref. [94] with permission.

Full size|PPT slide

Fig.17 Comparison between DBM simulation results and analytical solutions: temperature difference profile. Reproduced with permission from Ref. [94].

Full size|PPT slide

Fig.18 The results of Couette flow at steady state: (a) velocity profiles, the symbols are DSMC data and the lines denote DBM results; (b) viscous shear stress profiles, the symbols denote DSMC data, the dashed line represents the Lattice ES-BGK results, and the solid line represents the DBM results. Reproduced with permission from Ref. [93].

Full size|PPT slide

Fig.19 Pressure-driven flow. (a) Schematic, and (b) DBM simulation results of velocity contour under different Kn numbers. Reproduced with permission from Ref. [94].

Full size|PPT slide

Fig.20 (a) Profiles of pressure $p(x)$ along the central line under different Kn numbers. (b) Profiles of $u_x$ along the $y$ direction under different Kn numbers. Reproduced with permission from Ref. [94].

Full size|PPT slide

Fig.21 Inverse reduced mass flow rate under various Knudsen numbers, in which including the results of DBM simulation, NS simulation and experiment. Reproduced with permission from Ref. [94].

Full size|PPT slide

Fig.22 Simulation results of non-equilibrium cavity flow in steady state. (a) Temperature contour. (b) Heat flow streamline. (c) Horizontal velocity distribution on the vertical centerline. (d) Vertical velocity distribution on the vertical centerline. Reproduced with permission from Ref. [93].

Full size|PPT slide

Fig.23 Comparison of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ profile between DBM simulation results and analytical results. The three figures represent the cases of weak, medium, and strong TNE strength, respectively. The two symbols indicate the results of first-order and second-order DBM, respectively, with lines representing the analytical results at two TNE levels. Reproduced with permission from Ref. [77].

Full size|PPT slide

Fig.24 The DBM simulation results from different orders of DBM. Reproduced with permission from Ref. [77].

Full size|PPT slide

Fig.25 Comparison of ${\mathrm{\Delta }}_{2,xx}$ profile between DBM simulation results and analytical results. Each row of the pictures show corresponds to cases with weak, moderate, and strong TNE strength, and each column corresponds to the results of D2V13, D2V15, and D2V30, respectively. In the figures, the circles represent results from DBM at various levels, the green lines and green points are results calculated from the first theoretical analysis, and red lines denote results obtained from the second theoretical analysis. Reproduced with permission from Ref. [65].

Full size|PPT slide

Fig.26 Evolution of spike amplitude. The red points in the picture are experimental results from Ref. [102], and blue line represent DBM results. Reproduced with permission from Ref. [100].

Full size|PPT slide

Fig.27 Snapshots of schlieren images of the interaction between a shock wave and a heavy-cylindrical bubble. The odd rows represent experimental results from Ref. [34] with permission, and the even rows are DBM simulation results. Numbers in the picture represent the time in μs. Reproduced with permission from Ref. [103].

Full size|PPT slide

Fig.28 Evolution curves of the characteristic scale of shock−bubble interaction. The lines indicate DBM simulation results and the symbols are experimental results extracted from Ref. [34]. Reproduced with permission from Ref. [103].

Full size|PPT slide

Fig.29 Comparison between DBM results and NS results of pressure profile at the line $y=0.625\pi$ in simulation of Orszag−Tang magnetic turbulence. Reproduced with permission from Ref. [43].

Full size|PPT slide

Fig.30 (a) Profiles of the dimensionless mass flow rate with various rarefaction parameters under three different values of TMAC for the one-dimensional pressure driven flow through a microchannel. (b) Profiles of the dimensionless mass flow rate with rarefaction parameters for three different aspect ratios in the context of the two-dimensional pressure-driven flow through a microchannel. Reproduced with permission from Ref. [64].

Full size|PPT slide

4.1 Combustion, detonation and shock wave

The DBM study of complex flows started from combustion and detonation. The first work was completed by Bo Yan [81], a postdoctoral fellow of the first author’s research group, and then Chuandong Lin [82-89], Yudong Zhang [90], Yiming Shan [91], and some other graduate students of the research group participated in this direction. Currently, there is a relatively substantial body of literature on this subject from various research groups.

Detonation wave can be regarded as shock wave sustained by chemical reaction. The DBM in published literature is only concerned with the non-equilibrium induced by the small structure and/or fast flow behavior, so the description of chemical reaction uses the traditional phenomenological model. From the perspective of constructing new model, what needs to be verified is the capability to describe the structure and behavior of shock waves. Fig.13 shows the comparison of the shock wave structures obtained by DBM and DSMC. It should be pointed out that for the simulations of steady-state structures such as shock waves, the use of steady-state DBM is more effective.

4.2 Microscale flow

Micro-nano flow is a common phenomenon in micro-nano electromechanical systems. It is well-known that mass and heat transfer in micro-nano flows exhibit many unique behaviors, such as the appearance of Knudsen layer near the wall. Non-linear mass transport such as velocity slip, and non-Fourier heat transfer behavior such as temperature jump would come out within the Knudsen layer (as shown in Fig.14). There are two main reasons for these phenomena: (i) With the decrease in the transverse size of the channel, the average molecular distance is no longer negligible compared to the transverse characteristic scale of the flow behavior. (ii) The non-equilibrium driving mechanism in the fluid close to the wall is very different from that in internal fluid away from the wall, because the former carries the energy exchange (if no material exchange) between the system and the outside world. Therefore, in addition to considering the order of Kn numbers to be retained in the internal flow description, it is also necessary to consider the design of suitable kinetic boundary conditions [61, 92-94].

Couette flow, Poiseuille flow, and so on are classic problems in the study of slip phenomenon. For the convenience of discussion, the velocity Knudsen layer shown in Fig.14(a) and (b) and the temperature Knudsen layer shown in Fig.14(c) and (d) are discussed separately. The schematic diagram for the velocity Knudsen layer is shown in Fig.14(a). It depicts a scenario where two long parallel plates are filled with gas. The two plates move in the $y$ direction with velocities $u_L=-v$ and $u_R=v$, respectively, while velocities in the $x$ direction are zero. The temperature of the two plates is fixed at $T_0$. The following two cases would make the Kn number increase and the velocity slip phenomenon becomes more significant: (i) the distance between the plates decreases, (ii) the gas density between the two plates decreases. After entering the steady state, although nonlinear velocity distribution [BD segment in Fig.14(b)] and velocity slip appear in the Knudsen layer, the velocity distribution in the region away from the wall still changes linearly [AB segment in Fig.14(b)]. Extend the line segment AB to the wall, the intersection point is C, and the corresponding flow rate is $v_{\text{c}}$. At this time, the difference between the plate speed $v_w$ and $v_{\text{c}}$ is

(50)$v_{slip}=v_w-v_{\text{c}},$

which is defined as slip velocity. The distribution along $x$ direction of velocity difference $\mathrm{\Delta }u(x)$ between the curve BD and BC is called velocity Knudsen profile, which is an important parameter to describe the velocity Knudsen layer. The analytical solutions of $v_{slip}$ and $\mathrm{\Delta }u(x)$ can be found in Refs. [95, 96]. Considering the Maxwell-type boundary condition, the expression of slip velocity is $v_{slip}=k_s\lambda \sqrt{\pi }/2\cdot \text{d}v/\mathrm{d}x$, where $k_s$ is a coefficient related to the tangential momentum accommodation coefficients (TMAC, $\alpha$, which is used to measure the proportion of complete diffuse reflection) [96]. $\lambda$ represents the mean free path of molecules. In BGK collision model, ${\lambda }_{\mathrm{B}\mathrm{G}\mathrm{K}}=\mu /p\sqrt{8RT/\pi }$, and in hard sphere collision model, ${\lambda }_{\mathrm{H}\mathrm{S}}=4/5{\lambda }_{\mathrm{B}\mathrm{G}\mathrm{K}}$. Fig.15(a) shows the comparison of slip velocity normalized by $\text{d}v/\mathrm{d}x$ between DBM results and two types of analytical solutions at different Kn numbers [92]. In Fig.15(b), DBM simulation results for various TMAC ($\alpha$) are compared with analytical solutions. Fig.16 shows the curves of velocity difference profile $\mathrm{\Delta }u(x)$ for different Kn numbers, where $M$ is the number of discrete velocity directions [94]. The analytical solution of $\mathrm{\Delta }u(x)$ is $\mathrm{\Delta }u(x)=Y_0(\eta )\sqrt{2}\text{Kn}\cdot \text{d}{u}_y/\text{d}x$, where $Y_0$ is the velocity Knudsen-layer function [95] and $\eta =(1-x)/(\sqrt{2}\text{Kn})$. The DBM results are in good agreement with the analytical solutions, and it can be seen that DBM can accurately capture the velocity slip and the flow velocity behavior in the Knudsen layer.

The schematic diagram for the temperature Knudsen layer is depicted in Fig.14(c). Similar to that of plate shear flow, it involves two long parallel plates filled with gas, with zero velocity between the plates. The temperatures of the two plates are maintained at $T_L$ and $T_R$, respectively. After entering the steady state, although the nonlinear distribution of temperature [section BD in Fig.14(d)] and the temperature jump occurs in the Knudsen layer, the temperature distribution in the area away from the wall still changes linearly [section AB in Fig.14(d)]. Extending line segment AB to the wall, the intersection point is denoted as C, and the corresponding temperature is represented as $T_{\text{c}}$. At this point, the difference between the plate temperature $T_w$ and $T_{\text{c}}$ is

(51)$T_{slip}=T_w-T_{\text{c}},$

which is defined as the slip temperature. The temperature difference between the curves BD and BC, i.e., $\mathrm{\Delta }T(x)$, as a function of position, is known as the temperature Knudsen profile, which is an important parameter to describe the temperature Knudsen layer. The analytical solutions of $T_{slip}$ and $\mathrm{\Delta }T(x)$ can be found in Ref. [95]. Fig.17 displays the agreement of the temperature difference profiles $\mathrm{\Delta }T(x)$ between DBM results and analytical solutions, for cases with two different Kn numbers [94]. The analytical solution of temperature difference profiles is $\mathrm{\Delta }T(x)={\mathrm{\Theta }}_1(\eta )\sqrt{2}\text{Kn}\cdot \text{d}T/\text{d}x$, where ${\mathrm{\Theta }}_1(\eta )$ represents the temperature Knudsen-layer function. Fig.18(a) shows the comparison of velocity distribution between DBM results and DSMC results when Couette flow develops to a steady state, with two different Kn numbers [93]. Fig.18(b) shows the profile of the $xx$ component of shear stress, which contains the results of DBM, DSMC, and lattice ES-BGK under three Kn numbers. It can be seen that the DBM results are consistent with the DSMC results at several Kn numbers. However, at a high Kn number (Kn = 0.5), the DBM results are better than those of the lattice ES-BGK model.

In a pressure-driven Poiseuille flow, the wall is fixed and there is a pressure difference between the two ends of the pipe. The flow of fluid is driven by the pressure difference. Due to the action of the wall surface, the velocity of the fluid in the center of the pipeline is larger, gradually decreasing from the center to the wall, as shown in Fig.19(a). Fig.19(b) shows the $u_x$ contour for three different Kn numbers. Fig.20(a) illustrates the distribution of pressure $p$ along the $x$ direction at the center line. It can be seen that the pressure distribution is nonlinear and has negative curvature, which is consistent with the previous experimental results and analytical theory [94]. Fig.20(a) explains the phenomenon that appeared in Fig.19(b), showing that the more linear the pressure distribution, the smaller the difference in velocity distribution in the cross section. Fig.20(b) displays the velocity distribution $u_x(y)$ along the $y$ direction under different Kn numbers. It can be seen that as the Kn number increases, the deviation between the analytical solution based on the NS equations and DBM simulation results becomes more pronounced. This is physically reasonable because as the Kn increases, slip boundary analysis gradually becomes less applicable. Fig.21 shows the variation of inverse reduced mass flow rate with Kn number, including DBM simulation results, NS combined with slip boundary simulation results, and experimental results of two different gases ($\text{He}$ and N₂). It can be seen that when the Kn number is less than 0.1, the results of NS are in good agreement with the experimental results. With the increase of the Kn number, the difference between the results of NS and the experimental results becomes larger. However, the DBM result is still in good agreement with the experimental result when the Kn number approaches 0.5 [94].

4.3 Cavity flow

Two-dimensional square cavity flow is a common boundary-driven flow, which is also frequently encountered in non-equilibrium flow problems [93]. In a cavity with size $L\times L$ filled with gas, the upper wall of the cavity moves to the right at a constant horizontal velocity $U_w$, while the other three walls remain stationary. The temperature on the four walls is fixed at $T_0$. In a specific practical scenario, we choose argon as the gaseous medium in the square cavity. The mass of the argon atom is $\text{m}=6.63\times {10}^{-26}\mathrm{k}\mathrm{g}$, the temperature of the wall is $T_0=273\text{K}$, and the motion velocity of the upper wall is $U_w=50\mathrm{m}/\mathrm{s}$. DBM simulation uses dimensionless quantities. Fig.22(a) shows the temperature contour of the flow field when the square cavity flow develops to a steady state, and Fig.22(b) is the corresponding heat flux streamline. From Fig.22(a), it can be seen that the left side of the flow field is the low-temperature region and the right side is the high-temperature region. However, it is shown in Fig.22(b) that the heat flux is from the left side to the right side. This is inconsistent with Fourier’s thermal conductivity law but is reasonable. This is a typical strong thermodynamic non-equilibrium effect. According to the Chapman−Enskog multi-scale analysis, when the Knudsen number is small, only the first order term in Knudsen number is significant. In other words, when the degree of TNE is low, the linear response theory works. In this case, the heat flux is positively proportional to the negative temperature gradient. But with increasing the Knudsen number, the second and even higher order terms in the Knudsen number can no longer be ignored. In other words, with increasing the degree of TNE, the linear response theory gradually fails and the nonlinear responses begin to become more significant. The heat flux is no longer only related to the temperature gradient, but also influenced by velocity gradient and/or density gradient. When the negative contribution from density gradient and/or velocity gradient surpasses the contribution of the temperature gradient, the heat flux will flow from the low-temperature region to the high-temperature region. The DBM results in this case are confirmed by the DSMC predictions. Fig.22(c) and (d) demonstrate the comparison between DBM results and DSMC results of the velocity distribution, showing good agreement between the two.

4.4 Fluid collision

Collision of two fluids will lead to the increase of gradients of density, temperature, flow velocity, and pressure at the collision interface, resulting in strong non-equilibrium effects at the interface [77, 79]. Therefore, fluid collisions are often used to assess the ability of models to describe non-equilibrium effects. According to CE multi-scale analysis theory, the degree of nonlinearity of constitutive relation varies with the degree of non-equilibrium. Even for the same flow behavior, using different physical quantities to investigate, it may require different levels of accuracy in the DBM physics.

Fig.23 and Fig.24 show the distribution of the $xx$ component of ${\mathbf{\Delta }}_2^{\ast }$ (corresponding to the viscous stress tensor in the macroscopic fluid dynamics equations) for detecting the system’s non-equilibrium behavior [77]. The higher the degree of TNE (from the perspective of Kn number) considered in the modeling process, the corresponding constitutive relationship of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ in EHE is more accurate. The expressions for first-order and second-order non-equilibrium stresses can refer to Appendix B. It is clear that relaxation time $\tau$ and macroscopic gradients are the driving forces of non-equilibrium effects. For ${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$, the main driving factors are relaxation time $\tau$ and velocity gradient. But for ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$, it is also related to density gradient and temperature gradient. Quantities such as Kn number, relaxation time $\tau$, the spatial gradient and temporal change rate of macroscopic physical quantities, are all measures of TNE intensity/degree/strength from a certain perspective. The intensity of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ is also a measure of TNE strength from a particular perspective. Therefore, it is again stated that for complex flow studies, using a single perspective to define the non-equilibrium strength is easy to draw one-sided conclusions. This is the reason why DBM introduces a non-equilibrium intensity vector (with each component being a measure of TNE strength from a specific perspective) to carry out multi-perspectives cross-positioning for the non-equilibrium intensity.

Based on the above considerations, in order to demonstrate the ability of the DBM model to describe different degrees of non-equilibrium, three different cases are designed. If only see from the view of ${\mathrm{\Delta }}_{2,xx}^{\ast }$, cases 1, 2, and 3 correspond to weak, moderate, and strong non-equilibrium strengths, respectively. Among these cases, case 1 has an initial velocity of 0, so its strength of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ is the weakest. In case 2, the initial collision velocity is increased, thereby incorporating first-order non-equilibrium effects, resulting in a greater ${\mathrm{\Delta }}_{2,xx}^{\ast }$ strength compared to case 1. In case 3, we increase the relaxation time $\tau$ ten times compared with the $\tau$ in case 2, so that it produces a strength approximately ten times greater than that of case 2. However, the non-equilibrium strength of the system cannot be solely measured with a single indicator. For example, in case 1, its velocity gradient is small, so the density gradient and temperature gradient play a significant role, resulting in the relative value of ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$ to ${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ can not be negligible. Therefore, DBM further introduces the relative non-equilibrium strength quantity ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}/{\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ to more accurately describe the relative importance of different levels of non-equilibrium effects. From the perspective of relative non-equilibrium strength, the non-equilibrium strength in case 2 is smaller than that in cases 1 and 3.

In Fig.23, the physical accuracies of the used DBM are first order and second order (represented by circles with different colors), respectively. The analytic solution results correspond to first order [${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$] and second order [${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ + ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$] (represented by lines with different colors), respectively. It can be seen that in case 1, although the non-equilibrium strength from the view of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ is weak, the results of the first order DBM show significant deviations from both first-order and second-order analytical results. This is because in case 1, ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$ is non-negligible relative to ${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ [i.e., the relative non-equilibrium strength ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}/{\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ is relatively large]. The physical accuracy of second order non-equilibrium (from the perspective of Kn number) is not guaranteed in the first order DBM during the model construction process, so the first order DBM does not guarantee the accuracy of description of ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$. However, the results of second-order DBM match perfectly with the second-order analytical results. This is because in the construction process of second-order DBM, the accuracy of the second-order non-equilibrium description (from the Knudsen number perspective) is ensured. In case 2, although the non-equilibrium strength from the view of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ is increased, ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$ becomes negligible compared to ${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ [${\mathrm{\Delta }}_{2,xx}^{\ast (2)}/{\mathrm{\Delta }}_{2,xx}^{\ast (1)}$ is small]. As a result, the results of first-order DBM, second-order DBM, and first-order and second-order analytical results are all in agreement. In case 3, the non-equilibrium strength from the view of ${\mathrm{\Delta }}_{2,xx}^{\ast }$ is further increased, ${\mathrm{\Delta }}_{2,xx}^{\ast (2)}$ becomes non-negligible compared to ${\mathrm{\Delta }}_{2,xx}^{\ast (1)}$. In this case, the results of first-order DBM show deviations from the second-order analytical results, while the results of second-order DBM align well with the second-order analytical results. Fig.23 demonstrates one specific aspect of the complexity of non-equilibrium flow behavior.

The higher the physical accuracy of the DBM model, the more kinetic moments need to be retained in the modeling process, which in turn increases the computational cost of simulations. So, in the actual simulation study, how many orders of physical accuracy of DBM need to be considered? Meeting the needs, the simplest, the least cost, is the first choice. The required level of physical accuracy can be determined through the convergence of simulation results. The specific examples are shown in Fig.24, in which the above three working cases are simulated with first order to sixth order DBM (from the view of Kn number), respectively. It can be observed that, except for the results of the first order DBM, the rest of the DBM results are consistent. For cases 1 and 2, to obtain accurate viscous stresses, at least a second order model is required. However, for case 3, the second-order DBM results differ from the higher-order DBM results and the results from the third-order model and above overlap. Therefore, for case 3, to accurately describe the non-equilibrium behavior from the view of ${\mathrm{\Delta }}_{2,xx}^{\ast }$, at least a third order DBM is required.

4.5 Multiphase flow

Multiphase flows and phase change heat transfer phenomena are widely exist in the natural world and engineering applications [97-99]. Reference [65] further explores the ability of different orders of DBM (from the perspective of Kn number) to describe the non-equilibrium effect in liquid−vapor phase transition systems. Reference [65] analyzed the ability of the DBM model to describe different degrees of non-equilibrium from four basic perspectives, i,e., ${\mathrm{\Delta }}_{2,xx}$, ${\mathrm{\Delta }}_{3,xxx}$, ${\mathrm{\Delta }}_{3,1x}$, and ${\mathrm{\Delta }}_{4,2xx}$. Fig.25 displays the spatial distribution of typical thermo-dynamic non-equilibrium effects, i.e., ${\mathrm{\Delta }}_{2,xx}$, (the combination of hydrodynamic non-equilibrium and thermodynamic non-equilibrium) near the gas-liquid interface. If only look at the view from ${\mathrm{\Delta }}_{2,xx}$, each row of Fig.25 corresponds to weak, moderate, and strong non-equilibrium scenarios. Each column in the figure corresponds to three cases where the selection of discrete velocity is D2V13, D2V15, and D2V30, respectively. In terms of the retained order of the Knudsen number, the first two correspond to the NS level (retaining up to the first-order terms of the Knudsen number). However, D2V15 retains $M_{5,1}$ moment relations more than the D2V13 model, making the constitutive relations in the EHE more complete. The D2V30 case corresponds to the super-Burnett level (retaining up to the third-order terms of the Knudsen number). For comparison, the figure also provides analytical results for the first-order (${\mathrm{\Delta }}_{mn}^{(1)}$) and second-order (${\mathrm{\Delta }}_{mn}^{(2)}$) non-equilibrium effects. For the case of weak non-equilibrium (from the perspective of ${\mathrm{\Delta }}_{2,xx}$), as can be seen in the first row of Fig.25, there exists a deviation between the first-order analytical solution and the second-order analytical solution. It indicates that the second order non-equilibrium (${\mathrm{\Delta }}_{2,xx}^{(2)}$) cannot be ignored compared with the first order non-equilibrium (${\mathrm{\Delta }}_{2,xx}^{(1)}$), that is, the relative non-equilibrium strength ${\mathrm{\Delta }}_{2,xx}^{(2)}/{\mathrm{\Delta }}_{2,xx}^{(1)}$ is larger. Because the first order DBM does not guarantee the accuracy of the description of ${\mathrm{\Delta }}_{2,xx}^{(2)}$ in model construction, the results of D2V13 and D2V15 are both deviated from the analytical solutions. However, the D2V30 model which considers up to the third order non-equilibrium, provides numerical results that match well with the analytical results. For the case where the non-equilibrium strength is moderate (from the perspective of ${\mathrm{\Delta }}_{2,xx}$) and the relative non-equilibrium strength ${\mathrm{\Delta }}_{2,xx}^{(2)}/{\mathrm{\Delta }}_{2,xx}^{(1)}$ is small, as shown in the second row of Fig.25, all three models can effectively describe non-equilibrium. In the case of strong non-equilibrium (from the perspective of ${\mathrm{\Delta }}_{2,xx}$), as shown in the third row of Fig.25, the results from the D2V30 model which considers the third order non-equilibrium agree well with the analytic solutions, while the other two models deviate significantly. But relatively speaking, the results from D2V15 are better than that from D2V13. These results indicate that higher-order DBM models have a stronger ability to describe non-equilibrium effects. As discussed in Sections 4.4 and 4.5, the Kn number, relaxation time $\tau$, macroscopic quantity gradient, TNE quantity, and relative TNE quantity can all be used to describe the non-equilibrium strength, but the results obtained by relying on a single quantity and from a single perspective are obviously incomplete. It is more reasonable to introduce the non-equilibrium strength vector and describe the non-equilibrium strength from multiple perspectives.

4.6 Hydrodynamic instability

RM instability under re-shock condition is an important research topic in ICF field. In Ref. [100], the re-shock process of ${\text{SF}}_6$ and air was studied by using two-fluid DBM, and the corresponding amplitude evolution curve was shown in Fig.26. It can be seen that DBM can capture the amplitude evolution accurately. The RMI under same conditions were also simulated in Ref. [101]. It should be noted that in the literature [101], the numerical simulation uses Euler equation, so artificial viscosity has to be added, and in order to better conform to the experiment, the diffusion effect was also artificially corrected. In DBM simulations, such artificial interventions are not used.

4.7 Shock−bubble interaction

The interaction between a shock wave and a heavy-cylindrical bubble (shock−bubble interaction, SBI) is studied in Ref. [103]. Fig.27 shows the comparison of snapshots of schlieren images between experimental results and DBM simulation results. In Fig.27, the odd rows represent experimental results from Ref. [34], and the even rows are DBM simulation results. The schlieren images of DBM simulation are calculated from the density gradient formula, i.e., $|\mathrm{\nabla }\rho |/|\mathrm{\nabla }\rho {|}_{max}$, with $|\mathrm{\nabla }\rho |=\sqrt{(\mathrm{\partial }\rho /\mathrm{\partial }x{)}^2+(\mathrm{\partial }\rho /\mathrm{\partial }y{)}^2}$. The typical wave patterns and the bubble’s main characteristic structure are marked out in the figures. For example, “TS” represents the transmitted shock propagating downstream inside the bubble, which is generated after the incident shock impacts the bubble. The shock focusing occurs when two high-pressure regions meet, which in turn produces secondary transmitted shock waves (STS), Reflected Rarefaction Waves (RRW), and jet. These typical wave patterns during the SBI can all be captured by DBM. The evolution of bubble characteristic scales during the SBI process are demonstrated in Fig.28. The lines in the figure are DBM simulation results, and the symbols are the experimental results extracted from Ref. [34]. It can be seen that DBM can accurately describe the bubble deformation process in the SBI process.

4.8 Plasma system

The compressible Orszag−Tang (OT) magnetic turbulence problem was first proposed by Orsazg and Tang. It is widely used to validate the effectiveness of magnetohydrodynamics models due to the complex wave structures that emerge during the evolution of the flow field. Fig.29 shows the pressure distribution at the line $y=0.625\pi$ when the OT magnetic turbulence problem develops to steady state [43]. The black lines represent the DBM simulation results, and the red dots are the simulation results from Ref. [104]. DBM can accurately capture the pressure distribution in this problem.

4.9 Steady flow

Fig.30 presents a set of validation results for steady-state DBM simulations [64]. Specifically, Fig.30(a) and (b) show the curves of mass flow rate $Q_p$ versus the rarefaction coefficient $\delta$ in one- and two-dimensional pressure-driven flow problems, respectively. Fig.30(a) displays results for different TMAC $\alpha$ and Fig.30(b) shows the results for different aspect ratios $h/w$ of pipelines. The lines represent DBM simulation results, and the symbols represent data from previous literature. The results show that DBM can accurately describe the mass flow rate of micropipe flow with different TMAC and different aspect ratios.

5 Applied research

To study the unknown is the original intention and purpose of the new model. New models are often built because existing models do not meet requirements. A new model, if it does not exceed any of the physical functions of the earlier model, the necessity of its existence is greatly reduced! In addition to the introduction of some specific research works, the author also wants to convey a more important idea: as the degrees of discreteness and thermodynamic non-equilibrium increase, the complexity of system behavior rises sharply; it is incomplete to rely only on the macroscopic quantities in the NS model, and more physical quantities need to be introduced to describe the state and behavior of the system to ensure that the actual control ability of system behavior does not decline. Therefore, this paper presented below includes results obtained from multiple analysis methods and approaches. From the perspective of the ways of extracting complex physical field information, the following results are mainly obtained by three methods: macroscopic quantitative analysis, morphological feature analysis, and thermodynamic non-equilibrium effect/behavior analysis. Among these, the thermodynamic non-equilibrium effect/behavior analysis means the description method based on the non-conserved moment of $f-f^{(0)}$, as shown in Sections 3.1 and 3.3. It should be noted that the definition of non-equilibrium strength in the following results may not be the same. However they are all selected according to the specific system characteristics, and all describe the characteristics of the system behavior from the corresponding perspective. Morphological feature analysis refers to the morphological analysis method based on Minkowski measure [61]. In a $D$-dimensional space, a set of convex sets satisfying motion invariance and additivity can be completely described by $D+1$ Minkowski measures. For the two-dimensional case, the three Minkowski measures are the proportion $A$ of the area of the high-$\mathrm{\Theta }$ region, the length $L$ of the boundary between the high-$\mathrm{\Theta }$ and low-$\mathrm{\Theta }$ regions, and the Euler characteristic coefficient $\chi$ describing the degree of connectivity of the regions. Here $\mathrm{\Theta }$ can be any physical quantity, such as density, temperature, non-equilibrium quantity, etc. Different methods and different angles of the same method constitute a more complete description of the behavior of fluid systems.

Similarly, we first introduce the application of unsteady DBM and then introduce the application of steady DBM. In terms of physical results, the simulation results shown below are divided into two categories, one that can also be obtained by NS, and another that NS cannot provide. The former, such as the evolution characteristics of interface spikes in RT instability systems and the deformation process of bubbles during the shock-bubble interaction, will not be introduced here. The latter represents where DBM surpasses NS. These previously poorly understood discrete effects and non-equilibrium behavior characteristics contain a large number of physical functions to be developed. The physical problems that DBM is concerned with include, but are not limited to, hydrodynamic instability, phase separation/non-equilibrium phase change, bubble/droplet collision, fusion and fragmentation, shock and detonation, and plasma kinetic theory.

5.1 Detonation and shock wave

Entropy production rate is an important concern in many fields related to compression science, such as shock wave and detonation physics, ICF, and aerospace [105]. By detecting TNE features from different perspectives, we can study the main mechanisms causing entropy production in the system and their relative importance. In detonation problems involving chemical reactions, the entropy production rate contributes in three ways: chemical reactions, NOMF, and NOEF. As shown in Fig.31, the entropy production rate for detonation problems involving chemical reactions is investigated [90]. The results show that the entropy production rate ($\mathrm{\Delta }{S}_3$) caused by the chemical reaction is much larger than the other two ($\mathrm{\Delta }{S}_1$ and $\mathrm{\Delta }{S}_2$).

Fig.31 Comparison of three entropy production rates in fluid systems of detonation problems with chemical reactions. Reproduced with permission from Ref. [90].

Full size|PPT slide

In many cases, the initial shock waveform is not flat. In the case where the initial waveform can be approximated as a plane plus disturbance, such shock waves are often referred to as disturbance shock waves for convenience of description. In general, the waveform of such a disturbed shock gradually returns to a plane as it propagates through the material. Such a process often referred to as shock shaping. The RMI caused by the interaction between a plane shock wave and a disturbed interface has been extensively studied. In recent years, more attention has been paid to the RMI caused by interaction between a disturbed shock wave and plane interface [106-108]. At present, all the studies on shock shaping and RMI induced by the interaction between disturbed shock wave and plane interface in the literature are based on traditional fluid modeling. Shan et al. [109] are using DBM modeling and analysis method to investigate the kinetic behaviors, particularly those ignored by traditional fluid modeling.

5.2 Shock−bubble interaction

Below, we present examples from Refs. [103] and [110] to illustrate the specific applications of TNE quantities. These two works introduce a non-equilibrium intensity vector to study the non-equilibrium effect and behavior of the SBI process from multiple perspectives. Fig.32 shows the contour for two typical non-equilibrium quantities of two different fluid components during the SBI process [110]. In this figure, the first two rows are the spatial distribution of $|{\mathbf{\Delta }}_2^{\sigma \ast }|$, and the last two rows correspond to the spatial distribution of $|{\mathbf{\Delta }}_{3,1}^{\sigma \ast }|$. The odd rows are the distribution of component A and the even rows are the distribution of component B. Different columns in the figure correspond to different typical moments. The location and characteristics of shock waves (incident shock wave, transmitted shock wave, etc.) and rarefaction wave are clearly captured by the non-equilibrium quantity. Their spatial distribution is highly related to interfaces with strong macroscopic gradients, so they can be used to identify interface features. At the same time, these non-equilibrium quantities also have rich physical connotations in phase space. The positive and negative values of non-equilibrium quantity represent the direction that deviates from the equilibrium state, and the magnitude of a non-equilibrium quantity indicates the extent of deviation from equilibrium. The non-equilibrium strength of the flow field is analyzed from another perspective. Fig.33 shows the spatial profile of the average non-equilibrium strength along the $y$ direction obtained by summing and averaging the non-equilibrium quantity along the $x$ direction. Fig.33 also shows the viscous effects on the two types of average non-equilibrium strengths. It can be seen that the spatial profile of these non-equilibrium quantities shows interesting symmetries, in which the rich and complex non-equilibrium effects and mechanisms are far from being fully understood. Taking ${\overline{\mathrm{\Delta }}}_{2,\alpha \beta }^{A\ast }$ as an example, the profiles of ${\overline{\mathrm{\Delta }}}_{2,xx}^{A\ast }$ is symmetry about the line $y=0.6$, indicating that the upper and lower parts of the flow field deviate from equilibrium in the same way. From the perspective of ${\overline{\mathrm{\Delta }}}_{2,xy}^{A\ast }$, the profile is symmetric about the origin, suggesting that the upper and lower portions of the flow field deviate from equilibrium in opposite directions but with consistent magnitudes. The profiles of ${\overline{\mathrm{\Delta }}}_{2,xx}^{A\ast }$ and ${\overline{\mathrm{\Delta }}}_{2,yy}^{A\ast }$ are both symmetric about the line ${\overline{\mathrm{\Delta }}}_{2,\alpha \beta }^{A\ast }=0$, indicating that when looking at the non-equilibrium strength of the fluid field from these two different perspectives, yields results with consistent magnitudes but opposite directions. In other words, looking at the non-equilibrium strength of the system from different perspectives, the result may be the same or opposite. Fig.33 only displays results for a particular time step, while results for other time steps are equally important. It is because the non-equilibrium characteristics of the flow field are very complex and far from being fully understood, so multi-level and multi-perspective research is needed. Understanding the time evolution of non-equilibrium quantities is helpful in understanding the kinetics of SBI processes. Fig.34 shows the global TNE strength of the two fluid components, where $t_1$ corresponds to the moment when the incident shock wave just sweeps the bubble, and $t_2$ is the moment when the incident shock wave exits the flow field, and these two moments correspond to the extremum of the global TNE strength curve. After $t_2$, it can be seen that the global TNE strength curve is oscillatory, which is attributed to the influence of reflected shock waves. Hence, global TNE strength can provide a comprehensive reflection of certain flow field characteristics, serving as another perspective for understanding system properties. Different perspectives and different levels of description constitute a relatively complete understanding of the fluid system.

Fig.32 Contours of the two typical TNE quantities $|{\mathbf{\Delta }}_2^{\sigma \ast }|$ and $|{\mathbf{\Delta }}_{3,1}^{\sigma \ast }|$ at three different moments. Reproduced with permission from Ref. [110].

Full size|PPT slide

Fig.33 Profiles of average TNE strength along the $y$ direction. Reproduced with permission from Ref. [110].

Full size|PPT slide

Fig.34 Evolution curves of the quantity of global TNE strength. Reproduced with permission from Ref. [110].

Full size|PPT slide

5.3 Phase transition and phase separation

It is found that in the case of non-equilibrium phase transition and phase separation, as shown in Fig.35, the TNE strength $D$ gradually increases in the first stage of phase separation (Spin Decomposition stage, SD stage), and gradually decreases in the second stage (Domain Growth stage, DG stage). Therefore, its maximum point can be used as the first physical criterion to divide the two stages [111]. Similar characteristics exist in the evolution of morphological characteristic quantity $L$, so it can also serve as a physical criterion to divide the two stages. Fig.36 shows the evolution of logarithm of the characteristic length $\lg R$, non-equilibrium strength ${\overline{D}}_2^{\ast }$, boundary length $L$ and entropy production rate ${\dot{S}}_{sum}$ during phase transition. The changing trends of these quantities can be used to identify the two stages [112]. Specifically, the entropy production rate ${\dot{S}}_{sum}$ gradually increases in the first phase of phase separation (SD stage), and gradually decreases in the second phase (DG stage), so its maximum point can be used as the second physical criterion to divide the two stages.

Fig.35 Evolution curves of boundary length $L$ and non-equilibrium strength $D$ in non-isothermal phase separation process. Reproduced with permission from Ref. [111].

Full size|PPT slide

Fig.36 Evolution curves of several characteristic quantities in the phase separation process. Reproduced with permission from Ref. [112].

Full size|PPT slide

5.4 Phase transition and bubble coalescence

During the non-equilibrium phase transition and bubble coalescence, the characteristics of each component of average stress exhibit distinct stages [113]. Therefore, the quantity of ${\overline{\mathrm{\Delta }}}_{2,xx}^{\ast }$ inside the bubble can be used to calibrate the bubble merging features and divide the merging stages. As shown in Fig.37, the evolution of ${\overline{\mathrm{\Delta }}}_{2,xx}^{\ast }$ shows two characteristic instants: the moment $t_{umax}$ when the mean coalescence speed $\overline{u}$ gets the maximum (at this time the ratio of minor and major axes is about 0.5) and the moment when the ratio of minor and major axes becomes 1 for the first time. These two characteristic instants divide the process of bubble coalescence into three stages. As also shown in Fig.38, the average TNE strength ${\overline{D}}^{\ast }$, coalescence acceleration $\overline{a}$ and the slope of boundary length $\mathrm{d}L/\mathrm{d}t$ are highly correlated in the early stage, and their amplitudes reach their maximum values almost simultaneously. Therefore, their maximum point can be regarded as the end of the first stage of bubble coalescence and the beginning of the second stage. In another work, Ref. [114], it is shown that the temporal evolutions of the total TNE strength and the total entropy production rate can both provide physical criteria to distinguish the stages of droplet coalescence.

Fig.37 Evolution curves of average coalescence velocity $\overline{u}$ and ${\overline{\mathrm{\Delta }}}_{2,xx}^{\ast }$. Reproduced with permission from Ref. [113].

Full size|PPT slide

Fig.38 Evolution curves of the average strength ${\overline{D}}^{\ast }$, average coalescence velocity $\overline{a}$, the slope of boundary length $\mathrm{d}L/\mathrm{d}t$, and $\overline{(\mathrm{\nabla }\mathit{u}:\mathit{u}{)}^{0.5}}$. Reproduced with permission from Ref. [113].

Full size|PPT slide

5.5 Droplet collision

In droplet collision kinetics, the characteristic of the average non-organized momentum flow (NOMF) strength (${\overline{D}}_2^{\ast }$) can be used to distinguish the type of droplet collision and divide the collision stage [36]. As shown in Fig.39, in the collision−fusion case, in the late period of collision (after point D in the figure), ${\overline{D}}_2^{\ast }$ rapidly decreases and then oscillates towards to a steady state. However, in collision−separation case, ${\overline{D}}_2^{\ast }$ continues to exhibit significant oscillations during the descent.

Fig.39 Evolution curves of two TNE quantities (${\overline{D}}_2^{\ast }$ and ${\overline{D}}_3^{\ast }$) during the droplet collisions, where (a) represents the two droplets are fused together after the head-on collision, (b) indicates the two droplets are separated after collision. Reproduced with permission from Ref. [36].

Full size|PPT slide

5.6 RTI and anti-RTI: Effects of intermolecular potential and high compressibility

In the study of RTI with consideration of intermolecular potential effects, the surface tension has a stage-wise impact on perturbation amplitudes, bubble velocity, and two entropy generation rates [115]. The interface tension inhibits the evolution of RTI during the bubble acceleration stage. In the stage of bubble asymptotic velocity, the interface tension first promotes and then inhibits the evolution of RTI. In addition, as shown in Fig.40, the first maximum point of $\mathrm{d}L/\mathrm{d}t$ (green line in the figure) and the first maximum point of $\mathrm{d}{\dot{S}}_{\text{NOMF}}/\mathrm{d}t$ induced by NOMF (purple line in the figure) can both be used as physical criteria for the bubble velocity entering the asymptotic stage in RTI. Interface length $L$ is also one of the perspectives for the study/characterization of TNE behavior, and the exponential growth stage of $L$ corresponds to the bubble acceleration stage.

Fig.40 Evolution curves of change rate of entropy production rate $\mathrm{d}{\dot{S}}_{\text{NOMF}}/\mathrm{d}t$, change rate of boundary length $\mathrm{d}L/\mathrm{d}t$, and bubble velocity $V_{\text{bubble}}$. Reproduced with permission from Ref. [115].

Full size|PPT slide

RTI under constant acceleration has been studied extensively and a relatively mature physical understanding has been obtained, but the acceleration of the RTI occurring in ICF is variable. The research on RTI under variable acceleration is relatively few, and the physical understanding is insufficient. The reverse of acceleration will cause the anti-RTI mechanism and behavior, and Jia et al. [116] are conducting research on this issue. It has been shown that the magnitude of acceleration and the time of acceleration reversal have significant influences on the evolution of RTI, such as the location of the interface, the location and development trend of the spike bubble, the evolution of local and global TNE characteristics and effects. In the strongly compressible case, some new behaviors and mechanisms may also emerge in the evolution of the RTI system. For example, the density of the fluid from the upper layer which is originally higher may become lower dense than that of the surrounding originally lighter fluid, resulting in the exchange of identities of the light and heavy fluids, leading to anti-RTI mechanism and behavior. The emergence of these new mechanisms and behaviors invalidate earlier methods of determining the source of matter particles based on density. Therefore, Chen et al. [117] introduced tracer particles into the DBM model of van der Waals fluids to study the influence of flow compressibility on the evolution of single-mode RTI systems of van der Waals fluids from the perspective of statistical physics combined with morphological analysis.

5.7 KHI and coupled RT-KHI system

In order to solve the problem of analyzing various complex physical fields during the evolution of KHI, in the literature [24], Gan et al. combined tracking non-equilibrium behavior characteristics and morphological analysis techniques to carry out physical identification and tracking technology design of feature structures or patterns, and quantitatively characterize the width and development rate of KHI mixed layer. It is found that the width of the mixed layer, the non-equilibrium strength and the length of the interface boundary show a high spatiotemporal correlation. The viscosity, which is highly correlated with the NOMF, inhibits the development rate of KHI, prolongates the duration of linear phase, reduces the maximum disturbance kinetic energy, improves the overall and local thermal non-equilibrium strength, and expands the range of non-equilibrium. The thermal conduction effect, generally consistent with the viscous effect, inhibits the growth of structures with small wavelengths during the evolution of KHI, thus making it easier for the system to form large-scale structures. Different from the consistent inhibition of viscosity on KHI, the thermal conduction effect is first suppressed and then increased in phases when the heat conduction is strong. The inhibition is due to the heat conduction extending the interface width, reducing the macroscopic quantity gradient and the non-equilibrium driving force strength, and the enhancement is due to the broadened density transition layer making the KHI more easily absorb energy from the fluid on both sides, thus enhancing the KHI. The competition of these two effects accompanied the whole process of KHI evolution, the inhibitory effect is dominant in the early stage, and the enhancement effect is dominant in the late stage.

In the coupled RT-KHI system, the TNE quantities can serve as criteria for identifying the transition of dominant instability mechanisms [17, 118]. As shown in Fig.41, in the case where KHI dominates at an earlier time and the RTI dominates at a later time (corresponding to shear dominance and buoyancy dominance, respectively), the boundary length $L$ and NOEF strength show similar linear growth stages. The end of the linear growth stage of boundary length $L$ and NOEF strength $D_{3,1}$ can be used as a physical criterion for the transition from KHI dominance to RTI dominance [118].

Fig.41 Evolution curves of boundary length $L$ and strength of NOEF $D_{3,1}$ in coupled RT-KHI system. These two figures correspond to the cases with different initial shear velocity. Reproduced with permission from Ref. [118].

Full size|PPT slide

5.8 Shock wave/boundary layer interaction

There exist complex shock wave/boundary layer interaction phenomena in supersonic flow. The main sources of aerodynamic drag in scramjet engines are “friction drag” caused by boundary layer and “wave drag” caused by shock wave, both of which are directly related to entropy increase, so entropy increase is the key parameter to evaluate the aerodynamic drag. Song et al. [105] recently investigated this problem using the DBM modeling and analysis method. Their results show that for regular reflection, the non-equilibrium intensity of the reflected shock wave is stronger than that of the incident shock wave. For shock wave/laminar boundary layer interaction, entropy production rate caused by viscosity is dominant in shock wave, and entropy production rate caused by heat conduction is dominant in boundary layer. The intensity of the two entropy production rates increases with the increase of Mach number. The research results can provide theoretical guidance for evaluating the flow quality in the inlet.

5.9 Plasma shock wave

In plasma research, TNE characteristics can be used to physically distinguish plasma shock waves from neutral fluid shock waves. This distinction can help in the design of tracking techniques for the respective interfaces [119]. Fig.42 shows profiles of TNE quantity near the shock front in plasma shock wave under different Ma numbers. The left four figures are for ${\mathbf{\Delta }}_2^{\ast }$ and the right four figures are for ${\mathbf{\Delta }}_{4,2}^{\ast }$. In a neutral fluid shock wave, the TNE quantities ${\mathbf{\Delta }}_2^{\ast }$ and ${\mathbf{\Delta }}_{4,2}^{\ast }$ exhibit symmetry with their $xx$ and $yy$ components deviating from equilibrium in opposite directions (results for neutral fluids are not shown here). For plasma shock waves, the profiles of non-equilibrium quantity ${\mathbf{\Delta }}_2^{\ast }$ behave the same as that of neutral fluid. However, the profile of ${\mathbf{\Delta }}_{4,2}^{\ast }$ does not exhibit the symmetric deviation. Thus, these differences can be used to distinguish plasma shock waves from neutral fluid shock waves.

Fig.42 Profiles of ${\mathbf{\Delta }}_2^{\ast }$ and ${\mathbf{\Delta }}_{4,2}^{\ast }$ of plasma shock wave with different Ma numbers. Reproduced with permission from Ref. [119].

Full size|PPT slide

5.10 Plasma Orszag−Tang vortex

Through a DBM for plasma system, the kinetic behaviors of the Orszag−Tang vortex problem are investigated for the first time [43]. Fig.43(a) shows the contour of total TNE strength at $t=3$. It is observed that the total TNE strength is mainly distributed near the shock front. Fig.43(b) shows the evolution of four kinds of entropy production rates with time. In general, the entropy production rate is related to the difficulty of compression. From Fig.43(b), it is found that the compression difficulty is divided into two stages. Before $t=1.7$, the entropy production rate and compression difficulty increase with time. After $t=1.7$, the entropy production rate and compression difficulty decrease with time.

Fig.43 TNE effects of Orszag−Tang vortex problem. (a) Contour of total TNE strength at $t=3$. (b) Evolution of four kinds of entropy production rates with time. Reproduced with permission from Ref. [43].

Full size|PPT slide

5.11 Plasma RMI

The influence of magnetic field on KHI has been extensively studied since the 1960s, where transport characteristics such as viscosity and heat conduction are the focus of attention. The physical models used in previous studies are mainly macroscopic magnetohydrodynamic (MHD) models based on the continuum assumption. However, the MHD is valid only when the time- and length-scales of particle collisions are sufficiently small compared to those of the system behaviors under consideration. Rinderknecht et al. [40] pointed out that the kinetic behaviors occurring on the time- and length-scales of particle collisions may cause the hydrodynamic results to deviate from the experimental results. Song et al. [43] studied RMI kinetics using DBM modeling and analysis method. As a preliminary study, only the case of first order discrete/TNE effect is considered in this work. A few typical results are as follows: (i) As shown in Fig.44, before the RMI interface inversion, the magnetic field indirectly enhances the TNE strength ${\overline{D}}_T$ by inhibiting the interface development. After the interface inversion, the TNE strength ${\overline{D}}_T$ is significantly reduced. The minimum value of the global average TNE strength can be used as a physical criterion to determine the critical magnetic field strength under which the RMI interface amplitude could be reduced to 0. (ii) As shown in Fig.45 and Fig.46, the magnetic field has little effect on the entropy production rate contributed by the NOMF ${\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}$, but has a strong inhibition effect on the entropy production rate contributed by the NOEF ${\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}$. The minimum value of ${\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}$ can be used as a physical criterion to determine the critical magnetic field strength under which the RMI interface amplitude could be reduced to 0.

Fig.44 Evolution of global average TNE effects with different initial applied magnetic fields from $t=0$ to 500. Reproduced with permission from Ref. [43].

Full size|PPT slide

Fig.45 Evolution of entropy production rates from $t=0$ to 500 with magnetic fields ranging from 0.01 to 0.05. Reproduced with permission from Ref. [43].

Full size|PPT slide

Fig.46 Evolution of entropy production rates from $t=0$ to 500 with magnetic fields ranging from 0.10 to 0.30. Reproduced with permission from Ref. [43].

Full size|PPT slide

5.12 Plasma KHI

Song et al. [120] introduced tracer particles into DBM to study the influence of magnetic field on the transport and mixing characteristics of plasma KHI, as well as the influence of magnetic field on various kinetic behaviors. Preliminary results are as follows: (i) For the cases of weak magnetic fields, the magnetic field could inhibit the development of KHI at an early stage, but will lead to the formation of small-scale vortex structure and cause the KHI to enter the mixing stage earlier. The degree of mixing ${\overline{\chi }}_p$ decreases in the early stage but increases in the later stage. (ii) For the cases of strong magnetic fields, the development of KHI is significantly inhibited, and the mixing degree is greatly reduced. (iii) With the increase of Knudsen number, the TNE strength, entropy production rate and mixing degree all decrease.

The following two subsections continue to introduce the application of TNE effects in probing complex systems. Just as a fingerprint is a person’s identity. TNE feature is the identity feature of a system or behavior. The techniques described below can be used to study various complex flows, not just the types listed above.

5.13 Recovery of actual distribution function

The TNE feature of a flow state is its identity. The kinetic theory uses the distribution function $f$ to describe the system state, while the discrete distribution function $f_i$ itself does not have a clear physical meaning in LBM and unsteady DBM. Specifically, $f_i$ is by no means the actual distribution function of velocity ${\mathit{v}}_i$. What physically meaningful are the kinetic moments, such as ${\mathit{M}}_n^{\ast }\left(f_i\right)$ and ${\mathbf{\Delta }}_n^{\ast }={\mathit{M}}_n^{\ast }\left(f_i-f_i^{eq}\right)$. The equilibrium distribution function $f^{eq}$ can be obtained from the density $\rho$, velocity $\mathit{v}$ and temperature $T$. Therefore, DBM can use $f^{eq}$ and ${\mathbf{\Delta }}_n^{\ast }$ to recover the main features of the actual distribution function $f$ of the target region [121, 122]. Fig.47 shows the two-dimensional contours of the actual distribution function in velocity space ($v_r$,$v_{\theta }$) at the rarefaction wave, the material interface and the shock front, respectively, and Fig.48 displays the three-dimensional contours of actual distribution function at the corresponding positions. It can be seen that the distribution functions of different interfaces show different morphological appearances. DBM simulations and theoretical derivations allow for the representation of the morphological features of the real distribution functions. In literatures [122-124], an alternative method for recovering the main characteristics of the velocity distribution quantitatively is presented and developed through the macroscopic quantities and their spatial and temporal derivatives.

Fig.47 The sketch of the two-dimensional contours of actual distribution function in velocity space ($v_r$,$v_{\theta }$) at rarefaction front, material interface, and shock front, respectively. Reproduced with permission from Ref. [121].

Full size|PPT slide

Fig.48 The three-dimensional contour graph of the actual distribution function in velocity space ($v_r$,$v_{\theta }$) at rarefaction front, material interfaces, and shock fronts, respectively. Reproduced with permission from Ref. [121].

Full size|PPT slide

5.14 Identification of various interfaces

There generally more than one kinds of interfaces in complex flows. The TNE feature of an interface is its identity. By analyzing the TNE effect, the finer physical structure of the interface can be obtained, which can be used for the design of physical identification and tracking technology of different interfaces in the system. Fig.49 shows the characteristics of non-equilibrium quantities at the rarefaction wave, the material interface, and the shock front when the shock wave propagates from the heavy fluid to the light fluid. These distinct characteristics of these quantities can be used to discriminate different types of interfaces [121]. Fig.50 shows the contour of TNE quantity at the interface of the RTI system [125]. It is evident that TNE quantity is highly correlated with the macroscopic quantity gradient at the interface. Therefore, the distribution of TNE can well represent the characteristics of the interface and then can be used to identify the interface. As Fig.51 shows, the evolution of interface amplitude captured with the ${\mathrm{\Delta }}_{(3,1)y}^{\ast }$ matches the interface captured with the average temperature [126].

Fig.49 Profiles of four types of TNE quantities. Each column of the figure represent different TNE quantities, and each row indicate the case with various positions. Reproduced with permission from Ref. [121].

Full size|PPT slide

Fig.50 Contours of various TNE quantities at a typical moment during the RTI evolution. Reproduced with permission from Ref. [125].

Full size|PPT slide

Fig.51 Evolution curve of interface amplitude obtained from two kinds of tracking techniques. Reproduced with permission from Ref. [126].

Full size|PPT slide

6 Summary and prospect

The research content presented in this article is primarily motivated by two sets of scientific problems. These two sets of scientific problems can also be said to be “two dark clouds” faced in this cross-over field at that time. Each dark cloud is composed of a set of scientific problems. The first set of scientific questions are related to the LBM method: what exactly can LBM do as a system of method other than restore known governing equations (or as a revolutionary discrete format for known governing equations)? Can it be really beyond traditional fluid theory? If so, where and how is it beyond? If still not, in which way can it be hopeful to be truly beyond? The second set of scientific problems, independent of LBM and originated from the study of complex medium kinetics, correspond to the set mentioned in Section 2.1. Thinking about the first set of scientific problems becomes one of the methodological backgrounds and basis for the authors to seek solutions to the second set of scientific problems.

Compared to macroscopic (large-scale and slow-changing) behavior, the description of mesoscale behavior requires more physical quantities. The typical characteristics of mesoscale behavior are significant discrete effects and significant thermodynamic non-equilibrium effects. Therefore, before simulation, the ability to describe these discrete effects and non-equilibrium effects is the basic requirement for selecting/developing physical models. After simulation, it is the basic requirement of the analysis method to construct appropriate physical quantities to detect and identify these discrete effects and non-equilibrium effects, and then to give an intuitive geometric image correspondence. The discrete Boltzmann method is a physical model construction method and complex physical field analysis method developed based on the discrete Boltzmann equation for a series of fundamental scientific problems involved in the study of mesoscale behavior. It is the specific application and further development of coarse-grained description method, non-equilibrium behavior description method, and phase space description method in statistical physics under the framework of discrete Boltzmann equation. From a historical perspective, DBM was developed from the physical modeling branch of LBM, undergoing a process of taking, abandoning and adding. It aims at the two major problems of how to model and how to analyze in numerical experimental research. It no longer relies on the continuity assumption and near equilibrium approximation of traditional fluid modeling, and it abandons the “lattice gas” images of standard LBM. It adds detection, presentation, description, and analysis solutions based on phase space for discrete/TNE states and effects. Over time and according to the research needs, more information extraction techniques and complex physical field analysis schemes will be introduced.

Although the flow behavior targeted is complex, the DBM approach is inherently simple. As far as the modeling method is concerned, it only needs to check the kinetic moments of $f^{(0)}$ involved in the more accurate constitutive relationships according to the orders of tensors, and determine them as the kinetic moments that needs to preserve value, instead of deriving extremely complex extended hydrodynamic equations like the KMM approach. As far as the analysis methods are concerned, they are simply a concrete application and development of the basic principles of statistical physics. They have clear physical images, are easy to understand and convenient to use. The model equations of DBM include the system behavior evolution equation(s) and discrete velocity constraint equations, where the evolution equations, in addition to discrete Boltzmann equation, depending on the specific situation, may also include phase field evolution equation, chemical reaction evolution equation, electromagnetic field evolution equation, etc. With the model equations, if further with the concrete discrete scheme, numerical simulation can be carried out. With the model equations, if without specific discrete format, the Chapman−Enskog multi-scale analysis and subsequent analytical studies can be performed.

The “mesoscopic” characteristics of DBM can be summarized as follows: (i) In terms of scale, it connects micro and macro scales. (ii) In terms of functionality, it is beyond the NS but is weaker than molecular dynamics in describing the discrete behavior and non-equilibrium behavior. It should be pointed out that, compared with traditional fluid modeling, what the kinetic modeling brings are the improvement of understanding level, the broadening of application range, and the supplement of technology, but not a replacement. For the cases where the conventional fluid modeling is convenient, effective, and adequate, it remains the first choice for most practical applications. Because the associated theory, method, software and so on are more mature. The user’s own convenience is also an important factor in decision-making process.

The mesoscale behaviors promote the emergence of new technologies. The studies on mesoscale behavior have made gratifying progress, but there is still a long way to go. Fortunately, experimental studies on non-equilibrium flow are attracting increasing attention [127-134]. Due to various limitations of measurements for non-equilibrium effects, progress in experimental study has been relatively slower. This is precisely why numerical simulation studies have surged ahead of experiments, and it is also why numerical simulation should be accorded priority in many studies of non-equilibrium flow. Some practical applications of discrete/non-equilibrium effects have first been realized in numerical experimental studies. There will be more of these discoveries and applications over time. Numerical research provides cognitive basis for experimental research and engineering applications, which is the value of numerical simulation.

In addition to the relative discreteness caused by rarefaction or small structure and the fast flow in the case of no reaction or slow reaction, both the real gas effect at high temperature and fast chemical reaction can also cause the thermodynamic non-equilibrium flow. In addition to the further study of discrete and non-equilibrium effects, DBM modeling in the latter two cases and related non-equilibrium regulation techniques are meaningful research directions that need to be explored in the near future.

Astrophysics is one of the important application fields of fluid mechanics. Since many of the extreme conditions, behaviors, or physical images of high energy density physics techniques like ICF can be found in astrophysics, in addition to pure astrophysicists, scientists in fields related to high energy density physics techniques (such as ICF-related scientists) are often interested in astrophysics related research. Using fluid mechanics to model astrophysical processes has many specific applications, such as the evolution of stars, planets, accretion disks, and galaxies. Star behaviors are divided into atmospheric and internal convection regions, and accretion disk behaviors are divided into protostellar disks, black hole disks, binary disks, and so on. At present, the fluid mechanics simulations in astrophysics are mainly the NS simulation, Smoothed Particle Hydrodynamics (SPH) simulation and Particle In Cell (PIC) simulation. SPH and PIC simulations have a certain degree of kinetics simulation characteristics, but mainly focus on physics images of NS. The NS continuous description corresponds to the fact that we are far enough away from that celestial structure so that the field of view is large enough, and consequently the distance between the stars is negligible relative to the scale of our field of view. As our spacecraft’s observing instruments get closer and the scale of the field of view decreases, we will gradually observe larger and larger gaps and more inhomogeneity between the stars including different sizes of the stars, different translational and rotational speeds, etc. The existence of gravitation, etc., makes these asymmetries often imply non-equilibrium. Discreteness and non-equilibrium begin to be observed, which is the benefit of improved observation accuracy. The gradual approach of our spacecraft’s observing instruments corresponds to a gradual increase in the Knudsen number, because the scale corresponding to our observing field of view is gradually decreasing. DBM simulations can show changes in the view as our spacecraft’s observing instruments get closer. Of course, it is presented in the form of statistical features to the extent that DBM theory is valid.

Dust gas kinetics and granular flow kinetics have important engineering application background and prospect. “Particles” in dust gases and granular flows are “macroscopic” particles made up of a large number of molecules. Collisions between such particles are not elastic, but dissipative. Therefore, dust gas kinetics and granular flow kinetics are more complex than the usual complex fluids. The study of this kind of complex fluid system is also divided into macroscopic, mesoscopic and microscopic levels, facing the meso-scale dilemma of “insufficient physical function of traditional fluid modeling and impossibility of molecular dynamics simulation in terms of scale”. As long as the generalized BGK-like Boltzmann equation including particle collision dissipation effect can be obtained, the DBM modeling and analysis method introduced in this review can be generalized and applied. The theory of dust kinetics and granular flow kinetics, although more complex and beyond the discussion scope of the current review, has a history of nearly 30 years. Interested readers are referred to more specialized literature. Everything flows. As long as traditional fluid mechanics has applications, there will be kinetic behaviors that traditional fluid mechanics fails to well describe, and DBM can find a reasonable application. Although the former can also be studied using the long-established NS simulation, DBM can achieve a certain degree of cross-scale simulation under the same framework, thereby avoiding the problems that arise when different scale methods are interfaced.

The previous sections of this article do not cover the discussion of specific discrete format, because that is another field. Finally, just present a brief comment from the perspective of application. Any form of discretization of an originally continuous space will inevitably bring about symmetry breaking. In numerical simulation, symmetry breakings from different dimensions are coupled together to make effects. In the study of specific discrete format, not only the numerical accuracy but also the numerical stability should be guaranteed. In many practical applications, numerical stability is more important than numerical accuracy. The numerical stability depends not only on the discrete scheme, but also on the specific flow behavior. Because the traditional fluid modeling works in time−position space, the specific discrete format of time−position space has been deeply and systematically studied in the traditional computational fluid dynamics. With the emergence of kinetic methods using discrete particle velocities such as lattice gas method, LBM and DBM, the study of stable discrete format in time-position-velocity space has become an urgent research topic in computational mathematics. For more details on discrete format study, please refer to the more specialized literature.

References

[1] H. S. Tsien. Superaerodynamics, mechanics of rarefied gases. J. Aeronaut. Sci., 1946, 13(12): 653

[2] W. Thomson. Hydrokinetic solutions and observations. Lond. Edinb. Dublin Philos. Mag. J. Sci., 1871, 42(281): 362

[3] H. L. F. Helmholtz. On discontinuous movements of fluids. Monatsberichte der Königlichen Preussische Akademie der Wissenschaften zu Berlin, 1868, 36: 337

[4] R.Rayleigh, Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density, Proceedings of the London Mathematical Society s1‒14, pp 170–177 (1882)

[5] G. I. Taylor. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes (I). Proc. R. Soc. Lond. A, 1950, 201(1065): 192

[6] R. D. Richtmyer. Taylor instability in shock acceleration of compressible fluids. Commun. Pure Appl. Math., 1960, 13(2): 297

[7] E. E. Meshkov. Instability of the interface of two gases accelerated by a shock wave. Fluid Dyn., 1972, 4(5): 101

[8] Y. X. Liu, L. F. Wang, K. G. Zhao, Z. Y. Li, J. F. Wu, W. H. Ye, Y. J. Li. Thin-shell effects on nonlinear bubble evolution in the ablative Rayleigh–Taylor instability. Phys. Rev. Lett., 2022, 29(8): 082102

[9] Y. X. Liu, Z. Chen, L. F. Wang, Z. Y. Li, J. F. Wu, W. H. Ye, Y. J. Li. Dynamic of shock–bubble interactions and nonlinear evolution of ablative hydrodynamic instabilities initialed by capsule interior isolated defects. Phys. Plasmas, 2023, 30(4): 042302

[10] K. Lan. Dream fusion in octahedral spherical hohlraum. Matter Radiat. Extrem., 2022, 7(5): 055701

[11] Y. H. Chen, Z. C. Li, H. Cao, K. Q. Pan, S. W. Li, X. F. Xie, B. Deng, Q. Q. Wang, Z. R. Cao, L. F. Hou, X. S. Che, P. Yang, Y. J. Li, X. A. He, T. Xu, Y. G. Liu, Y. L. Li, X. M. Liu, H. J. Zhang, W. Zhang, B. L. Jiang, J. Xie, W. Zhou, X. X. Huang, W. Y. Huo, G. L. Ren, K. Li, X. D. Hang, S. Li, C. L. Zhai, J. Liu, S. Y. Zou, Y. K. Ding, K. Lan. Determination of laser entrance hole size for ignition-scale octahedral spherical hohlraums. Matter Radiat. Extrem., 2022, 7(6): 065901

[12] K. Lan, Y. Dong, J. Wu, Z. Li, Y. Chen, H. Cao, L. Hao, S. Li, G. Ren, W. Jiang, C. Yin, C. Sun, Z. Chen, T. Huang, X. Xie, S. Li, W. Miao, X. Hu, Q. Tang, Z. Song, J. Chen, Y. Xiao, X. Che, B. Deng, Q. Wang, K. Deng, Z. Cao, X. Peng, X. Liu, X. He, J. Yan, Y. Pu, S. Tu, Y. Yuan, B. Yu, F. Wang, J. Yang, S. Jiang, L. Gao, J. Xie, W. Zhang, Y. Liu, Z. Zhang, H. Zhang, Z. He, K. Du, L. Wang, X. Chen, W. Zhou, X. Huang, H. Guo, K. Zheng, Q. Zhu, W. Zheng, W. Y. Huo, X. Hang, K. Li, C. Zhai, H. Xie, L. Li, J. Liu, Y. Ding, W. Zhang. First inertial confinement fusion implosion experiment in octahedral spherical Hohlraum. Phys. Rev. Lett., 2021, 127(24): 245001

[13] X. M. Qiao, K. Lan. Novel target designs to mitigate hydrodynamic instabilities growth in inertial confinement fusion. Phys. Rev. Lett., 2021, 126(18): 185001

[14] Y. B. Gan, A. G. Xu, G. C. Zhang, Y. J. Li. Lattice Boltzmann study on Kelvin–Helmholtz instability: Roles of velocity and density gradients. Phys. Rev. E, 2011, 83(5): 056704

[15] C. D. Lin, A. G. Xu, G. C. Zhang, K. H. Luo, Y. J. Li. Discrete Boltzmann modeling of Rayleigh–Taylor instability in two-component compressible flows. Phys. Rev. E, 2017, 96(5): 053305

[16] F. Chen, A. G. Xu, G. C. Zhang. Collaboration and competition between Richtmyer–Meshkov instability and Rayleigh–Taylor instability. Phys. Fluids, 2018, 30(10): 102105

[17] F. Chen, A. G. Xu, Y. D. Zhang, Y. B. Gan, B. B. Liu, S. Wang. Effects of the initial perturbations on the Rayleigh–Taylor–Kelvin–Helmholtz instability system. Front. Phys., 2022, 17(3): 33505

[18] L. Chen, H. L. Lai, C. D. Lin, D. M. Li. Specific heat ratio effects of compressible Rayleigh–Taylor instability studied by discrete Boltzmann method. Front. Phys., 2021, 16(5): 52500

[19] G. Zhang, A. G. Xu, D. J. Zhang, Y. J. Li, H. L. Lai, X. M. Hu. Delineation of the flow and mixing induced by Rayleigh–Taylor instability through tracers. Phys. Fluids, 2021, 33(7): 076105

[20] Y. F. Li, H. L. Lai, C. D. Lin, D. M. Li. Influence of the tangential velocity on the compressible Kelvin–Helmholtz instability with non-equilibrium effects. Front. Phys., 2022, 17(6): 63500

[21] H. Lai, C. Lin, Y. Gan, D. Li, L. Chen. The influences of acceleration on compressible Rayleigh–Taylor instability with non-equilibrium effects. Comput. Fluids, 2023, 266(15): 106037

[22] C. D. Lin, K. H. Luo, Y. B. Gan, Z. P. Liu. Kinetic simulation of nonequilibrium Kelvin–Helmholtz instability. Commum. Theor. Phys., 2019, 71(1): 132

[23] C. D. Lin, K. H. Luo, A. G. Xu, Y. B. Gan, H. L. Lai. Multiple-relaxation-time discrete Boltzmann modeling of multicomponent mixture with non-equilibrium effects. Phys. Rev. E, 2021, 103(1): 013305

[24] Y. B. Gan, A. G. Xu, G. C. Zhang, C. D. Lin, H. L. Lai, Z. P. Liu. Nonequilibrium and morphological characterizations of Kelvin–Helmholtz instability in compressible flows. Front. Phys., 2019, 14(4): 43602

[25] W.F. ChenW. W. Zhao, Moment Equations and Numerical Methods for Rarefied Gas Flows, Beijing: Science Press, 2017 (in Chinese)

[26] Z. H. Li, A. P. Peng, H. X. Zhang, X. G. Deng. Numerical study on the gas-kinetic high-order schemes for solving Boltzmann model equation. Sci. China Phys. Mech. Astron., 2011, 54(9): 1687

[27] G.KarniadaskisA.BeskokN.Aluru, Microflows and Nanoflows: Fundamentals and Simulation, New York: Springer-Verlag, 2005

[28] A. Keerthi, A. K. Geim, A. Janardanan, A. P. Rooney, A. Esfandiar, S. Hu, S. A. Dar, I. V. Grigorieva, S. J. Haigh, F. C. Wang, B. Radha. Ballistic molecular transport through two-dimensional channels. Nature, 2018, 558(7710): 420

[29] G. López Quesada, G. Tatsios, D. Valougeorgis, M. Rojas-Cárdenas, L. Baldas, C. Barrot, S. Colin. Design guidelines for thermally driven micropumps of different architectures based on target applications via kinetic modeling and simulations. Micromachines (Basel), 2019, 10(4): 249

[30] N. Kavokine, R. R. Netz, L. Bocquet. Fluids at the nanoscale: From continuum to subcontinuum transport. Annu. Rev. Fluid Mech., 2021, 53(1): 377

[31] X. Jiang, G. A. Siamas, K. Jagus, T. G. Karayiannis. Physical modelling and advanced simulations of gas–liquid two-phase jet flows in atomization and sprays. Prog. Eenerg. combust., 2010, 36(2): 131

[32] J. C. Ding, T. Si, M. J. Chen, Z. G. Zhai, X. Y. Lu, X. S. Luo. On the interaction of a planar shock with a three-dimensional light gas cylinder. J. Fluid Mech., 2017, 828: 289

[33] Y. Liang, Z. G. Zhai, X. S. Luo. Interaction of strong converging shock wave with SF₆ gas bubble. Sci. China Phys. Mech. Astron., 2018, 61(6): 064711

[34] J. C. Ding, Y. Liang, M. J. Chen, Z. G. Zhai, T. Si, X. S. Luo. Interaction of planar shock wave with three-dimensional heavy cylindrical bubble. Phys. Fluids, 2018, 30(10): 106109

[35] X. S. Luo, M. Li, J. C. Ding, Z. G. Zhai, T. Si. Nonlinear behaviour of convergent Richtmyer–Meshkov instability. J. Fluid Mech., 2019, 877: 130

[36] Y. D. Zhang, A. G. Xu, J. J. Qiu, H. T. Wei, Z. H. Wei. Kinetic modeling of multiphase flow based on simplified Enskog equation. Front. Phys., 2020, 15(6): 62503

[37] Z.L. JiangH. H. Teng, Gaseous Detonation Physics and Its Universal Framework Theory, Singapore: Springer, 2022

[38] J.P. WangS. B. Yao, Principle and Technology of Continuous Detonation Engine, Beijing: Science Press, 2018 (in Chinese)

[39] J. Xiao, P. Liu, C. X. Wang, G. W. Yang. External field-assisted laser ablation in liquid: An efficient strategy for nanocrystal synthesis and nanostructure assembly. Prog. Mater. Sci., 2017, 87: 140

[40] H. G. Rinderknecht, P. A. Amendt, S. C. Wilks, G. Collins. Kinetic physics in ICF: Present understanding and future directions. Plasma Phys. Contr. Fusion, 2018, 60(6): 064001

[41] F. Vidal, J. P. Matte, M. Casanova, O. Larroche. Ion kinetic simulations of the formation and propagation of a planar collisional shock wave in a plasma. Phys. Fluids B Plasma Phys., 1993, 5(9): 3182

[42] D. Bond, V. Wheatley, R. Samtaney, D. I. Pullin. Richtmyer–Meshkov instability of a thermal interface in a two-fluid plasma. J. Fluid Mech., 2017, 833: 332

[43] J.H. SongA. G. XuL.MiaoF.ChenZ.P. Liu L.F. WangN. F. WangX.Hou, Plasma kinetics: Discrete Boltzmann modelling and Richtmyer–Meshkov instability, Phys. Fluids 36, 016107(2023)

[44] P. L. Yao, H. B. Cai, X. X. Yan, W. S. Zhang, B. Du, J. M. Tian, E. H. Zhang, X. W. Wang, S. P. Zhu. Kinetic study of transverse electron-scale interface instability in relativistic shear flows. Matter Radiat. Extrem., 2020, 5(5): 054403

[45] H. B. Cai, X. X. Yan, P. L. Yao, S. P. Zhu. Hybrid fluid–particle modeling of shock-driven hydrodynamic instabilities in a plasma. Matter Radiat. Extrem., 2021, 6(3): 035901

[46] R. K. Agarwal, K. Y. Yun, R. Balakrishnan. Beyond Navier–Stokes: Burnett equations for flows in the continuum-transition regime. Phys. Fluids, 2001, 13(10): 3061

[47] P.K. KunduI. M. CohenD.Dowling, Fluid Mechanics, 4th Ed., Ch. 13, pp 537–601, Oxford: Elsevier Academic Press, 2008

[48] H. J. Zhou, Y. Zhang, Z. C. Ouyang. Bending and base-stacking interactions in double-stranded DNA. Phys. Rev. Lett., 1999, 82(22): 4560

[49] X. Y. Zhang, S. Boccaletti, S. G. Guan, Z. H. Liu. Explosive synchronization in adaptive and multilayer networks. Phys. Rev. Lett., 2015, 114(3): 038701

[50] G. H. Tang, C. Bi, Y. Zhao, W. Q. Tao. Thermal transport in nano-porous insulation of aerogel: Factors, models and outlook. Energy, 2015, 90: 701

[51] Z. C. Zong, S. C. Deng, Y. J. Qin, X. Wan, J. H. Zhan, D. K. Ma, N. Yang. Enhancing the interfacial thermal conductance of Si/PVDF by strengthening atomic couplings. Nanoscale, 2023, 15(40): 16472

[52] L. N. Yang, B. S. Yang, B. W. Li. Enhancing interfacial thermal conductance of an amorphous interface by optimizing the interfacial mass distribution. Phys. Rev. B, 2023, 108(16): 165303

[53] H. Zhao. Identifying diffusion processes in one-dimensional lattices in thermal equilibrium. Phys. Rev. Lett., 2006, 96(14): 140602

[54] Y. Li, X. Y. Shen, Z. H. Wu, J. Y. Huang, Y. X. Chen, Y. S. Ni, J. P. Huang. Temperature-dependent transformation thermotics: From switchable thermal cloaks to macroscopic thermal diodes. Phys. Rev. Lett., 2015, 115(19): 195503

[55] Z. Wang, W. C. Fu, Y. Zhang, H. Zhao. Wave-turbulence origin of the instability of anderson localization against many-body interactions. Phys. Rev. Lett., 2020, 124(18): 186401

[56] N. B. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, B. W. Li. Phononics: Manipulating heat flow with electronic analogs and beyond. Rev. Mod. Phys., 2012, 84(3): 1045

[57] L. Wang, D. H. He, B. B. Hu. Heat conduction in a three-dimensional momentum-conserving anharmonic lattice. Phys. Rev. Lett., 2010, 105(16): 160601

[58] L. Zhao, C. L. Wang, J. Liu, B. H. Wen, Y. S. Tu, Z. W. Wang, H. P. Fang. Reversible state transition in nanoconfined aqueous solutions. Phys. Rev. Lett., 2014, 112(7): 078301

[59] I. Maasilta, A. J. Minnich. Heat under the microscope. Phys. Today, 2014, 67(8): 27

[60] S.G. Chen, Nonequilibrium Statistical Mechanics, Beijing: Science Press, 2010 (in Chinese)

[61] A.G. XuY. D. Zhang, Complex Media Kinetics, Beijing: Science Press, 2022 (in Chinese)

[62] J.H. LiG. WeiW.WangN.YangX.H. Liu L.M. WangX. F. HeX.W. WangJ.W. WangM.Kwauk, From Multiscale Modeling to Meso-Science: A Chemical Engineering Perspective, Berlin: Springer-Verlag, 2013

[63] W. L. Huang, J. H. Li, P. P. Edwards. Mesoscience: exploring the common principle at mesoscales. Natl. Sci. Rev., 2018, 5(3): 321

[64] Y. D. Zhang, X. Wu, B. B. Nie, A. G. Xu, F. Chen, R. Wei. Lagrangian steady-state discrete Boltzmann model for non-equilibrium flows at micro–nanoscale. Phys. Fluids, 2023, 35(9): 092008

[65] Y. B. Gan, A. G. Xu, H. L. Lai, W. Li, G. L. Sun, S. Succi. Discrete Boltzmann multi-scale modelling of non-equilibrium multiphase flows. J. Fluid Mech., 2022, 951: A8

[66] C. L. Tien. Molecular and microscale transport phenomena: A report on the 2nd US Japan Joint Seminar, Santa Barbara, California, 7−10 August, 1996. Microscale Thermophys. Eng., 1997, 1(1): 71

[67] B.M. McCoy, Advanced Statistical Mechanics, Oxford: Oxford university press, 2010

[68] T. D. Lee, C. N. Yang. Statistical theory of equations of state and phase transitions. II. Lattice gas and Ising model. Phys. Rev., 1952, 87(3): 410

[69] S.Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford: Oxford University Press, 2001

[70] Y.L. HeY. WangQ.Li, Lattice Boltzmann Method: Theory and Applications, Beijing: Science Press, 2009 (in Chinese)

[71] Z.L. GuoC. Shu, Lattice Boltzmann Method and its Applications in Engineering, Beijing: World Scientific Publishing, 2013

[72] H.B. HuangM. C. SukopX.Y. Lu, Multiphase Lattice Boltzmann Methods: Theory and Application, John Wiley & Sons, 2015

[73] A. G. Xu, G. C. Zhang, Y. B. Gan, F. Chen, X. J. Yu. Lattice Boltzmann modeling and simulation of compressible flows. Front. Phys., 2012, 7(5): 582

[74] A. G. Xu, C. D. Lin, G. C. Zhang, Y. J. Li. Multiple relaxation-time lattice Boltzmann kinetic model for combustion. Phys. Rev. E, 2015, 91(4): 043306

[75] A.G. XuG. C. ZhangY.D. ZhangY.B. Gan, Discrete Boltzmann modeling of nonequilibrium effects in multiphase flow, Presentation at the 31st International Symposium on Rarefied Gas Dynamics, see also FLOWS: Physics & Beyond, 1001 (2018)

[76] A.G. XuJ. H. SongF.ChenK.XieY.J. Ying, Modeling and analysis methods for complex fields based on phase space, Chinese Journal of Computational Physics 38(6), 631 (2021) (in Chinese)

[77] D. J. Zhang, A. G. Xu, Y. D. Zhang, Y. B. Gan, Y. J. Li. Discrete Boltzmann modeling of high-speed compressible flows with various depths of non-equilibrium. Phys. Fluids, 2022, 34(8): 086104

[78] Y.D. Zhang, Modeling and research of non-equilibrium flows and multi-phase flows: Based on discrete Boltzmann method, Nanjing: Nanjing University of Science & Technology, 2019 (in Chinese)

[79] Y. B. Gan, A. G. Xu, G. C. Zhang, Y. D. Zhang, S. Succi. Discrete Boltzmann trans-scale modeling of highspeed compressible flows. Phys. Rev. E, 2018, 97(5): 053312

[80] L. Wu, J. M. Reese, Y. H. Zhang. Solving the Boltzmann equation deterministically by the fast spectral method: Application to gas microflows. J. Fluid Mech., 2014, 746: 53

[81] B. Yan, A. G. Xu, G. C. Zhang, Y. J. Ying, H. Li. Lattice Boltzmann model for combustion and detonation. Front. Phys., 2013, 8: 94

[82] C. D. Lin, A. G. Xu, G. C. Zhang, Y. J. Li. Double distribution-function discrete Boltzmann model for combustion. Combust. Flame, 2016, 164: 137

[83] C. D. Lin, K. H. Luo. Mesoscopic simulation of nonequilibrium detonation with discrete Boltzmann method. Combust. Flame, 2018, 198: 356

[84] C. D. Lin, A. G. Xu, G. C. Zhang, Y. J. Li. Polar coordinate lattice Boltzmann kinetic modeling of detonation phenomena. Commum. Theor. Phys., 2014, 62(5): 737

[85] C.D. LinK. H. LuoL.L. FeiS.Succi, A multicomponent discrete Boltzmann model for nonequilibrium reactive fows, Commum. Theor. Phys. 7(1), 14580 (2017)

[86] C.D. LinK. H. Luo, MRT discrete Boltzmann method for compressible exothermic reactive flows, Comput. Fluids 166, 176 (2018)

[87] C. D. Lin, K. H. Luo. Discrete Boltzmann modeling of unsteady reactive flows with nonequilibrium effects. Phys. Rev. E, 2019, 99(1): 012142

[88] C. D. Lin, K. H. Luo. Kinetic simulation of unsteady detonation with thermodynamic nonequilibrium effects. Combust. Explos. Shock Waves, 2020, 56(4): 435

[89] C. D. Lin, X. L. Su, Y. D. Zhang. Hydrodynamic and thermodynamic nonequilibrium effects around shock waves: Based on a discrete Boltzmann method. Entropy (Basel), 2020, 22(12): 1397

[90] Y. D. Zhang, A. G. Xu, G. C. Zhang, C. M. Zhu, C. D. Lin. Kinetic modeling of detonation and effects of negative temperature coefficient. Combust. Flame, 2016, 173: 483

[91] Y. M. Shan, A. G. Xu, Y. D. Zhang, L. F. Wang, F. Chen. Discrete Boltzmann modeling of detonation: Based on the Shakhov model. Proc. Inst. Mech. Eng., C J. Mech. Eng. Sci., 2023, 237(11): 2517

[92] Y. D. Zhang, A. G. Xu, G. C. Zhang, Z. H. Chen. Discrete Boltzmann method with Maxwell-type boundary condition for slip flow. Commum. Theor. Phys., 2018, 69(1): 77

[93] Y. D. Zhang, A. G. Xu, G. C. Zhang, Z. H. Chen, P. Wang. Discrete Boltzmann method for non-equilibrium flows: Based on Shakhov model. Comput. Phys. Commun., 2019, 238: 50

[94] Y. Zhang, A. Xu, F. Chen, C. Lin, Z. H. Wei. Non-equilibrium characteristics of mass and heat transfers in the slip flow. AIP Adv., 2022, 12(3): 035347

[95] Y.Sone, Molecular Gas Dynamics: Theory, Techniques, and Applications, Springer Science & Business Media, 2007

[96] Y.Onishi, A rarefied gas flow over a flat wall, Bull. Univ. Osaka Prefect. Ser. A Eng. Nat. Sci. 22(2), 91 (1974)

[97] X. S. Chen, V. Dohm, N. Schultka. Order-parameter distribution function of finite O(n) symmetric systems. Phys. Rev. Lett., 1996, 77(17): 3641

[98] A. J. Wagner, J. M. Yeomans. Breakdown of scale invariance in the coarsening of phase-separating binary fluids. Phys. Rev. Lett., 1998, 80(7): 1429

[99] M. R. Swift, W. R. Osborn, J. M. Yeomans. Lattice Boltzmann simulation of nonideal fluids. Phys. Rev. Lett., 1995, 75(5): 830

[100] Y. M. Shan, A. G. Xu, L. F. Wang, Y. D. Zhang. Nonequilibrium kinetics effects in Richtmyer–Meshkov instability and reshock processes. Commum. Theor. Phys., 2023, 75(11): 115601

[101] M. Latini, O. Schilling, W. S. Don. High-resolution simulations and modeling of reshocked single-mode Richtmyer–Meshkov instability: Comparison to experimental data and to amplitude growth model predictions. Phys. Fluids, 2007, 19(2): 024104

[102] B.D. CollinsJ.W. Jacobs, PLIF flow visualization and measurements of the Richtmyer–Meshkov instability of an air/SF6 interface, J. Fluid Mech. 464, 113 (2002)

[103] D. J. Zhang, A. G. Xu, J. H. Song, Y. B. Gan, Y. D. Zhang, Y. J. Li. Specific-heat ratio effects on the interaction between shock wave and heavy-cylindrical bubble: Based on discrete Boltzmann method. Comput. Fluids, 2023, 265: 106021

[104] G.S. JiangC. L. Wu, A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics, J. Comput. Phys. 150(2), 561 (1999)

[105] J.H. SongA. G. XuL.MiaoY.G. LiaoF.W. Liang F.TianM. Q. NieN.F. Wang, Entropy increase characteristics of shock wave/plate laminar boundary layer interaction, Acta Aeronautica et Astronautica Sinica 44(21), 528520 (2023) (in Chinese)

[106] S. F. Liao, W. B. Zhang, H. Chen, L. Y. Zou, J. H. Liu, X. X. Zheng. Atwood number effects on the instability of a uniform interface driven by a perturbed shock wave. Phys. Rev. E, 2019, 99(1): 013103

[107] L. Y. Zou, M. Al-Marouf, W. Cheng, R. Samtaney, J. C. Ding, X. S. Luo. Richtmyer–Meshkov instability of an unperturbed interface subjected to a diffracted convergent shock. J. Fluid Mech., 2019, 879: 448

[108] L. Y. Zou, Q. Wu, X. Z. Li. Research progress of general Richtmyer–Meshkov instability. Sci. Sin. Phys. Mech. Astron., 2020, 50: 104702

[109] Y.M. ShanA. G. XuY.D. ZhangL.F. Wang Wall-heating phenomena in shock wave physics: Physical or artificial? (in preparation)

[110] D. J. Zhang, A. G. Xu, Y. B. Gan, Y. D. Zhang, J. H. Song, Y. J. Li. Viscous effects on morphological and thermodynamic non-equilibrium characterizations of shock-bubble interaction. Phys. Fluids, 2023, 35(10): 106113

[111] Y. B. Gan, A. G. Xu, G. C. Zhang, S. Succi. Discrete Boltzmann modeling of multiphase flows: Hydrodynamic and thermodynamic non-equilibrium effects. Soft Matter, 2015, 11(26): 5336

[112] Y. D. Zhang, A. G. Xu, G. C. Zhang, Y. B. Gan, Z. H. Chen, S. Succi. Entropy production in thermal phase separation: A kinetic-theory approach. Soft Matter, 2019, 15(10): 2245

[113] G. L. Sun, Y. B. Gan, A. G. Xu, Y. D. Zhang, Q. F. Shi. Thermodynamic nonequilibrium effects in bubble coalescence: A discrete Boltzmann study. Phys. Rev. E, 2022, 106(3): 035101

[114] G.L. SunY. B. GanA.G. XuQ.F. Shi, Droplet coalescence kinetics: Thermodynamic nonequilibrium effects and entropy production mechanism, arXiv: 2311.06546 (2023), Phys. Fluids (2024) (in press)

[115] J. Chen, A. G. Xu, D. W. Chen, Y. D. Zhang, Z. H. Chen. Discrete Boltzmann modeling of Rayleigh–Taylor instability: Effects of interfacial tension, viscosity, and heat conductivity. Phys. Rev. E, 2022, 106(1): 015102

[116] Y.Q. JiaA. G. XuL.F. WangD.W. Chen, Study on Rayleigh−Taylor instability under variable acceleration: Based on discrete Boltzmann method, Beijing: Symposium on Interfacial Instability and Multi-media Turbulence, 2023 (in Chinese)

[117] J.ChenA. G. XuD.W. ChenY.D. ZhangZ.H. Chen, Kinetics of RT instability in van der Waals fluid: The influence of compressibility (in preparation)

[118] F. Chen, A. G. Xu, Y. D. Zhang, Q. K. Zeng. Morphological and non-equilibrium analysis of coupled Rayleigh–Taylor–Kelvin–Helmholtz instability. Phys. Fluids, 2020, 32(10): 104111

[119] Z. P. Liu, J. H. Song, A. G. Xu, Y. D. Zhang, K. Xie. Discrete Boltzmann modeling of plasma shock wave. Proc. Inst. Mech. Eng., C J. Mech. Eng. Sci., 2023, 237(11): 2532

[120] J.H. SongL. MiaoA.G. XuF.ChenL.Li X.Hou, Plasma kinetics: Transport and mixing characteristics induced by Kelvin–Helmholtz instability (in preparation)

[121] C. D. Lin, A. G. Xu, G. C. Zhang, Y. J. Li, S. Succi. Polar-coordinate lattice Boltzmann modeling of compressible flows. Phys. Rev. E, 2014, 89(1): 013307

[122] Y. D. Zhang, A. G. Xu, G. C. Zhang, Z. H. Chen, P. Wang. Discrete ellipsoidal statistical BGK model and Burnett equations. Front. Phys., 2018, 13(3): 135101

[123] X. L. Su, C. D. Lin. Nonequilibrium effects of reactive flow based on gas kinetic theory. Commum. Theor. Phys., 2022, 74(3): 035604

[124] X. L. Su, C. D. Lin. Unsteady detonation with thermodynamic nonequilibrium effect based on the kinetic theory. Commum. Theor. Phys., 2023, 75(7): 075601

[125] H.L. Lai, Modeling and Simulation of Compressible Rayleigh-Taylor Instability by Discrete Boltzmann, Post-doctoral research report of Beijing Institute of Applied Physics and Computational Mathematics, 2015 (in Chinese)

[126] H. L. Lai, A. G. Xu, G. C. Zhang, Y. B. Gan, Y. J. Ying, S. Succi. Nonequilibrium thermohydrodynamic effects on the Rayleigh-Taylor instability in compressible flows. Phys. Rev. E, 2016, 94(2): 023106

[127] A.G. Chen, Study of shock detachment distance in rarefied flow field, Hefei: The 2nd National Conference on Shock Wave and Shock Tube, 2022 (in Chinese)

[128] W.M. Bao, Future aerospace technology development and aerodynamics challenges, Tianjin: The 2nd Aerodynamics Congress, 2023 (in Chinese)

[129] G.S. Zhu, Several encountered aerodynamic problems and progress in hypersonic flight, Tianjin: The 2nd Aerodynamics Congress, 2023 (in Chinese)

[130] B.C. Ai, Research on some problems of aerodynamics of space vehicle, Tianjin: The 2nd Aerodynamics Congress, 2023 (in Chinese)

[131] A.G. ChenJ. WangZ.H. LiZ.Q. LiY.Tian Z.Y. Long, Velocity measurement investigation of rarefied flow field by pulse electron beam fluorescence technique, Phys. Gases 6(5), 67 (2021) (in Chinese)

[132] W.ChenH. Y. HuL.Wang, Characterization method and measurement technique of thermodynamic non-equilibrium characteristics of high enthalpy flow field, Dalian: The 17th National Conference on Physics and Mechanics, 2023 (in Chinese)

[133] P.M. Danehy, Molecular based hypersonic nonequilibrium flow optical diagnosis (in Chinese, translated by Qifeng Chen) (to be published)

[134] P.M. DanehyB. F. BathelC.T. Johansen., . Molecular Based Hypersonic Non-Equilibrium Flow Optical Diagnosis, AIAA Progress Series, 2013

Declarations

The authors declare that they have no competing interests and there are no conflicts.

Acknowledgements

The authors thank Academician Shigang Chen for his academic guidance and training for so many years, thank Prof. Xing Wei et al. for their discussions on fluid mechanical simulation of astrophysical processes, and thank Academician Jinghai Li and Prof. Limin Wang et al. for their enlightening comments, discussions and suggestions. Dr. Jiahui Song organized some materials for the current review article. This work was supported by the National Natural Science Foundation of China (Grant Nos. 12172061 and 11875001), the Opening Project of State Key Laboratory of Explosion Science and Technology (Beijing Institute of Technology) (Grant No. KFJJ23-02M), the Foundation of National Key Laboratory of Shock Wave and Detonation Physics (Grant No. JCKYS2023212003), the Foundation of National Key Laboratory of Computational Physics, Hebei Outstanding Youth Science Foundation (Grant No. A2023409003), Hebei Natural Science Foundation (Grant No. A2021409001), the Central Guidance on Local Science and Technology Development Fund of Hebei Province (Grant No. 226Z7601G), and the “Three, Three and Three Talent Project” of Hebei Province (Grant No. A202105005).

RIGHTS & PERMISSIONS

2024 The Authors