Entropy stable modeling of nonisothermal multicomponent diffuseinterface twophase flows with realistic equations of state^{†}^{†}thanks: This work is supported by KAUST research fund to the Computational Transport Phenomena Laboratory at KAUST.
Abstract
In this paper, we consider mathematical modeling and numerical simulation of nonisothermal compressible multicomponent diffuseinterface twophase 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 convexconcave 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. Multicomponent twophase flow; Nonisothermal flow; Entropy stability; Convex splitting.
AMS subject classifications. 65N12; 76T10; 49S05
1 Introduction
Various nonisothermal multicomponent twophase 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 nonisothermal twophase 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 diffuseinterface models for multiphase flows have been developed in the literature. The pioneering work is that the densitygradient 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 phasefield models for immiscible and incompressible twophase 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 multicomponent twophase flows with partial miscibility and realistic equations of state (e.g. PengRobinson equation of state [54]) are intensively studied in recent years [57, 31, 32, 19, 35, 39, 55]. 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 [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 multicomponent twophase 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 diffuseinterface models for the singlecomponent 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 singlecomponent 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 liquidvapor flows is rigorously derived using the thermodynamical laws. The nonisothermal diffuseinterface 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 PengRobinson equation of state [54] extensively employed in petroleum and chemical industries due to its accuracy and consistency for numerous realistic gasliquid fluids including N, CO, and hydrocarbons. It is noted that different from the singlecomponent fluids, a reference velocity, such as massaverage velocity, molaraverage velocity and so on, usually needs to be selected for multicomponent 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 multicomponent 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 nonisothermal multicomponent diffuseinterface twophase model based on the thermodynamical laws and realistic equations of state (e.g. PengRobinson 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 twophase 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 nonnegativity 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, easytoimplement 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 entropystable numerical scheme has been developed fundamentally on the concept of functional entropy variables. Recently, in [37], using the convexconcave splitting of Helmholtz free energy density, the authors propose an entropy stable numerical method for the singlecomponent fluids with the PengRobinson equation of state, and rigorously prove that the first and second laws of thermodynamics are preserved by this method. But only singlecomponent flows are considered in these existing methods. In this paper, we focus on the numerical schemes preserving the laws of thermodynamics for the general multicomponent 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 [59, 63, 18, 6, 24]. For isothermal singlecomponent and multicomponent diffuseinterface models with PengRobinson 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 nonisothermal singlecomponent fluids is also developed in [37]. But the convex splitting approach is never studied yet for nonisothermal multicomponent fluids with PengRobinson equation of state. In this paper, we will develop the convexconcave splitting of Helmholtz free energy density based on PengRobinson 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 convexconcave splitting schemes usually result in the energydissipation feature for isothermal systems, whereas in this work we will prove that the convexconcave splitting of Helmholtz free energy density with respect to molar densities and temperature leads to the entropy stability of a numerical scheme for the nonisothermal systems. Another great challenge in designing the entropystable 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 nonisothermal multicomponent diffuseinterface twophase 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 convexconcave 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 Nonisothermal multicomponent twophase flow model
In this section, we present the system of equations modeling the nonisothermal multicomponent twophase 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
\hb@xt@.01(2.1a)  
\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
\hb@xt@.01(2.2) 
where represents the variational derivative. The entropy density, denoted by , can be defined as
\hb@xt@.01(2.3) 
We denote the internal energy density by . The internal energy, entropy, and temperature have the following relation
\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
\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
\hb@xt@.01(2.6) 
where . We get the formulation of the general pressure from (2.1), (2.5) and (2.6) as
\hb@xt@.01(2.7)  
where
We denote by the molar weight of component , and we further define the mass density of the mixture as
\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
\hb@xt@.01(2.9) 
where is a reference velocity and is the diffusion flux of component .
Let us define the stress tensor
\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
\hb@xt@.01(2.11) 
where . Thanks to the following relation (which will be proved in Subsection 3.4)
\hb@xt@.01(2.12) 
we obtain a new formulation of the momentum balance equation as
\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
\hb@xt@.01(2.14) 
where is the heat flux and
\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
\hb@xt@.01(2.16a)  
\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 semidefinite.
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.

The model has thermodynamicallyconsistent unified formulations for the general reference velocity and the mass diffusion and heat fluxes.

The model can characterize the compressibility, partial miscibility and heat transfer between different phases through the diffuse interfaces.

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 multicomponent mixture. Moreover, and have the tightly dependent relations as shown below.

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 semidefinite. 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
\hb@xt@.01(2.17) 
To determine the diffusion flux precisely, we give two typical choices of , which correspond to molaraverage velocity and massaverage velocity respectively.
 (J1)

In the first case, the diffusion mobility parameters are taken as
\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 molaraverage velocity.
 (J2)

In the second case, we take
\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 massaverage 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 timedependent volume , we define the internal energy and kinetic energy () within as
\hb@xt@.01(3.1) 
The gravitational potential energy has the form
\hb@xt@.01(3.2) 
We recall the first law of thermodynamics
\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
Cauchy’s relation between face force and the stress tensor of component gives , and as a result,
\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
\hb@xt@.01(3.5) 
where represents the general heat flux.
Applying the Reynolds transport theorem and the Gauss divergence theorem, we deduce that
\hb@xt@.01(3.6)  
where we have also used the relation . We can also derive that
\hb@xt@.01(3.7)  
On the other hand, the following overall mass balance equation can be obtained from the component mass balance equations (2.9)
\hb@xt@.01(3.8) 
Substituting (3.8) into (3.7) yields
\hb@xt@.01(3.9)  
where
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 nonnegative terms. For this purpose, we need to derive the transport equation of Helmholtz free energy density to reduce (3.11).
We define , where and
Firstly, using the component mass balance equations, we derive the transport equation of as
\hb@xt@.01(3.12)  
For the gradient contribution of Helmholtz free energy density, we can derive
\hb@xt@.01(3.13)  
and
\hb@xt@.01(3.14) 
Combining (3.12)(3.2), we deduce the transport equation of the Helmholtz free energy as
where we have also used the formulation of the general pressure . Using the identity (3.26), we deduce that
\hb@xt@.01(3.16) 
Substituting (3.2) into (3.2), we obtain
\hb@xt@.01(3.17)  
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 nonconservative terms in (3.2) shall be nonnegative. The fourth term on the righthand 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 nonconservative term on the righthand side of (3.2) disappear. In order to find the detailed formulation of the stress tensor, we consider the third term on the righthand side of (3.2), which is a nonconservative term and thus shall be nonnegative. 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 righthand side of (3.2) is always nonnegative
\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
\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