Implicit-Explicit Schemes for Compressible Cahn–Hilliard–Navier–Stokes Equations

According to Kynch’s theory (see [20]) for sedimentation of homogeneous monodisperse suspensions, consisting of solid spherical particles of the same diameter and density immersed in a viscous fluid, two interfaces form in the settling process: a descending interface between the clear liquid and the initial homogeneous mixture and an ascending interface between the maximally concentrated mixture and the initial homogeneous mixture.

In [28], it is observed that, after several days, a colloidal monodisperse suspension of polystyrene particles forms a layered structure when sedimenting, a fact that contradicts Kynch’s theory. A spinodal decomposition, governed by the Cahn–Hilliard equation [9], is then conjectured as the underlying mechanism that explains this phenomenon.

The Cahn–Hilliard equation models the evolution of a mixture of immiscible fluids by letting an interface between them be diffuse, with some small, but positive, thickness. This is achieved by considering an expression for the Ginzburg–Landau free energy with terms that depend on the gradient of the difference of mass-fractions of the components of the mixture, in such a way that the Cahn–Hilliard equation is a gradient descent equation for the free energy, in contrast to other alternative models, as Volume Of Fluid methods [17] and Level Set methods [23, 24], that do not take energy into account.

The Cahn–Hilliard equation cannot explain, by itself, the cited layering phenomenon, since it does not take into account the gravitational force. The incorporation of dynamics into the evolution of binary fluids can be achieved by using Lagrangian techniques [1, 6, 22]. This leads to systems of equations that couple the evolution of individual mass-fractions, bulk momenta and, possibly, energy.

Among these possibilities, and although (quasi) incompressible versions of these equations might be more suitable for explaining the cited layering phenomenon, in this paper we focus on the approach in [1]. Ignoring temperature, this yields a system of equations, the isentropic compressible Cahn–Hilliard–Navier–Stokes equations, which are a system of fourth-order partial differential equations.

The difficulties in the numerical treatment of Cahn–Hilliard–Navier–Stokes equations stem from the fact that their second, third and fourth order differential terms in spatial coordinates have discretizations with large eigenvalues, due to their inherent differential nature, e.g., the finite difference discretization of a fourth-order derivative with respect to x will have \(\varDelta x^4\) in the denominator, where \(\varDelta x\) is the grid size in the x coordinate.

When solving a system of ODE \(z'=f(z)\), with Jacobian matrices \(f'(z)\) with negative eigenvalues with large absolute values (stiffness), implicit methods can perform larger time steps, whereas an explicit method is typically constrained to time steps proportional to \(\frac{1}{|\lambda |}\), where \(|\lambda |\) is the largest absolute value of eigenvalues of \(f'(z)\). This situation does not hold if \(f'(z)\) has positive eigenvalues, i.e., an implicit solver is forced also to perform small time steps.

In some situations where \(f(z)=g'_e(z)-g'_c(z)\), with an expansive term \(g'_e(z)\) and a contractive term \(g'_c(z)\), where \(g_c, g_e\) are strictly convex functions, an Implicit-Explicit scheme that treats \(g'_c(z)\) implicitly and \(g'_e(z)\) explicitly gives a gradient-stable method, as shown in [13].

Since the Cahn–Hilliard equation has diffusive (contractive) and anti-diffusive (expansive) terms and the viscous terms in the Cauchy stress tensor of the Navier–Stokes equations are also diffusive, efficient time-stepping schemes for the Cahn–Hilliard–Navier–Stokes equations should treat implicitly high order terms that are contractive and explicitly the rest, specially if they do not contribute to stiffness.

Effective numerical schemes for the Cahn–Hilliard equation can be traced back to [8, 11,12,13]. Numerical methods for incompressible and quasi-incompressible models have the additional difficulty of having to solve an equation for the pressure related to the divergence of the velocity (this being 0 in the case of incompressible models) (see [10, 15, 18, 19, 21, 26, 31]). Numerical methods for compressible models do not have this burden but need anyway to cope with the issues inherent to convective terms. In [16] relaxation techniques are used for this task.

The contribution of this paper is the proposal of second order Implicit-Explicit Runge–Kutta schemes for the isentropic compressible Cahn–Hilliard–Navier–Stokes equations with stability restriction on the time step dictated only by the convective terms, so that it is proportional to the spatial grid size. Some features of the scheme are that only linear systems have to be solved at each stage of the Implicit-Explicit Runge–Kutta method and those are efficiently solved by the V-cycle multigrid algorithm with Gauss-Seidel smoothing.

The paper is organized as follows: in section 2 the compressible isentropic Cahn–Hilliard–Navier–Stokes equations are introduced; in Sect. 3 implicit-explicit Runge–Kutta finite difference numerical schemes for the two-dimensional equations are proposed; in Sect. 4 we perform some numerical experiments to assess the efficiency of our proposals; finally, in Sect. 5 we draw some conclusions and give some perspectives for future research.

2.1 Model Description

The models in this exposition are based on [1]. We denote by \(c_i\), the mass concentration of species \(i=1,2\), by \(c=c_1-c_2\), by \(\rho \) the density of the mixture and by \(\varvec{v} \) its bulk velocity (we use boldface for vector variables).

