Hydrodynamic equations for dense fluid mixtures with multistep interaction between particles

Generalization of the kinetic equation for a dense fluid with a multistep interaction potential to the case of a dense mixture is presented. Derivation of hydrodynamic equations is performed starting with the kinetic equation for the one-particle distribution function and the transport equation for the potential energy density. Expressions for momentum and heat fluxes are obtained. The transport equation for the kinetic energy density is shown to contain a new term of the “source” type, which describes fast exchange processes between kinetic and potential energies of the system.


Introduction
Construction of kinetic equations for dense fluids remains one of the main problems in the kinetic theory.Besides the difficulties characteristic of moderate and high densities and connected with the increasing role of interparticle correlations, that go beyond the pair collision approximation and take higher-order collisions into account, there appears another circumstance.When the average interparticle distance decreases, the interaction energy contribution to the total energy becomes appreciable.Apart from the short-range repulsion, the processes occurring between particles at distances of molecular attraction become important.Thus, a problem arises concerning the construction of the kinetic equations, which, through the corresponding collision integral could explicitly include microscopic processes due to the long-range attractive part of the interaction potential.
Neither the standard Enskog theory (SET) [1] nor its modifications (MET) [2,3] and (RET) [4,5] solve the problem as far as they are valid for hard spheres.The Enskog-type kinetic theories [6,7] suggested for a "hard spheres + smooth tail" potential treat the long-range interactions by a non-dissipative mean-field term.But it indirectly effects the transport coefficients through the pair distribution function and the contributions to thermodynamic quantities.
The square-well (SW) potential is the simplest one for which in the case of dense fluid the intermolecular attraction can be irreversibly included into the kinetic equation.For the first time this was performed in the so-called DRS theory [8].Though the transport coefficients were calculated, still substantial shortcomings consisted in disregarding the energy conservation law as well as in the absence of the H theorem.These deficiencies were overcome in the revised version RDRS [9].The H theorem was proved and the new kinetic theory was shown to give the correct equilibrium solution.The theory got its continuation in [10], where it was shown that the DRS and RDRS bulk viscosity coefficients differ from each other.Later, the spectrum of the linearized kinetic equation of RDRS was investigated and the new mode of the exchange type was discovered in [11].
The multistep (MS) interaction potential is more general and can better reproduce the course of the realistic interaction both at short and at long intermolecular distances.The corresponding kinetic equation for the one-component fluid was obtained in [12] based on the Zubarev's nonequilibrium statistical operator method starting with the Liouville equation.The H theorem was proved [12] and the transport coefficients for one-component fluid Ar along the gas-liquid saturation curve were calculated together with temperature and density dependencies [13,14].
Estimation of the effect of the soft attractive part on the transport coefficients in the high density region remains rather important.Two ways, i.e., indirect and direct, of dealing with the soft interaction can be distinguished.For the first one, the kinetic variational theory III, the stochastic theory, and the renormalized Kirkwood theory have been analyzed in a recent short review [15].Their capability of describing the self-diffusion coefficient as well as the bulk and shear viscosities for the one-component Lennard-Jones (LJ) fluid was investigated.In the semi-phenomenological approach of Morioka [16] a simple method was formulated for the shear viscosity calculation for fluids with a soft attraction, based on the combination of the Enskog hard sphere result and the evaluation of the average cross-section of the momentum transfer for the continuous potential.
For the second case, considerations involve time-correlation functions and memory-kernel methods starting with the Green-Kubo formulas rather than with an appropriate kinetic equation.In [17] the results for the self-diffusion coefficient of a dense LJ fluid are obtained in the binary collision approximation for the friction memory kernel.An approximate treatment of pair dynamics for the SW potential described by a Smoluchowski diffusion-like equation is used in evaluating the self-diffusion and shear viscosity coefficients [18].The deviation from the Stokes-Einstein law is observed if higher-particle time correlations are taken into account.
The SW interaction is applied to the calculation of the high-temperature correction to transport coefficients for the moderate density fluid with a hard core [19].Finally, both the square-well and square-barrier potentials have been used to mimic cohesive forces and soft material effect for inelastic particles in granular matter systems.In order to investigate the interplay between these forces of different nature, the 1D dense inelastic hard rod gas was studied by the kinetic theory and computer simulations [20].
Here, we consider a generalization of the kinetic equation for a system of particles with the multistep interaction to the case of dense fluid mixture and are concerned with the issue of deriving the hydrodynamic equations.In order to get an adequate kinetic description, the transport equation for the interaction energy density needs also to be considered as it was done in [9][10][11].However, in these works the proper attention was not paid to derivation of the kinetic energy density transport equation starting with the corresponding kinetic equation.The many-component version is of great interest since it permits to investigate diffusion and thermal diffusion processes as well as to elucidate the effect of the long-range interaction on the related transport coefficients, which is lacking in both SW and MS potentials.
The paper is organized as follows.In section 2 the kinetic equation is presented.In section 3 transport equations for the partial mass densities, average velocity, and kinetic energy density are derived.Section 4 is devoted to the potential energy density.In the last section we give some conclusions.

