Entropy stable modeling of non-isothermal multi-component diffuse-interface two-phase flows with realistic equations of stateThis work is supported by KAUST research fund to the Computational Transport Phenomena Laboratory at KAUST.

# Entropy stable modeling of non-isothermal multi-component diffuse-interface two-phase flows with realistic equations of state††thanks: This work is supported by KAUST research fund to the Computational Transport Phenomena Laboratory at KAUST.

Jisheng Kou School of Mathematics and Statistics, Hubei Engineering University, Xiaogan 432000, Hubei, China.    Shuyu Sun Corresponding author. Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia. Email: shuyu.sun@kaust.edu.sa.
###### Abstract

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.

Key words. Multi-component two-phase flow; Non-isothermal flow; Entropy stability; Convex splitting.

AMS subject classifications. 65N12; 76T10; 49S05

## 1 Introduction

Various non-isothermal multi-component two-phase flows are ubiquitous in nature and industry, and thus their research carries broad and far-reaching 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 gas-liquid 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, [9, 17, 1, 5, 28, 2, 23, 14, 58, 40, 15, 65, 64, 29, 7, 12] 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 [57, 31, 32, 19, 35, 39, 55]. The multi-component 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 [26, 27, 45, 46, 33, 36, 49, 20] and prediction of surface tension [47, 30, 38, 20], but the fluid motion is never considered in these applications. The models of compositional fluid flows in porous media, for example, [56, 48, 25, 44], 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 non-isothermal 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 [42, 53, 8, 61, 10, 43, 52] 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 [43, 22]. 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, CO, 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 velocities 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, 50, 51, 43, 42], 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 phase-field models [59, 63, 18, 6, 24]. 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 [57, 19, 35, 36, 55, 41], 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 energy-dissipation 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.

### 2.1 Notations and thermodynamical relations

We consider a fluid mixture composed of components in a variable temperature field. The temperature is denoted by . Let denote the molar density of the th component, and then we denote the molar density vector by .

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 , is expressed as

 f(n,T)=fb(n,T)+f∇(n,T), \hb@xt@.01(2.1a) f∇(n,T)=12M∑i,j=1cij(n,T)∇ni⋅∇nj, \hb@xt@.01(2.1b)

where stands for the Helmholtz free energy density of a bulk fluid and is the cross influence parameter generally relying on molar densities and temperature.

The chemical potential of component is defined as

 μi=(δf(n,T)δni)T,n1,⋯,ni−1,ni+1,⋯,nM,  i=1,⋯,M, \hb@xt@.01(2.2)

where represents the variational derivative. The entropy density, denoted by , can be defined as

 s=−(δfδT)n. \hb@xt@.01(2.3)

We denote the internal energy density by . The internal energy, entropy, and temperature have the following relation

 ϑ=f+sT. \hb@xt@.01(2.4)

The following relation between the pressure, Helmholtz free energy and chemical potential [35] holds for the bulk and inhomogeneous fluids

 p=M∑i=1μini−f, \hb@xt@.01(2.5)

which allows us to define the general thermodynamical pressure.

According to (2.2), the chemical potential of component can be deduced from (2.1) as

 μi=μbi−M∑j=1∇⋅(cij∇nj)+12N∑j,k=1∂cjk∂ni∇nj⋅∇nk, \hb@xt@.01(2.6)

where . We get the formulation of the general pressure from (2.1), (2.5) and (2.6) as

 p = pb−M∑i,j=1ni∇⋅(cij∇nj)+12N∑i,j,k=1ni∂cjk∂ni∇nj⋅∇nk \hb@xt@.01(2.7) −12M∑i,j=1cij∇ni⋅∇nj,

where

We denote by the molar weight of component , and we further define the mass density of the mixture as

 ρ=M∑i=1niMw,i. \hb@xt@.01(2.8)

Let be the absolute value of the gravity acceleration and is the height referred to a given reference platform.

### 2.2 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 is written as

 ∂ni∂t+∇⋅(uni)+∇⋅Ji=0, \hb@xt@.01(2.9)

where is a reference velocity and is the diffusion flux of component .

Let us define the stress tensor

 σ=pI−τ,  τ(u)=(λ∇⋅u)I+ηε(u),   ε(u)=∇u+∇uT, \hb@xt@.01(2.10)

where is the identity tensor, is the volumetric viscosity, is the shear viscosity and . We assume that and . The momentum balance equation is

 ρ(∂u∂t+u⋅∇u)+M∑i=1Mw,iJi⋅∇u=−∇⋅σ−M∑i,j=1∇⋅(cij∇ni⊗∇nj)+ρg, \hb@xt@.01(2.11)