We denote by \(\Omega \) the open set in \(\mathbb {R}^3\) that is filled by the fluids and by \(\varepsilon \) a parameter related to the thickness of the diffuse interface of the fluid mixture. The Ginzburg–Landau free energy in some region \(V\subseteq \Omega \) of the immiscible compressible two-phase fluid is

where \(\psi (c) = \frac{1}{4}(c^2-1)^2\) is a double-well potential function and \(f_{e}\) is the specific Helmholtz free energy of an equivalent one-phase fluid.

The isentropic compressible Cahn–Hilliard–Navier–Stokes equations with gravitation are the following equations:

where \(\mathop {\text {div}}\) is the divergence operator with respect to \(\varvec{x}\in {\mathbb {R}}^3\), the first equation is the continuity equation for the mixture, the second is the equation for conservation of bulk momenta. In these equations

is the stress tensor, \(p(\rho , c)=\rho ^2\frac{\partial f(\rho , c)}{\partial \rho }\) is the fluid pressure, \(\nu (c), \lambda (c)>0\) are the viscosity coefficients, \(\varvec{g}\) is the gravitational acceleration, and

is the chemical potential.

The \((\rho , \rho \varvec{v})\)-subsystem, with \(\varepsilon =0\), form the compressible isentropic Navier–Stokes equations. The equation for \(\rho c\), for constant \(\rho \) (which may be assumed to be 1) and \(\varvec{v}=0\), is the Cahn–Hilliard equation (see [9]):

These equations are supplemented by initial conditions \(\rho _0, \varvec{v}_0, c_0\) and the boundary conditions

where \(\varvec{n}\) is the outward normal vector to the boundary.

In [1] it is proved that these equations admit weak solutions, with renormalization of \(\rho \) in the sense of Di Perna and Lions, in any interval \([0, T], T>0\), provided \(\gamma > \frac{3}{2}\), \(0\le \rho _0\in L^{\gamma }(\Omega ), \rho _0|v_0|^2 \in L^1(\Omega ), c_0\in H^1(\Omega )\).

We henceforth consider \(\nu (c), \lambda (c)\) constant and \(p=p(\rho )=C_{p}\rho ^{\gamma }\), for a positive constant \(C_{p}\) and the adiabatic constant \(\gamma > 1.5\), which corresponds to \(f_e(\rho )={C_p}\frac{\rho ^{\gamma -1}}{\gamma -1}\). Therefore, the equation for the conservation of bulk momenta can be rewritten as:

As expected, the \(\rho \) and \(q=\rho c\) variables are conserved since the respective associated fluxes,

vanish at the boundary due to (4).

The two-dimensional version of these equations, for \(\varvec{v}=(v_1, v_2)\), is:

where \(\varDelta w=w_{xx}+w_{yy}\), for any \(\mathcal {C}^2\) function \(w=w(x,y)\), and the gravity acceleration \(\varvec{g}=(0,g)\) acts along the y coordinate.

We also consider the one-dimensional version of these equations:

where gravity g acts along the x coordinate.

2.2 Spinodal Decomposition

The Cahn–Hilliard equation was proposed in [9] to model the separation of a homogeneous mixture of two incompressible fluids, the first of them stable with respect to the presence of small quantities of the second one, and this one unstable with respect to the presence of small quantities of the first one. The boundary of the unstable region in (c, T, p)-space, T, p being temperature and pressure, respectively, is given by the equation \(\frac{\partial ^2 G}{\partial c^2}=0\), where G(c, T, p) is the Gibbs free energy density of the fluid, and is usually named the spinodal.

To analyze the spinodal decomposition (see [11] for further details), we consider the linearization of (3) about a constant state \(c_0\) in the spinodal region \((-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}})\), which satisfies \(\psi ''(c_0)<0\):

for assumedly small u, with \(\int _{\Omega } u=0\). Notice then that

therefore the linearized Cahn–Hilliard equation:

is deduced from the Cahn–Hilliard equation (3).

By separation of variables, functions of the form

satisfy the homogeneous Neumann boundary conditions and \(\int _{\Omega } u(\varvec{x}, t)d\varvec{x} =0\). Since

u is solution of (7) if

whose solutions are

The linearized equation (7) will therefore develop instabilities when \(t\rightarrow \infty \) if the coefficient of t in the exponent of (9) is \(>0\), i.e., provided

for some \( k_i \ge 1\), and this will be so if \(\varepsilon \pi ^2+\psi ''(c_0) < 0\). But these instabilities, triggered by the linearized equation, will grow until some point, when the nonlinear character of the Cahn–Hilliard equations makes the linearization no longer valid, and will stop, since energy decreases with time as shown below, by using differentiation under the integral, integration by parts and Neumann boundary conditions for c:

for \(F(c)=E(1, c)\), where E is given in (1).

3.1 Spatial Semidiscretization

In this section, our purpose is to design finite difference numerical methods for the efficient approximation of solutions of the two-dimensional equations in Section 2. For this, we consider \(\Omega =(0,1)^2\) and the equispaced computational grid given by the \(M^2\) nodes \(\varvec{x}_{i,j}=((i-\frac{1}{2})h, (j-\frac{1}{2})h), i,j=1,\dots ,M\), where \(h=\frac{1}{M}\) and we denote by (x, y) the spatial variable \(\varvec{x}\).