Kinetic equation
We will consider a many-component system of classical particles with M components.The particles interact by pair central forces as it holds for simple fluids.A realistic potential such as the Lennard-Jones type is approximated by a multistep function.It consists of the hard core part and a set of repulsive and attractive walls of finite heights (figure 1).The following parameters fully define the geometry of the MS potential: σ ij0 , σ r ijl , σ a ijl determine allocations of the hard-core, repulsive, and attractive walls, K r ij is the number of the repulsive walls (not including the hard core), K a ij is the number of the attractive walls.Parameters r ijl , a ijl denote the values of the MSP between walls, while characterize the heights of the walls and are introduced in such a way that ∆ r ijl , ∆ a ijl > 0. The plateaux are numbered in the direction of r ij increasing: from the hard core for the repulsive part and from the first attractive wall for the attractive one.For the bottom of the well we use two designations r ij,K r ij +1 and a ij0 .Besides, The distinctive feature of the interactions between two particles at the walls is that they are instantaneous: the interaction time tends to zero τ int → 0+.Due to this, ternary and higher-order processes at the walls are much less probable and we can restrict ourselves in the kinetic equation to the pair collision approximation.The kinetic equation for the one-component case was obtained in [12] based on the theoretical scheme for constructing the kinetic theory of dense fluid within the framework of the non-equilibrium statistical operator method by Zubarev.We generalize the theory to the mixture and write the corresponding kinetic equation for the one-particle distribution function of species i in the form: where The collision integral I E+MSP i consists of two contributions corresponding to the structure of the MS potential: where f ij 2 are the two-particle distribution functions.The Enskog type collision integral I E ij describes pair interactions at the hard core, while I MSP ij contributes due to the processes at repulsive and attractive walls of finite heights.
where v i , v j and v i , v j are pre-and postcollision velocities of particles i, j at the hard core contact, • σ, the unit vector σ is used to determine the mutual arrangement of the two particles and in the direct collision (v i , v j ) → (v i , v j ) it is directed from the center of particle j to the center of particle i, The contribution I MSP ij in equation ( 3) reads: where the term under summation describes the pair (ij)-processes at the wall of type q ("r" or "a") with number l.There are three types of such processes depending on the interplay between the value of the kinetic energy of relative motion and the height of the potential wall {q, l} as well as on the character of the relative motion (approaching or mutual removal): • p = ⊕ descending and acceleration at the wall (entering): each of two particles receives an additional momentum along the line of centers; the two momenta are directed one against another (for a wall of "a"-type) and one from another (for an "r"-type wall); the interaction energy decreases; • p = ascending and deceleration at the wall (escape): each particle loses a lot of momentum along the line of centers; the interaction energy increases; • p = ⊗ reflection from the wall (bound state): the kinetic energy of relative motion along the line of centers is not enough to overcome the potential barrier; both the kinetic and potential energies do not change.
The most inner sum p in equation ( 6) incorporates the described three types of processes.Together with the hard core collision, thereafter denoted as (E), these pair processes can be divided into two subsets correspondingly to either the two-particle both kinetic and interaction energies conserve separately or the exchange between them occurs: nonexchange (E, ⊗) and exchange processes (⊕, ).
For convenience we will ascribe formal numerical values to the introduced symbolic values of parameters q and p: This permits to present the collision integrals I qp ijl for all p and q in the following compact form.The reflection processes are described with and corresponds to the descending or ascending processes if the appropriate value of p is chosen.In these expressions numerical values of p and q are used to determine the position of particle j for the function f ij 2 , the argument of the θ function in equation ( 9), and to fix the type of the limiting value (right or left) for the function f ij 2 (..) −q in expression (8), and f ij 2 (..) ∓qp in expression (9).In all other cases the parameters q and p are symbols used for designation.Owing to the discontinuity of the pair potential each f ij 2 is also discontinuous in the coordinate space.Its right (+) and left (−) limiting values in equations ( 4), ( 8), (9) mean: The other symbols introduced in the above expressions mean: σ q ijl = σ q ijl σ, v i , v j are velocities of particles i and j just after the reflection from the q-type wall with number l which are to be determined from the same relations (5) as for the hard core collisions; v qp il , v qp jl are velocities just after the process of type p at wall {q, l} and can be determined from the pair collision rules for processes of the exchange type (p = ⊕, ): where v q ijl = (2∆ q ijl /µ ij ) 1/2 is the wall height in velocity units, µ ij = m i m j /(m i + m j ) is the reduced mass.
To a certain degree the instantaneous processes at the walls correspond to noncompleted scattering trajectories in the region of smooth behaviour of the intermolecular potential of the real fluid.To a greater extent this is associated with the sloping attractive part.Instead of dealing with these complicated contributions from the noncompleted trajectories we effectively replace them with the instantaneous processes at the walls.From this viewpoint the MS potential dominates over the square-well potential since it treats the long-range attraction effects more accurately.
As the collision integrals contain such pair processes, when during the typical time interval of changing of the one-particle distribution function, the two particles i and j are allowed to interact only at one wall, while successive interactions at two or more neighbouring walls which are described with more complicated terms [21,22] are out of consideration, the following condition must be satisfied [13,14]: where l free is the mean free path, ∆σ is the minimum separation between two walls ∆σ = min |σ q ij,l+1 − σ q ijl |, i, j = 1 ÷ M, q = r, a, l = 1 ÷ K q ij .The density restriction appears from the estimation for l free given, for example, for the one-component dense fluid in [13,14] and results in the possibility of applying the introduced kinetic equation to high densities only.The estimation (13) remains valid for the case of mixture as well.
If the number of walls increases and at the same time separations between them decrease, then the MSP approximates the realistic potential better and better.Eventually, condition ( 12) is violated and we cannot proceed further but are forced to seek an optimal solution to the choice of the MS potential parameters.Nevertheless, if we formally consider such a limiting behaviour of the theory for φ MSP ij → φ tail ij , where φ tail ij is, for example, the "hard spheres + smooth tail" interaction, then as shown in [14], the irreversible MSP collision integral is reduced to the well-known term of the mean-field type [6].Works [21,22] present the analysis of successive processes of the "attractive wall -hard core -attractive wall" type and some others for the SW case as well as clears up their effects on transport coefficients of moderately dense gas.

