AN ENERGY DISSIPATIVE SEMI-DISCRETE FINITE-DIFFERENCE METHOD ON STAGGERED MESHES FOR THE 3D COMPRESSIBLE ISOTHERMAL NAVIER–STOKES–CAHN–HILLIARD EQUATIONS

. We consider the initial-boundary value problem for the 3D regularized compressible isothermal Navier–Stokes–Cahn–Hilliard equations descri- bing ﬂows of a two-component two-phase mixture taking into account capillary eﬀects. We construct a new numerical semi-discrete ﬁnite-diﬀerence method using staggered meshes for the main unknown functions. The method allows one to improve qualitatively the computational ﬂow dynamics by eliminating the so-called parasitic currents and keeping the component concentration inside the physically reasonable range (0 , 1). This is achieved, ﬁrst, by discretizing the non-divergent potential form of terms responsible for the capillary eﬀects and establishing the dissipativity of the discrete full energy. Second, a logarithmic (or the Flory–Huggins potential) form for the non-convex bulk free energy is used. The regularization of equations is accomplished to increase essentially the time step of the explicit discretization in time. We include 3D numerical results for two typical problems that conﬁrm the theoretical predictions.


1.
Introduction. Flows of multiphase, in particular, two-phase fluids with a significant role of capillary effects are widespread in nature and technology. Therefore, the development and improvement of mathematical models describing such flows as well as the corresponding algorithms for modeling their computational dynamics are important problems. To model two-phase flows, various methods are applied such as the level set method, the volume of fluid method, methods based on the lattice Boltzmann equations and the phase-field type methods (also called diffuse interface methods).
In the last methods, to describe the spatial distribution of phases, a special order parameter C is introduced, and it is supposed that some thermodynamic potential that determines the system state, most often, the Helmholtz free energy, depends both on C (in a non-convex way) and ∇C [2,27,29]. This makes it possible to take into consideration, in particular, the surface tension. Importantly, the interface itself is represented by not a surface but rather a layer of small non-zero thickness within which in fact the interfacial forces act. The concept of phase is not explicitly introduced at all, but due to the mentioned type of free energy, subdomains occupied by separate almost unmixed components are formed in the flow, and they are interpreted as phases.
An example of a phase-field type model is the 3D Navier-Stokes-Cahn-Hilliard (NSCH) system of quasilinear equations of the 4th order in space which describes the two-phase two-component flows of viscous compressible isothermal fluids, where C is the mass concentration of a component [2,27]. Due to the complexity of this system, the compressible case is not well represented in the literature, although compressibility (or at least quasi-compressibility [18,21]) is necessary to describe correctly a number of effects.
In this work, we use namely this NSCH system (taking into account the potential body force), with an additional so-called quasi-hydrodynamic (QHD) regularization [4,8]. This regularization, together with the more complex quasi-gas dynamic (kinetic) one associated with it, has proved itself well in the numerical solution of various problems [11,16,30]. In them, thermodynamically consistent terms with a small parameter τ > 0 are introduced into the equations. They can be considered as numerical regularizers allowing one to apply conditionally stable explicit in time and central finite-difference approximations in space which are convenient and promising in modern parallel computations. A linearized dissipativity analysis of such schemes for 1D and multi-D one-component barotropic gas equations has recently been accomplished in [37,38]. Though all the theoretical results in this paper are valid uniformly in τ ≥ 0, numerical experiments demonstrate that, in the absence of regularization (i.e., for τ = 0), the time step has to be reduced by several orders of magnitude to get adequate solutions, and thus the simulation time increases unacceptably. For certain classes of flows, the regularizing terms can also provide improved agreement with the experiment [11,16]. Other approaches to regularization are also known including recent ones in [20,31].
To ensure adequate computational flow dynamics, it is important to construct discretizations of the flow equations inheriting their basic properties including the fulfillment of the basic conservation laws, the symmetry and/or positive definiteness of some important terms and the correct approximation of equilibrium solutions, etc. They usually provide qualitatively correct computational dynamics even for relatively rough meshes. On the contrary, simple discretizations of phase-field models describing two-phase flows with capillary effects are often characterized by the presence of the so-called parasitic (or spurious) currents; they are vortex-like velocity fields in the vicinity of the interphase boundary which do not disappear when the system tends to an equilibrium state. This purely numerical effect is able to distort computational flow dynamics significantly. The parasitic currents can occur also for other models allowing for the interphase boundary resolution and capillary effects [1,12,22,28].
In the context of phase-field models, it was noted previously in [24,25] that in order to eliminate this effect, it is essential to establish the correct balance between the kinetic energy (associated with the surface tension) and free energy at the discrete level. In particular, the fulfillment of this property in combination with the symmetry and positive definiteness of the discrete Navier-Stokes viscosity tensor allows one to guarantee the total energy dissipativity for spatial finite-difference discretizations. To ensure this balance (and, therefore, the total energy dissipativity), it is proposed to rewrite the term in the momentum equation responsible for capillary effects, in a non-divergent potential form. Strictly speaking, this leads to non-conservative in momentum discretizations. However, in practice when solving a number of problems, the total energy dissipativity property turns out to be more important, and also the conservativeness in momentum holds with a rather high accuracy. For the regularized one-component compressible Navier-Stokes equations, a similar approach was implemented in [36].
In this paper, we construct a new spatial finite-difference discretization of the initial-boundary value problem for the regularized NSCH equations which has the total energy dissipativity property. This is achieved using a preliminary representation of capillary forces in the momentum equation in the potential form. Moreover, for the main sought functions, we apply staggered meshes (C-mesh according to the Arakawa classification [3]) which is implemented firstly for the regularizations of considered type. For the same equations, alternative discretizations on non-staggered meshes for all the main sought functions having the total energy dissipativity property have quite recently been constructed in the simplified 1D [5] and 3D periodic [7] statements. It turned out that the discretizations of the Navier-Stokes viscosity tensor and regularizing terms are simpler in the case of the staggered main meshes.
The second element in the paper allowing one to achieve the qualitative improvement in computational flow dynamics is the involvement of a logarithmic, or the Flory-Huggins [23], form for the non-convex term in the Helmholtz free energy which forbids the C values outside its physical range (0, 1). Despite the fact that this form is physically well-motivated, the simpler 4th order polynomial form (actually the simplest approximation of the logarithmic one) is common in the literature but often leads to appearing of the C values outside of (0, 1) and thereby, strictly speaking, to non-physical solutions [13,17,32,33].
In addition, advantages of the constructed method are the fulfillment of the discrete total mass and component mass conservation laws (their importance has recently been discussed in [18,21]) and the property that the method is well-balanced meaning it is able to reproduce the equilibrium solutions adequately.
The constructed method is discretized explicitly in time and practically implemented and tested on two typical problems such as the equilibrium state of a spherical droplet (for various radii and the model parameters) and the 3D spinodal decomposition which can be considered as a spontaneous unmixing of the initially homogeneous mixture perturbed randomly. For the first problem, a good agreement is obtained between the observed pressure difference inside and outside the droplet and similar one according to the well-known Young-Laplace law. Also the kinetic energy of the system tends to zero and thus the possible parasitic currents disappear as time grows (in contrast to the method based on the divergent "stress form" of capillary forces in the momentum equation). In all computations, the nonincrease of the total energy in time is observed. The presented numerical simulation results illustrate and confirm the theoretical conclusions.
Notice that ultimately our research focuses largely on computations in complex 3D spatial domains having the voxel representation such as the pore space of rock samples [6,9,15]. This forms a part of digital rock physics technologies and is of considerable interest for modern industrial applications [9,15].