where . Thanks to the following relation (which will be proved in Sub-section 3.4)

 ∇p+M∑i,j=1∇⋅(cij∇ni⊗∇nj)=M∑i=1ni∇μi+s∇T, \hb@xt@.01(2.12)

we obtain a new formulation of the momentum balance equation as

 ρ(∂u∂t+u⋅∇u)+M∑i=1Mw,iJi⋅∇u=−M∑i=1ni∇μi−s∇T+∇⋅τ(u)+ρg, \hb@xt@.01(2.13)

which indicates that the fluid motion is driven by the gradients of chemical potentials and temperature.

We denote the total energy density . The total energy conservation equation is expressed as

 ∂et∂t+∇⋅(uet+σ⋅u)=−∇⋅(q−π), \hb@xt@.01(2.14)

where is the heat flux and

 π=M∑i,j=1cij∂ni∂t∇nj−12M∑i=1Mw,i|u|2Ji. \hb@xt@.01(2.15)

The mass diffusion fluxes in (2.9) and the heat flux in (2.14) are assumed to be linearly related to and as

 Ji=−M∑j=1Li,j∇μj+Mw,jghT+Li,M+1∇1T, \hb@xt@.01(2.16a) q=−M∑j=1LM+1,j∇μj+Mw,jghT+LM+1,M+1∇1T. \hb@xt@.01(2.16b)

The mobility matrix 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 and the diffusion flux for a multi-component mixture. Moreover, and have the tightly dependent relations as shown below.