We denote by

and aim to approximate (5) by a spatial semi-discretization consisting of \(4M^2\) ordinary differential equations

for \(4M^2\) unknowns \(u_{k,i,j}(t)\in {\mathbb {R}}\), which are approximations of \( u_k(\varvec{x}_{i,j}, t)\). We define \(M\times M\) matrices \(\varrho , V_1, V_2, C\) as follows: for \(i,j=1,\dots ,M\)

Here, as in the rest of this section, we drop the dependence of U on t to obtain the cited spatial semidiscretization. The spatial semidiscretization is split as follows:

for the following approximations:

3.1.1 The Operator \(\mathcal {L}_1\)

The approximation of the only nonzero source term \({\mathcal {L}}_1(U)_{3,i,j} \approx \rho (x_{i,j},t) g\) is taken pointwise, i.e.,

3.1.2 The Operator \(\mathcal {C}\)

The convective term \({\mathcal {C}}\) is obtained through finite differences of numerical fluxes obtained by WENO5 reconstructions [3, 4] on Global Lax–Friedrichs flux splittings [27], which is fifth-order accurate for finite difference schemes, based on point values. Other schemes for systems of hyperbolic conservation laws, as MUSCL, could be used as well, see [29] and references therein, but we choose this one for its relative ease of implementation while achieving high resolution.

3.1.3 Basic Finite Difference Operators

We define finite difference operators for functions on \(M\times M\) grids, for fixed \(h>0\), to approximate first-order derivatives along the x dimension:

and the shift operator

and likewise in the y direction.

The properties of these operators are the following:

1. 1.

   \(D_{x}^{1*}f_{i,j}\) is a second-order accurate approximation for \(f_{x}(x_{i-\frac{1}{2},j})\) when \(f_{i,j}=f(x_{i,j})\) and \(f\in {\mathcal {C}}^3\) with \(f(x_{0,j})=f(x_{M,j})=0\), which is used to approximate pure double derivatives.

2. 2.

   \(D_{x}^1f_{i,j}\) is a second-order accurate approximation for \(f_{x}(x_{i+\frac{1}{2},j})\) when \(f_{i,j}=f(x_{i,j})\) and \(f\in {\mathcal {C}}^3\) with \(f_{x}(x_{M+\frac{1}{2},j})=0\), which is used to approximate pure double derivatives. These two operators are related as \(D_{x}^{1*}=-(D_{x}^1)^T\).

3. 3.

   \(D_{x}f_{i,j}\) is a second-order accurate approximation for \(f_{x}(x_{i,j})\) for \(1<i,j<M\) and first-order accurate otherwise, when \(f_{i,j}=f(x_{i,j})\) and \(f\in {\mathcal {C}}^3\).

4. 4.

   \(D_{x}^*f_{i,j}\) is a second-order accurate approximation for \(f_{x}(x_{i,j})\) for \(1<i,j<M\) or \(j=1,M\) and first-order accurate otherwise, when \(f_{i,j}=f(x_{i,j})\) and \(f\in {\mathcal {C}}^3\) with \(f(x_{i,\frac{1}{2}})=f(x_{i,M+\frac{1}{2}})=0\).

For the sake of notation, for functions f, g on \(M\times M\) grids we denote

3.1.4 The Operator \(\mathcal {L}_2\)

To approximate the terms \(\mathcal {L}_2(U)_{k,i,j}, k=2,3,\) that involve derivatives of c in the conservation of momenta, we consider the following expressions:

As an instance, for \(C=(c_{i,j})\):

Symmetrically:

These approximations take the boundary conditions (4) on c into account and are second-order accurate at interior points \(1<i,j<M\) and first-order accurate at the rest of the points.

3.1.5 The Operator \(\mathcal {L}_4\)

To approximate the terms \(\mathcal {L}_4(U)_{k,i,j}, k=2,3,\) that involve derivatives of \(\varvec{v}\) in the conservation of momenta, we consider the following finite difference approximations

for any \(\mathcal {C}^4\) function z on [0, 1] such that \(z(0)=z(1)=0\). This approximation is second-order accurate for \(1<i<M\) and first-order accurate otherwise.

These approximations lead to the \(M\times M\) blocks \({\mathcal {L}}_4(U)_2, {\mathcal {L}}_4(U)_3\)

where D, E are the matrices associated to the finite difference operators in (14) and (16), respectively:

(17)

The matrices in (17) fail to be symmetric due to the boundary conditions. This could be circumvented with staggered grids for the velocity, but we do not consider this possibility in this paper.

3.1.6 The Operator \(\mathcal {L}_3\)