VLADISLAV BALASHOV AND ALEXANDER ZLOTNIK
The paper is organized as follows. In Section 2, the regularized Navier-Stokes-Cahn-Hilliard equations are written out, and the conservation laws and the total energy balance equation as well as equations for equilibrium solutions are given for them. The logarithmic (or Flory-Huggins) potential is discussed. In Section 3 the mesh notation is introduced, a semi-discrete finite-difference method is constructed, and the discrete total mass and component mass conservation laws as well as the equations for the mesh equilibrium solutions are derived. In Section 4, we first state and prove the symmetry and positive definiteness of the constructed discrete Navier-Stokes viscosity tensor and then the main result on the discrete total energy dissipativity. Section 5 is devoted to numerical experiments with the fully discrete method using the explicit Euler approximation in time. Numerical studies of two typical 3D problems concerning the equilibrium state for a single droplet and spinodal decomposition are accomplished including demonstration of non-increasing of the discrete total energy and the absence of parasitic currents.

2.
Regularized Navier-Stokes-Cahn-Hilliard equations. We rely on the Navier-Stokes-Cahn-Hilliard (NSCH) system of equations describing isothermal viscous compressible flow of two-phase two-component mixture with surface effects [2,27]. Its regularized version [4,8] consists of the total mass, momentum and component mass balance equations where ρ > 0 is the total density, u = (u 1 , u 2 , u 3 ) is the velocity of both mixture components, 0 < C < 1 is the mass concentration of one of two considered components all depending on (x, t), with x = (x 1 , x 2 , x 3 ) ∈ Ω := (0, X 1 ) × (0, X 2 ) × (0, X 3 ) and t ≥ 0. Hereafter ∂ t = ∂/∂t and ∂ i = ∂/∂x i , the divergence of a tensor is taken with respect to its first index, symbols ⊗ and · correspond to the tensor and inner products of vectors respectively. In addition, Φ(x) is the potential of a stationary body force and M = M (C) > 0 is the mobility coefficient. The regularized mass flux J = (J 1 , J 2 , J 3 ) is given by the formulas It includes the regularizing momentum m (rather than a regularizing velocity following [7,35] but in contrast to the standard approach) with the regularization parameter τ = τ (ρ, C) > 0. The tensor Π = Π N S − Q + Π τ consists of the Navier-Stokes viscous stress tensor Π N S , the capillary stress tensor Q and the regularizing tensor Π τ : where ∇u = {∂ i u j } 3 i,j=1 , η = η(ρ, C) > 0 and ζ = ζ(ρ, C) ≥ 0 are the dynamic and bulk viscosity coefficients as well as I is the unit tensor.
The regularizing terms in the equations are physically motivated [4,11,16,30]. In numerical computations, they provide (conditional) stability of explicit central finite-difference discretizations of the equations.
It remains to define the pressure p and the (generalized) chemical potential µ. We set the Helmholtz free energy of the mixture as follows where the constant λ 1 > 0 is a model parameter, c si andρ i > 0 are the sound speed and reference density of ith component, i = 1, 2.
In the present work we are mainly interested in the logarithmic form of Ψ sep (C), which is also known as the Flory-Huggins potential [23] Ψ (The theoretical constructions are valid for general Ψ sep (C) but numerical simulations are accomplished namely for its form (7).) The parameters ω 1 > 0 and ω 2 > 0 depend on the mixture temperature and the component properties; hereafter we assume that they are constant. If ω 1 ≥ ω 2 , then Ψ sep (C) ≥ 0, and thus Ψ sep (C) is globally convex on (0, 1).
If ω 2 > ω 1 , then Ψ sep (C) is non-convex and has two minima that is of great concern in simulation of surface effects. In this case, it is possible to construct the line tangent to the graph ofΨ 0 (C) at two points C − and C + , see Fig. 1; herẽ Ψ 0 (C) := Ψ 0 (ρ, C), whereρ > 0 is fixed. The values C − and C + are the mixture equilibrium concentrations [10,26]. The values C + and C − are solutions of the following non-linear system [26] Hereafter The more popular choice for the non-convex part of free energy is Ψ p4 = A ψ C 2 (1− C) 2 which can be regarded as a polynomial approximation of Ψ ∆ (C) := Ψ sep (C) − Ψ sep (C + ). The function Ψ p4 has minima at C = 0, 1 but has no singularities as C → +0 and C → 1 − 0. This easily allows C to leave the range (0, 1) [13,17,32,33].