4. For the special case of taken such that , the terms in (2.11) and in (3.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 is symmetric, positive definite or positive semi-definite. We take for , and furthermore, we take , where is the Fourier thermal conductivity coefficient of the mixture. As a result, we obtain the heat flux and the diffusion flux of component as

 q=−K∇T,   Ji=−M∑j=1Lij∇μj+Mw,jghT,  i=1,⋯,M. \hb@xt@.01(2.17)

To determine the diffusion flux precisely, we give two typical choices of , which correspond to molar-average velocity and mass-average velocity respectively.

(J1)

In the first case, the diffusion mobility parameters are taken as

 Lii=M∑j=1DijninjnR,     Lij=−DijninjnR,  j≠i,  1≤i,j≤M, \hb@xt@.01(2.18)

where stands for the universal gas constant and the mole diffusion coefficients satisfy and for . In this case, we have , which means that the reference velocity is the molar-average velocity.

(J2)

In the second case, we take

 Lii=M∑j=1DijniρjMw,iρR,     Lij=−DijninjρR,  j≠i,  1≤i,j≤M, \hb@xt@.01(2.19)

where are the mass diffusion coefficients satisfying and for . In this case, it holds that , so the reference velocity becomes the mass-average velocity.

It is easy to prove that the above choices of 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 , we define the internal energy and kinetic energy () within as

 U=∫V(t)ϑdV,   E=12∫V(t)ρ|u|2dV. \hb@xt@.01(3.1)

The gravitational potential energy has the form

 H=∫V(t)ρghdV. \hb@xt@.01(3.2)

We recall the first law of thermodynamics

 d(U+E+H)dt=d_Wdt+d_Qdt, \hb@xt@.01(3.3)

where is the work done by the face force , and is the heat transfer from external environment of . The work done by is expressed as

 d_Wdt=∫∂V(t)Ft⋅uds.

Cauchy’s relation between face force and the stress tensor of component gives , and as a result,

 d_Wdt=−∫∂V(t)(σ⋅ν)⋅uds=−∫V(t)∇⋅(σ⋅u)dV, \hb@xt@.01(3.4)

where is the unit normal vector towards the outside of . 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

 d_Qdt=−∫∂V(t)ϕq⋅νds=−∫V(t)∇⋅ϕqdV, \hb@xt@.01(3.5)

where represents the general heat flux.

Applying the Reynolds transport theorem and the Gauss divergence theorem, we deduce that

 dUdt = ∫V(t)∂ϑ∂tdV+∫V(t)∇⋅(uϑ)dV \hb@xt@.01(3.6) = ∫V(t)(∂f∂t+∇⋅(uf))dV+∫V(t)T(∂s∂t+∇⋅(us))dV +∫V(t)s(∂T∂t+u⋅∇T)dV,

where we have also used the relation . We can also derive that

 dEdt = 12∫V(t)∂(ρu⋅u)∂tdV+12∫V(t)∇⋅(u(ρu⋅u))dV \hb@xt@.01(3.7) = ∫V(t)(ρu⋅∂u∂t+12u⋅u∂ρ∂t)dV +12∫V(t)((ρu⋅u)∇⋅u+(u⋅u)u⋅∇ρ+ρu⋅∇(u⋅u))dV = ∫V(t)(ρu⋅∂u∂t+12u⋅u∂ρ∂t)dV +12∫V(t)((u⋅u)∇⋅(ρu)+2ρu⋅(u⋅∇u))dV = ∫V(t)ρu⋅(∂u∂t+u⋅∇u)dV+12∫V(t)u⋅u(∂ρ∂t+∇⋅(ρu))dV.

On the other hand, the following overall mass balance equation can be obtained from the component mass balance equations (2.9)

 ∂ρ∂t+∇⋅(ρu)+M∑i=1Mw,i∇⋅Ji=0. \hb@xt@.01(3.8)

Substituting (3.8) into (3.7) yields

 dEdt = ∫V(t)u⋅(ρDuDt+M∑i=1Mw,iJi⋅∇u)dV \hb@xt@.01(3.9) −12∫V(t)M∑i=1Mw,i∇⋅(|u|2Ji)dV,

where

Substituting (3.4), (3.5), (3.6), (3.9) into (3.3), and taking into account the arbitrariness of , we obtain

 +u⋅(ρDuDt+M∑i=1Mw,iJi⋅∇u)−12M∑i=1Mw,i∇⋅(|u|2Ji) =−∇⋅ϕq−∇⋅(σ⋅u), \hb@xt@.01(3.10)

where is the material derivative as

 DTDt=∂T∂t+u⋅∇T.

Furthermore, we rewrite (3.1) as

 ∂s∂t+∇⋅(us) = −1T∇⋅ϕq−1TσT:∇u−sTDTDt \hb@xt@.01(3.11) −1T(∂(f+ρgh)∂t+∇⋅(u(f+ρgh))) +12TM∑i=1Mw,i∇⋅(|u|2Ji) −uT⋅(ρDuDt+M∑i=1Mw,iJi⋅∇u+∇⋅σ).

### 3.2 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).

We define , where and

 γ∇=(δf∇δT)n=12M∑i,j=1∂cij∂T∇ni⋅∇nj.

Firstly, using the component mass balance equations, we derive the transport equation of as

 ∂fb∂t+∇⋅(fbu) = M∑i=1μbi∂ni∂t+γb∂T∂t+fb∇⋅u+u⋅∇fb \hb@xt@.01(3.12) = −M∑i=1μbi(∇⋅(niu)+∇⋅Ji)+γb∂T∂t +fb∇⋅u+M∑i=1u⋅μbi∇ni+u⋅γb∇T = −M∑i=1(μbini∇⋅u+u⋅μbi∇ni+μbi∇⋅Ji) +γbDTDt+fb∇⋅u+M∑i=1u⋅μbi∇ni = −pb∇⋅u−M∑i=1μbi∇⋅Ji+γbDTDt.

For the gradient contribution of Helmholtz free energy density, we can derive

 ∂f∇∂t = 12∂(∑Mi,j=1cij∇ni⋅∇nj)∂t \hb@xt@.01(3.13) = 12M∑i,j=1∂cij∂T∂T∂t∇ni⋅∇nj+12N∑i,j,k=1∂cjk∂ni∂ni∂t∇nj⋅∇nk +M∑i,j=1cij∇∂ni∂t⋅∇nj = γ∇∂T∂t−12N∑i,j,k=1∂cjk∂ni(∇⋅(niu)+∇⋅Ji)∇nj⋅∇nk −M∑i,j=1cij∇nj⋅∇(∇⋅(niu)+∇⋅Ji) = γ∇∂T∂t−M∑i,j=1∇⋅((∇⋅(uni))cij∇nj) −12N∑i,j,k=1∂cjk∂ni(ni∇⋅u+u⋅∇ni+∇⋅Ji)∇nj⋅∇nk +M∑i,j=1∇⋅(uni)∇⋅(cij∇nj)−M∑i,j=1cij∇nj⋅∇(∇⋅Ji) = γ∇∂T∂t−M∑i,j=1∇⋅((∇⋅(uni))cij∇nj) −12N∑i,j,k=1∂cjk∂ni(ni∇⋅u+u⋅∇ni+∇⋅Ji)∇nj⋅∇nk +M∑i,j=1ni(∇⋅u)∇⋅(cij∇nj)+M∑i,j=1(u⋅∇ni)∇⋅(cij∇nj) −M∑i,j=1∇⋅((∇⋅Ji)cij∇nj)+M∑i,j=1(∇⋅Ji)∇⋅(cij∇nj),

and

 ∇⋅(f∇u) =12∇⋅(uM∑i,j=1cij∇ni⋅∇nj) =12(M∑i,j=1cij∇ni⋅∇nj)∇⋅u+12u⋅∇(M∑i,j=1cij∇ni⋅∇nj). \hb@xt@.01(3.14)

Combining (3.12)-(3.2), we deduce the transport equation of the Helmholtz free energy as

 ∂f∂t+∇⋅(fu) = ∂fb∂t+∇⋅(fbu)+∂f∇∂t+∇⋅(f∇u) = γbDTDt−pb∇⋅u−M∑i=1μbi∇⋅Ji +γ∇∂T∂t−M∑i,j=1∇⋅((∇⋅(uni))cij∇nj) −12N∑i,j,k=1∂cjk∂ni(ni∇⋅u+u⋅∇ni+∇⋅Ji)∇nj⋅∇nk +M∑i,j=1(∇⋅u)ni∇⋅(cij∇nj)+M∑i,j=1(u⋅∇ni)∇⋅(cij∇nj) −M∑i,j=1∇⋅((∇⋅Ji)cij∇nj)+M∑i,j=1(∇⋅Ji)∇⋅(cij∇nj) +12(M∑i,j=1cij∇ni⋅∇nj)∇⋅u+12u⋅∇(M∑i,j=1cij∇ni⋅∇nj) = −p ∇⋅u+γbDTDt−M∑i=1μi∇⋅Ji+γ∇∂T∂t +M∑i,j=1∇⋅(∂ni∂tcij∇nj)−12N∑i,j,k=1(u⋅∇ni)∂cjk∂ni∇nj⋅∇nk +M∑i,j=1(u⋅∇ni)∇⋅(cij∇nj)+12u⋅∇(M∑i,j=1cij∇ni⋅∇nj),

where we have also used the formulation of the general pressure . Using the identity (3.26), we deduce that

 M∑i,j=1(∇ni)∇⋅(cij∇nj)+12∇(M∑i,j=1cij∇ni⋅∇nj) +12M∑i,j=1cij∇(∇ni⋅∇nj) =γ∇∇T+12N∑i,j,k=1(∇ni)∂cjk∂ni∇nj⋅∇nk+M∑i,j=1∇⋅cij(∇ni⊗∇nj). \hb@xt@.01(3.16)

Substituting (3.2) into (3.2), we obtain

 ∂f∂t+∇⋅(fu) = −p ∇⋅u+γDTDt−M∑i=1μi∇⋅Ji+M∑i,j=1∇⋅(∂ni∂tcij∇nj) \hb@xt@.01(3.17) +u⋅(M∑i,j=1∇⋅cij(∇ni⊗∇nj)).

The transport equation of the gravity potential energy can be easily derived as

 ∂ρgh∂t+∇⋅(ρghu) = gh(∂ρ∂t+∇⋅(ρu))+uρg∇h \hb@xt@.01(3.18) = −M∑i=1Mw,igh∇⋅Ji−ρg⋅u.

Substituting (3.17) and (3.18) into (3.11), we reduce the entropy balance equation as

 ∂s∂t+∇⋅(us)=−∇⋅(ϕqT)+ϕq⋅∇1T−1T(σ−pI):∇u−s+γTDTDt +12M∑i=1Mw,i∇⋅(1T|u|2Ji)−12M∑i=1Mw,i|u|2Ji⋅∇1T +M∑i=1∇⋅(μi+Mw,ighTJi)−M∑i,j=1∇⋅(1T∂ni∂tcij∇nj) +M∑i,j=1(cij∂ni∂t∇nj)⋅∇1T−M∑i=1Ji⋅∇μi+Mw,ighT −uT⋅(ρDuDt+M∑i=1Mw,iJi⋅∇u+∇⋅σ+M∑i,j=1∇⋅cij(∇ni⊗∇nj)−ρg). \hb@xt@.01(3.19)

### 3.3 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.2) shall be non-negative. The fourth term on the right-hand side of (3.2) will vanish due to the relation 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.2) disappear. In order to find the detailed formulation of the stress tensor, we consider the third term on the right-hand side of (3.2), 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 . For the realistic viscous fluids, in terms of Newtonian fluid theory, the Cauchy stress tensor given in (2.10) shall be included in the total stress, and consequently, we get , which ensures that the third term on the right-hand side of (3.2) is always non-negative

 −1T(σ−pI):∇u=1Tτ(u):∇u=λT|∇⋅u|2+η2T|ε(u)|2≥0. \hb@xt@.01(3.20)

So far the momentum balance equation reaches its complete form (2.11) with the stress tensor (2.10).

W define the heat flux

 q = ϕq−12M∑i=1Mw,i|u|2Ji+M∑i,j=1cij∂ni∂t∇nj. \hb@xt@.01(3.21)

The fluxes and may rely on and according to Curie’s Principle. Onsager’s principle suggests that and have a linear dependent relationship with and