Kinetic modeling of multiphase flow based on simplified Enskog equation
Yudong Zhang, Aiguo Xu, Jingjiang Qiu, Hongtao Wei, Zung-Hang Wei
aa r X i v : . [ phy s i c s . c o m p - ph ] S e p Kinetic modeling of multiphase flow based on simplified Enskogequation
Yudong Zhang , Aiguo Xu , ∗ , Jingjiang Qiu , Hongtao Wei , Zung-Hang Wei †
1, School of Mechanics and Safety Engineering,Zhengzhou University, Zhengzhou 450001, China2, Laboratory of Computational Physics,Institute of Applied Physics and Computational Mathematics,P. O. Box 8009-26, Beijing 100088, China3, Center for Applied Physics and Technology,MOE Key Center for High Energy Density Physics Simulations,College of Engineering, Peking University, Beijing 100871, China (Dated: September 8, 2020) ∗ Corresponding author. E-mail: Xu [email protected] † Corresponding author. E-mail: [email protected] bstract A new kinetic model for multiphase flow was presented under the framework of the discreteBoltzmann method (DBM). Significantly different from the previous DBM, a bottom-up approachwas adopted in this model. The effects of molecular size and repulsion potential were describedby the Enskog collision model; the attraction potential was obtained through the mean-field ap-proximation method. The molecular interactions, which result in the non-ideal equation of stateand surface tension, were directly introduced as an external force term. Several typical benchmarkproblems, including Couette flow, two-phase coexistence curve, the Laplace law, phase separation,and the collision of two droplets, were simulated to verify the model. Especially, for two types ofdroplet collisions, the strengths of two non-equilibrium effects, ¯ D ∗ and ¯ D ∗ , defined through thesecond and third order non-conserved kinetic moments of ( f − f eq ), are comparatively investigated,where f ( f eq ) is the (equilibrium) distribution function. It is interesting to find that during thecollision process, ¯ D ∗ is always significantly larger than ¯ D ∗ , ¯ D ∗ can be used to identify the differentstages of the collision process and to distinguish different types of collisions. The modeling methodcan be directly extended to a higher-order model for the case where the non-equilibrium effect isstrong, and the linear constitutive law of viscous stress is no longer valid. PACS numbers: 05.20.Dd, 51.10.+y, 47.11.-j
Keywords: multiphase flow, discrete Boltzmann method, Enskog equation, non-equilibrium characteristics . INTRODUCTION Multiphase flow is ubiquitous in nature and industrial processes, such as the formationof raindrops, fabrication of pills, oil and gas exploitation, micro-fluidic chip[1–4], etc. Un-derstanding the hydrodynamic and thermodynamic characteristics of multi-phase flow is ofgreat significance for fundamental research and engineering applications.Generally, multiphase flow refers to the fluid flow consisting of two or more phases. It ismore complicated than a single-phase flow process, due to phase interfaces and the interac-tion between the different phases [5]. The multiphase flow process is often accompanied byabundant interfacial dynamic behaviors, including interface generation, movement, defor-mation, fragmentation, and fusion, etc. These complex interfacial behaviors are challengingto measure experimentally and very difficult, if not impossible, to obtain analytic solutions.As computation power increased in recent decades, numerical simulation has graduallybecome a critical alternative approach to study multiphase flow problems [6–8]. In gen-eral, the numerical simulation methods of multiphase flow can be divide into three levels:macroscopic, microscopic, and mesoscopic methods. The macroscopic approaches are mainlybased on the numerical solutions of the hydrodynamic (Navier-Stokes, NS) equations com-bined with interface capture schemes [6]. Although those methods have been successfullyapplied to the large-scale multiphase flow problems, they fail to describe the micro-nanoscale interface structure because the equations based on the continuum assumption may notbe applicable in the micro-nano scale [8, 9]. In addition, the fluctuations are significantwith the decreasing of system size and can even perceivably affect the flow characteristicsat the micro-nano scale [10]. The surface tension and viscosities might need to be furtherconsidered with the effect of fluctuation. The microscopic model mainly refers to the molec-ular dynamic (MD) method, which is more fundamental and accurate [11]. However, thecomputational burden of MD simulation is often unaffordable; consequently, the time andspace scales it can simulate are limited. As a bridge between macroscopic and microscopic,the mesoscopic model has developed rapidly over the past three decades [12–16]. It has beensuccessfully applied in a wide variety of multiphase problems, including wetting, droplet colli-sion and splashing, boiling and evaporation, phase separation, and hydrodynamic instability[6, 17–19], etc.The mesoscopic model for multiphase flow, also known as the kinetic model, is mainly3ased on the Boltzmann equation. The vast majority, if not all, of the mesoscopic multi-phase flow models in the literature, are developed from the lattice Boltzmann model (LBM),including the color-gradient model [20], the Shan-Chen model [21], the free energy model[22, 23], and the phase-field model [24–27], etc. It is undeniable that they are very convenientin modeling multiphase flow processes, and have achieved great success in the research ofvarious multiphase flow problems [17]. However, it should be noted that the standard LBM,on which most multiphase flow models base, is just a solver of the incompressible Navier-Stokes (NS) equations. Although it may be more convenient to deal with the interparticleforces and the complicated boundaries [13], it can not provide kinetic information beyondthe NS equations.Based on the thermal finite different LBM of Watari and Tsutahara [28], Gonnella etal. first proposed the thermal multiphase model [29], in which the nonideal gases and thesurface tension effects are introduced by an extra term I ki . This multiphase flow modelcan correctly reproduce the transport equations of nonideal fluid established by Onuki [30],through a Chapman-Enskog multiscale expansion. Gan et al.[31] further developed a thermalmultiphase model with negligible spurious velocities by introducing the Windowed FastFourier Transform scheme when calculating the spatial derivatives in the convection termand the external force term. Then morphological characteristics of phase separation in thethermal systems were systematically studied [32].In the early studies, the non-equilibrium behaviors in various complex flow were all thosedescribed by hydrodynamic equations. For the convenience of description, we refer thosebehaviors to as Hydrodynamic Non-Equilibrium (HNE) behaviors. Besides the HNE, themost relevant Thermodynamic Non-Equilibrium (TNE) practices in various complex flowsystems are attracting more and more attention with time[33–36].It is understandable that when the TNE is very weak, the loss of considering TNE is notmeaningful. But when the TNE is strong, the situation will be significantly different. Forexample, it has been shown that TNE’s existence directly affects the density, temperature,and pressure, as well as the magnitude and direction of flow velocity. Without consideringTNE, the density, flow velocity, temperature, and pressure given will have a significantdeviation [37]. The existence of TNE is the underlying cause of the appearance of heat flowand viscous stress. If insufficiently considered (considered only the linear response part), theamplitude of heat flow and viscous stress obtained may be too large[38].4hen the strength of TNE beyond some threshold value, it may change the directions ofheat flow and viscous stress. If the TNE is not sufficiently considered (considered only thelinear response part), the obtained viscous stress and heat flow maybe, even in the wrongdirections [9]. In addition, in the phase separation system, both the mean TNE strength[39]and the entropy production rate, one of the quantities describing the TNE effects[40], in-crease with time in the spinnodal decomposition stage and decrease with time in the domaingrowth stage. Therefore, both the peak value points of the mean TNE strength and the en-tropy production rate can work as physical criteria to discriminate the two stages[35, 39, 40].In the system with combustion, the TNE behaviors help to better understand the physicalstructures of the von Neumann peak and various non-equilibrium detonation[34, 37, 41–47]. In the system with Rayleigh-Taylor instability, the TNE behaviors around interfaceshave been used to physically identify and distinguish various interfaces and design relevantinterface-tracking schemes. With increasing the compressibility, more observable TNE ki-netic modes appear for given observation precision[36, 48]. The correlation between themean density nonuniformity and mean TNE strength is almost 1. The correlation betweenthe mean temperature nonuniformity and mean Non-Organized Energy Flux (NOEF) isalmost 1. The correlation between the mean flow velocity nonuniformity and the meanNon-Organized Momentum Flux (NOMF) is also high, but generally less than 1 [36, 49].The TNE effect helps to understand the effect of system dispersion (described by Knudsennumber) on its kinetic behaviors[50]. In the system with Kelvin-Helmholtz instability, viasome defined TNE quantity, for example, the heat flux intensity, the density interface andtemperature interface can both be observed so that the material mixing and energy mix-ing in the Kelvin-Helmholtz instability evolution can be investigated simultaneously[9, 51].The TNE behaviors were used to understand the mixing entropy in multi-component flowsbetter[52].Since the traditional hydrodynamic model is incapable of capturing enough TNE behav-iors, the works mentioned above on both HNE and TNE resorted to the recently proposeddiscrete Boltzmann method (DBM)[33, 36]. Compared with the standard LBM, the modelsystem’s evolution does not resort to the simple “propagation + collision” scenario of virtualparticle which is inherited from the lattice-gas-automaton model. Any suitable numericalscheme can be used to solve the discrete Boltzmann equation(DBE) according to the specificcase. The non-conserved kinetic moments of ( f i − f eqi ) are used to describe how the system5eviates from its thermodynamic equilibrium state and to fetch information on the TNE be-haviors, where f i ( f eqi ) is the discrete (equilibrium) distribution function and i is the indexof discrete velocity. The kinetic moment relations to be based in the construction processof DBM are dependent on the research perspective and the non-equilibrium extent whichthe DBM aims to describe[36]. Besides by theory, results of DBM have been confirmed andsupplemented by results of molecular dynamics[53–55], direct simulation Monte Carlo[9, 56]and experiment[46].The existing multiphase flow models based on DBM is a semi-kinetic and semi-top-downapproach. The discrete Boltzmann equation is modified by an extra term to be consistentwith the correct hydrodynamic equations of nonideal fluid. The specific forms of the co-efficients (expressed as the gradient of macroscopic quantities) in the external force termare determined by the Chapman-Enskog multiscale expansion [29, 57]. Therefore, in themodeling of multiphase flow DBM, the hydrodynamic equations of nonideal fluid need tobe known first. This partly limits the ability of DBM to study the non-equilibrium complexflow, because the hydrodynamic equations are often unknown when the non-equilibriumstrength is strong. In this paper, this restriction is addressed by developing a bottom-upkinetic modeling method of multiphase flow under the framework of DBM.The rest of this paper is organized as follows. In Sec. II, the simplified Enskog equationis derived, and the multiphase model based on the discrete Enskog equation is presented.Several benchmark tests of multiphase flow are simulated to verify the new model, and thenon-equilibrium characteristics during the collision of two droplets are analyzed in Sec.III.Section IV concludes the present paper. II. MODELS AND METHODSA. Enskog equation and its simplification
According to molecular kinetic theory, the evolution equation of molecular velocity dis-tribution function reads ∂f∂t + v · ∇ f + a · ∇ v f = (cid:18) ∂f∂t (cid:19) c , (1)where f = f ( r , v , t ) is the molecular velocity distribution function, r and v represent theposition space and velocity space coordinates, respectively, a is the acceleration generated6y the total extra force. (cid:0) ∂f∂t (cid:1) c denotes the rate of change in distribution function due tothe collisions between molecules. Given that the probability of two molecules colliding ismuch higher than that of three or more molecules colliding simultaneously, the hypothesisof two-body collision is reasonable.Based on the elastic molecular collision model, when the molecular volume effect is ig-nored, the collision term in the Boltzmann equation can be obtained as [58] (cid:18) ∂f∂t (cid:19) c = Z ∞−∞ Z π ( f ∗ f ∗ − f f ) v r σd Ω d v , (2)where f ∗ and f ∗ represent the post-collision molecular velocity distribution function withvelocity v ∗ and v ∗ , respectively, f and f represent the pre-collision molecular velocitydistribution function with velocity v and v , respectively. v r = | v − v | is the value of therelative velocity, which remains unchanged pre to and after the collision. σ and Ω denotethe differential collision cross-section and solid angle, respectively.However, this assumption is not appropriate for dense gases or liquid. With increasingthe number density, compared with the mean distance between neighboring molecules, thesize of a gas molecule is no longer negligible anymore, the volume effects of the moleculeshould be taken into account, the Boltzmann collision operator should be replaced by theEnskog collision operator which reads [59] (cid:0) ∂f∂t (cid:1) E = R ∞−∞ R π (cid:2) χ (cid:0) r + d ˆe r (cid:1) f ∗ ( r ) f ∗ ( r + d ˆe r ) − χ (cid:0) r − d ˆe r (cid:1) f ( r ) f ( r − d ˆe r ) (cid:3) v r σd Ω d v , (3)where d is the diameter of the hard-sphere molecules and χ represents the collision prob-ability correction considering the molecular volume effect. By using the Taylor expansionand keeping to the first derivative term, one can get that χ ( r + d ˆe r ) = χ ( r ) + d ∇ χ · ˆe r χ ( r − d ˆe r ) = χ ( r ) − d ∇ χ · ˆe r f ∗ ( r + d ˆe r ) = f ∗ ( r ) + d ∇ f ∗ · ˆe r f ( r − d ˆe r ) = f ( r ) − d ∇ f · ˆe r (cid:0) ∂f∂t (cid:1) E = χ R ∞−∞ R π ( f ∗ f ∗ − f f ) v r σd Ω d v + d χ R ∞−∞ R π ( f ∗ ∇ f ∗ + f ∇ f ) · ˆe r v r σd Ω d v + d R ∞−∞ R π ∇ χ · ˆe r ( f ∗ f ∗ − f f ) v r σd Ω d v (4)where χ , f ∗ , and f are all at r . If f in the latter two terms are approximated by the localequilibrium distribution function f eq , f eq = ρ πRT ) / exp[ − ( v − u ) RT ] . (5)Then the Enskog collision operator becomes (cid:0) ∂f∂t (cid:1) E = χ R ∞−∞ R π ( f ∗ f ∗ − f f ) v r σd Ω d v − f eq bρχ n ( v − u ) · h ρ ∇ ρ + T ∇ T (cid:16) ( v − u ) RT − (cid:17)i + RT ( v − u ) ( v − u ) : ∇ u + (cid:16) ( v − u ) RT − (cid:17) ∇ · u o − f eq bρ ( v − u ) · ∇ χ, (6)where bρ = πnd .For isothermal systems, the collision term can be reduced to (cid:0) ∂f∂t (cid:1) E = χ R ∞−∞ R π ( f ∗ f ∗ − f f ) v r σd Ω d v − f eq bρχ n ( v − u ) · ρ ∇ ρ + RT ( v − u ) ( v − u ) : ∇ u + h ( v − u ) RT − i ∇ · u o − f eq bρ ( v − u ) · ∇ χ. (7)In addition, the last two terms in the second row of Eq. (7) satisfy the following relationships: Z f eq bρχ ( RT ( v − u )( v − u ) : ∇ u + 35 " ( v − u ) RT − ∇ · u ) d v = 0 , (8) Z f eq bρχ ( RT ( v − u )( v − u ) : ∇ u + 35 " ( v − u ) RT − ∇ · u ) v d v = 0 . (9)The sum of these two terms have no effect on the evolution of density and momentum [60],so they can be ignored in the modeling framework of this paper. As a result, the Enskogcollision term can be further simplified as (cid:0) ∂f∂t (cid:1) E = χ R ∞−∞ R π ( f ∗ f ∗ − f f ) v r σd Ω d v − f eq bρχ ( v − u ) · ρ ∇ ρ − f eq bρ ( v − u ) · ∇ χ. (10)8imilar to the modeling of DBM, the first term on the right-hand side can be replacedby the BGK collision operator, consequently Eq. (7) becomes (cid:18) ∂f∂t (cid:19) E = − τ ( f − f eq ) − f eq bρχ ( v − u ) · ∇ ln (cid:0) ρ χ (cid:1) . (11)Now we can get the simplified Enskog model for the approximately incompressible fluidsystems when the temperature is uniform, which reads ∂f∂t + v · ∇ f + a · ∇ v f = − τ ( f − f eq ) − f eq bρχ ( v − u ) · ∇ ln (cid:0) ρ χ (cid:1) . (12)Similar to previous studies, in this work we consider the case where the force term a · ∇ v f can be approximated by a · ∇ v f = a · ∇ v f eq = − F · ( v − u ) ρRT f eq , (13)where m is the mass of the molecule. As a result, the simplified Enskog equation becomes ∂f∂t + v · ∇ f − F · ( v − u ) ρRT f eq = − τ ( f − f eq ) − ( v − u ) · ∇ ( bρ RT χ ) ρRT f eq . (14)The last term on the right side of the equation represents the repulsion interaction ofmolecules. B. Multiphase flow model based on Enskog equation
Now the repulsion interaction of molecules is introduced by using the Enskog collision op-erator instead of the Boltzmann collision operator, then the key to modeling multiphase flowis to incorporate the intermolecular attraction. According to the average filed approxima-tion, the mutual attraction between molecules can be regarded as an average force potential[24], which reads Φ( r ) = Z r >d φ attr ( r ) ρ ( r ) d r , (15)where r = | r − r | is the distance between r and r , φ attr ( r ) denotes the attractionpotential. Expanding ρ ( r ) at r and ignoring terms higher than the second order, Φ( r )can be approximated by Φ( r ) = − aρ − K ∇ ρ, (16)where the first term mainly affects the equation of state (EOS) with a = − R r>d φ attr ( r ) d r and the second term contributes to the surface tension. The coefficient of surface tension is K = − R r>d r φ attr ( r ) d r , and it is assumed to be a constant in this work.9ased on the expression of Φ( r ) in Eq. (16), the total external force of molecules at r can be calculated as F = − ρ ∇ Φ = ∇ ( aρ ) + ρ ∇ ( K ∇ ρ ) . (17)Substituting the expression of F into Eq. (14) gives ∂f∂t + v · ∇ f − F ′ · ( v − u ) ρRT f eq = − τ ( f − f eq ) , (18)where F ′ = [ −∇ ψ + ρ ∇ ( K ∇ ρ )] with ψ = bρ RT χ − aρ . There are two terms in F ′ , thefirst term is related to the EOS and the second one corresponds to the surface tension.According to Chapman-Enskog multi-scale expansion, the hydrodynamic equations formultiphase flow can be obtained from Eq. (18) as follows: ∂ρ∂t + ∇ · ( ρ u ) = 0 , (19) ∂ ( ρ u ) ∂t + ∇ · ( ρ uu + P I + Π ) − ρ ∇ ( K ∇ ρ ) = 0 , (20)where P represents the EOS of real gas in the form of P = ρRT (1 + bρχ ) − aρ . (21)When χ = − bρ/ − bρ/ , the Carnahan-Starling EOS can be obtained as P c = ρRT η + η − η (1 − η ) − aρ , (22)with η = bρ . The van der Waals (VDW) EOS P v = ρRT − bρ − aρ , (23)can also be derived from Eq. (21) when χ = − bρ . The Π in Eq. (20) represents the viscousstress which has an expression as Π = µ (cid:18) ∇ u + ∇ T u − D ∇ · u (cid:19) , (24)where D denotes spatial dimension and µ is the coefficient of dynamic viscosity, µ = τ ρRT .The last term, ρ ∇ ( K ∇ ρ ), on the left-hand side of Eq. (20) denotes the surface tension F s ,which is equivalent to the following expression [19] F s = ∇ (cid:20) ( K ∇ ρ · ∇ ρ + Kρ ∇ ρ ) I − K ∇ ρ ∇ ρ (cid:21) . (25)10 . Discrete Enksog model for multiphase flow After the previous derivation, the kinetic equation of multiphase flow is obtained as Eq.(18). Correspondingly, the evolution of discrete Enskog equation for multiphase flow reads ∂f ki ∂t + v ki · ∇ f ki − ( v ki − u ) · F ′ ρRT f eqki = − τ ( f ki − f eqki ) , (26)where f ki = f ki ( r , t ) is the distribution function of the discrete velocity v ki and the subscript“ ki ” indicates the index of discrete velocity. f eqki is the discrete local equilibrium distributionfunction; its general expression of any order has been given in the previous literature [9].In the previous Chapman-Enskog multi-scale expansion precess, only the 0th to 3rd kineticmoments are needed to recover Eqs. (19) and (20). Therefore, for computational efficiency,the discrete local equilibrium distribution function is approximated by the 3rd-order Hermitepolynomial which reads f eqki = ρF k " (1 − u T ) + 11! (1 − u T ) v ki · u T + 12! ( v ki · u ) T + 13! ( v ki · u ) T , (27)where F k is the weight coefficient which can be expressed as F = 24 T − c + c ) T + c c Tc ( c − c ) ( c − c ) , (28) F = 24 T − c + c ) T + c c T c ( c − c ) ( c − c ) , (29) F = 24 T − c + c ) T + c c T c ( c − c ) ( c − c ) , (30)and F = 1 − F + F + F ) . (31)where c , c , and c are the values of three groups of discrete velocities. The scheme of thediscrete velocity model is shown in Fig. 1. The value of c k does not affect the final results aslong as the calculation can be performed stably. It should be noted that the new model canbe easily extended to higher-order cases only if the discrete local equilibrium distributionfunction satisfies more kinetic moment relations. Higher-order terms need to be kept in theHermite polynomials of Eq. (27) to achieve this aim. Accordingly, more discrete velocitiesare required to satisfy the higher-order isotropy of the discrete model. More details of thehigher-order model are available in Ref. [9].11 IG. 1: Schematic diagram of discrete velocity model.
Once the weight coefficients F k are calculated, the discrete local equilibrium distributionfunction f eqki in Eq. (27) is also known, then the discrete distribution function in t + △ t ( f t + △ tki ) can be solved from f tki by the following evolution equation f t +∆ tki = 1 τ (cid:0) f tki − f eqki (cid:1) − v ki · ∇ f tki + ( v ki − u ) · F ′ ρRT f eqki . (32)It involves the calculation of spatial derivatives in both the convection term ∇ f tki and theexternal force term F ′ . Various types of different schemes, such as finite-difference, Nine-Point Stencils (NPS), Nonoscillatory Non-free-parameter and Dissipative (NND), and FastFourier Transform (FFT) schemes, etc. can be used depending on the specific problems.In addition, from Eq. (32) we see that f ki and its corresponding local equilibrium distri-bution function f eqki are both known at a certain time. Consequently, the non-equilibriumquantities ∆ ∗ and ∆ ∗ can be calculated as ∆ ∗ = X ki ( f ki − f eqki )( v ki − u )( v ki − u ) , (33) ∆ ∗ = X ki ( f ki − f eqki )( v ki − u )( v ki − u )( v ki − u ) . (34)Correspondingly, the quantities of non-equilibrium strength, D ∗ and D ∗ , can be obtained,which are defined as D ∗ = q(cid:12)(cid:12) ∆ ∗ ,xx (cid:12)(cid:12) + (cid:12)(cid:12) ∆ ∗ ,xy (cid:12)(cid:12) + (cid:12)(cid:12) ∆ ∗ ,yy (cid:12)(cid:12) (35)and D ∗ = q(cid:12)(cid:12) ∆ ∗ ,xxx (cid:12)(cid:12) + (cid:12)(cid:12) ∆ ∗ ,xxy (cid:12)(cid:12) + (cid:12)(cid:12) ∆ ∗ ,xyy (cid:12)(cid:12) + (cid:12)(cid:12) ∆ ∗ ,yyy (cid:12)(cid:12) , (36)12espectively. In previous studies [9], it has been found that some quantities of non-equilibrium strength may perform better than individual non-equilibrium components indescribing some characteristics of fluid interfaces. In this work, the role of these non-equilibrium quantities in multiphase flow will be further investigated. III. SIMULATION AND DISCUSSIONA. Couette flow
As the first test, the Couette flow is simulated to examine the viscous stress calculatedby the discrete model with the discrete velocity model shown in Fig. 1. Two parallel plates,filled with fluid between them, are placed in the y direction. The left plate is stationarywhile the right plate moves at a constant speed U y . The fluid is driven by the plate dueto the viscous stress. This flow process is simulated by using the discrete model withoutconsidering the non-ideal gas effects and surface tension.The speed of the moving plate is U y = 0 . T = 1 .
0. The computational meshes are N x × N y = 100 × dx = dy = 0 .
002 and time step dt = 0 . τ . The initial velocity in the flow filed is zero.The non-slip boundary condition is used on both left and right boundary. The first-orderforward difference is used for time discretization and the NND scheme is used for spatialdiscretization. The simulation results are shown in Fig. 2.Figure 2 (a) gives the profiles of velocity U y along the x direction at several differenttimes with a fixed relaxation time τ = 0 . U y withseveral different relaxation time at t = 1 .
0. The analytical solutions denoted by the solidlines are also plotted for comparison. The simulation results are in good agreement with theanalytical solutions, which verifies the accuracy of the viscous stress calculated by the newdiscrete model.
B. Two-phase coexistence curve
The two-phase coexistence problem is one of the most typical multiphase flow problemsthat can be described by a single-component multiphase flow model. According to the13
IG. 2: Verification of viscosity stress calculated by the new model. (a) Velocity profiles of Couetteflow at several different times, (b) Velocity profiles of Couette flow with several different viscositycoefficients represented by relaxation time τ . The symbols are DBM results while the solid linesare analytical solution. VDW EOS, below the critical temperature T c , the system allows for two coexisting fluidwith different densities under the same pressure P . The high density ( ρ l ) corresponds tothe liquid phase, while the lower density ( ρ v ) to the vapor phase. Seen from the EOS shownin Fig. 3 (a), there are many groups of coexistence points at different pressures, as long asthe value of pressure is between A and B ; the liquid phase corresponds to the AB segmentwhile the vapor phase to the CD segment. However, only one set of liquid-vapor pointssatisfy both the mechanical equilibrium and thermodynamic equilibrium (more specifically,the chemical potential equilibrium) at a given temperature T . This set of liquid-vaporpoints can be calculated by the Maxwell equal area rule.The liquid-vapor coexistence points at several different temperatures are calculated byusing the new multiphase flow model. The computation meshes are N x × N y = 250 × dx = dy = 4 × − . The time step and the relaxation time are dt = 5 × − and τ = 0 .
02, respectively. The VDW EOS in Eq. (23) is adopted with a = 9 / b = 1 /
3. The coefficient of surface tension is K = 2 × − . The first-orderforward difference is used for time discretization. The NND scheme is adopted to calculatethe spatial derivative in the convection term while the NPS scheme is used to calculate thespatial derivatives in the force term F ′ . Periodic boundary conditions are adopted on both14 IG. 3: (a) P versus 1 /ρ for the equation of state of van der Waals, (b) Verification of non-idealequation of state by two-phase coexistence curve. left and right boundaries. The initial conditions are set as follows: ( ρ, T, u x , u y ) L = ( ρ v , . , , ρ, T, u x , u y ) M = ( ρ l , . , , ρ, T, u x , u y ) R = ( ρ v , . , ,
0) (37)where the subscript “L”,“M”, and “R”indicate the regions 0 < x ≤ N x / N x / < x ≤ N x / N x / < x ≤ N x , respectively. ρ v and ρ l are the theoretical vapor and liquid densitiesat T = 0 .
99. According to the VDW EOS combined with the equal area rule we get that ρ v = 0 . ρ l = 1 . .
01 and wait for the system to achieve another equilibrium state. The aboveprocess is repeated until the temperature drops to T = 0 .
85. The vapor and liquid densitiescalculated by the new model at different temperatures are plotted in Fig. 3 (b). The symbolsare simulated by the new multiphase flow model while the solid line represents the analyticsolution. The simulation results are in excellent agreement with the theoretical solutions,which verified the accuracy of the new multiphase model in calculating the equilibrium phasetransition.To further validate the effect of surface tension, the transition curves between the liquidand vapor phases with different coefficients of surface tension, i.e., K = 2 × − , K =15 IG. 4: The density interface between the liquid and vapor phases at T = 0 .
95. The symbols arecalculated by DBM and the solid lines are theoretical solutions. × − , and K = 1 × − , at T = 0 .
95 are calculated and the results are shown in Fig. 4.The theoretical solutions are also plotted for comparison. The theoretical solution of phaseinterfaces given by the modified VDW theory reads [61] x − x = − p a/K Z ρ ∗ ( x ) ρ ∗ ( x ) p Φ ∗ ( ρ ∗ ) − Φ ∗ ( ρ ∗ l ) dρ ∗ (38)with Φ ∗ ( ρ ∗ ) = ρ ∗ ξ − ρ ∗ T ∗ [ln(1 /ρ ∗ −
1) + 1] − ( ρ ∗ ) (39)and ξ = T ∗ ln(1 /ρ ∗ s − − ρ ∗ s T ∗ − ρ ∗ s + 2 ρ ∗ s (40)where ρ ∗ and T ∗ denotes the reduced density and temperature, ρ ∗ = bρ and T ∗ = bT /a , ρ ∗ s represents the reduced density of the liquid or vapor corresponding to the equilibrium stateat T ∗ .The theoretical solutions are shown as solid lines and the DBM results are represented bysymbols. It shows that the density profiles of the interface simulated by the new multiphasemodel agree well with the theoretical solutions, which verifies the accuracy of the surfacetension of the new model. 16 IG. 5: Simulation results of circular droplets suspended in its vapor: (a) density contour mapand (b) pressure contour map.
C. Droplet suspension
To further validate the effect of surface tension, a circular droplet suspended in its vaporis simulated. According to the Laplace law, under a fixed surface tension coefficient, thepressure difference ∆ P inside and outside the droplet is inversely proportional to the radiusof the droplet R , i.e., ∆ P ∼ /R .Several droplets with various radii are simulated under two different coefficients of surfacetension K = 2 × − and K = 2 × − , respectively. The computational meshes are N x × N y = 128 ×
128 with the spatial steps dx = dy = 0 .
02 and time step dt = 0 . R = 0 . P inside and outside the circular droplet is detected during thecalculation. The profile of ∆ P with time is shown in Fig. 6(a) and it eventually tends toa stable value. The calculation continues until ∆ P is almost unchanged, then the pressuredifference ∆ P and the corresponding radius R is recorded. Several circular droplets withdifferent radii are simulated. Figure 6(b) shows the relation between the ∆ P and 1 /R , inwhich symbols are simulation results and the solid lines are linear fitting. Clearly, ∆ P is17 IG. 6: Verification of surface tension by the Laplace’s law: (a) evolution of pressure difference∆ P between the inside and the outside the droplet with time t and (b) the relationship between∆ P and the radius of droplet R . proportional to 1 /R , which agrees well with the Laplace law. The accuracy of the surfacetension calculated by the new multiphase model is well verified from this simulation. D. Phase separation
In this part, a dynamic phase separation process is simulated by using the new multiphaseflow model. The initial conditions are set as( ρ, T, u x , u y ) = (1 . δ, . , ,
0) (41)where δ is a random noise with amplitude 0.01, which provides a starting point for phaseseparation. In the phase separation process, the temperature is fixed at T = 0 .
92, which isequivalent to that the system contacts with a large heat source with a temperature T = 0 . N x × N y = 300 ×
300 with spatial steps ∆ x = ∆ y = 0 . t = 1 × − and τ = 0 .
01, respectively. Periodicboundary conditions are adopted on all boundaries and corners. The first-order forwarddifference is used for time discretization. The NND and NPS schemes are used to calculatethe spatial derivatives in the convection term and the force term. Figure 7 shows the densitycontour maps at several different times. The simulation results are consistent with those inprevious literature[22, 40, 57]. 18
IG. 7: Density contour maps in the process of isothermal phase separation at (a) t = 0 .
5, (b) t = 1 .
0, (c) t = 2 .
0, (d) t = 5 .
0, (e) t = 9 .
0, and (f) t = 14 .
0, respectively.
E. Droplets collision
Droplet collision is another typical multiphase flow problem and has a broad applicationbackground in industrial production. It involves complex interfacial dynamics, includinginterface movement, fusion, and separation, etc. In this part, we simulate the head-oncollision between two droplets using the new multiphase flow model. Two droplets withthe opposite speeds are suspended in the same horizontal position. The droplet on the leftmoves to the right at speed U , while the droplet on the right moves to the left at the samespeed. There are two situations after the collision depending on the initial speed U . If U is much small, the two droplets eventually fuse together, whereas if U is large enough, theyseparate again along the vertical direction.Two cases with U = 0 . U = 0 .
5, respectively, are simulated by using the newmultiphase flow model. The computational meshes are N x × N y = 60 ×
120 with the spatialstep ∆ x = ∆ y = 0 .
02. The time step is ∆ t = 1 × − and the relaxation time is τ = 2 × − .The coefficient of surface tension is set as K = 2 × − and the temperature is fixed at T = 0 .
92. The boundary conditions and the discrete schemes are the same as the previoussimulation of phase separation. The snapshots of the collision processes for two cases areshown in Figs. 8 and 9, respectively. Figure 8 shows the case with U = 0 .
2, in which two19
IG. 8: Head-on collision of two droplets which are fused together after the collision.FIG. 9: Head-on collision of two droplets which are separated after the collision. droplets fuse together after the collision and figure 9 shows the case with U = 0 . D ∗ and D ∗ , areshown in Fig. 10, in which the left and right subgraphs correspond to the two kinds ofdroplets collision with U = 0 . U = 0 .
5, respectively. The spatial mean values ¯ D ∗ and¯ D ∗ are used to represent the total mean non-equilibrium strength, i.e.¯ D * m = 1 l x l y Z Z D * m dxdy, m = 2 , l x l y represents the total area of the simulation region.Figure 10 (a) shows that at least five stages of the collision can be identified from thenon-equilibrium strength ¯ D ∗ . In the first stage, two droplets come close to each other,and ¯ D ∗ gradually decreases until it reaches point “A” in the ¯ D ∗ ( t ) curve. In the secondstage (AB segment on the profile of ¯ D ∗ ), two droplets contact and the interface begins tofuse, and ¯ D ∗ increases with time. After the phase interface fusion is completed, the thirdstage (BC segment) begins, and the merged droplet elongate vertically, accompanied by adecrease of ¯ D ∗ . Because the kinetic energy after collision is not enough to overcome theconstraint of surface tension, the merged droplet can not be separated again. As a result,in the fourth stage (CD segment), the droplet contracts vertically and the overall trendof the non-equilibrium strength is upward. In the fifth stage (after the “D” point in the¯ D ∗ ( t ) curve), the droplet lengthens and contracts alternately in the horizontal and verticaldirections, and finally tends to a stable state. Correspondingly, the ¯ D ∗ shows the trend ofoscillation attenuation and finally tends to a stable value. From the profile of ¯ D ∗ , only thefirst three stages can be well identified, which indicates that the roles of non-equilibriumstrength of different orders are much different.For the second type of droplet collision with U = 0 .
5, there are also five stages as shownin Fig. 10 (b). The profile of ¯ D ∗ in the first four stages is similar to the case with U = 0 . D ∗ continues to decrease until the separated droplets leave thesimulation region at t = 5. The profiles after t = 5 is meaningless because the droplet hasleft the simulation region. So these two types of collisions can be well distinguished by thelast stage of ¯ D ∗ . Likewise, from the profile of ¯ D ∗ , only the first three stages can be wellidentified. As a result, it is difficult to distinguish these two types of collision cases from theprofile of ¯ D ∗ .The evolution of the non-equilibrium strength with time can be explained by the changeof the velocity gradient in the flow field. Because ∆ ∗ corresponds to the viscous stress Π IG. 10: Profiles of two kinds of non-equilibrium strength in the collision process, where (a) isfor head-on collision of two droplets which are fused together after the collision, (b) is for head-oncollision of two droplets which are separated after the collision. It is clear that ¯ D can be usedto characterize the stage behaviors of the collision process and ¯ D shows significantly differentbehaviors in the last stage ( after the “D” point in the ¯ D ∗ ( t ) curve ). in the hydrodynamic equations, which can be expressed as Eq. (24) in the NS level. In thisapproximation, D ∗ can be expressed as D ∗ ≈ (2 µ ) / (cid:20) ∇ u : ∇ u + 2( ∂u x ∂y ∂u y ∂x − ∂u x ∂x ∂u y ∂y ) (cid:21) / . (43)We further find that the last two terms on Eq. (43) have little effect on the overall evolu-tionary characteristics. For the sake of convenience, we will only use ( ∇ u : ∇ u ) / in thefollowing analysis. Figure 11 shows the profiles of ¯ D ∗ and the spatial average of ( ∇ u : ∇ u ) / .The definition of ( ∇ u : ∇ u ) / is the same as that of ¯ D ∗ m in Eq. (42). Figures 11 (a) and (b)correspond to the first and second types of collision processes, respectively. It shows that¯ D ∗ and ( ∇ u : ∇ u ) / have similar evolutionary characteristics in both Fig. 11 (a) and (b),which is consistent with our previous analysis.For the first type of droplet collision, the collision process is roughly divided into fivestages as shown in Fig. 11 (a). In the first stage, the two droplets approach each otherindependently. Part of the kinetic energy of the droplets is transferred to the flow field ofvapor phase, leading to a decrease of the system’s velocity gradient whole. In the secondstage, two droplets collide and fuse. During the fusion process, the interface contracts rapidly22ue to the work done by surface tension, which results in the increase of velocity gradient inthe flow field. After the droplets are completely fused, the new droplet begins to stretch inthe y direction in the third stage. In this process, the surface tension does the negative work,causing the interface to slow down and the velocity gradient to decrease. In the fourth stage,because the kinetic energy of the interface motion is not enough to overcome the constraintof surface tension, the interface begins to contract in reverse when it stretches to the longestin the y direction. The surface tension does positive work, so the interface accelerates andthe velocity gradient increases. Soon the interface shrinks to its minimum in the y directionand begins to stretch again in the x direction. The velocity gradient decreases in the processof stretching and increases in the process of shrinking. Finally, in the fifth stage, the dropletis alternately stretched and contracted in the x and y direction, but the maximum positionit can reach gradually shortens until it reaches a stable state.For the second type of droplet collision, the first stages are the same as those of the firsttype of collision. However, in the fourth stage, the fused droplet rupture and separate again.The interface of the two separated droplets shrinks rapidly, leading to an increase of velocitygradient. After that, in the fifth stage, the single droplet decelerates independently, so thevelocity gradient decreases until the droplets leave the simulation area. ∆ ∗ is a higher-ordernon-equilibrium characteristic quantity related to heat flux but its corresponding macro-scopic physical significance is not yet clear. Because the temperature is set as a constant inthe simulation, the heat flow in the flow field is sharply limited, so the amplitude is relativelysmall and its evolution characteristics are not visible. IV. CONCLUSION
In this work, under the framework of the discrete Boltzmann method, a new kinetic mul-tiphase flow model based on a discrete Enskog equation was proposed. Compared with theprevious multiphase flow discrete Boltzmann model, a bottom-up modeling approach wasadopted in the new model. The repulsion potential was introduced by using the Enskog col-lision operator instead of the Boltzmann collision operator, and the attraction potential isincorporated through average filed approximation. The effect of total intermolecular forcesis ultimately taken into account via the force term of the DBM. The discrete local equilib-rium distribution function is solved by the Hermite polynomials. The third order Hermite23
IG. 11: Profiles of ¯ D ∗ and ( ∇ u : ∇ u ) / in the collision process, where (a) is for head-on collisionof two droplets which are fused together after the collision, (b) is for head-on collision of twodroplets which are separated after the collision. polynomial is adopted in order to improve the computational efficiency as much as possible,although it is straightforward to take the higher-order ones. Several benchmarks, includingthe Couette flow, two-phase coexistence curve, droplet suspension, phase separation, anddroplet collisions are simulated, through which the accuracy of the new multiphase flowmodel is verified. In addition, for two kinds of droplet collisions, the non-equilibrium char-acteristics are comparatively investigated via two TNE strength quantities, ¯ D ∗ and ¯ D ∗ . Itis found that during the collision process, ¯ D ∗ is always significantly larger than ¯ D ∗ , ¯ D ∗ can24e used to identify the different stages of the collision process and to distinguish differenttypes of collisions. Acknowledgements
Many thanks to the anonymous referees for their helpful comments and suggestion.This work is supported by the China Postdoctoral Science Foundation under Grant No.2019M662521, National Natural Science Foundation of China under Grant No. 11772064,CAEP Foundation (under Grant No. CX2019033) and the opening project of State KeyLaboratory of Explosion Science and Technology (Beijing Institute of Technology) underGrant No. KFJJ19-01 M. [1] Chen Y, Xie Q, Sari A, Bardy P, and Saeedi A, Oil/water/rock wettability: Influencing factorsand implications for low salinity water flooding in carbonate reservoirs, Fuel 215, 171-177(2018).[2] Chen Y and Deng Z, Hydrodynamics of a droplet passing through a microfluidic T-junction.J. Fluid Mech. 819, 401-434 (2017).[3] Tice J, Song H, Lyon A, and Ismagilov R, Formation of Droplets and Mixing in MultiphaseMicrofluidics at Low Values of the Reynolds and the Capillary Numbers, Langmuir 19(22),9127-9133 (2003).[4] Gunther A and Jensen K, Multiphase microfluidics: from flow characteristics to chemical andmaterials synthesis, Lab on A Chip 6(12), 1487 (2006).[5] Christopher E. Brennen, Fundamentals of Multiphase Flow, Cambridge: Cambridge Univer-sity Press, 2005.[6] Saurel R and Pantano C, Diffuse-Interface Capturing Methods for Compressible Two-PhaseFlows, Annu. Rev. Fluid Mech. 50(1), 105-130 (2018).[7] Frezzotti A, Barbante P, and Gibelli L, Direct simulation Monte Carlo applications to liquid-vapor flows, Phys. Fluids 31(6), 062103 (2019).[8] Worner M, Numerical modeling of multiphase flows in microfluidics and micro process engi-neering: a review of methods and applications, Microfluid. Nanofluid. 12(6), 841-886 (2012).
9] Zhang Y, Xu A, Zhang G, Chen Z, and Wang P, Discrete Boltzmann method for non-equilibrium flows: based on Shakhov model, Comput. Phys. Commun. 238, 50-65 (2019).[10] Moseler M and Landman U. Formation, Stability, and Breakup of Nanojets. Science, 289,1165 (2000).[11] Zhan S, Su Y, Jin Z, Zhang M, Wang W, Hao Y, and Li L, Study of liquid-liquid two-phaseflow in hydrophilic nanochannels by molecular simulations and theoretical modeling, Chem.Eng. J. 395, 125053 (2020).[12] Wolfram S, Cellular automaton fluids 1: Basic theory, J. Stat. Phys. 45(3-4), 471-526 (1986).[13] Chen S and Doolen G, Lattice Boltzmann Method for fluid flows, Annu. Rev. Fluid Mech.30(1), 329-364 (1998).[14] Succi S, The lattice Boltzmann equation: for fluid dynamics and beyond. Oxford: Oxforduniversity press, 2001.[15] He X and Doolen G D, Thermodynamic Foundations of Kinetic Theory and Lattice BoltzmannModels for Multiphase Flows, J. Stat. Phys. 107(1-2), 309-328 (2002).[16] Qin R, Mesoscopic interparticle potentials in the lattice Boltzmann equation for multiphasefluids, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 73(6), 066703 (2006).[17] Li Q , Luo K, Kang Q, He Y, Chen Q and Liu Q, Lattice Boltzmann methods for multiphaseflow and phase-change heat transfer, Prog. Energy Combust. Sci. 52, 62-105 (2016).[18] Qin R, Thermodynamic properties of phase separation in shear flow, Comput. Fluids 117,11-16 (2015).[19] Timm K, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, and Viggen E, The Lattice Boltz-mann Method - Principles and Practice, Springer, 2017.[20] Grunau D, Chen S, and Eggert K, A lattice Boltzmann model for multiphase fluid flows, Phys.Fluids 5(10), 2557-2562 (1993).[21] Shan X and Chen H, Lattice Boltzmann model for simulating flows with multiple phases andcomponents, Phys. Rev. E 47(3), 1815-1819 (1993).[22] Swift M R, Osborn W R, and Yeomans J M, Lattice Boltzmann simulation of nonideal fluids,Phys. Rev. Lett. 75(5), 830-833 (1995).[23] Xu A, Gonnella G, and Lamura A, Phase-separating binary fluids under oscillatory shear,Phys. Rev. E 67, 056105 (2003).[24] He X, Chen S, and Zhang R, A Lattice Boltzmann Scheme for Incompressible Multiphase Flow nd Its Application in Simulation of Rayleigh-Taylor Instability, J. Comput. Phys. 152(2),642-663 (1999).[25] Liang H, Li Q, Shi B, and Chai Z, Lattice boltzmann simulation of three-dimensional Rayleigh-Taylor instability, Phys. Rev. E 93, 033113 (2016).[26] Wang H, Yuan X, Liang H, Chai Z, and Shi B, A brief review of the phase-field-based latticeBoltzmann method for multiphase flows, Capillarity 2 (3), 33-52 (2019).[27] Sun D, A discrete kinetic scheme to model anisotropic liquid-solid phase transitions, Appl.Math. Lett. 103, 106222 (2020).[28] Watari M and Tsutahara M, Two-dimensional thermal model of the finite-difference latticeBoltzmann method with high spatial isotropy, Phys. Rev. E 67(3), 036306 (2003).[29] Gonnella G, Lamura A, and Sofonea V, Lattice Boltzmann simulation of thermal nonidealfluids, Phys. Rev. E 76(3), 036703 (2007).[30] Onuki A, Dynamic van der Waals Theory of Two-Phase Fluids in Heat Flow, Phys. Rev. Lett.94(5), 054501 (2015).[31] Gan Y, Xu A, Zhang G, and Li Y, FFT-LB Modeling of Thermal Liquid-Vapor System,Commun. Theor. Phys. 57(4), 681-694 (2012).[32] Gan Y, Xu A, Zhang G, and Li Y, Phase separation in thermal systems: A lattice Boltzmannstudy and morphological characterization, Phys. Rev. E 84(4), 046715 (2011).[33] Xu A, Zhang G, Gan Y, Chen F, and Yu X, Lattice Boltzmann modeling and simulation ofcompressible flows, Front. Phys. 7(5), 582-600 (2012).[34] Xu A, Zhang G, and Ying Y, Progess of discrete Boltzmann modeling and simulation ofcombustion system, Acta. Phys. Sin. 64, 184701 (2015).[35] Xu A, Zhang G, and Gan Y, Progress in studies on discrete Boltzmann modeling of phaseseparation process, Mech. Eng. 38, 361-374 (2016).[36] Xu A, Zhang G, Zhang Y. Discrete Boltzmann modeling of compressible flows. in: G. Z.Kyzas, A.C. Mitropoulos (Eds.), Kinetic Theory, InTech, Rijeka, 2018, Ch. 02.[37] Lin C and Luo K, Discrete Boltzmann modeling of unsteady reactive flows with nonequilibriumeffects, Phys. Rev. E, 99: 012142 (2019).[38] Gan Y, Xu A, Zhang G, Zhang Y, and Succi S. Discrete Boltzmann trans-scale modeling ofhigh-speed compressible flows, Phys. Rev. E 97(5), 053312 (2018).[39] Gan Y, Xu A, Zhang G and Succi S, Discrete Boltzmann modeling of multiphase flows: ydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11, 5336 (2015).[40] Zhang Y, Xu A, Zhang G, Gan Y, Chen Z and Succi S, Entropy production in thermal phaseseparation: a kinetic-theory approach, Soft Matter 15, 2245-2259 (2019).[41] Yan B, Xu A, Zhang G, Ying Y, and Li H, Lattice Boltzmann model for combustion anddetonation, Front. Phys. 8, 94-110 (2013).[42] Xu A, Lin C, Zhang G, and Li Y, Multiple-relaxation-time lattice Boltzmann kinetic modelfor combustion, Phys. Rev. E 91, 043306 (2015).[43] Lin C, Xu A, Zhang G, and Li Y, Double-distribution-function discrete Boltzmann model forcombustion, Combust. Flame 164, 137-151 (2016).[44] Zhang Y, Xu A, Zhang G, Zhu C, and Lin C, Kinetic modeling of detonation and effects ofnegative temperature coefficient, Combust. Flame 173, 483-492 (2016).[45] Lin C and Luo K, MRT discrete Boltzmann method for compressible exothermic reactiveflows, Comput. Fluids 166, 176-183 (2018).[46] Lin C, Luo K, Fei L, and Succi S, A multi-component discrete Boltzmann model for nonequi-librium reactive flows, Sci. Rep. 7, 14580 (2017).[47] Xu A, Zhang G, Zhang Y, Wang P, and Ying Y, Discrete Boltzmann model for implosion andexplosion related compressible flow with spherical symmetry, Front. Phys. 13, 135102 (2018).[48] Lai H, Xu A, Zhang G, Gan Y, Ying Y, and Succi S, Non-equilibrium thermohydrodynamiceffects on the Rayleigh-Taylor instability incompressible flow, Phys. Rev. E. 94(2), 023106(2016).[49] Chen F, Xu A, and Zhang G, Viscosity, heat conductivity, and Prandtl number effects in theRayleigh Taylor Instability, Front. Phys. 11(6), 114703 (2016).[50] Ye H, Lai H, Li D, Gan Y, Lin C, Chen L, and Xu A, Knudsen number effects on two-dimensional Rayleigh-Taylor instability in compressible fluid: based on a discrete Boltzmannmethod, Entropy 22, 500 (2020).[51] Gan Y, Xu A, Zhang G, Lai H, and Liu Z, Nonequilibrium and morphological characterizationsof Kelvin-Helmholtz instability in compressible flows, Front. Phys. 14(4), 43602 (2019).[52] Lin C, Xu A, Zhang G, Luo K, and Li Y, Discrete Boltzmann modeling of Rayleigh-Taylorinstability in two-component compressible flows, Phys. Rev. E 96, 053305 (2017).[53] Liu H, Kang W, Zhang Q, Zhang Y, Duan H, and He X, Molecular dynamics simulations ofmicroscopic structure of ultra strong shock waves in dense helium, Front. Phys. 11, 1152062016).[54] Liu H, Zhang Y, Kang W, Zhang P, Duan H, and He X, Molecular dynamics simulation ofstrong shock waves propagating in dense deuterium, taking into consideration effects of excitedelectrons, Phys. Rev. E 95, 023201 (2017).[55] Liu H, Kang W, Duan H, Zhang P, and He X, Recent progresses on numerical investigationsof microscopic structure of strong shock waves in fluid, Sci. Sinica Phys. Mech. Astron. 47,070003 (2017).[56] Meng J, Zhang Y, Hadjiconstantinou N, Radtke G, and Shan X, Lattice ellipsoidal statisticalBGK model for thermal non-equilibrium flows, J. Fluid Mech. 718, 347-370 (2013).[57] Gan Y, Xu A, Zhang G, and Succi S, Discrete Boltzmann modeling of multiphase flows:hydrodynamic and thermodynamic non-equilibrium effects. Soft Matter 11(26), 5336-5345(2015).[58] Shen Q, Rarefied gas dynamics: fundamentals, simulations and micro flows, Springer, 2005.[59] Chapman S, Cowling T, and Burnett D, The mathematical theory of non-uniform gases: an ac-count of the kinetic theory of viscosity, thermal conduction and diffusion in gases, Cambridge:Cambridge university press, 1990.[60] Guo Z and Zheng C, Theory and Applications of Lattice Boltzmann Method, Beijing: Sciencepress, 2008.[61] Bongiorno V, Davis H T, Modified Van der Waals theory of fluid interfaces, Physical ReviewA, 12(5), 2213-2224 (1975).[62] Huang H, Sukop M and Lu X. Multiphase Lattice Boltzmann Methods: Theory and Applica-tion. John Wiley & Sons, Inc, 2015.