Hydrodynamic level of description
Local densities of the conserved quantities, mass, momentum, and total energy form the basis of the hydrodynamic level.In turn, these are equivalent to the partial mass densities, ρ i (r, t), the average mass velocity, V(r, t), and the internal energy density, e(r, t).The latter is equal to the total energy e(r, t) minus the convective part and has two contributions: e(r, t) = e k (r, t) + e p (r, t), (14) where e k is the kinetic energy density in the local reference system, e p is the density of the potential energy of interaction.
Define total densities ρ(r, t), ρ(r, t)V(r, t), e k (r, t) through the corresponding partial densities, which are first moments of f i (r, v i , t): The transport equations for the introduced variables can be obtained in a general form by differentiating definitions (15) with t and using equation (2): where is a notation of averaging.Each term in the left-hand side of equation ( 16) can be immediately averaged for each ψ i and we receive for total quantities: where After summing over all species the right-hand side (r.h.s.) of equation ( 16) reads: Further transformations of this expression having, for example, for where f ij 2 denotes f ij 2 (r, v i , r + σ ij0 , v j ) + (see equation ( 4)), consist in reducing to a divergence of some flux: Here ψ i denotes the corresponding postcollision value of ψ i .The same should be done for all the processes at walls in equation (19).Quantities under the divergence are related to the fluxes of momentum and kinetic energy.The described transformations are carried out due to the well-known symmetrization procedure with the use of the momentum and energy conservation laws.
For the unexchange processes all the details of derivation are similar to the case of the Enskog kinetic equation, e.g.see [6], and are presented elsewhere [23].The final results for the r.h.s. of equation ( 19) for the hard-core and reflection processes are: For the exchange processes, the symmetrization scheme has some peculiarities related to the appearance of the source term, Appendix A. The final result reads: where the source term s k is of the form: Expressions for all the contributions to P and q are given in Appendix B, equations ( 49), ( 50), (51).
Collecting the contributions to the r.h.s., given by equations ( 22) and ( 23) we finally get the system of transport equations for the one-particle macroscopic quantities: where J md i = ρ i (V i − V), P MSP = P ⊗ + P ex , q MSP = q ⊗ + q ex .The appearance of the source term s k in the last transport equation was not stressed in any of the previous articles dealing with both the SW interaction [8][9][10][11]24] and the MS potential [12][13][14].It is related to accounting for the processes at the walls of finite height in an irreversible way and describes the exchange rate between the kinetic and potential energy densities.Therefore, s k has no analogue both in the case of the Enskog equation including its mean-field extensions [6,7], and in the case of the Boltzmann kinetic equation.In the next section the relation of this term to the transport equation for e p is considered.