The term \((\psi '(c) )_{xx}=(\psi ''(c)c_x)_x\) needs special care, since it is not negatively definite, due to \(\psi ''(c)=3c^2-1\) changing sign in \((-1,1)\). Since only negatively definite (right-hand side) terms should be treated implicitly, we consider the splitting \(\psi '=\phi _{+}+\phi _{-}\), with

as in [13, 30], where it is shown that time-stepping schemes for the Cahn–Hilliard equations that are built by considering the term containing \(\phi _+\) as implicitly treated and the term containing \(\phi _-\) as explicitly treated are unconditionally stable.

For \(\chi =\phi _{\pm }\), taking into account that \((\chi (c))_{z}=\chi '(c)c_z\), \(z=x,y\), so \(\chi (c)\) satisfies Neumann boundary conditions, we have the following second-order accurate approximations:

These yield approximations

where we denote by \(\mathcal {M}_{\pm }(C)\) the tensor built from the values of \(\chi '=\phi '_{\pm }\) that appear in (18).

It can be seen that the boundary conditions

are equivalent to

for \(\xi =\frac{1}{\rho }\varDelta c\).

If f(x, y) satisfies the Neumann boundary condition \(\nabla f(x,y)\cdot \varvec{n}(x,y)=0\) and \(f_{i,j}=f(x_{i,j})\) then \((\varDelta _h f)_{i,j}\), \(\varDelta _h=\varDelta _{x,h}+\varDelta _{y,h}\), yields a second-order accurate approximation of \(\varDelta f(x_{i,j})\) for \(f\in {\mathcal {C}}^4\), where

Therefore we get the second-order accurate approximation

which, together with (19), yields the second-order accurate approximation

where \(D(v), v\in {\mathbb {R}}^{M\times M},\) is the diagonal operator on \(M\times M\) matrices given by

3.2 Vector Implementation

An \(M\times M\) A matrix can be cast into a \(M^2\) (column) vector J(A) as follows:

i.e., J(A) is formed by stacking the columns of A.

In order not to clutter notation, we henceforth identify by X the vector J(X) for \(X=\varrho , V_1, V_2, C, \mathcal {C}(U), \mathcal {L}_k(U), k=1,2,3,4\). For \(M^2\) vectors a, b we denote by \(a*b\) the \(M^2\) vector such that

Using this notation, we reorder the 4 \(M\times M\) matrices functions that we have used in the previous section into a \(4M^2\) (column) vector function U(t) as follows

and, correspondingly:

so that (10) is equivalent to

with U(t) and \(\mathcal {L}(U(t))\) being \(4M^2\) vectors.

With this notation \({\mathcal {L}}_4(U)\) is given by

where \(I_M\) is the \(M\times M\) identity matrix, \(\otimes \) is the Kronecker product, the matrices D, E have been defined in (17) and \(\varvec{O}_{M^2}\) denotes the \(M^2\times M^2\) null matrix.

Likewise, \(\mathcal {L}_3(U)\) is given by

where \(\varvec{O}_{M}\) denotes the \(M^2\) null vector and \(\mathcal {M}_{+}(C),\varDelta _h\) are the \(M^2\times M^2\) matrices corresponding to the homonymous tensors in Subsection 3.1.6.

Remark 1

The numerical schemes for the one-dimensional case are obtained in a straightforward manner from those for the two-dimensional case.

3.3 Implicit-Explicit Schemes

An explicit Runge–Kutta solver for (21), as the Explicit Euler method

applied to obtain a fully discrete scheme would require \(\varDelta t \propto \varDelta x^4\) for stability, which would yield a prohibitively expensive numerical scheme.

Instead, we use the technique of doubling variables and partitioned Runge-Kutta schemes [5, 25] to obtain Linearly Implicit Explicit (LIMEX) schemes. We denote the variables that are to be treated explicitly with a tilde and we define

where

i.e., the coefficients of the matrix \(\mathcal {M}_+(\widetilde{C})\) are computed from the explicit data \(\widetilde{C}\), whereas all the term corresponding to the negative part \(\phi _{-}\) of the splitting \(\psi '=\phi _++\phi _{-}\) is computed from \(\widetilde{C}\).

From (20) and (11), the requirement \(\widetilde{\mathcal {L}}(U, U)={\mathcal {L}}(U)\) is met.

Now we have that the IVP

is equivalent to

A partitioned Runge–Kutta scheme, in which there are two different s stages Butcher tableaus, one explicit and one (diagonally) implicit

can be applied to (23). It can be seen (cf. [14, 25]) that if both Butcher tableaus yield second-order accurate schemes and \(\beta =\widetilde{\beta }\), then the resulting partitioned Runge-Kutta scheme is second-order accurate and that third-order accuracy is also achievable.

The s intermediate stages in the partitioned Runge–Kutta method are computed for \(i=1,\dots ,s\) as follows:

since \(\widetilde{\beta }_{j}=\beta _j\), \(\forall j\) if \(U^{n}=\widetilde{U}^{n}\), so there is no need of doubling variables, which we henceforth assume.

The definitive time-stepping expression is the following:

where

In this paper we consider Stiffly Accurate Runge–Kutta solvers, i.e., the last row of the \(\alpha \) matrix coincides with \(\beta ^T\). Specifically, we consider the following Butcher tableaus

The DIRKSA scheme is the only 2-stages second order stiffly accurate DIRK method with \(\alpha _{ij}\ge 0\). The order for EE-IE is 1 and it is 2 for *-DIRKSA.

3.4 Solution of the Nonlinear Systems

One needs to solve

for \(U^{(i)}\), where

As we shall see, although \({\mathcal {L}}_2, \widetilde{{\mathcal {L}}}_3, {\mathcal {L}}_4\) are not linear, only linear systems for \(V_k^{(i)}, C^{(i)}\) have to be solved.

For the first variable, with block superscript notation for the operators and \({\mathcal {K}}\) variables, we get:

so \(\varrho ^{(i)}\) is explicitly computable.

For the fourth variable \(Q^{(i)}\), since \(\varrho ^{(i)}\) is already known, this system can be cast for the \(C^{(i)}\) variables:

which is equivalent to

If \(\varrho ^{i}_k>0\,\forall k\), then the matrix of this system is symmetric and positive definite, since it is the sum of a diagonal positive matrix and two symmetric and positive semidefinite matrices.

For the second and third variables, one needs to solve:

Since \({\mathcal {L}}_1^j(U^{(i)})\) and \({\mathcal {L}}_2^j(U^{(i)}), j=2,3,\) do not depend on \(V_1^{(i)},V_2^{(i)}\), they can be computed from previous steps and there only remains the following equation to be solved:

If \(\varrho ^{(i)}_k>0\,\forall k\), then the matrix of this system should be close to symmetric and positive definite, since the matrix

is a discretization of the self-adjoint elliptic operator

under the boundary conditions (4).

3.5 Linear Solvers

We have used the multigrid V-cycle algorithm with 4 pre- and post- Gauss-Seidel smoothings and direct solution when the size of the projected systems is \(\le 4\) (see [7]) for the solution of systems (25) and (26).

The matrix in system (25) is symmetric and positive, so there is the possibility of using the conjugate gradient method, using the following approximation

as preconditioner, for \(\varDelta _h\) can be efficiently diagonalized by discrete cosine transforms.

We next analyze the condition number of the preconditioned matrix:

We drop the superindex \(^{(i)}\) for simplicity. The matrix \(\mathcal {M}_{+}(\widetilde{C})\) can be expressed as

Since \(\phi _+'=2\) and

we get the following:

Therefore, for \(0\ne z\in {\mathbb {R}}^{M^2}\):

therefore the condition of the preconditioned matrix is bounded above by

which is close to 1 if \(\varrho \) is nearly constant. Therefore, it is expected to be a good preconditioner in this case.

3.6 Time-Step Selection

The time-step stability restrictions of the purely convective subsystem is

where CFL\(^*\) is some number and the maximum of the characteristic speeds, \(\text {cs}\), is computed, at each time step, as

The scheme is not ensured to be bound preserving, i.e., it might happen that density become negative or the c-variable be outside \([-1,1]\). Purely convective models might develop vacuum regions and the high-order WENO reconstructions might create overshoots near high gradients. Coping with these difficulties is certainly challenging and out of the scope of this paper.

In our case, there is no guarantee that the solution of (25) be in \([-1,1]\). We have used in our simulation the strategy of decreasing CFL\(^*\) (by a factor of 1/2 in our implementation) when |c| reaches some threshold (1.5 in the experiments) and increasing it (by a factor of 2 in our implementation) until a prescribed maximum value CFL otherwise. The rationale behind this strategy is that decreasing \(\varDelta t\) asymptotically ensures bound preservation.

Although we do not have a proof that this strategy always succeeds, the simulations have been successful, in the sense that |c| has been kept always below the threshold, at the expense of local reductions of the parameter CFL\(^*\).

The objectives of the experiments in this section are the following:

1. 1.

   Showing that the order of the global errors in some experiments coincides with the expected design order of the scheme used to obtain them.

2. 2.

   Showing that some IMEX schemes can perform time steps \(\varDelta t\) with the same stability restrictions as the purely convective subsystem, see (27).

3. 3.

   Testing the behavior of different issues for the algorithms, such as conservation, number of iterations for the linear solvers, etc.

In all numerical experiments, the adiabatic constant \(\gamma \) has been set to 5/3 and the constant \(C_p\) to 1, unless otherwise stated.

All the results have been obtained with a C++ implementation, using the GNU C++ compiler with optimizations -O3 and running in a single core of an AMD EPYC 7282 3.0 GHz CPU. The matrices of systems (25) and (26) are stored by diagonals.

4.1 One-Dimensional Experiments

4.1.1 Test 1

The purpose of this test is to show that the Explicit Euler solver needs that \(\varDelta t\) be proportional to \(\varDelta x^4\) and that our proposal is stable for \(\varDelta t\) selected by (27) with \(CFL^*=1\) for one-dimensional simulations.

We consider the following initial condition for a one-dimensional test: (\(c_0\) is in the unstable region \((-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}})\))

with parameters \(g=-10\), \(\nu _*=2\nu +\lambda =2\), \(\varepsilon =10^{-4}\). The Explicit Euler scheme (22) blows up for \(M=4000\), and \(\varDelta t=\varDelta x^3\) for \(t\approx { 10^{-10}}\) , whereas, with the same setup and \(\varDelta t = \varDelta x^4\), the code has taken 640000 time steps to get to a simulation time of \(t= 10^{-9}\) without blow up.

The EE-IE and *-DIRKSA blow up for \(M=100\) and \(\varDelta t\) computed by (27) for \(\text {CFL}^*=1.1\) for \(t\approx 10^{-1}\), whereas they do not for \(\text {CFL}^*=1\) and \(M=10000\).

4.1.2 Test 2, 3 and

The purpose of these tests is to have a comparison with tests taken from [16] and the effect of using \(\gamma < 1.5\). The setup of the test requires some changes to the equations, due to a different scaling of the Ginzburg–Landau free energy. Specifically, we obtain approximate solutions for the following system of equations:

for \(\nu ^*=1\), \(\varepsilon =10^{-3}\), with \(M=100\), simulated time \(T=100\) and the following initial conditions

1. 1.

   Test 2 (Example 3 in [16]), \(\gamma =1.5\):

   $$\begin{aligned} \rho _0(x)&=0.9+\frac{1}{1000} \cos (2 \pi (x-1)), \\ v_0(x)&= \frac{1}{2} e^{-50 (x-0.5)^2},\\ c_0(x)&=\tanh \left( \frac{x-0.5}{0.1 \sqrt{5}}\right) . \end{aligned}$$
2. 2.

   Test 3 (Example 4 in [16]), \(\gamma =1.5\):

   $$\begin{aligned} \rho _0(x)&=0.9+\frac{1}{1000} \cos (2 \pi (x-1)), \\ v_0(x)&= \frac{1}{2} e^{-50 (x-0.5)^2},\\ c_0(x)&=0.75+\frac{1}{5}\tanh \left( \frac{x-0.5}{0.1 \sqrt{5}}\right) . \end{aligned}$$
3. 3.

   Test 4 (Example 3 in [16]), \(\gamma =1.4\):

   $$\begin{aligned} \rho _0(x)&=0.9+\frac{1}{1000} \cos (2 \pi (x-1)), \\ v_0(x)&= \frac{1}{2} e^{-50 (x-0.5)^2},\\ c_0(x)&=\tanh \left( \frac{x-0.5}{0.1 \sqrt{5}}\right) . \end{aligned}$$

Notice that \(c_0\) does not satisfy the boundary conditions, which are, in this case:

The results of the comparisons of Test 2 and Test 3 are displayed in Fig. 1 and show some qualitative similarity with those from [16]. We cannot explain the quantitative difference between our results and those from the cited reference.

Tests 2 and 3: Test 2 corresponds to the Example 3 in [16] and displays phase separation and Test 3 corresponds to the Example 4 in [16] and does not display phase separation

Full size image

The results of the comparisons of Test 2 and Test 4 are displayed in Fig. 2 and do not show much difference between using \(\gamma =1.5\) and \(\gamma =1.4\).

Test 4: comparison of Example 3 in [16] with \(\gamma =1.5\) and \(\gamma =1.4\)

Full size image

4.2 Two-Dimensional Tests

4.2.1 Order Test

This test aims to assess that the *-DIRKSA method achieves second-order accuracy in the global errors. For this purpose, we add a forcing term to the equations so that the solution is prescribed. Specifically, the solution in this case is

Notice that these functions satisfy the boundary conditions (4). The equations that are solved are the following:

The parameters that have been used for this test are the following:

For these tests, we have used \(\varDelta t\) given by (27) with CFL=0.4.

For \(M\times M\) grids, with \(M=2^{l}, l=3,\dots ,8\), the global errors for the approximations \(u_{k,i,j}^n\) obtained by the *-DIRKSA method for \(t_n=T=0.01\), are computed as

and are displayed in Table 1, where it can be observed the convergence of the experimental orders of convergence \(\text {EOC}_{M}=\log _2\left( e_{M}/e_{2M}\right) \) towards 2, which is compatible with the scheme being second-order accurate. The analogous experiment is performed for the EE-EI scheme resulting in experimental orders of convergence that decrease away from 2.

4.2.2 Tests 1, 2 and 3

For the following two-dimensional tests, we have used \(\varDelta t\) given by (27) with CFL=0.4, which is a safe setup for simulations with the corresponding explicit schemes for the convective part only (isentropic Euler equations).

We consider the following tests:

• Test 1 (\(c_0\) is in the unstable region \((-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}})\)):

  $$\begin{aligned}&\rho _0(x, y)=0.1\cos (2 \pi x) \cos (\pi y)+1.25\\&\varvec{v}_{0}(x, y) =( \sin (\pi x) \sin (\pi y), \sin (\pi x) \sin (2\pi y))\\&c_0(x, y)=0.1\cos ( \pi x) \cos (\pi y) \end{aligned}$$