VLADISLAV BALASHOV AND ALEXANDER ZLOTNIK
Other polynomial approximations can be obtained by applying the Taylor expansion of Ψ sep (C) [14,19].
The functions p and µ are introduced using Ψ 0 (ρ, C): It should be noted that the nonlinear terms with ∇C appearing in Q in combination with the non-convexity in C of Ψ 0 are responsible for the surface effects. If τ = 0, then m = 0, and system (1)-(3) is reduced to the standard NSCH system for viscous compressible two-component isothermal mixture flows.
We supplement system (1)-(3) with the initial conditions where ρ 0 , u 0 and C 0 are the given initial functions. Let n be the outward normal to ∂Ω. Under the boundary conditions J · n| ∂Ω = 0 and ∂ n µ| ∂Ω = 0, equations (1) and (3) imply the following total mass and component mass conservation laws For the considered system, the following important result is also valid [8].
Consequently, the total energy is non-increasing: ∂ t E ≤ 0 for t ≥ 0.
Finally, recall the equilibrium solutions ρ = ρ S (x) > 0, u = 0 and 0 < C = C S (x) < 1 of the considered system. Similarly to [7] it can be checked that under the boundary condition ∂ n µ| ∂Ω = 0 they satisfy the following system of equations Notice that the above regularization does not influence these equations.
In the very particular case of the plane interface with ρ = const, C = C(x 1 ), x 1 ∈ R and Φ = 0, according to [10,26] the parameters λ 1 , ω 1 and ω 2 are connected with the surface tension coefficient by the following equalities where The interface thickness can be defined using the slope of the line tangent to the graph of C(x) at x 0 = C −1 (1/2), see also Fig. 4a below and [10]: Assuming ω 2 ∝ ω 1 we conclude that σ ∝ √ ω 1 λ 1 and l ∝ λ 1 /ω 1 .