Interaction energy density and the closure relation
The density of the potential energy of interaction e p (r, t) is expressed by means of the twoparticle distribution functions f ij 2 : Its transport equation can be derived from the second equation of the BBGKY hierarchy with the appropriate choice of the collision operators which correspond to the MS potential.It reads: where q p (r, t) is the peculiarly molecular or "heat" flux of the potential energy density.The source s p (r, t) in the r.h.s.like describes the change of e p due to collisions.However, only the exchange-type processes p = ⊕, contribute.s p (r, t) is identical (with the opposite sign) to its kinetic counterpart (24): The result for s p is also obtainable using heuristic ideas about the number of collisions [23] and is to be considered as a direct generalization to the MSP case of that used in [9][10][11] in which, however, a connection to the transport equation for e k was not mentioned.Due to the observed property the transport equation for the total energy density does not contain any source term as it should be for a density of a conserved quantity.
In order to close the equations for f i and e p a closure relation should be supplied.This can be done in the spirit of [9] (also see [10,11,24]) using the approximation which ignores pair correlations in the velocity space: where g ij 2 is a functional of the local number densities n k (r, t) and the reciprocal quasi-temperature β(r, t) (also called the inverse potential energy temperature) so that g ij 2 has the same cluster expansion (n-vertex, f -bond) as the equilibrium counterpart.But in nonequilibrium case n k (r, t) replaces each n k and β ij (r i , r j , t) = [β(r i , t) + β(r j , t)]/2 replaces 1/k B T at each bond.β(r, t) is a Lagrange multiplier conjugated to the potential energy density [9,11,24] and is treated in the theory using the transport equation for e p (r, t).The functional g ij 2 is discontinuous at each point of discontinuity of the MSP and obeys: Thus, the equations for f i and e p become closed.The closure problem in the kinetic theory of fluids with attraction was analyzed in detail in [7] within the framework of the approach of maximization of entropy subject to some constraints (see also [6,25]).Some useful analysis is given in [26,27].
In approximation (30) the expression for q p reduces to where is the number density flux.Under assumption (30) the heat flux q p is of intrinsically diffusive nature and vanishes for a one-component fluid q p | M =1 = 0. Perhaps, this result can be used as another way of varifying the applicability of assumption (30) specifically for the systems with the MS potential, when the comparison of theoretical and computer simulation results is made.
If we subtract the convective term (ρV 2 )/2 from equation (25) and combine the result with equation ( 27), the equation for the internal energy density is received: with the extra contribution q p contrary to the total stress tensor.The right-hand side term describes the rate of internal energy change due to the processes of the viscous molecular friction ( † designates the transpose of ∇V).