• Test 2 (\(c_0\) is in the stable region):

  $$\begin{aligned}&\rho _0(x, y)=0.1\cos (2 \pi x) \cos (\pi y)+1.25\\&\varvec{v}_{0}(x, y) =( \sin (\pi x) \sin (\pi y), \sin (\pi x) \sin (2\pi y))\\&c_0(x, y)=0.75+0.1\cos ( \pi x) \cos (\pi y) \end{aligned}$$

• Test 3: \(\rho =1, \varvec{v}_0=0,\) \(c_0\) has been set to a uniform random sample of 0 mean and \(10^{-10}\) standard deviation.

Notice that these functions satisfy the boundary conditions (4) (for test 3, almost within roundoff error).

In Fig. 3, we show the time evolution of the conservation errors for \(\rho \) and \(q=\rho c\) for Test 3, and parameters \(\nu =10^{-3}, \lambda =10^{-4}, {\varepsilon }=10^{-4}, M=256\), with linear systems solved by the V-cycle multigrid algorithm with relative decrease of residual of \(10^{-12}\) as stopping criterion. Specifically, we approximate

Conservation errors for test 3

Full size image

In Fig. 4 we show the time evolution of the CFL\(^*\) parameter, according to subsection 3.6, for a case for Test 1, with parameters \(\nu =10^{-3},\lambda =10^{-4}, \varepsilon =10^{-5}, M=256\), with the largest number of variations in CFL\(^*\), according to subsection 3.6, in all test cases.