3.
A semi-discrete method on staggered meshes.
be an auxiliary (the adjoint) mesh and ω * hk = {x k(l−1/2) } N k l=1 be its internal part. Denote by H(ω) the space of functions defined on a mesh ω.
We define the averages and difference quotients for v ∈ H(ω hk ) and y ∈ H(ω * hk ) We introduce the inner products respectively in H(ω hk ), H(ω * hk ) and H(ω hk ) Similarly to [36], we will apply the following formula for v ∈ H(ω hk ) and y ∈ H(ω * hk ). In particular, if (s * k y) 0 = (s * k y) N k = 0, then the last two terms on the right disappear; moreover, one has We need to define the collections of 3D meshes and corresponding inner products and (y, z) * ≡ (y, z) 1 * ,2 * ,3 * = (((yz, 1) 1 * , 1) 2 * , 1) 3 * . Essentially, for a permutation (k, l, m) of (1, 2, 3), below we use the mesh ω hk,l,m and the inner product (v, w)k ,l,m meaning that the mesh ω hĩ is taken in ith direction, for i = k, l, m. We also need some similar 2D meshes ω hl,m and inner products (v, w)l ,m .

VLADISLAV BALASHOV AND ALEXANDER ZLOTNIK
Location of nodes of ω hk,l * ,m * (thick dots) and ω hk * ,l,m (red crosses), where u k and Π lk , l = k, are respectively defined.