Summary
We have presented the kinetic equation for dense fluid mixtures with the multistep potential of interaction.Together with the transport equation for the potential energy density it generalizes the corresponding one-component case [12][13][14] as well as the kinetic equation for SW particles [8,9], if only one attractive wall is retained.If all the wall heights are set equal to zero, we reproduce the kinetic equation of RET for mixtures [4].Each of the transport equations for e k and e p contains the source-type term, s k and s p , which are identical with each other but differ in the sign and describe the energy exchange phenomena in a non-equilibrium state but should vanish in equilibrium.
Further, we denote variables v qp il , v qp jl by v a , v b .Then ψ i (r, v i ) which means "the value of the molecular quantity before the process {qlp}" (e.g.before descending at the wall: p = +1), has, after denoting in terms of the new variables v a , v b , the opposite sense -"the value after the process {ql p}" (respectively, after ascending at the wall: p = −1).The designation must be changed as follows: ψ i (r, v i ) → ψ q p il (r, v q p al ).In terms of the new variables equation (37) becomes: The received expression substantially differs in arguments of the functions θ and f ij 2 from R qp0 ijl in equation ( 34) and should be combined with that from the process of p type at the same wall, that is to say with R q p0 ijl .This is the very circumstance that makes different the two symmetrization procedures for the nonexchange and exchange processes, arising due to their different symmetry properties.In an nonexchange process (E or ⊗) the inversion of the relative velocity normal component occurs while for an exchange one this component keeps its direction.This yields: Having in mind this special feature we must rearrange the direct (0) and inverse ( * ) contributions from processes p = ±1 in accordance with the following rule: Now, carry out the symmetrization with respect to indices i, j for the expression which is obtained from equation (39) by the replacement of p with p.To this end we perform the following manipulations: a) permute indices i → ← j and corresponding arguments of f ji 2 ; b) change the designations of velocities v a → ← v b ; c) choose σ = −σ as a new variable and denote it again by σ.Then, the term with permuted ji indices reads: Here f ij 2 describes such a configuration in which position r is occupied by particle j with velocity v b while particle i with velocity v a lies in the position r + qpσ q ijl .The conservation laws with corresponding spatial coordinates read: where respectively for the mass, momentum, and energy.Substitute equation (43) into equation ( 42) and take half of the sum of the obtained result and equation (41): × [ψ qp il (r, v qp al ) − ψ i (r, v a )]f ij 2 (r, v a , r − qpσ q ijl , v b ) qp − [ψ qp il (r + qpσ q ijl , v qp al ) − ψ i (r + qpσ q ijl , v a )]f ij 2 (r + qpσ q ijl , v a , r, v b ) qp + 1 2 (−∆ψ qp ijl )(σ q ijl ) 2 dv a dv b dσ v baσ θ(v baσ + p−1 2 v q ijl )f ij 2 (r + qpσ q ijl , v a , r, v b ) qp .(45) The first term in the r.h.s.can be written as a divergence, if we present the difference of the products within the curly brackets as The resulting expression reads: [ij+ji]l = −∇• 1 2 qp(σ q ijl ) 3 dv i dv j dσ v jiσ θ(v jiσ + p−1 2 v q ijl )σ × f ij 2 (r + λqpσ q ijl , v i , r + (λ − 1)qpσ q ijl , v j ) qp + 1 2 (−∆ψ qp ijl )(σ q ijl ) 2 dv i dv j dσ v jiσ θ(v jiσ + p−1 2 v q ijl )f ij 2 (r, v i , r − qpσ q ijl , v j ) qp , (48) where in the second term we have returned to the previous presentation of f ij 2 , carrying out the inverse manipulations at indices, velocities, and σ.Just this term forms the source s k .The summation ij ql p=±1 gives the final results (23).

B. Contributions to fluxes
Here, expressions for the nonexchange (E, ⊗) and exchange (ex) contributions to the stress tensor and the heat flux are given which emerge after the symmetrization and appear for the first time in (22), (23).