Test 1. Left: Time evolution of the CFL parameter; Right: Time evolution of \(\min c\), \(\max c\)

Full size image

We show in Figs. 5, 6, 7, 8, 9, 10, 11, 12, 13 and 14 some snapshots of the results obtained for all the tests with *-DIRKSA, \(M=256\), \(g=-10\), \(\nu =10^{-3}, \lambda =10^{-4}, \varepsilon =10^{-4}\), which correspond to flows with Reynolds number roughly in the range \([10^2, 10^3]\).

In Fig. 5, it can be seen that the c-component in the initial condition for Test 1 lies entirely within the spinodal region \((-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}})\). Therefore, at the early stages of the simulation, at \(T=0.1\), separation occurs forming complex patterns. Meanwhile, it can be appreciated from the pictures corresponding to the \(\rho \)-variable that gravity is acting so that density increases at the bottom boundary, \(y=0\). This trend continues in Fig. 6, where it can be appreciated that the maximal density continues increasing, a sharp upgoing front develops for all the variables, but c, where the diffuse interface experiments many topological changes, with growing regions, a phenomenon named nucleation. At the final stages of the simulation, as seen in Fig. 7, the maximal density decreases, nucleation continues and the bulk flow enters into a seemingly turbulent regime.

Lastly, in Fig. 8 the results for \(T=1.0\) and \(M=64\) and \(M=128\) are displayed in order to be compared with the corresponding result in Fig. 7. As empirically expected, there are larger differences in the c variable than in the others.