3.2.
A semi-discrete finite-difference method and some its properties. We discretize in space the system of equations (1)-(3) taking into account (13) on the above staggered meshes and construct the following discrete total mass, momentum and component mass balance equations where t ≥ 0 and k = 1, 2, 3. Hereafter div h J = δ i J i and (k, l, m) is a permutation of (1, 2, 3). Also the summation over repeated indices i and j (and only over them) is assumed, and δ (kj) is the Kronecker symbol. Notice that in the sums like (1 − δ (kj) )a kj , the terms a kk may not be defined at all. The main sought functions ρ > 0 and 0 < C < 1 are defined onω * h whereas u k is defined on ω hk,l * ,m * , k = 1, 2, 3, see Fig. 2a-2c; also Φ is a function given onω * h . The quantities J, Π, G h and µ, see (4)-(6), (9) and (13), are discretized as follows The construction of equations (20)- (22) and formulas (23)-(29) is non-trivial and, in particular, the choice of the averaged quantities becomes more clear below in the course of proving the energy dissipativity for the method. To this end, we have exploited some ideas from [36] where energy dissipative spatial discretizations on the non-staggered main mesh (i.e., for the main sought functions) were considered.
We supplement the above equations by the initial and boundary conditions for all k, l ∈ {1, 2, 3}, k = l, where once again (k, l, m) is a permutation of (1, 2, 3). For such k, l, the first boundary condition (31) implies The boundary conditions (32) and (33) mean that u k is odd in x l with respect to x l = 0, X l as well as C and µ are even in x k with respect to x k = 0, X k . In addition, we will assume that ρ has the same parity property as C and µ. Moreover, actually it is convenient to define C and ρ on the enlarged meshω l=−1 , k = 1, 2, 3, conserving their same parity properties. This ensures that formulas (23)- (29) are correct on the mentioned meshes.
The derived equations are quite similar to the equilibrium differential equations (14)- (15). This means that the constructed discretization is able to reproduce the equilibrium solutions correctly, i.e. it is well-balanced. Moreover, the result holds uniformly in variable τ k ≥ 0, k = 1, 2, 3, in contrast to the case of a non-staggered main mesh [7]. 4. The energy dissipativity of the semi-discrete method. We begin with the properties of the Navier-Stokes terms (26)- (27).
where u k and v k are defined on ω hk,l * ,m * , k = 1, 2, 3. The following self-adjointness and positive definiteness properties of the Navier-Stokes terms hold in the latter property, the condition δ * i u k | x k =0,X k = 0 for all i, k ∈ {1, 2, 3}, i = k is assumed to be valid.
Proof. Rearranging the terms and changing the order of summation, we easily obtain thus (38) is valid. Here (i, j, k) is a permutation of (1, 2, 3). Next taking v = u, extracting the summands for x j = 0, X j in the third term on the right in (40) and using the assumed boundary condition in it as well as symmetrizing the rest, we obtain that clearly implies (39).
Note that this lemma is more general and its proof is much simpler than the corresponding result in the case of the non-staggered main mesh [36]. Now we state and prove the main theoretical result of the paper.
the following dissipativity property holds The 2nd-4th terms on the left are clearly non-negative, in particular, see (39). Consequently, the discrete total energy is non-increasing: ∂ t E h ≤ 0 for t ≥ 0.
Proof. The proof comprises three steps. 1. We first take the inner product in H(ω * h ) of the mass balance equation (20) Owing to formulas (19) and the boundary conditions (31) we have Owing to the definition G h = Ψ 1ρ (ρ, C) + 1 2 λ 1 |∇ h C| 2 and then the first formula (19) together with the second boundary condition (31) as well as the mass balance equation (20) we can write Using the formulas δ * k s k = s * k δ k , 1 2 δ k (v 2 ) = (s k v)δ k v and 1 2 δ * k (y 2 ) = (s * k y)δ * k y and then the second formula (19) together with the second boundary condition (31) and property (34) we also get Thus we can rewrite (42) as follows 2. Now we take the inner product in H(ω hk,l * ,m * ) of the momentum balance equation (21) and u k , replace k with i and sum up the results in i = 1, 2, 3. Using the formula (19) two times together with the first boundary condition (31) and (32) as well as formula (18) together with (32) we obtain Also owing to the second formula (19) together with the first boundary condition (31) we have as well as owing to the first property (34) and the same formula and boundary condition we get Consequently adding equalities (42) and (44) as well as taking into account the two last formulas and the definitions J k = (s * k ρ)u k − m k and (38), after reduction of terms, we derive 3. Finally we turn to the component mass balance equation (22). Applying the formula δ k (J k s * k C) = (δ k J k )C + s k (J k δ * k C) and the mass balance equation (20) we transform its left-hand side to the non-divergent form We take the inner product in H(ω * h ) of the equation and µ, see (29), apply both formulas (19) together with the boundary conditions (33) and get Applying sequentially formulas (19) together with the first boundary condition (33) we can write down the chain of transformations Summing up equality (45) and the last one, using the formula Ψ 1ρ ∂ t ρ + Ψ 1C ∂ t C = ∂ t Ψ 1 , reducing some terms and recalling definition (23)-(24) of m i , we complete the derivation of (41).

