{'fraction_non_alphanumeric': 0.03585224183810696, 'fraction_numerical': 0.01505050216029071, 'mean_word_length': 4.7445759368836296, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entropy stable modeling of non-isothermal multi-component diffuse-interface two-phase flows with realistic equations of state

In this paper, we consider mathematical modeling and numerical simulation of non-isothermal compressible multi-component diffuse-interface two-phase flows with realistic equations of state. A general model with general reference velocity is derived rigorously through thermodynamical laws and Onsager's reciprocal principle, and it is capable of characterizing compressibility and partial miscibility between multiple fluids. We prove a novel relation among the pressure, temperature and chemical potentials, which results in a new formulation of the momentum conservation equation indicating that the gradients of chemical potentials and temperature become the primary driving force of the fluid motion except for the external forces. A key challenge in numerical simulation is to develop entropy stable numerical schemes preserving the laws of thermodynamics. Based on the convex-concave splitting of Helmholtz free energy density with respect to molar densities and temperature, we propose an entropy stable numerical method, which solves the total energy balance equation directly, and thus, naturally satisfies the first law of thermodynamics. Unconditional entropy stability (the second law of thermodynamics) of the proposed method is proved by estimating the variations of Helmholtz free energy and kinetic energy with time steps. Numerical results validate the proposed method.

1. Introduction. Various non-isothermal multi-component two-phase flows are ubiquitous in nature and industry, and thus their research carries broad and farreaching significance. An industrial example is the phase transition of hydrocarbon mixtures in the reservoir; at specified thermodynamical conditions, a hydrocarbon mixture may split into gas and liquid (oil) to stay in an equilibrium state; when the thermal enhanced oil recovery is employed, intentionally introduced heat disrupts the equilibrium states and vaporizes part of the oil, thereby changing the physical properties such that oil flows more freely through the reservoir [11,37]. Another example is the utilization of supercritical fluids as solvents in chemical analysis and synthesis [43]. In the natural world, many common phenomena, such as boiling, evaporation and condensation, are also related to physical properties and motion of non-isothermal two-phase flows [42].
For realistic fluids, the interfaces between multiple fluids always exist and play a very important role in the mass and energy transfer between different phases. Partial miscibility of multiple fluids, a common phenomenon displayed by realistic fluids in experiments and practical processes, takes place through the interfaces. Moreover, capillarity effect, a significant mechanism of flows in porous media, is also caused by the anisotropic attractive force of molecules on the interfaces [30]. To describe the gasliquid interfaces, the diffuse-interface models for multiphase flows have been developed in the literature. The pioneering work is that the density-gradient contribution on the interfaces is introduced by van der Waals in the energy density (see [51] and the references therein), and on the base of it, the Korteweg stress formulation is induced by composition gradients (see [42,51] and the references therein). Various phase-field models for immiscible and incompressible two-phase flows have been developed and simulated in the literature, [1,2,5,7,9,12,14,15,17,23,28,29,40,58,64,65] for instance.
Modeling and simulation of compressible multi-component two-phase flows with partial miscibility and realistic equations of state (e.g. Peng-Robinson equation of state [54]) are intensively studied in recent years [19,31,32,35,39,55,57]. The multicomponent models with realistic equations of state are traditionally applied for simulation of many problems in chemical and petroleum engineering, for example, the phase equilibria calculations [20,26,27,33,36,45,46,49] and prediction of surface tension [20,30,38,47], but the fluid motion is never considered in these applications. The models of compositional fluid flows in porous media, for example, [25,44,48,56], describe the fluid motion through Darcy's law, but using the sharp interface. A general diffuse interface model for compressible multi-component two-phase flows with partial miscibility are developed in [34,35] based on the thermodynamic laws and realistic equations of state. It uses molar densities as the primal state variables, and takes a general thermodynamic pressure as a function of the molar density and temperature, thereby never suffering from the difficulty of constructing the pressure equation. However, the temperature field in this model is assumed to be homogeneous and constant.
There exist many situations, such as boiling, evaporation, condensation and thermal enhanced oil recovery, in which phase transitions and fluid motions are highly influenced by an inhomogeneous and variable temperature field. In [3], the nonisothermal diffuse-interface models for the single-component and binary fluids are developed by including gradient contributions in the internal energy. In [50,51], Onuki generalizes the van der Waals theory for the single-component fluids by including gradient contributions in both the internal energy and the entropy. Subsequently, improvements and applications of the model in [50,51] are investigated in [8,10,42,43,52,53,61] and the other literature; especially, in [42] a continuum mechanics modeling framework for liquid-vapor flows is rigorously derived using the thermodynamical laws. The non-isothermal diffuse-interface models are extended to the compressible binary fluids in [22,43]. The aforementioned research works are done on the basis of the van der Waals equation of state, but rarely concerning the other equations of state, for example, the Peng-Robinson equation of state [54] extensively employed in petroleum and chemical industries due to its accuracy and consistency for numerous realistic gas-liquid fluids including N 2 , CO 2 , and hydrocarbons. It is noted that different from the single-component fluids, a reference velocity, such as mass-average velocity, molar-average velocity and so on, usually needs to be selected for multi-component fluids [13], so the models are expected to be compatible with general reference velocity. However, up to now, the rigorous generalization of the aforementioned models to multi-component flows with general reference velocity is still an open problem.
In this paper, we will generalize the aforementioned model; more precisely, we will derive a general non-isothermal multi-component diffuse-interface two-phase model based on the thermodynamical laws and realistic equations of state (e.g. Peng-Robinson equation of state) with mathematical rigors. A significant feature of the general model is that it has a set of unified formulations for general reference veloci-ties and related mass diffusion fluxes. Moreover, a general thermodynamic pressure, which is a function of the molar density and temperature, is used and consequently, it is free of constructing the pressure equation.
The entropy balance equation plays a fundamental role in the derivations of diffuse interface two-phase flow models [3,42,43,50,51], by which we can apply the entropy production principle (the second law of thermodynamics) to derive the forms of thermodynamical fluxes, including the stress tensor, the mass diffusion and heat transfer fluxes. Different from the existing derivation approach using the Gibbs relation and other thermodynamical relations, we derive the entropy balance equation directly from the total energy balance equation (the first law of thermodynamics) based on general mass balance equations involving a general reference velocity. The transport equation of Helmholtz free energy density is derived to further reduce the entropy balance equation into a form composed of conservative terms and entropy production terms, from which we derive the forms of thermodynamical fluxes using the non-negativity principle of entropy production and Onsager's reciprocal principle. Consequently, the derived general model satisfies the first and second laws of thermodynamics and Onsager's reciprocal principle, thereby ensuring the thermodynamical consistency.
In the momentum conservation equations of the existing models, the pressure and surface tension are formulated as the primary driving force of the fluid motion except for the external forces. In this paper, we prove a novel relation among the pressure, temperature and chemical potentials, which leads to a new formulation of the momentum conservation equation indicating that the gradients of chemical potentials and temperature become the primary driving force except for the external forces. It will be shown that on the basis of the new formulation of momentum conservation equation, we can conveniently design efficient, easy-to-implement and entropy stable numerical algorithms.
A key and challenging issue in numerical simulations is how to design an algorithm preserving the first and second laws of thermodynamics obeyed by the model. To our best knowledge, there are only a few works regarding such algorithms in the literature. In [42], a provably entropy-stable numerical scheme has been developed fundamentally on the concept of functional entropy variables. Recently, in [37], using the convex-concave splitting of Helmholtz free energy density, the authors propose an entropy stable numerical method for the single-component fluids with the Peng-Robinson equation of state, and rigorously prove that the first and second laws of thermodynamics are preserved by this method. But only single-component flows are considered in these existing methods. In this paper, we focus on the numerical schemes preserving the laws of thermodynamics for the general multi-component flow model.
The proposed numerical scheme will be based on the convex splitting approach, which is first proposed in [16,18] and has been popularly employed in various phasefield models [6,18,24,59,63]. For isothermal single-component and multi-component diffuse-interface models with Peng-Robinson equation of state, the convex splitting schemes have been intensively studied recently [19,35,36,41,55,57], and very recently, a convex splitting scheme for the non-isothermal single-component fluids is also developed in [37]. But the convex splitting approach is never studied yet for non-isothermal multi-component fluids with Peng-Robinson equation of state. In this paper, we will develop the convex-concave splitting of Helmholtz free energy density based on Peng-Robinson equation of state; in particular, we will prove that this Helmholtz free energy density is always concave with respect to temperature.
In the algorithm proposed in [37], the internal energy equation is solved instead of the total energy balance equation. The proposed numerical algorithm in this paper will solve the total energy balance equation directly for ease of preserving the first law of thermodynamics. The key issue becomes how to gain entropy stability; in other words, the method shall be designed to satisfy the second law of thermodynamics. It is well known that the convex-concave splitting schemes usually result in the energydissipation feature for isothermal systems, whereas in this work we will prove that the convex-concave splitting of Helmholtz free energy density with respect to molar densities and temperature leads to the entropy stability of a numerical scheme for the non-isothermal systems. Another great challenge in designing the entropy-stable scheme is the very tightly coupling relationship among molar density, energy (temperature) and velocity. Such relationship will be well treated through very careful mathematical and physical observations. The rest of this paper is organized as follows. In Section 2, we will introduce the mathematical model for non-isothermal multi-component diffuse-interface two-phase flows, and rigorous derivations for this model will be provided in Section 3. In Section 4, we will propose an entropy stable numerical scheme based on the convex-concave splitting of Helmholtz free energy density and prove its entropy stability. Numerical results will be provided in Section 5 to validate the proposed method. Finally, some concluding remarks are given in Section 6.
2. Non-isothermal multi-component two-phase flow model. In this section, we present the system of equations modeling the non-isothermal multi-component two-phase flows with general reference velocity, which is composed of the mass balance equations, the momentum balance equation and total energy balance equation. Moreover, a new formulation of the momentum balance equation is derived through the relationship between the pressure, chemical potentials and temperature.

Notations and thermodynamical relations.
We consider a fluid mixture composed of M components in a variable temperature field. The temperature is denoted by T . Let n i denote the molar density of the ith component, and then we denote the molar density vector by n = [n 1 , n 2 , · · · , n M ] T .
The diffuse interfaces, which always exist between two phases in a realistic fluid, play an extremely important role in the heat and mass transfer between multiple phases. The key feature of diffuse interface models is to introduce the local density gradient contribution in the energy density of inhomogeneous fluids. The general Helmholtz free energy density, denoted by f , is expressed as where f b stands for the Helmholtz free energy density of a bulk fluid and c ij is the cross influence parameter generally relying on molar densities and temperature. The chemical potential of component i is defined as where δ δni represents the variational derivative. The entropy density, denoted by s, can be defined as We denote the internal energy density by ϑ. The internal energy, entropy, and temperature have the following relation The following relation between the pressure, Helmholtz free energy and chemical potential [35] holds for the bulk and inhomogeneous fluids which allows us to define the general thermodynamical pressure. According to (2.2), the chemical potential of component i can be deduced from (2.1) as . We get the formulation of the general pressure from (2.1), (2.5) and (2.6) as We denote by M w,i the molar weight of component i, and we further define the mass density of the mixture as Let g be the absolute value of the gravity acceleration and h is the height referred to a given reference platform.

Model equations.
We now state the system of model equations, which basically consists of the mass balance equations, the momentum balance equation and total energy conservation equation. First, the mass balance equation for component i is written as where u is a reference velocity and J i is the diffusion flux of component i. Let us define the stress tensor where I is the identity tensor, ξ is the volumetric viscosity, η is the shear viscosity and λ = ξ − 2 3 η. We assume that η > 0 and λ > 0. The momentum balance equation is where g = −g∇h. Thanks to the following relation (which will be proved in Sub- we obtain a new formulation of the momentum balance equation as which indicates that the fluid motion is driven by the gradients of chemical potentials and temperature. We denote the total energy density e t = ϑ + 1 2 ρ|u| 2 + ρgh. The total energy conservation equation is expressed as where q is the heat flux and The mass diffusion fluxes in (2.9) and the heat flux in (2.14) are assumed to be linearly related to ∇ µi T and ∇ 1 T as The mobility matrix L = (L i,j ) M+1 i,j=1 shall be symmetric in terms of Onsager's reciprocal principle. Moreover, the second law of thermodynamics requires that it shall be positive definite or positive semi-definite.
For the boundary conditions, we assume that all boundary terms will vanish when integrating by parts is performed; for example, we can use homogeneous Neumann boundary conditions or periodic boundary conditions.
We have the following comments on the above model. 1. The model has thermodynamically-consistent unified formulations for the general reference velocity and the mass diffusion and heat fluxes. 2. The model can characterize the compressibility, partial miscibility and heat transfer between different phases through the diffuse interfaces. 3. It is different from the case of the pure substance that there exist many choices of the reference velocity u and the diffusion flux J i for a multi-component mixture. Moreover, u and J i have the tightly dependent relations as shown below. 4. For the special case of J i taken such that .21) will vanish. However, in general cases, these terms are essential to ensure the thermodynamical consistency. 5. Due to the use of the general thermodynamic pressure, it is not necessary to construct the pressure equation, indeed, the pressure can be explicitly calculated from (2.7) if the molar density and temperature are available. We now present a few choices of the mass diffusion fluxes and heat flux. In general, these flux formulations shall be chosen such that the mobility matrix L is symmetric, positive definite or positive semi-definite. We take L i,M+1 = 0 for 1 ≤ i ≤ M , and furthermore, we take L M+1,M+1 = KT 2 , where K is the Fourier thermal conductivity coefficient of the mixture. As a result, we obtain the heat flux and the diffusion flux of component i as To determine the diffusion flux J i precisely, we give two typical choices of L ij (1 ≤ i, j ≤ M ), which correspond to molar-average velocity and mass-average velocity respectively.
(J1) In the first case, the diffusion mobility parameters are taken as where R stands for the universal gas constant and the mole diffusion coefficients D ij satisfy D ii = 0 and D ij = D ji > 0 for i = j. In this case, we have M i=1 J i = 0, which means that the reference velocity is the molar-average velocity.
(J2) In the second case, we take where D ij are the mass diffusion coefficients satisfying D ii = 0 and D ij = D ji > 0 for i = j. In this case, it holds that M i=1 M w,i J i = 0, so the reference velocity becomes the mass-average velocity. It is easy to prove that the above choices of L is symmetric positive semidefinite, and consequently, they obey Onsager's reciprocal principle [13] and the second law of thermodynamics.
3. Derivations of the model. In this section, we show the rigorous derivations of the model equations given in Section 2. The component mass balance equations (2.9) are assumed to hold with general reference velocity, but the mass diffusion fluxes are undetermined. The formulations of the mass diffusion fluxes, the momentum balance equation and total energy conservation equation will be derived using the laws of thermodynamics and Onsager's reciprocal principle.
3.1. Primary thermodynamical equations. In a time-dependent volume V (t), we define the internal energy and kinetic energy (E) within V (t) as The gravitational potential energy has the form We recall the first law of thermodynamics where W is the work done by the face force F t , and Q is the heat transfer from external environment of V (t). The work done by F t is expressed as Cauchy's relation between face force F t and the stress tensor σ of component i gives F t = −σ · ν, and as a result, where ν is the unit normal vector towards the outside of V (t). We note that the other external forces are ignored in this work, but the model derivations can be easily extended to the cases in the presence of additional external forces. The heat flux is expressed as where φ q represents the general heat flux. Applying the Reynolds transport theorem and the Gauss divergence theorem, we deduce that where we have also used the relation ϑ = f + T s. We can also derive that On the other hand, the following overall mass balance equation can be obtained from the component mass balance equations (2.9) where Du Dt = ∂u ∂t + u · ∇u. Substituting (3.4), (3.5), (3.6), (3.9) into (3.3), and taking into account the arbitrariness of V (t), we obtain where DT Dt is the material derivative as Furthermore, we rewrite (3.10) as (3.11)

Entropy equation.
In order to derive the model equations using the second law of thermodynamics, we will deduce an entropy equation, which consists of conservative and non-negative terms. For this purpose, we need to derive the transport equation of Helmholtz free energy density to reduce (3.11).
Firstly, using the component mass balance equations, we derive the transport equation of f b as For the gradient contribution of Helmholtz free energy density, we can derive and Combining (3.12)-(3.14), we deduce the transport equation of the Helmholtz free energy f as where we have also used the formulation of the general pressure p. Using the identity (3.26), we deduce that (∇n i · ∇n j ) ∇c ij Substituting (3.16) into (3.15), we obtain The transport equation of the gravity potential energy can be easily derived as (3.19)

Physical principles of the model derivations.
The second law of thermodynamics states that the total entropy can never decrease over time for an isolated system, so the non-conservative terms in (3.19) shall be non-negative. The fourth term on the right-hand side of (3.19) will vanish due to the relation s = −γ given in (2.3). Thanks to Galilean invariance, we must formulate the momentum conservation equation as the form of the equation (2.11) with an undetermined stress tensor, thereby making the last non-conservative term on the right-hand side of (3.19) disappear. In order to find the detailed formulation of the stress tensor, we consider the third term on the right-hand side of (3.19), which is a non-conservative term and thus shall be non-negative. This term shall vanish for the reversible process, and in this case, the total stress only contains the reversible part, i.e., the pressure pI. For the realistic viscous fluids, in terms of Newtonian fluid theory, the Cauchy stress tensor τ (u) given in (2.10) shall be included in the total stress, and consequently, we get σ = pI − τ (u), which ensures that the third term on the right-hand side of (3.19) is always non-negative So far the momentum balance equation reaches its complete form (2.11) with the stress tensor (2.10). W define the heat flux The fluxes q and J i may rely on ∇ µi T and ∇ 1 T according to Curie's Principle. Onsager's principle suggests that q and J i have a linear dependent relationship with ∇ µi T and ∇ 1 T and moreover, this linear relationship shall be symmetric. So there exists the linear, symmetric relationship given in (2.16). To obey the second law of thermodynamics, the sum of two related non-conservative terms shall be non-negative, i.e., (3.22) which requires that the mobility matrix L given in (2.16) must be positive definite or positive semi-definite. We rewrite (3.10) as which leads to the total energy conservation equation (2.14) taking into account (3.21). Finally, the entropy equation (3.19) is reduced into the following form The above model equations are established such that the entropy production (non-conservative) terms in (3.24) are all non-negative, thereby obeying the second law of thermodynamics.
3.4. The proof of the relation (2.12). Using the formulations of the pressure and chemical potentials, we can deduce where we have also used the following identity Consequently, we obtain the relation (2.12) taking into account s = −γ.

Entropy stable numerical method.
We have already shown that the proposed model obeys the first and second laws of thermodynamics. In this section, we focus on developing the semi-implicit time marching scheme preserving these laws of thermodynamics. The energy conservation equation will be solved in the proposed scheme, and thus, the first law of thermodynamics can be naturally satisfied. The key effort put in designing the efficient numerical scheme is how to preserve the second law of thermodynamics, which states that the total entropy can never decrease over time for an isolated system. Here, the entropy stability of a numerical scheme means that it obeys the second law of thermodynamics at the time-discrete level.
We will develop an entropy stable numerical method based on the use of an intermediate velocity and the convex-concave splitting of Helmholtz free energy density with respect to molar density and temperature. The convex-concave splitting approaches usually leading to the energy-dissipation feature for numerical simulation of the system with a constant temperature will be proved to result in the entropy stability of a numerical scheme for the non-isothermal systems. A great challenge in designing the entropy-stable scheme is the very tightly coupling relationship among molar density, energy (temperature) and velocity, the treatment of which requires very careful mathematical and physical observations. 4.1. Convex-concave splitting of Helmholtz free energy density. We analyze the convex-concave splitting of Helmholtz free energy density based on Peng-Robinson equation of state, which is widely used in the oil reservoir and chemical engineering. The gradient term of Helmholtz free energy density is always convex with respect to molar densities. For the bulk Helmholtz free energy density for the multi-component mixture, there are a few approaches proposed in [19,35,36] to get a strict convex-concave splitting scheme. Here, we follow the approach proposed in [35] to introduce a stabilization term which has a diagonal positive definite Hessian matrix, and thus is convex with respect to molar densities. Let us introduce a stabilization parameter θ ≥ 0, and then we reformulate the bulk Helmholtz free energy density as a sum of the convex part f convex b and the concave part  15 We now turn to consider the convex-concave property of Helmholtz free energy density with respect to the temperature. It is noted that ψ p i in the formulation of f b is the molar heat capacity of the ideal gas at the constant pressure. We recall the following thermodynamical relation for the ideal gas where R is the universal gas constant and ψ v i > 0 is the molar heat capacity of component i at the constant volume for the ideal gas.
Lemma 4.1. The second derivative of the bulk Helmholtz free energy density f b with respective to the temperature satisfies and thus f b is concave with respect to the temperature. Proof. Let n = M i=1 n i . Fixing molar densities, we calculate the second derivative of f b with respect to T as where we have used the relation (4.5). The first term on the right-hand side of (4.7) is negative. We note that the parameters a where T r,i = T /T c,i , T c,i is the critical temperature and m i is the parameter given in (A.2). It is obtained that A ij > 0 by the physical definitions of parameters, and thus we have a ′′ (T ) > 0 taking into account k ij < 1. Consequently, (4.6) is reached. For the gradient term of Helmholtz free energy density f ∇ , we take the constant influence parameters c ij in numerical tests as usual practices in the dynamical van der Waals model. This choice makes the derivatives of f ∇ with respect to the temperature zero when molar densities are fixed.
It is easy to prove that Helmholtz free energy density based on van der Waals equation of state has the property (4.6), and can also be split into the sum of a convex part and a concave part with respect to molar density. In the following, we assume that the influence parameters c ij may rely on the temperature but independent of molar densities and the convex-concave splitting can be reached for the Helmholtz free energy density considered in this paper; more precisely, f b can be written as the form (4.2a) and ∂ 2 f (n,T ) ∂T 2 n ≤ 0 always holds. 16 4.2. Semi-implicit time scheme. We now contruct the semi-implicit time marching scheme. We divide the time interval T = (0, t f ], where t f > 0 into N subintervals T k = (t k , t k+1 ], where t 0 = 0 and t N = t f . The time step size is denoted as δt k = t k+1 − t k . For a scalar function v(t) or a vector function v(t), we denote by v k or v k its approximation at the time t k .
For chemical potentials, we treat molar densities implicitly in their convex parts, while we apply the explicit treatment for their concave parts; meanwhile, the temperature is always implicit in both parts. More precisely, the time discrete chemical potential of component i at the (k + 1)th time step is expressed as where c k+1 ij = c ij (T k+1 ). The relation parameters in the mass diffusion fluxes and the heat flux may rely on molar densities and temperature, and here, we take the following schemes (4.9b) The above parameter matrix (L i,j ) M i,j=1 still ensures the symmetry and semi-positive definite property. We introduce an intermediate velocity as Based on the formulations of u k ⋆ and J k+1 i , we propose the following semi-implicit scheme for the mass balance equation of component i The momentum balance equation is discretized as where τ k+1 = λ k ∇ · u k+1 I + η k ε(u k+1 ). (4.13) The total energy at the kth time step is expressed as (4.14) We also introduce the intermediate total energy as The pressure at the (k + 1)th time step is calculated as The total energy conservation equation has the following discrete form 4.3. Unconditional entropy stability. We will derive an equality of entropy to prove that the proposed numerical scheme is unconditionally entropy stable.
The entropy equality will be deduced from the discrete total energy conservation equation through reducing the kinetic energies and Helmholtz free energies, and for this purpose, we need first to prove the following lemmas on the variations of the kinetic energies and Helmholtz free energies.
Lemma 4.2. The variation of kinetic energies between the kth and (k + 1)th time steps is estimated as Proof. For the kinetic energy at the (k + 1)th time step and the intermediate kinetic energy, we derive their difference (4.20) now we deduce the variation of gravity potential energy with time steps Using (4.19), (4.25) and (4.30), and taking into account we reduce the kinetic energies, Helmholtz free energies and gravity potential energies from the total energy equation to derive that Furthermore, we reduce (4.32) as which leads to (4.28). We now analyze the terms in the right-hand side of (4.28). The first three terms are the conservative terms, which do not produce the entropy, while the rest terms are the entropy production terms, which shall be non-negative according to the second law of thermodynamics. The non-negativity of the fourth and fifth terms is true due to the choice principles of q k+1 and J k+1 i . The last term is obviously non-negative. Consequently, the proposed scheme obeys the second law of thermodynamics. For the isolated system, the homogeneous Neumann boundary conditions are applied (that is, all terms vanish on the boundary when integration by parts is applied), and thus the inequality (4.29) can be obtained by integrating (4.28) over the domain Ω.

Numerical tests.
In this section, we will apply the proposed method to simulate non-isothermal multi-component two-phase flow problems. A binary mixture composed of methane (C 1 ) and pentane (C 5 ) is located in a square domain Ω with the length 20 nm. To validate the proposed method, we simulate the dynamics of the isolated systems, showing the entropy increase with time steps as proved in (4.29). The mass-average velocity is employed with the diffusion mobility formulations given by (2.19) taking the mass-diffusion coefficients D 12 = D 21 = 10 −8 m 2 /s. The volumetric viscosity and the shear viscosity are taken equal to zero. The heat diffusion coefficient is taken as K = n×10 −3 J/(m·s·K), where n = M i=1 n i . The initial liquid-phase molar densities are n L C1 = 6.8663 kmol/m 3 and n L C5 = 4.7915 kmol/m 3 respectively, while the initial gas-phase molar densities are n G C1 = 7.4302 kmol/m 3 and n G C5 = 0.6736 kmol/m 3 respectively. The gravity effect is ignored. The parameter θ = 0 is taken in (4.2). The rest physical parameters can be found in Appendix. We use a uniform     rectangular mesh with 40 × 40 elements. We employ the cell-centered finite difference method and the upwind scheme to discretize the mass balance equations and total energy conservation equation, and the finite volume method on the staggered mesh [62] for the momentum balance equation. These spatial discretization methods can be equivalent to the special mixed finite element methods with specified quadrature rules [4,21]. In this example, a square-shaped droplet is located in the center of the domain at the initial time. The initial molar density distributions of C 1 and C 5 are depicted in Figures 5.2(a) and 5.3(a) respectively. The initial temperature is uniformly taken equal to 310 K. The time step size is taken as 10 −12 s, and 60 time steps are simulated.
The total entropy profiles with time steps are depicted in Figure 5.1(a), while  Figure 5.5, we show the fluid motion including the velocity field and magnitudes of both velocity components. It is clearly observed that the initially square-shaped droplet is gradually changing to a circle, and subsequently, the mixture vapor is partially cooled and condenses into the liquid phase due to the lower temperature in the droplet region.   6. Conclusions. A general model obeying the thermodynamic laws and Onsager's reciprocal relations and characterizing compressibility and partial miscibility between multiple fluids has been developed for non-isothermal multi-component diffuse-interface two-phase flows with realistic equations of state. This model is formulised by a set of nonlinear and coupling equations, including the component mass balance equations, the momentum conservation equation and the total energy balance equation, which is unified for general reference velocities. We have proved an important relation between pressure, chemical potential and temperature, and on the basis of it, we obtain a new formulation of the momentum conservation equation, which indicates that the gradients of chemical potentials and temperature become the primary driving force of the fluid motion except for the external forces. For numerical simulations, we have proposed an efficient, entropy stable numerical method based on the convex-concave splitting of Helmholtz free energy density with respect to molar densities and temperature. We have proved unconditional entropy stability of the proposed method by deriving the variations of Helmholtz free energy and kinetic energy with time steps. Numerical results are provided to validate the proposed method. The above formulations are referred to [60]. We list some physical parameters of the substances in Tables C.1, C.2 and C.3.