Results for Test 1. Left: Initial condition, with c-variable inside spinodal region; Right: Results for \(T=0.1\), where density increases at the bottom and separation is clearly visible in the c-variable

Full size image

Results for Test 1, \(T=0.3\) (left) and \(T=0.5\) (right) where it can be seen that density continues increasing at the bottom, forming an upgoing front, and nucleation is beginning, as seen in the c-variable

Full size image

Results for Test1, \(T=0.7\) (left) and \(T=1.0\) (right) where it can be seen in the velocity that vorticity has developed and nucleation is increasing, as seen in the c-variable

Full size image

Results for Test1, \(T=1.0\) \(M=64\) (left) and \(T=1.0\) \(M=128\) (right)

Full size image

In Fig. 9, it can be seen that the c-component in the initial condition for Test 2 lies above the spinodal region. Therefore, as seen in the pictures for the c-variable in the snapshots shown in this Figure and also in Fig. 10, the fluid remains almost homogeneous, with c tending to 3/4 in the whole domain, which corresponds to mass fractions \(c_1=\frac{7}{8}, c_{2}=\frac{1}{8}\), which are exactly the initial proportions of the individual species. This fact means that the rest of the equations behave like a uniform fluid governed by the compressible Navier–Stokes equations under gravitation, with Reynolds number high enough for a seemingly turbulent regime.

Results for Test 2. Left: Initial condition, with c-variable above the spinodal region; Right: Results for \(T=0.3\), where it can be seen that density continues increasing at the bottom, forming an upgoing front, while the c-variable is converging towards 0.75

Full size image

Results for Test 2, \(T=0.6\) (left) and \(T=1.0\) (right) where it can be seen in the velocity that vorticity has developed, while the c-variable has almost fully converged to 0.75

Full size image

In Fig. 11, it can be seen that the c-component in the initial condition for Test 3 is almost 0, thus lies entirely within the spinodal region \((-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}})\). Therefore, at the early stages of the simulation, at \(T=0.01\), a typical spinodal decomposition begins appearing in the form of a medium-frequency pattern, corresponding to a solution as in (8)-(9), for \(k_1, k_2\in \mathbb {N}\) that minimize

corresponding to the expression in (9), for \(c_0=0\), taking into account that \(\psi ''(c_0)=-1\).

The rest of the simulation can be seen in Figs. 12, 13, 14, where the spinodal decomposition continues until nucleation. Density increases at the bottom, but no clear turbulence is appreciated, maybe due to the simulation having been carried only until \(T=0.29\).