5.
Numerical study of the fully discrete scheme. We apply the explicit Euler method for the approximation in time of the above introduced spatial discretization. The section is devoted to numerical study of this fully discrete scheme. In all the tests below, we study the two-component mixture with the identical state equations and, therefore, with c s1 = c s2 = c s . The following constants are used η = 5 · 10 −4 Pa · s, ζ = 0, M (C) ≡ M = 10 −8 kg · s/m 3 and c s = 200 m/s. We take the initial density ρ 0 = 2.5 kg/m 3 and velocity u 0 = 0.
The logarithmic form (7) of Ψ sep (C) is used, and we set ω 2 = 2ω 1 that after solving system (8) gives C − ≈ 0.02125 and C + ≈ 0.97875. The parameters λ 1 and ω 1 and the initial concentration C 0 are taken depending on a specific test problem.  We take the regularization parameter τ = 0.5h/c s associating it with the spatial mesh. But we do not apply a regularization of η contrary to usual QHD-approach.
For our computations, the high-performance computer "Zhores" installed at Skolkovo Institute of Science and Technology [34] was used. 5.1. Stationary single droplet. The first test problem deals with numerical study of equilibrium state for a single droplet with a radius R. The initial concentration is set as where x c = 0.5(X, X, X) is the center of Ω. Here β = 10 4 √ 5 m −1 , and ε c = 0.01 defines the "offset" of C 0 (x) values from 0 and 1. For R = 0.25X, C 0 (x) in the domain x 1,2 ∈ [0.5X, 0.87X], x 3 = 0.5X is represented in Fig. 3a.
We consider λ 1 = λ * 1 /k * and ω 1 = ω * 1 /k * , where λ * 1 = 2 · 10 −6 J · m 2 /kg and ω * 1 = 2 · 10 3 J/kg, with k * = 1, 2, 10 and call these cases respectively (I), (II) and (III). According to formulas (16)-(17) they correspond to the same interface thickness l ≈ 3.75h but different surface tension coefficients: σ I = 0.08831, σ II = 0.04416 and σ III = 0.008831 in J/m 2 .  Note that the form of C 0 corresponds to the interface thickness ≈ 26h (cf. Fig. 3a). In Fig. 3a-3d, clearly l is reduced drastically in time. In Fig. 4a-4b, distributions of C and ρ along the segment with x 1 ∈ [0.5X, X] and x 2 = x 3 = 0.5X are presented at t = 20 · 10 3 ∆t. From Fig. 4a apparently one has the observable thickness ≈ 4h that is close to l ≈ 3.75h predicted by 1D theory. Also here we emphasize that since the logarithmic potential is used, C remains within the physically reasonable range (0, 1). However, its maximum value C max exceeds C + and its minimum value C min stays slightly greater then C − , see subfigures in Fig. 4a. This phenomenon is similar to the so-called "bulk shift" effect [17,33] meaning that the values of C are "shifted" with respect to its equilibrium values. But thus if the latter values are 0 and 1 as for the popular polynomial potential, then C leaves (0, 1) that makes the computational dynamics false (non-physical) even if it does not destroy computations.
In Fig. 4b, ρ falls abruptly in the vicinity of interface. We can observe this effect since the compressible isothermal case is considered (see also [5]).
Also we observe that ρ in > ρ out , where ρ in and ρ out are the average densities inside and outside the droplet. Due to the state equation p = c 2 s ρ the pressure difference occurs. It is well known that it is related to the droplet radius by the Young-Laplace law: and the state equation implies ∆p = c 2 s (ρ in − ρ out ). Here σ L is the Young-Laplace surface tension coefficient. It can be different from (16) since that is obtained based on the planar approximation. Looking ahead we notice that nevertheless these coefficients are close to each other.
In (47), R a stands for the observable value of droplet radius calculated using the droplet volume 4 3 R 3 a = h 3 N C>0.5 , where N C>0.5 denotes the number of mesh cells with C > 0.5. We define ρ in and ρ out as follows where the sums are taken over the mesh cells such that |∇ h C| <ε c /h and respectively C > 0.5 or C < 0.5 (to identify the cells inside and outside the droplet) as well as N in and N out are the respective numbers of such cells; hereε c = 10 −2 is chosen. The condition on |∇ h C| allows to distinguish between the cells located in the bulk phases and within the interface (where |∇ h C| is large); see Fig. 4a in this respect. Note that these formulas are more robust than those used by us earlier [7].
In Fig. 5 (left), evolution of σ L (t) is presented in the cases (I)-(III) obtained in simulation of the droplet equilibrium state with R = 0.18X. In each case, σ L finally tends to the values predicted by the 1D theory. However, the time moments at which σ L almost attains these values vary significantly. According to Fig. 5 (right) presenting σ L (t) in the case (III) for other radii, this time moment (t ≈ 1.5·10 −3 s) is almost independent of R but it depends heavily on the initial interface thickness. To verify the latter claim, we have performed an additional simulation with R = 0.375X and the parameters like in the case (III) but with β = 10 5 m −1 , which implies the interface thickness ≈ 6h. The result is also presented in Fig. 5, and in this situation the surface tension attains σ III considerably faster. This peculiarity is important for simulation practice: the initial interface thickness should correspond to the parameter values to provide the prescribed surface tension as quickly as possible. In its turn, the correct surface tension coefficient is necessary to ensure the correct computational flow dynamics. Similar simulations have been performed for other radii R = 0.18X, 0.25X, 0.3X and 0.375X in (46) as well. In Fig. 6, the dependence of observable pressure difference ∆p on 1/R a is given. The average values of σ L presented there have been   Figure 9. Evolution of C min and C max (for spinodal decomposition).
calculated according to the line slopes. They are close to σ I , σ II and σ III . We note that R a slightly differs from R, i.e. |R − R a | h, in some cases. In Fig. 7a and 7b, evolution of the mean kinetic energyĒ kin := ρ, 1 2 s i (u 2 i ) * /|Ω| and the total energy E h is demonstrated for the droplet with R = 0.25X in the cases (I)-(III); here |Ω| = X 3 is the volume of Ω. ClearlyĒ kin oscillates and goes to zero whereas E h is non-increasing as time grows. In Fig. 7b,Ẽ is used for convenience of presentation. The same behavior has been obtained for other R. Note that C 0 is the same in these cases; however, the initial total energy E h | t=0 is different (see Fig. 7b) since the model parameters are different too.
In Fig. 8, we also compare the behaviorĒ kin (t) for the scheme constructed in this paper (scheme A) and one based on the standard divergent form of capillary forces from [6] (scheme B). The stationary droplet for R = 0.25X and the parameters in the case (III) are considered. Until t ≈ 6 · 10 −6 the simulation results are almost identical. But, for scheme B,Ē kin (t) is increasing rather than decreasing at

Figure 11
large times that is related to the presence of the so-called "parasitic (or spurious) currents". This known unpleasant numerical phenomenon manifests itself in the form of vortex-like velocity fields that appear in the vicinity of interface and do not diminish when the system goes to its equilibrium state. On the contrary, scheme A demonstrates rapid decreasing ofĒ kin (t) thus eliminating the parasitic currents. For the regularization parameter, the formula τ = α * h/c s with α * = 0.5 has been used above as in [5][6][7][8]. Let us analyze the choice of other values of α * for scheme A. To reduce the total amount of computations, we accomplish this for the 2D version of the problem posed at the beginning of this section in the case of R = 0.25X and the set of parameters (I).
First, let α * = 0, i.e., there is no regularization. The attempt to use the same time step ∆t = 10 −8 s as well as its fivefold smaller value ∆t = 2 · 10 −9 s, leads us to numerical instability and rapid collapse of the solution (in the latter case it happens slightly later). This can be seen in Fig. 10 and 11a where the graphs ofĒ kin (t) (here in the 2D case |Ω| = X 2 ) and E h (t) −Ẽ withẼ = 9.16711642 · 10 −2 J are given. The unstable computations are characterized by a sharp increase in both E h (t) and E kin (t). With a further decrease of ∆t down to 10 −9 s, the computations become already stable and allows us to get the equilibrium state (up to machine precision), see again Fig. 10 and 11a. Thus, here for the explicit Euler discretization in t, under the regularization one can use 10 times larger step ∆t than in its absence. In stable calculations, the solutions obtained for α * = 0 and 0.5 are close enough, see Fig. 10 and 11a. Their difference decreases further as the system tends to its equilibrium state as it follows from Fig. 11b and comparison of solutions at time t = 10 −3 s: the difference between them in the uniform grid norm is 9.4 · 10 −9 for ρ (the relative difference is meant for this function) and 1.3 · 10 −7 for C. This behavior of solutions at large t is consistent with the fact that the contribution of the regularizing momentum m disappears as the system tends to its equilibrium state.
Let us also discuss the choice of some α * ≥ 0.5. Fig. 12a shows the corresponding graphs ofĒ kin (t). For α * = 1 and ∆t = 10 −8 s the computations are stable. For α * = 2 and ∆t = 10 −8 s or 7 · 10 −9 s the computations are unstable and collapse rather quickly (the former case is not shown in the figure), whereas for ∆t = 5·10 −9 s the stability takes place. Thus, a noticeable increase in α * affects negatively the value of ∆t that provides stability of computations; this is also known both from the similar practical experience [16] and the theoretical analysis of linearized stability [37,38].
Moreover, as one can see from Fig. 12b (see also Fig. 10), the increase in α * leads to a decrease in the oscillation amplitude ofĒ kin (t). However, further computations for α * = 2 with twice and four times less h and correspondingly less ∆t lead to an increase in the amplitude. This means that the presence of oscillations ofĒ kin (t) is physically natural during the formation of the equilibrium phase interface in the problem and also confirms the existing experience [5][6][7][8] that α * = 0.5 provides a more adequate numerical solution (cf. also the above case α * = 0).
It was also verified that in all the stable computations performed, E h (t) is nonincreasing, and its graphs are very close to those shown in Fig. 11a.