Results for Test 3. Left: Initial condition, where it can be seen that the c-component is in the range \((-5\, 10^{-11}, 5\, 10^{-11})\), thus lies entirely within the spinodal region \((-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}})\); Right: for \(T=0.01\) onset of a spinodal decomposition begins appearing with an amplitude around \(10^{-8}\)

Full size image

Results for Test 3, where it can be seen that density increases at the bottom and the spinodal decomposition continues its development with an amplitude around \(4\, 10^{-5}\) for \(T=0.02\) (left) and 0.2 for \(T=0.03\) (right)

Full size image

Results for Test 3, where it can be seen that density continues increasing at the bottom and the spinodal decomposition has given a fully developed separation pattern for \(T=0.04\) (left), which has continued towards the onset of nucleation for \(T=0.14\) (right)

Full size image

Results for Test 3, where nucleation is visible in the c-variable near (0.9, 0.5), with an ever-shrinking region of the second fluid disappearing

Full size image

Results for Test 3 for times \(T=0.4,0.6,0.8,1.0\)

Full size image

4.2.3 Test 4

This test is taken from [10] and simulates the evolution dynamics of two merging bubbles, which are initially tangent. Specifically, the initial condition is

The parameters in our model are \(\varepsilon =10^{-3}, C_p=1000, \gamma =5/3\). We have tried two viscosity parameters \(\lambda =\nu =0.1\) and \(\lambda =\nu =1\).

There are substantial differences between the results reported here and those in the cited reference. This may be explained because our model is compressible and that in the cited reference is incompressible, besides the additional terms in the Cauchy stress tensor being radically different. The most prominent difference is that there seems to be a factor of 2500 between the simulated times reported here and there. Other issue is that there are little differences when changing the viscosity parameter.

The 0 level set of the c variable at times \(T=0,0.002, 0.008,0.04,0.08, 0.12\) (which are 2500 times smaller than those times in the cited reference \(T=0,5,20,\) 100, 200, 300) is displayed in Fig. 16. In Fig. 17 the velocity field and the pressure fluctuations about there mean value are displayed for \(T=0.002, 0.008,\) 0.04, 0.12

Results for Test 4: 0-level set of c for \(T=0,0.002, 0.008,0.04,0.08, 0.12\)

Full size image

Results for Test 4: the velocity field and the pressure fluctuations about there mean value for \(T=0.002, 0.008,0.04, 0.12\)

Full size image

4.3 Linear Solvers

The multigrid solver mentioned in subsection 3.5 for the solution of (26) has shown a very satisfactory performance throughout all tests, requiring almost always only one iteration to achieve convergence to double precision digits.

The performance of the solvers for the solution of (25) is more complex. In Tables 2, 3 and 4 we show the average number of multigrid iterations to solve (25) for tests 1, 2, 3, respectively, with \(\nu =10^{-2},10^{-3},10^{-4}\) and \(\varepsilon =10^{-3},10^{-4},10^{-5}\), for a stopping criterion based on a relative decrease of the residual by a factor of \(10^{-6}\).

It can be deduced from these tables that the number of iterations grows slowly with M for all cases.

In Table 5 we show a comparison between the multigrid and the preconditioned conjugate gradient solvers in subsection 3.5. The results have been obtained with *-DIRKSA, \(M=16,32,64,128,256\), \(g=-10\), \(\nu =10^{-2}, \lambda =10^{-3}, \varepsilon =10^{-4}\). It can be deduced that the preconditioned conjugate gradient solver uses less CPU time than the multigrid solver, although the latter takes fewer iterations than the former.

In this paper, we propose efficient linearly implicit-explicit schemes for the two-dimensional compressible isentropic Cahn–Hilliard–Navier–Stokes equations. Some tests are performed to show that they achieve second-order accuracy under time-step stability restrictions dictated only by the convective part of the equations.

If the pressure law is not stiff, then the proposal yields an efficient time-stepping method to solve the isentropic Cahn–Hilliard–Navier–Stokes equations. For large pressure coefficient \(C_p\), the time steps have to be smaller, so including the pressure term in the implicitly treated terms might allow again larger time steps, dictated only by the velocity (akin to low-Mach-number solvers [2]).

As future research, we plan to extend these techniques to other, stiffer pressure laws, and to a three-dimensional setting with Galerkin techniques. We also plan the extension of these techniques to quasi-incompressible models (see [22]).

A crucial part of the algorithms is the iterative linear solvers used for solving the system related to the Cahn–Hilliard subequation. We plan to explore the possibility of using the multigrid solver as preconditioner for the conjugate gradient solver.

Third-order accuracy can be obtained by using fourth-order accurate finite-differences for the spatial discretization and third-order accurate partitioned Runge–Kutta schemes, at the expense of larger stencils and four Runge–Kutta stages.

The code can be parallelized by decomposing the computational domain in smaller rectangles and letting a processor do the job in each rectangle. Communication would be needed for data associated to neighboring cells in different rectangles. More global communication would be needed in case of using multigrid solvers.

The use of bound-preserving high-order reconstructions for discretizing the convective part is planned as a possible mechanism for reducing the regions where \(|c|>1\). Exploring situations and conditions for the solution of the Cahn–Hilliard subsystem (25) being in \([-1, 1]\) will be also undertaken.