Spinodal decomposition.
The last test deals with numerical study of spinodal decomposition which can be regarded as a spontaneous "unmixing" of the initially perturbed homogeneous mixture. Now we take C 0 (x) as the uniformly distributed random variable in [0.5 − 0.05, 0.5 + 0.05]. The parameters λ 1 , ω 1 and ω 2 are set as in the case (III). In Fig. 9, evolution of C max and C min is presented. At the very beginning of the process, C tends to 0.5 but afterwards C max increases and C min decreases rapidly. In zoomed parts of the plots, they both leave the interval (C − , C + ) but still stay within the physical range (0, 1). In Fig. 13a, evolution ofĒ kin is shown. After the first time stepĒ kin increases considerably which is presumably since C 0 is a discontinuous, oscillating from node to node, function. ThenĒ kin decreases: the system is relaxed (cf. the corresponding time interval in Fig. 9). AfterwardsĒ kin increases; moreover, this begins almost simultaneously with increasing of C max (or decreasing of C min ). After the mixture decomposes,Ē kin is "decreasing" non-monotonically.
In Fig. 13b, the presented evolution of total energy E h is non-increasing. Interestingly, E h is almost constant as long as the mixture stays almost homogeneous, but as soon as the system becomes considerably heterogeneous (C max ≈ 0.8 and C min ≈ 0.2, see Fig. 9), E h begins to decrease rapidly.
In Fig. 14a-14d, evolution of isosurface C = 0.5 is given. Initially the system is quite fine-grained. As it becomes unmixed, the structures turn coarser and coarser. Notice the perpendicular position of the interface with respect to ∂Ω that correctly corresponds to the used first boundary condition (33). 6. Conclusion. In this paper, a new spatial finite-difference discretization of the 3D compressible isothermal Navier-Stokes-Cahn-Hilliard system of equations has been constructed. The system describes flows of a two-component two-phase mixture taking into account capillary effects. The equations were pre-regularized according to the quasi-hydrodynamic approach, and the initial-boundary value problem for them has been considered. Unlike previous similar papers, staggered meshes have been used in the discretization for the main unknown functions (the total density, velocity of both mixture components and the mass concentration of a component). This has allowed us to simplify significantly the discretization of regularizing terms which in turn has simplified the implementation of the method. The main property of the constructed discretization is its total energy dissipativity, and the corresponding theorem has been proved. The fulfillment of this property is closely related to solving the well-known problem of eliminating the so-called spurious currents in numerical solutions. In addition, the discretization is well-balanced, i.e., (a) t = 2 · 10 5 ∆t = 2 · 10 −3 s (b) t = 3 · 10 5 ∆t = 3 · 10 −3 s (c) t = 6 · 10 5 ∆t = 6 · 10 −3 s (d) t = 8 · 10 5 ∆t = 8 · 10 −3 s Figure 14. Isosurfaces C = 0.5 at different time moments (for spinodal decomposition).
it is able to reproduce the equilibrium solutions adequately. Importantly, both of these properties remain valid in the absence of regularization. A numerical study of the spatial discretization combined with the explicit Euler method in time has been carried out as well. This is the simplest way to implement the discretization and, at the same time, sufficient to evaluate its quality. To this end, the problem of the 3D equilibrium state of a spherical droplet and its compliance with the classical Young-Laplace law as well as the problem of 3D spinodal decomposition have been considered. In all cases, non-increasing of the total energy of the system and tending of its kinetic energy to 0 (thereby, the absence of spurious currents) at large times have been observed. In addition, the logarithmic form of Helmholtz free energy has been applied in the simulations preventing the concentration values from going beyond the physically meaningful interval (0, 1) throughout the entire computation, which is often the case when using some other its forms.
The question of choosing the regularization parameter τ > 0 has been raised. In computations associated with the former problem, for an adequate choice of τ > 0, the admissible time step has turned out to be 10 times more compared to the case τ = 0, i.e., in the absence of regularization. At the same time, overestimating the value of τ has lead to negative effects such as reducing the time step to ensure stability of computations and excessive smoothing of numerical solutions.