A Robust and Accurate Adaptive Approximation Method for a Diffuse-Interface Model of Binary-Fluid Flows

We present an adaptive simulation framework for binary-fluid flows, based on the Abels-Garcke-Gr\"un Navier-Stokes-Cahn-Hilliard (AGG NSCH) diffuse-interface model. The adaptive-refinement procedure is guided by a two-level hierarchical a-posteriori error estimate, and it effectively resolves the spatial multiscale behavior of the diffuse-interface model. To improve the robustness of the solution procedure and avoid severe time-step restrictions for small-interface thicknesses, we introduce an $\varepsilon$-continuation procedure, in which the diffuse interface thickness ($\varepsilon$) are enlarged on coarse meshes, and the mobility is scaled accordingly. To further accelerate the computations and improve robustness, we apply a modified Backward Euler scheme in the initial stages of the adaptive-refinement procedure in each time step, and a Crank--Nicolson scheme in the final stages of the refinement procedure. To enhance the robustness of the nonlinear solution procedure, we introduce a partitioned solution procedure for the linear tangent problems in Newton's method, based on a decomposition of the NSCH system into its NS and CH subsystems. We conduct a systematic investigation of the conditioning of the monolithic NSCH tangent matrix and of its NS and CH subsystems for a representative 2D model problem. To illustrate the properties of the presented adaptive simulation framework, we present numerical results for a 2D oscillating water droplet suspended in air, and we validate the obtained results versus those of a corresponding sharp-interface model.


Introduction
Binary-fluid flows in which the two fluid components are separated by a molecular transition layer are omnipresent in science and engineering. An important high-tech application of (typically, aqueous-aereous) binary fluids pertains to inkjet printing, where functionalized minuscule liquid droplets are deposited on a substrate. Inkjet printing is in fact the most widespread technological application of microfluidics [36]. Inkjet-printing technology exhibits exponential growth in a similar manner as Moore's law in semi-conductor technology [36], and to sustain this technological development, it is imperative to enable the deposition of increasingly many, increasingly small droplets per unit time. Modeling and simulation of binary-fluid flows in inkjet printing are essential to support the necessary technological advancements. Due to the small spatial and temporal length scales of the jetting process, in combination with the complexity of the involved microfluidic processes, such as wetting on fluid-solid interfaces, interface coalescence and break-up, phase transition, surfactant transport and dynamic surface-energy phenomena, etc., advanced numerical simulation techniques form an indispensable research instrument, complementary to experimental investigations.
Mathematical-physical models for binary-fluid flows can generally be subsumed under two categories, viz. sharp-interface and diffuse-interface models. In sharp-interface models, the transition layer between the two fluid components is modeled as a manifold of co-dimension one (e.g. a 2D surface in a 3D ambient space), at which the interaction of the two fluid components is accounted for by kinematic and dynamic (possibly extended by thermal and thermodynamic) interface conditions. In diffuse-interface models, the transition layer is represented as a thinbut-finite layer, in which the two components are mixed in a proportion that varies continuously and monotonously between the two pure species across the transition layer. Diffuse-interface models are significantly more versatile than their sharp-interface counterpart, in the sense that diffuse-interface models can intrinsically provide a description of phenomena that cannot be intrinsically accounted for in sharp-interface models, such as topological changes of the fluid-fluid interface due to coalescence or break-up [37], and wetting of solid substrates [55,44,32]. Diffuse-interface models therefore provide a modeling paradigm par excellence for complex microfluidic processes, e.g. those occurring in inkjet printing. In situations where both models are valid, the sharp-interface model can generally be obtained in the so-called sharpinterface limit, i.e. in the limit as the diffuse-interface thickness passes to zero [3,37,56]. It is to be noted, however, that contemporary theory on sharp-interface limits is still incomplete, and that for instance scaling of model parameters is presently unresolved. By virtue of their intrinsic representation of wetting phenomena, diffuse-interface models have recently also been propounded as a compelling modeling option for binary fluids in elasto-capillary fluid-solid interaction [50,47,51,5,6,16].
Diffuse-interface models for two immiscible incompressible fluid species are generally described by the Navier-Stokes-Cahn-Hilliard (NSCH) equations. The NSCH system is presently a class of models rather than a single model, as various variants exist. Historically, one of the first instances of the NSCH model is the so-called "model-H", presented by Hohenberg and Halperin [30] in the late 1970s. This model implicitly assumes that the densities of the two fluid species match, which leads to a solenoidal (incompressible) mixture-velocity field. Since its introduction, the Hohenberg-Halperin model has also often been used in non-matching density scenarios. However, in this case the model can violate thermodynamic consistency. An NSCH model that allows for non-matching densities was proposed in the late 1990s by Lowengrub and Truskinovsky [37]. This model is based on a barycentric mixture velocity and the mass fraction features as order parameter (phase variable). The resulting NSCH model is quasi-incompressible, i.e. the mixture velocity in this model is generally non-solenoidal. A new quasi-incompressible NSCH system in which the volume fraction acts as order parameter, was introduced by Shokrpour et al. in [46]. This model was developed in conjunction with a linearly-implicit, thermodynamically consistent time integration scheme. A deficiency of the quasi-incompressible NSCH models in [37,46] is that these models are not generally compatible with the underlying Navier-Stokes equations in pure-species regimes: both models insist that the pressure is harmonic if the phase variable is uniform. This limitation can be avoided by introducing a degenerate mobility, but this comes at the expense of severe complications in analytical and numerical treatments of the model. Abels, Garcke and Grün presented a thermodynamically consistent NSCH model for binary-fluids with non-matching densities, based on volume-averaged mixture velocity [4]. By virtue of the volume-averaged mixture-velocity premise, the mixture velocity in the AGG NSCH model is solenoidal. Mass conservation and thermodynamic consistency are inbuilt in the AGG NSCH model via an auxiliary relative mass flux in the balance laws for mass and linear momentum.
Despite the significant progress that has been made in the formulation and analysis of NSCH models for binary-fluid flows, and in the development of numerical techniques for these models, many outstanding challenges remain. One of the main challenges in numerical approximation of the NSCH system, pertains to maintaining robustness and efficiency of the simulations for realistic settings of the model parameters, e.g. for binary fluids comprising water and air, due to: 1. spatial multiscale behavior : In most applications, the thickness of the fluid-fluid interface is many orders of magnitude smaller than other relevant length scales in the problem under consideration, e.g. the radius of a droplet or of a liquid jet. For instance, x-ray measurements [15] and molecular-dynamics simulations [49] have conveyed that the thickness of a liquid-vapor interface is typically a fewÅ, while the radius of a pico-liter droplet, which are the smallest droplets used in current technologies, is approximately 10 µm. This corresponds to a scale difference of more than 10 4 , which is beyond the resolution capabilities of standard (non-adaptive) approximation methods. For many applications it is unnecessary to adhere to the actual value of the interface thickness, and it suffices that the interface thickness is sufficiently small relative to other characteristic length scales. A ratio of 1:100 is typically indicated as a guideline [56]. Nonetheless, the resulting resolution requirements are generally prohibitive if the problem under consideration exhibits additional multiscale behavior, e.g. due to the occurrence of satellite droplets [38]. 2. temporal multiscale behavior : Typically, the diffusion process at the diffuse interface carries a time scale that is significantly smaller than other characteristic time scales in the problem under consideration, e.g. the period of oscillation of a droplet. In fact, it is conjectured that in certain cases the diffusive time scale needs to pass to zero in the limit as the diffuse-interface thickness passes to zero [3,56]. Efficiency of a time-integration scheme for NSCH systems therefore requires that such a scheme is uniformly stable with respect to the interface-diffusion time scale. Recently, there has been much progress in the area of unconditionally thermodynamically stable time integration schemes for NSCH systems; see, e.g., [46,22] and references therein. These methods however tend to be very dissipative, and their accuracy is inadequate to efficiently solve problems with weakly dissipative dynamics, such as oscillations of low-viscosity droplets. 3. ε-conditional stability: A further complication related to time integration schemes for NSCH systems, is that the admissible time step is generally restricted by a CFL type condition in relation to the diffuse-interface thickness, i.e. the admissible time step is limited by ε/v, where ε indicates the diffuse-interface thickness parameter, and v the (normal) velocity of the interface (see Section 4.1 below). This restriction may emerge from conditional stability of the time integration scheme, or from potential non-robustness of the iterative solution procedure for implicit time-integration methods. This implies that to retain robustness, increasingly small time steps are required as the diffuse-interface thickness decreases, while the (macroscopic) time scale related to the interface motion is invariant. Such a CFL-type condition leads to prohibitive computational expenses in the sharp-interface limit ε → +0. 4. heterogenous parameters: In practice, the components of the binary fluid can exhibit significantly different properties, leading to large parameter variations in the NSCH system. For a binary fluid composed of water and air, for instance, the density and viscosity ratios are approximately 1000:1 and 100:1, respectively. This can cause ill-conditioning in the algebraic systems emerging from discrete approximations of the NSCH system. Moreover, in combination with non-compliance of the phase variable, large density and viscosity contrasts can lead to inadmissible negative densities and viscosities of the mixture. Non-compliance refers to the fact that the phase variable can assume values outside of the admissible codomain. The phase variable generally represents a difference of fractions (or a fraction), and its codomain is restricted to the interval [−1, 1] (resp. [0, 1]). The codomain of the phase variable must however be extended, thus including inadmissible values, to ensure existence of a solution to the NSCH system; see Ref. [26]. With standard constitutive relations for the density and viscosity, these inadmissible values of the phase variable can yield negative mixture densities or viscosities. 5. incompatibility and spurious velocities: Equilibrium solutions of the NSCH system are generally characterized by a particular profile of the phase variable across the diffuse interface (typically, tanh(·)) and vanishing velocity. However, as the particular equilibrium interface profile is generally incompatible with the discrete approximation space, discretization errors in the phase variable occur, which manifest as discretization errors in the interface stresses (in particular, the Korteweg stress tensor), and induce discretization errors in the velocity field. These discretization errors in the velocity field are commonly referred to as spurious velocities or parasitic currents [41]. These spurious velocities scale inversely proportional with the viscosity, and can become excessively large if the viscosity in one of the fluid components is small. 6. ill-conditioning of linear systems: Implicit time-integration methods for NSCH systems generally require the solution of a sequence of nonlinear algebraic systems of equations. The solution of each of these nonlinear problems by Newton's method, in turn, requires the solution of a sequence of linear algebraic problems, corresponding to the linear tangent problems in Newton's method. The linear-algebraic systems emerging from the NSCH equations are generally hard, owing to the intrinsic asymmetry of these systems, the aggregation of disparate subsystems (NS and CH), the occurrence of small parameters, and the occurrence of heterogeneous coefficients that may differ by orders of magnitude; see item 4 above.
In this work, we present an adaptive simulation framework for binary fluids based on the AGG NSCH model, which effectively addresses the aforementioned issues. To resolve the spatial multiscale behavior related to thin diffuse interfaces, we consider adaptive approximations based on high-regularity truncated hierarchical B-splines [31,19]. In particular, we utilize C 0 continuous cubic spline bases for all variables except the pressure field, which is approximated by C 0 quadratic splines. The velocity-pressure pair thus corresponds to a member of the Taylor-Hood family. The adaptive-refinement procedure is guided by a two-level hierarchical aposteriori error estimate. The error estimate is used in a support-based marking-and-refinement procedure [34]. The adaptive procedure efficiently controls discretization errors and, hence, intrinsically suppresses spurious velocities. To enhance the robustness of the Newton solution process throughout the iterative-refinement procedure, we introduce an ε-continuation approach in which the diffuse-interface parameter (ε) and mobility (m) are adapted to the resolution of the mesh. This continuation procedure moreover serves to bypass the severe time-step restriction that would otherwise occur for small ε (see item 3 above). To determine the appropriate scaling of m := m ε with respect to ε, we derive the characteristic diffusive time scale corresponding to the Cahn-Hilliard equations, which governs the rate of contraction (resp. expansion) of the diffuse interface caused by decreasing (resp. increasing) of ε. To further accelerate the computations and improve robustness, we apply a Backward Euler scheme with a second-order contractive/expansive splitting of the double-well potential [54] in the initial stages of the adaptive-refinement procedure in each time step, and a Crank-Nicolson scheme in the ultimate stages to restore second-order accuracy in time. Instability due to non-compliance of the phase variable in combination with large density and viscosity contrasts, is circumvented by applying a soft-clipped interpolation of the mixture density [26,14] and an Arrhenius interpolation of the mixture viscosity [50,8]. To avoid non-robustness of the Newton solution procedure due to ill-conditioning of the linear tangent problems, we introduce a partitioned solution procedure in each time step and each iteration of the adaptive-refinement process, where we make use of the inherent composition of the NSCH system of NS and CH subsystems. In addition, we present a systematic investigation of the conditioning of the monolithic NSCH tangent matrix and of its NS and CH subsystems, based on numerical computations of an oscillating liquid droplet in an ambient fluid in 2D. For this test case, we also examine the convergence behavior of preconditioned GMRES, regarding two different preconditioners, corresponding to two different partitions of the NSCH system into its NS and CH subsystems. To illustrate the properties of the presented adaptive simulation framework, we finally present numerical results for a 2D oscillating droplet with parameter values representative of a binary fluid composed of water and air, in which the interface-thickness parameter is approximately 70× smaller than the minimal radius of curvature of the interface, and we validate the AGG NSCH model by comparing the obtained results with those of a corresponding sharp-interface model. This paper is structured as follows: In Section 2 we present the strong and weak forms of the diffuse-interface AGG NSCH model, together with the applied constitutive relations. In Section 3 we introduce the approximation setting and the spatial and temporal discretization schemes, followed by a presentation of the error-estimation procedure and the adaptiverefinement process. Section 4 considers the Newton solution procedure for the nonlinear algebraic systems emerging from the discretization, and introduces the ε-continuation process and the partitioned solution method for the linear tangent problems. Section 5 presents numerical experiments. Finally, in Section 6 we conclude with a summary and a discussion.

Problem statement
We consider a binary fluid consisting of two immiscible, incompressible species separated by a finite layer which is comprised of a mixture of those species. The setting for the binary fluid under consideration is a liquid droplet (l) submerged in an ambient fluid (a); see the illustration in Figure 1. To describe the motion of the binary fluid, we consider the incompressible Navier-Stokes-Cahn-Hilliard (NSCH) model as presented by Abels, Garcke and Grün in [4]. This NSCH model is commonly referred as the AGG model, after its originators. It is noteworthy that alternate NSCH models exist, which can generally be classified into incompressible models [30,20,45] and quasi-incompressible models [46,37]. Our choice for the AGG model is motivated by its thermodynamic consistency for non-matching species densities and the convenience of implementation related to the incompressibility condition. It is also notable that the AGG model consistently reduces to the underlying single-species Navier-Stokes equations in the pure species setting. Well-posedness of the AGG model in various settings has recenty been established in [24,1,2]. Let us also note the application of the AGG NSCH model in elastocapillarity in [50].

NSCH model
To define a setting for the NSCH model, we regard an open time interval (0, t fin ) ⊆ R >0 and a spatial domain corresponding to a simply connected time-independent subset Ω ⊂ R d (d = 2, 3). The two fluid species l and a are both present according to the volume fractions ϕ l , ϕ a : (0, t fin ) × Ω → [0, 1] such that ϕ l + ϕ a = 1. The order parameter ϕ, identifying the presence of each species, is defined as from which the identities ϕ l = (1 + ϕ)/2 and ϕ a = (1 − ϕ)/2 follow. Thus, ϕ = 1 and ϕ = −1 correspond to pure l-species and a-species, respectively, while ϕ ∈ (−1, 1) corresponds to a mixture of the species. The mixture density ρ : (0, t fin ) × Ω → R >0 , where ρ l , ρ a > 0 denote the specific densities of the two species, is related to the volume fractions by ρ = ρ l ϕ l + ρ a ϕ a , and thus In the AGG NSCH model, the mixture velocity u : (0, t fin )×Ω → R d is volume-averaged, which leads to a solenoidal mixture velocity. The equation of motion contains an auxiliary mass flux. 6 The AGG NSCH diffuse-interface model [4] is given by: Here, is the auxiliary relative mass flux [4], Ψ (ϕ) = 1 4 (ϕ 2 − 1) 2 is the standard polynomial doublewell potential for the energy density of the mixture, and the variables p : (0, t fin ) × Ω → R and µ : (0, t fin ) × Ω → R represent the fluid pressure and chemical potential, respectively. We assume the mobility coefficient, m > 0, to be constant. We note that a degenerate mobility as function of the phase indicator ϕ, i.e. m := m(ϕ) with m : [−1, 1] → R ≥0 such that m vanishes at and only at ϕ = ±1, is preferable from a modeling perspective [9], as it eliminates the Ostwald-ripening effect in pure species. However, the introduction of a degenerate mobility complicates numerical-approximation procedures [9], and we therefore opt to retain a nondegenerate constant mobility. The fluid stress consists of a viscous component, τ , a capillary component, ζ, and the isotropic component, pI. The viscous stress tensor is defined by the constitutive relation where ε := 1 2 (∇u+(∇u) T ) is the symmetric strain-rate tensor and tr (ε) = ∇·u is its trace. We note that for the AGG NSCH model, the second term in (5) can be omitted on account of (3) 2 . We have however retained this term (with Λ = −2/d in accordance with Stokes' hypothesis; see [51]) and the term does not vanish for the weakly-solenoidal Taylor-Hood finite-element approximations considered in this work. The mixture viscosity η := η (ϕ) is defined according to the Arrhenius relation Λ η = ρlma ρaml is the intrinsic liquid-ambient volume ratio, where m l and m a are the molar densities corresponding to species l and a, respectively; see also [50]. Furthermore, the constitutive relation for the capillary stress tensor is given by The parameter σ is a rescaling of the liquid-ambient surface tension σ la according to 2 √ 2σ = 3σ la . The model parameter ε > 0 dictates the transverse length scale (or thickness) of the diffuse interface between the two species. It is noteworthy that in (3) the separate Navier-Stokes and Cahn-Hilliard equations are coupled through the relative mass flux, capillary stress tensor, and advection term ∇ · (ϕu).
Equation (1) insists that the codomain of ϕ is restricted to [−1, 1], to enable an unambiguous physical interpretation. We must however desist from this restriction on the codomain of the 7 order parameter to ensure existence of a solution to the NSCH system; see Ref. [26]. This implies that ϕ can assume values outside the interval [−1, 1], and accordingly the mixture density in (2) can assume negative values. To avoid the occurrence of negative densities, we extend the definition of the mixture density to all of R via a so-called soft-clipped linear interpolation. As the assignment of the labels l and a to the fluid species is in principle arbitrary, without loss of generality we assume ρ l > ρ a . Defining the regularization parameter δ = ρ a / (ρ l − ρ a ), we consider the following extended density function: see also [26,14]. The mixture-density definition (7) coincides with (2) in the interval [−1, 1], but is smoothly extended in such a manner that negative densities are prevented. A general treatment of auxiliary conditions for the NSCH system (3) is beyond the scope of this work, and we restrict ourselves here to a consideration of the initial and boundary conditions as applied in the numerical experiments in Section 5. The considered test cases pertain to a liquid droplet suspended in an ambient fluid in R 2 ; see the illustration in Figure 1. We consider configurations that display bilateral symmetry along two axes and, accordingly, only part of the domain is modeled and the symmetry is accounted for by corresponding boundary conditions. The boundary Γ := ∂Ω of the domain Ω is partitioned into two complementary disjoint open subsets, corresponding to the symmetry boundary Γ sym ⊂ Γ and the external boundary Γ ext := int(Γ \ Γ sym ). On the symmetry boundaries, we impose a free-slip condition for the velocity field and homogeneous Neumann conditions for the order parameter and the chemical potential: Let us note that the homogeneous Neumann condition on the order parameter ϕ implies that intersections of the diffuse interface with the symmetry boundary occur at a 90 • contact angle. On the lateral boundary Γ ext , we impose a no-slip condition, along with the same homogeneous Neumann boundary conditions for µ and ϕ: We assume that the domain is large enough to render the effect of the boundary conditions imposed on Γ ext on the droplet dynamics negligible compared to those induced by the viscous and capillary forces. For the combination of boundary conditions in (8)-(9), the pressure variable p is only determined up to a constant. We therefore in addition impose the condition that p vanishes on average, i.e. p := Ω p dx = 0.
NSCH equations such as (3) intrinsically constitute mixture models, describing the evolution of a binary fluid composed of a mixture of two immiscible constituents. However, NSCH models are usually applied as diffuse-interface models for binary fluids, by regarding only phasesegregated configurations. Conforming to this perspective, we consider initial conditions that correspond to a liquid volume set in an ambient background fluid, separated by a phasetransition with a prescribed profile. Specifically, denoting the liquid volume by Ω l ⊂ Ω and the corresponding liquid-ambient interface by Γ la := ∂Ω l \ ∂Ω, the initial condition for the order parameter is given by: where d ± (x, Γ la ) represents the signed (Euclidian) distance from x to Γ la , positive for x ∈ Ω l . The function s → tanh(s/ √ 2ε) represents an equilibrium solution for the phase field of the Cahn-Hilliard equations in one spatial dimension (see, e.g. [17], [51, §8.4.3]) and, accordingly, the phase field (10) is meta-stable if the radius of curvature of ∂Ω l is sufficiently large relative to ε. The initial conditions for (3) are completed by a specification of the mixture velocity. We restrict ourselves in this work to homogeneous initial conditions, implying that the binary-fluid is initially stationary and that its dynamics are induced by the unbalance of the capillary and pressure forces corresponding to the initial configuration of the liquid volume encoded in (10).

Weak formulation of the NSCH system
In this section we present a weak formulation of the NSCH problem (3). This weak formulation also forms the basis for the considered adaptive finite-element approximation. We will not attempt to specify the function spaces for the weak formulation precisely, and instead proceed formally. We define the test space for the weak formulation of (3) as: The auxiliary condition p = 0 will be accounted for in the weak formulation by means of the Lagrange-multiplier formalism. Accordingly, the test space is extended by a constant function, which we identify with an element of R. The trial space is a linear space of U-valued functions on the considered time interval (0, t fin ), denoted by L ((0, t fin ) ; U). The trial space for the Lagrange multiplier is denoted by L ((0, t fin ); R). The aggregated NSCH problem can then be condensed into the following weak formulation: The semi-linear forms a 0 , a 1 , a 2 : U × U → R are given by and the linear form l : U → R contains the exogenous data. To disambiguate the tensor notation, we note that components (with summation on repeated indices). The semi-linear form a 0 contains the two time-derivative terms. The first three terms of a 1 , together with the second term in a 2 , correspond to the equation of motion in (3) 1 . The fourth and fifth terms of a 1 represent the convective Cahn-Hilliard equation (3) 3 . Let us note that ∇ · (ϕu) has been replaced by u · ∇ϕ, which is a consistent modification for solenoidal u. The remaining terms in a 1 define the chemical potential in (3) 4 . The first term in a 2 accounts for the solenoidality of the volume-weighted mixture velocity (3) 2 .

Error estimation and adaptive approximation
In most applications of the NSCH equations it is necessary to select the transverse length scale of the diffuse interface, ε, significantly smaller than other characteristic length scales in the problem under consideration, e.g. the radius of a liquid droplet or the size of the domain. To account for the resulting spatial multi-scale behavior of the NSCH equations, we consider adaptive approximations of (12). Specifically, we consider adaptive approximations based on Truncated Hierarchical B-splines [23] (THB-splines) analogous to [50]. In this work, we restrict ourselves however to conventional C 0 continuous approximation, so that the obtained results in Section 5 extend immediately to conventional finite-element approximations. The adaptiverefinement procedure follows the standard recursive SEMR (Solve → Estimate → Mark → Refine) process [10,21]. To enhance robustness and efficiency of the adaptive-approximation procedure, we apply a first order Backward Euler time-integration scheme on coarse meshes and a second-order Crank-Nicolson scheme on the finest meshes.

Approximation of the NSCH system
We regard a discretization of the NSCH system (12) with respect to the spatial dependence based on THB splines. The details of the THB-splines are omitted here, and the reader is referred to [50, §3.1] for an overview. We define the hierarchical approximation space as: where H k,α L is the piece-wise polynomial THB spline space of order k and regularity α. The subscript L = 0, 1, . . . , L max denotes the refinement level, and the spaces U L are hierarchical in the sense that U L ⊂ U L+1 . The operator T # is the pullback of the transformation T : Ω → (0, 1) d from the physical domain to the unit square. Note that with the indicated choice of k and α, we have C 0 continuous spline bases of degree 3 for all but the (adjusted) pressure field p, the latter's degree being 1 lower. Accordingly, the velocity-pressure pair is a member of the Taylor-Hood family.
To enhance the stability of the discrete approximation, the convective term in the Navier-Stokes equations in (13b) is replaced by a skew-symmetric formulation according to Ref. [35], by adding the following semi-linear form to the formulation: By means of the product rule and Gauss' theorem, one can establish that s according to (14) forms a partition of zero for solenoidal u. Augmenting the weak formulation with s causes the volumetric terms associated with convective transport, except for the term corresponding to ∇ρ, to cancel if v is replaced by u, irrespective of solenoidality of u: In pure-species regions (i.e. regions where ∇ρ = 0), this modification eliminates the artificial energy production that could otherwise occur on account of inexact solenoidality of the Taylor-Hood approximation of the velocity.
To discretize the weak form in time, we introduce a partition of the time interval (0, t fin ) into subintervals of length τ n = 2 −l τ max for some l := l n ∈ Z ≥0 and fixed τ max . The timestep-reduction factor l is determined a-posteriori, based on the convergence of the Newton procedure; see Section 4.1 below. To enhance the stability of the formulation, on coarse meshes we apply a first order Backward Euler (BE) scheme to discretize the equations in the time dependence, with a second order contractive-expansive splitting of the double-well potential with stabilization [54]: Here, s + ands are given by s + replaces the derivative of the double-well potential by a second order contractive-expansive splitting, while the artificial-diffusion terms acts as stabilization, where the coefficient β needs to satisfy conditions as stated in [54, Thm.1] for unconditional energy stability. In this work, we select β := 2m, which allows some leniency if |ϕ| > 1. We use Π to denote a generic projection operator. In particular, in (16), the projection ΠU n−1 Lmax (incl. Πϕ n−1 Lmax ) corresponds to an L 2 -projection of the super-refined approximation (see Section 3.2) at the maximum level of refinement in the previous time step to the level-L approximation space in time step n; see Remark 1 below for further elaboration. On fine meshes in the sequence of adaptive refinements, we apply a second order Crank-Nicolson (CN) scheme to discretize the equations in the temporal dependence: It is to be noted that in (17), the bilinear form a 2 according to (13c) is treated fully implicitly, as are the terms related to the constraint on the average pressure, α p , and the corresponding constraint on the test space, χ q . Such a fully implicit treatment of constraints and corresponding Lagrange multipliers is standard in the CN scheme; see, e.g., [33]. Note also that in (16) and (17), terms involving the approximation in the previous time step, i.e. with index n − 1, are evaluated at the super-refinement of the maximum adaptive-refinement level.
Remark 1. Both (16) and (17) contain terms which depend on the approximation obtained in the previous time step at the super-refinement of the maximum level of refinement, U n−1 Lmax . The assembly of such terms carries a relatively high computational cost, even if the operative approximation space, U n L , is still coarse, i.e. if L < L max . In the right member of the CN scheme (17), U n−1 Lmax appears in the bilinear form a 0 , and in the semi-linear forms a 1 and s. In a standard BE scheme (i.e. according to (16) without projections), U n−1 Lmax only appears in the bilinear form a 0 and in the stabilization terms s + ands. Because the semi-linear forms a 1 and, to a lesser degree, s are significantly more complex than a 0 ,s and s + , on coarse meshes the computational expense of the CN scheme is considerably higher than that of a BE scheme. From the perspective that for L < L max each approximation U n L merely serves to guide the adaptive algorithm to construct a suitably refined approximation space at level L + 1 (and to provide an initial estimate to the Newton iteration procedure to solve for U n L+1 ; see Section 4.1), we opt to use the BE scheme (16) at the coarsest refinement levels L = 0, 1, 2, . . . , L be . On the subsequent refinement levels L be + 1, . . . , L max , we switch to the CN approximation (17) to recover second-order temporal accuracy. Application of the BE scheme (with an expansive-contractive decomposition of the double-well potential with stabilization) on coarse meshes, carries the additional advantage that this scheme is generally more stable, which enhances the robustness and efficiency of the Newton solution procedure on the coarse meshes; see Section 4.1.
To further reduce the computational complexity of the BE scheme at the coarse levels, instead of using U n−1 Lmax directly, we use its L 2 -projection ΠU n−1 Lmax onto U n L . The projection needs to be evaluated only once at each level of refinement, and it generally provides a sufficiently accurate approximation for the a 0 and s + terms in (16) for the purpose of guiding the adaptive-refinement procedure. We retain ϕ n−1 Lmax in the right-hand side of (16) because this term involves the gradient of the phase field instead of the phase field in algebraic form, and it needs to be evaluated only once. Especially in the nonlinear term s + in (16), the use of the projection Πϕ n−1 Lmax instead of ϕ n−1 Lmax is beneficial, because the terms that emanate from s + in the nonlinear solution procedure (see Section 4.1) require frequent re-evaluation.
Remark 2. In accordance with [33], in practice we observe that the CN scheme provides significantly better accuracy than the BE scheme. The inaccuracy of the BE scheme is actually prohibitive, in the sense that the time step which is required by the BE scheme to adequately resolve the dynamics of the interface motion, is generally many orders of magnitude smaller than the characteristic time scale of the motion itself. The aforementioned switch between BE and CN within the adaptive refinement procedure however relies on the assumption that within a single time step, the BE approximation is sufficiently accurate to guide the refinement process on the coarse levels, i.e. the diffuse-interface features corresponding to the CN approximation on level L be + 1 should be adequately resolved by the super-refinement U n Lbe . Even if the overall accuracy of the BE scheme is too limited, it can generally serve this purpose.
The systems (16) and (17) are considered for consecutive time steps n = 1, 2, . . .. For n = 1, the approximation U 0 Lmax that appears in the right members of (16) and (17), is replaced by the initial conditions. With reference to Section 2.1, we note that the initial-boundaryvalue problem for the NSCH equations is generally furnished with initial conditions for the order parameter ϕ and for the mixture velocity u. Because the right-hand side of the BE scheme (16) only depends on a 0 and s + and, hence, only on ϕ and u, the prescribed initial conditions suffice to initialize the BE scheme. The CN scheme (17) depends on a 0 , a 1 and s and, therefore, on ϕ, u and µ. To start the CN scheme, the initial chemical potential must thus be specified in addition to the initial order parameter and mixture velocity. The chemical potential can be consistently initialized based on the initial datum for the order parameter, ϕ 0 , by invoking its definition (3) 4 . In practice, we accomplish this by an L 2 (Ω) projection onto a highly resolved hierarchical-refined approximation space,Ṽ , according to: It is to be noted that the gradient of ϕ 0 that appears in the right member of (18) cannot generally be computed exactly for initial data of the form (10). However, for the considered test cases, we have access to very accurate series expansions of ϕ 0 and, hence, its gradient.

Error estimation and adaptive spatial refinement
The diffuse interface introduces features into the solution with a length scale ε that is generally significantly smaller than other characteristic length scales in the problem under consideration. To efficiently resolve the various features of the solution, including those in the vicinity of the diffuse interface, within each time step in the time-integration process we construct a sequence of adaptively refined meshes, following the standard recursive SEMR process. The adaptive-refinement procedure is directed by a hierarchical a-posteriori error estimate. The error-estimation-and-refinement procedure is described in detail in [50], but its main elements are repeated here for completeness.
We employ a so-called two-level hierarchical error estimate wherein, for each field, the approximation in U L is compared to the approximation in its super-refinement U L to estimate the error. More precisely, for each time step n and refinement level L, we define the super-refinement U n L of the approximation space U n L by uniform subdivision of each knot span (element) of U n L into 2 d equal parts. We denote by U n L ∈ U n L the super-refined approximation obtained from (16) or (17) (depending on the refinement level L). Considering appropriate norms for each of the fields f ∈ F := {u, p, ϕ, µ}, we can then estimate the error in U n L as: The applied marking strategy is based on support-wise error indicators [34,50,52,39]. Specifically, denoting by b f L any basis function of the super-refined approximation space for field f , we define the error indicators: denotes the open support of the basis function b f L . By virtue of the overlap in the basis functions, the sum of the support-wise error indicators yields an upper bound to the error estimate (19). For each field f ∈ F, we then mark a minimal set of basis functions such that the sum of the corresponding error indicators exceeds a prescribed fraction a f ∈ [0, 1) of the total sum: Note that we may set a f = 0 for some f ∈ F, in which case the fields f in question are not considered in the adaptive-refinement process. In the refinement step, all basis functions b f L of which the support intersects with the support of the marked set will be replaced by their corresponding hierarchical refinements; see [50] for details. In practice, some marginal auxiliary refinements are introduced in a so-called completion step to retain a regular structure.

Solution methods
The nonlinear systems (16) and (17) translate into systems of nonlinear algebraic equations for the coefficients relative to the THB-spline basis functions. Devising a robust solution procedure for the nonlinear algebraic systems emerging throughout the adaptive-refinement procedure and time-integration process, is nontrivial. We consider a solution procedure for the nonlinear problems based on a Newton method with line search. To enhance the robustness of the Newton procedure, we introduce a continuation approach in which the diffuse-interface parameter ε and mobility m are adapted to the resolution of the mesh. This continuation procedure also serves to bypass the severe time-step restriction that could otherwise occur for small ε. Furthermore, a time-step-reduction procedure is implemented as a fail-safe to the Newton solution method. The latter also aids in retaining robustness of the solution procedure in the presence of the wide range of dynamical features that generally occurs in the evolution of the NSCH equations for small ε. Motivated by the severe ill-conditioning that can occur in the linear tangent problems in Newton's method for the coupled NSCH system for small ε (see Section 5.1), we regard both a monolithic (i.e. fully coupled) solution procedure and an iterative solution procedure based on a partitioning of the NSCH system into its composing NS and CH subsystems.

Newton solution procedure with ε-continuation
The nonlinear systems (16) and (17) give rise to a system of nonlinear algebraic equations of the generic form R n L (U n L ) = 0, where U n L contains the coefficients of U n L ∈ U n L relative to the applied THB-spline bases. Below, we will generally suppress the affixes n and L, unless these are relevant for the exposition. Let us note that the super-refinements in the a-posteriori errorestimation process are essentially of the same form, and the discussion below applies mutatis mutandis to these super-refinements. The basic Newton procedure with line search can be summarized as follows: given an initial estimate U 0 , repeat for k = 0, 1, 2, . . .: where A := A(U k ) denotes the derivative of the residual R. The step size (or relaxation factor) S is determined from minimizing a cubic interpolation of the residual norm between U k and U k+1 , based on the residual and the tangent matrix at these approximations; see [53] for details on the implementation. The iteration (21) terminates when the norm of the residual drops below a certain prescribed tolerance. The robustness and efficiency of the Newton procedure (21) depend sensitively on the error in the initial estimate, on the nonlinearity of the residual function, and on the properties of the tangent matrix, notably, its condition number in the vicinity of the solution. The Newton procedure is generally non-robust for the nonlinear problems that emerge for small values of the diffuse-interface parameter ε on coarse meshes. Specifically, the under-resolved thin diffuseinterface in the first few iterations of the adaptive algorithm can cause divergence of the Newton procedure, or may lead to prohibitively small steps in the line search.
An additional problem in the solution of NSCH systems, pertains to the fact that the solution from a previous time step fails to represent an accurate approximation to the solution in the time step under consideration if the interface motion within a time step is large in comparison with the transverse length scale, ε. This problem is in fact universal, independent of the adaptiverefinement process. The adaptive-refinement process however enables much smaller values of ε than would generally be feasible on uniform meshes and, hence, it prominently exposes the problem. To illustrate this issue, we consider the uniform translation of a 1D diffuse interface described by the usual tangent-hyperbolic profile, with length scale ε: where v > 0 denotes the translation velocity of the diffuse interface; see the illustration in Figure 2. Considering a time step τ , it holds that: where the approximation holds to reasonable accuracy in the interval 0 ≤ vτ /ε ≤ 1. Equation (23) conveys that in order for ϕ(0, x) to provide a suitable approximation to ϕ(τ, x), it is necessary that vτ is sufficiently small compared to ε. This implies that if the diffuse-interface thickness ε is reduced, the time step in the time-integration procedure has to be decreased proportionally. Otherwise, the initial approximation provided to the Newton procedure is too Figure 2: A 1D diffuse interface at times t = 0, τ given by Equation (22), uniformly translated over a distance vτ with length scales ε (a) and ε 0 (b). The absolute differences of the functions before and after translation are presented in (c) and (d), respectively.
far off from the actual solution, leading to divergence of the Newton process or excessively small steps in the line-search procedure.
To improve the robustness of the Newton procedure on coarse meshes and to avoid the severe time-step restriction that can emerge for small ε, we introduce a so-called ε-continuation process, viz. we enlarge the parameter ε for the first iterations of the adaptive refinement process. As opposed to [50], in conjunction with the interface-thickness parameter ε, we rescale the mobility parameter m, so that the continuation effectively pertains to (ε, m). This rescaling of the mobility is imperative to retain robustness if significant interface displacements occur within a time step. To implement the ε-continuation procedure, we select K ≤ L max and set ε 0 := 2 K ε, such that the diffuse interface with thickness corresponding to ε 0 , is adequately resolved by the baseline approximation space U n 0 . Correspondingly, we introduce a mobility parameter m 0 = C τ 2 3K m, with C τ a sufficiently large constant such that the diffuse interface is appropriately equilibrated within a time step; see Remark 4 below. Note that, accordingly, the constant C τ depends on the time step in the time-integration process. In iteration L = 0, 1, . . . , K − 1 of the refinement process, we solve problem (16) or (17) with interface-thickness parameter ε L = 2 −L ε 0 and mobility m L = 2 −3L m 0 , i.e. the interface thickness is decreased during the adaptive refinement and the mobility is scaled correspondingly. For L = K, . . . , L max , we solve (16) or (17) with the original parameters, ε L = ε and m L = m. For all refinement levels except the coarsest, the initial approximation for the Newton procedure is obtained from a projection of the super-refinement at the previous refinement level i.e., symbolically, (U n L ) 0 = ΠU n L−1 . The projection Π is based on a Stokes projector for velocity to ensure solenoidality of the initial approximation, and an L 2 (Ω) projection for the order parameter and the chemical potential. The initial approximation for the super-refinements at all levels is obtained via the same projection of the approximation itself, i.e. (U n L ) 0 = ΠU n L . The initial approximation for the Newton process at the coarsest level is provided by a projection of the coarsest-level solution of the previous time step, i.e. (U n 0 ) 0 = ΠU n−1 0 . By virtue of the ε-continuation process, the interface is adequately resolved on all meshes that occur in the sequence of adaptive refinements. Moreover, the error in the initial approximation for each Newton process, except at the coarsest level, is essentially independent of the motion of the interface, and it depends only on the continuation error, i.e. the error induced by the transition from length scale ε L−1 to ε L . At the coarsest level, the initial error does depend on the motion of the interface but, with reference to (23), the error is only proportional to vτ /ε 0 . The time step restriction can therefore be effectively diminished by selecting the diffuse-interface thickness on the coarsest mesh, ε 0 , sufficiently large.
Remark 3. In the transition from refinement K to K + 1, the mobility decreases by the factor C τ . Our experience is that the robustness of the solution procedure is insensitive to C τ . However, alternate, more gradual continuation schemes can be devised if necessary.

Remark 4.
On all levels of the adaptive-refinement procedure and, accordingly, for all values of ε L , the phase field at the previous time step that appears in the right-hand side of (16) and (17) corresponds to the original interface-thickness parameter, ε. This implies that in the initial stages of the ε-continuation process, where ε L > ε, the diffuse interface essentially has to expand from length scale ε to length scale ε L within a single time step, τ . To accommodate this expansion, the mobility parameter, m, has to be chosen sufficiently large. The characteristic diffusive time scale of the Cahn-Hilliard equations, T diff , is given by see Appendix A. To maintain the same interfacial expansion rate throughout the different stages of the ε-continuation process, we therefore introduce a scaling of the mobility according to m L ∝ ε 3 L . The diffuse interface is generally suitably equilibrated for τ ≥ T diff ; see Appendix B.

Direct solution method
Central to the Newton iteration procedure, is the solution of the linear tangent problems in (21) 1 . It appears that the properties of tangent problems corresponding to the NSCH equations have not been systematically investigated previously. An important objective of this work is therefore to elucidate the properties of these tangent problems, by considering the condition number of the tangent matrix A appearing in (21) 1 and of its corresponding submatrices; see Section 5.1. Foregoing a detailed investigation, one can however anticipate that inversion of the tangent matrix A is generally hard, based on the characteristic features of the NSCH system, viz., the intrinsic asymmetry of the system, the aggregation of disparate subsystems (NS and CH), the occurrence of small parameters (notably, ε and m), and the occurrence of heterogeneous coefficients that may differ by orders of magnitude (notably, the species densities ρ l and ρ a , and species viscosities η l and η a ).
Considering the complexity of solving (21) 1 , a robust and general solution procedure is required. Direct solvers are generally regarded as the most robust solution methods for problems of the form (21) 1 . The robustness of direct solvers is however not uniform, and depends on underlying strategies such as reordering, rescaling and pivoting. In this work, we apply the PARDISO sparse direct solver [43,42,12], which is among the most robust direct solvers available; see, e.g., [25]. Specifically, we use either the Intel oneMKL PARDISO solver [18], or the PARDISO 6.0 solver [7,13,11].

Partitioned solution method
Motivated by the complexity of solving the coupled problem (21) 1 directly, we also consider a partitioned solution strategy, in which the tangent matrix is subdivided into parts corresponding to the NS and CH subsystems. Let us note that the Lagrange multipliers which have been introduced to impose the auxiliary condition p = 0 (see Section 2.2) are included in the NS subsystem. We opt to introduce the partitioning of the NSCH system at the level of the linear tangent problems (21) 1 , instead of at the level of the (nonlinear) problem formulation. The partitioning strategy can hence be conceived of as a preconditioning approach for the linear tangent problems.
To specify the partitioned solution strategy, we note that in correspondence with the composition of the NSCH system (3) of the NS and CH subsystems, we can introduce a decomposition of the coefficients, U := (U ns , U ch ) and, similarly, of the residual, R := (R ns , R ch ). Each linear tangent problem (21) 1 in the Newton procedure can then be decomposed as: In the setting of (24), the partitioned solution strategy consists of omitting one of the offdiagonal blocks, A ns-ch or A ch-ns , in the actual solution step. The resulting system is then block upper or lower triangular, and inversion only requires inversion of the subsystem matrices, A ns-ns and A ch-ch . Instead of using the partitioned method directly as an iterative solution method, we apply it as a preconditioner in a GMRES method [40] for solving (21) 1 . In particular, instead of solving (21) 1 , we solve P A δ = −P R by means of GMRES, where P ∈ {P upp , P low } and if the lower triangular A ch-ns block is omitted and, correspondingly, the upper-triangular part is retained, and if the upper triangular A ns-ch block is omitted. In the considered partitioned solution procedure, the subsystem solves corresponding to the inverse matrices A −1 ns-ns and A −1 ch-ch are again performed by means of the PARDISO sparse direct solver. The principal advantage of the partitioned solution method relative to the direct solution method, is that the subsystem matrices A ns-ns and A ch-ch are much better conditioned and more amenable to rescaling than the aggregated tangent matrix A; see Section 5.1.

Remark 5.
The partitioned solution methods and preconditioners associated with (25) and (26) are fundamentally different. Essentially, in the upper-triangular preconditioner (25) the effect of the transport velocity u in the CH equations is treated explicitly and the effect of the order parameter ϕ and the chemical potential µ (via the auxiliary mass flux J in (4) and capillary stress tensor ζ in (6)) in the NS equations is treated implicitly, while for the lower-triangular preconditioner (26) these effects are treated vice versa. The terms ∇ · (u ⊗ J ) and ∇ · ζ in (3) 1 contain higher-order derivatives of ϕ and µ, while u appears only in algebraic (i.e. non-differentiated) form in (3) 3 and as a second-order derivative in (3) 1 via the divergence of the viscous-stress tensor, ∇ · τ . This suggests that P −1 upp is a much better approximation of the principal part of the NSCH system than P −1 low and, accordingly, that P upp is much more effective as a preconditioner than P low . This conjecture is confirmed by the numerical experiments in Section 5.2.
Let us also note that the A ns-ch block that is retained in the upper-triangular preconditioner (and discarded in the lower-triangular preconditioner) carries the fluid-fluid surface tension, via the capillary tensor ζ. In the upper-triangular preconditioner, surface tension is therefore treated implicitly in the solution procedure, while in the lower-triangular preconditioner it is treated explicitly.

Numerical experiments
To establish the main properties of the AGG NSCH model (3), the adaptive approximation method as outlined in Section 3, and the solution procedures presented in Section 4, we conduct numerical experiments for an oscillating liquid droplet suspended in ambient fluid. We restrict ourselves here to a two-dimensional setting. We consider an initially ellipsoidal droplet; cf. Figure 1. We exploit the symmetry of the configuration along the principal axes of the ellipse and accordingly consider only a quarter of the liquid-ambient domain. Specifically, the considered computational domain is a square with sides of 50 µm, i.e. Ω := (0, 50) 2 µm 2 , and prescribe symmetry conditions according to (8) on The complementary part of the domain boundary is referred to as an external boundary, and is denoted by Γ ext := ∂Ω \ Γ sym . On this part of the boundary, we impose boundary conditions according to (9). The initial configuration of the droplet corresponds to an ellipse with foci (±10 µm, 0 µm): The ellipsoidal initial configuration of the liquid droplet is represented by an initial condition conforming to (10) for the order parameter. For the mixture-velocity field, we impose the homogeneous initial condition (11), with Γ la corresponding to Ω l in (28). Let us note that the initial configuration does not correspond to an eigenmode of the (linearized) oscillating droplet: for infinitesimal radial perturbations, an ellipsoidal shape formally corresponds to the second mode of oscillation, but the corresponding velocity field does not vanish at any moment. An overview of the model-parameter values for the various considered test cases is presented in Table 1. Section 5.1 considers the condition number of the linear (sub)systems that occur in the tangent problems in Newton's method for the coupled NSCH system, and its dependence on relevant system parameters. In Section 5.2 we analyse the convergence behavior of the GMRES method with the P upp and P low preconditioners, corresponding to the two variants of the partitioned solution method discussed in Section 4.3. Finally, in Section 5.3, we present a simulation with a thin diffuse interface and high spatiotemporal resolution, and validate the AGG NSCH model by comparing the obtained results with those of a corresponding sharp -interface model.

Conditioning of the tangent matrices
To assess the complexity of solving the linear tangent problems in (21) by means of the direct solution method in Section 4.2 and the partitioned solution methods in Section 4.3, we investigate the dependence of the 1-norm condition number of the tangent matrix A and of the submatrices A ns-ns and A ch-ch on the parameters of the NSCH model and its discretization. The condition number is an important characteristic in assessing the complexity of linear systems but, evidently, it provides only partial information. To account for the fact that most modern direct solution methods apply an equilibration (reordering and rescaling) of the matrix before entering the solution procedure, we consider both the condition number of the original tangent matrices, and of corresponding equilibrated matrices. In this work, we apply the PARDISO sparse direct solver [43,42,12]. Because the equilibrated systems as generated and used by the PARDISO solver are not accessible externally, we instead consider the systems obtained from MATLAB's equilibrate function to assess the effect of equilibration. In view of the size of the linear systems, instead of computing the condition numbers directly, we utilize MATLAB's built-in condest function, which is based on the algorithms described in [28] and [29], to compute a repeatable 1 lower bound on the 1-norm condition number.
The analysis is performed on the aforementioned oscillating-droplet test case. We vary the interface thickness parameter ε and (minimal) mesh width h such that ε/h ≥ 1/2, thus ensuring adequate resolution in the vicinity of the diffuse interface. In addition, we limit ourselves to ε ≥ 2 −10 × 10 2 µm, to restrict the size of the linear systems to a feasible range. We regard configurations with parameter values as presented for test cases 1A-G in Table 1. For each test case, the entries represented by the * symbol indicate parameter(s) to be varied. We restrict our presentation to results for the Crank-Nicolson time discretization according to (17). However, the Backward Euler approximation (16) yields similar behavior (results not displayed). The presented condition-number estimates pertain to the tangent matrices in the last Newton iteration at refinement level L max in the first time step. We have verified that the presented condition-number estimates are representative, by probing several linear systems in later time steps and during the Newton iterations, for various cases (results not displayed). Table 2 presents an overview of the estimated lower bound on the 1-norm condition num-bers for test case 1A, exploring the (ε, h)-parameter dependence. In this test case, we consider uniform meshes comprising (5 × 5) 2 k elements (k = 0, 1, 2, 3). The number of DOFs corresponding to the finest mesh exceeds 65×10 3 . Preliminary tests for k = 4 conveyed that the condition-number estimates are unreliable for k > 3 (results not displayed). Let us note that the empty entries in Table 2 correspond to (ε, h) combinations for which the diffuse interface is not adequately resolved. These entries therefore have no significance. The results in Table 2 convey that the condition number of the aggregated NSCH system is generally significantly larger than that of the NS and CH subsystems separately, especially for small ε and h. Considering the h-dependence of the condition numbers of the original matrices, one can observe that the condition numbers of the CH and NS subsystem are essentially independent of h, while the condition number of the aggregated NSCH system deteriorates under h-refinement.
Considering the effect of equilibration, we note that equilibration is generally very effective, reducing the condition number by orders of magnitude. However, the results indicate that the effectiveness of equilibration decreases under h-refinement. For the CH and NS subsystems, the condition number of the equilibrated matrix increases, while the condition number of the original matrices is essentially independent of h. For the aggregated NSCH system, the ratio of the condition number of the original matrix to the condition number of the equilibrated matrix generally decreases under h-refinement. In general, the condition number of the equilibrated NSCH system deteriorates more strongly under h-refinement than that of the NS and CH subsystems. Considering the ε-dependence of the condition numbers, we observe that the condition number of the aggregated NSCH system deteriorates monotonously as ε decreases. This dependence is however weak and, approximately, cond(B) = O(ε 1/2 ) (as ε → +0). The condition number of the CH subsystem varies weakly and non-monotonously in relation to ε.
The results in Table 2 convey that the condition number of the CH subsystem varies at most by a factor of 2 across the full range of considered values of ε. The estimated condition number of the NS subsystem is essentially independent of ε, in accordance with the fact that the NS subsystem does not depend on ε in the matched-density scenario considered in test case 1A. For sufficiently small h, equilibration leads to monotonous decay of the condition number of the aggregated NSCH system as ε decreases. Also in this case, the dependence of the condition number on ε is weak. For the CH subsystem, equilibration leads to monotonous decay of the condition number as ε decreases. Evidently, the condition number of the equilibrated NS subsystem is essentially independent of ε, analogous to the condition number of the original NS subsystem.
In actual computations, the mesh width h must be reduced if the interface-thickness parameter ε is decreased, in order to ensure an appropriate resolution of the diffuse interface. The blueshaded diagonal entries of Table 2 correspond to combinations of (ε, h) for which ε/h = 2 −3 ×10 and the diffuse-interface is marginally resolved. It must be noted, however, that in a typical application of the NSCH system, one would select the mobility m proportional to ε 3 , such that the diffuse time scale is fixed; see test case 1G and Appendix A. These variations in m also affect the condition number; see e.g. Figure 3 below. Therefore, we restrict ourselves here to a basic characterization of the behavior of the equilibrated condition number under simultaneous ε and h refinement. From Table 2, one can infer that the condition number of the equilibrated NSCH tangent matrix scales approximately as ε −5/2 under simultaneous ε and h ∝ ε refinement, and the condition numbers of the equilibrated NS and CH subsystems scale approximately as ε −2 . Table 3 presents the estimated lower bound on the 1-norm condition numbers for test case 1B in Table 1. Test case 1B pertains to an adaptively refined approximation, in which  Table 1. Parameters of interest are the interface thickness ε and mesh width h. No local mesh refinement has been applied. The labels 'orig' and 'equi' indicate the condition numbers of the matrices before and after equilibration, respectively. The monolithic system is represented by 'full'. The Cahn-Hilliard and Navier-Stokes subsystems are indicated by 'CH' and 'NS', respectively. The blue shaded entries correspond to (ε, h) combinations for which ε/h = 2 −3 ×10, and the diffuse interface is marginally resolved. the interface-thickness parameter ε and the number of refinement levels L max are varied. The coarsest mesh comprises 10×10 square elements with mesh width h 0 = 5 µm. For completeness, we mention that we apply K = L max continuation steps, and set (here and throughout) the refinement factor to a f = 0.95 for f ∈ {ϕ, u} and a f = 0 for f ∈ {p, µ}; see Section 3. Table 3 reports the condition numbers corresponding to the super-refinement of the maximum refinement level. To relate the adaptive-approximation results in Table 3 to the uniform-mesh results in Table 2, we note that for L max = 1 the minimal element size is 2 −2 × 5 µm, which coincides with the element size in the right-most column of Table 2. Comparing the condition numbers in left-most column of Table 3 to those in the right-most column of Table 2 (for identical values of ε), we observe that the adaptive approximations generally exhibits moderately worse condition numbers for the original tangent matrices, and moderately better condition numbers for the equilibrated matrices. The condition numbers of the aggregated NSCH system, and of the NS and CH subsystems, generally deteriorate as the number of refinement levels L max increases. In addition, the effectiveness of equilibration diminishes as L max increases, in the sense that the ratio of the condition number of the original matrix to the condition number of the equilibrated matrix decreases. The condition numbers of the equilibrated NSCH system and the NS and CH subsystems increase approximately proportional to 2 7 2 Lmax , 2 Lmax , and 2 7 2 Lmax , respectively. As opposed to the case of uniform refinements, the estimated condition number for the NS subsystem exhibits a (weak) dependence on ε. This is caused by the fact that the adaptively-refined meshes implicitly depend on ε. Considering the ε-dependence of the equilibrated condition numbers, one can observe that the condition numbers of the equilibrated NSCH system and of the equilibrated CH subsystem display a weak monotonous decay as ε decreases, approximately proportional to ε 1/2 , while the condition number of the equilibrated NS system is essentially independent of ε.
Regarding the scenario in which ε is reduced and L max is simultaneously increased to ensure proper resolution of the diffuse interface, as indicated by the blue-shaded entries in Table 3, one  Table 1. Parameters of interest are the interface thickness ε, and the number of refinements L max . The labels 'orig' and 'equi' indicate the condition numbers of the matrices before and after equilibration, respectively. The monolithic system is represented by 'full'. The Cahn-Hilliard and Navier-Stokes subsystems are indicated by 'CH' and 'NS', respectively. The blue shaded entries correspond to combinations of (ε, L max ) for which the ratio ε/h of ε to the minimal element size of the super-refinement h = 2 −(Lmax+1) h 0 equals 2 −3 ×10, and the diffuse interface is marginally resolved. The yellow highlighted entry corresponds to settings that are repeated in all test cases 1B-F. can infer that the condition number of the equilibrated NSCH tangent matrix scales approximately as ε −3 , the condition number of the equilibrated NS subsystem scales approximately as ε −1 , and the condition number of the equilibrated CH subsystem scales approximately as ε −3 under simultaneous ε and h ∝ ε refinement, asymptotically as ε → +0. These proportionalities are consistent with the previously established behavior of the equilibrated condition numbers with respect to ε and L max separately.
In test case 1C we vary the mobility m to determine its effect on the conditioning of the tangent matrices. Other parameters for this test case are listed in Table 1. Figure 3 plots the condition numbers of the aggregated NSCH system and of the NS and CH subsystems, and of their equilibrated counterparts, versus the mobility parameter. For m ≤ 4 × 10 −13 m d s kg −1 , the diffusive time scale associated with the interface-equilibration process exceeds the period of oscillation of the droplet; see Section 5.3 and Appendix A. Such small values of m therefore have limited significance in this setting. On the other hand, for m ≥ 10 −7 m d s kg −1 , the diffusive interface-equilibration process corresponding to the Cahn-Hilliard equations is excessively dissipative, leading to strong artificial damping of the droplet oscillation. Such large values of m are therefore also irrelevant. From Figure 3 one can observe that the condition numbers of the original NSCH tangent matrix and of the CH subsystem are essentially independent of m. The condition number of the original NS subsystem is completely independent of m, as is the condition number of the equilibrated NS subsystem, in agreement with the fact that NS subsystem does not contain m. Equilibration causes the condition numbers of the NSCH system and of the CH subsystem to reduce as m decreases. In particular, Figure 3 indicates that for 4×10 −13 m d s kg −1 ≤ m ≤ 10 −7 m d s kg −1 , the condition numbers of the equilibrated NSCH system and CH subsystem are approximately proportional to m. This proportionality ceases for even smaller values of the mobility but, as explained previously, such small values of m have limited relevance.  Table 1. Parameter of interest is the mobility coefficient m. Matrix naming and computational application are as before in Tables 2 and 3. The vertical line corresponds to settings that are repeated in all test cases 1B-F.
Test case 1D is concerned with the dependence of the condition numbers on the time-step size τ ; see Table 1. The corresponding results are presented in Figure 4. One can observe that the condition number of the original aggregated NSCH tangent matrix is essentially independent of τ . The condition number of the original NS subsystem increases as the time step is reduced, approximately proportional to τ −2 . For the CH subsystem, the condition number increases approximately proportional to τ −1 under time-step refinement. Equilibration effectively improves the conditioning of the NSCH tangent matrix and of the NS and CH subsystems. The condition numbers of the equilibrated NSCH tangent matrix and of the equilibrated CH subsystem decrease approximately proportional to τ . The condition number of the equilibrated NS subsystem increases approximately as τ −1/2 as τ decreases.
Test case 1E pertains to the dependence of the condition numbers on variations in the density. The results of test case 1E are presented in Table 4. We consider liquid and ambient densities in the physically relevant range 10 0 kg m −d ≤ ρ l , ρ a ≤ 10 4 kg m −d . The results in Table 4 indicate that the condition numbers of the aggregated NSCH tangent matrix and of the NS and CH subsystems generally depend only weakly on ρ l and ρ a , both for the original matrices and for their equilibrated counterparts. Because the densities do not appear explicitly in the CH subsystem, one would expect that the condition number of the CH subsystem is entirely independent of ρ l and ρ a . Nonetheless, the condition numbers of both the original and the equilibrated CH subsystem vary by a factor of approximately 10 across the considered range of (ρ l , ρ a ). These variations are indirect, induced by e.g. differences in the adaptively-refined meshes.  Table 1. Parameter of interest is the time step size τ . Matrix naming and computational application are as before in Tables 2 and 3, and Figure 3. The vertical line corresponds to settings that are repeated in test cases 1B-F.   Test case 1F serves to assess the sensitivity of the condition numbers to variations in the liquid and ambient viscosities. We consider liquid and ambient viscosities ranging between 10 −5 kg m 2−d s −1 and 10 −1 kg m 2−d s −1 . The lower limit resembles the viscosity of air at atmospheric conditions, while the upper limit resembles the viscosity of low-viscosity oil. Other parameters are as listed in Table 1. We mention that the convergence of the Newton procedure with ε-continuation (see Section 4.1) is robust for the considered range of viscosities, but the robustness deteriorates for viscosity values that are significantly lower than the lower bound 10 −5 kg m 2−d s −1 . In addition, the magnitude of discretization errors in the form of spurious currents near the interface, increases as the viscosity diminishes. The condition number of the aggregated NSCH system does not display a significant dependence on the viscosities, neither for the original system nor for the equilibrated system. The condition number of the CH subsystem is independent of the viscosities, in accordance with the fact that η l and η a do not appear in the CH subsystem. The condition number of the NS subsystem depends weakly on the viscosities for η l and η a between 10 −5 kg m 2−d s −1 and 10 −2 kg m 2−d s −1 . However, for larger values of η l and η a , the condition number displays a sharp increase. Equilibration renders the condition number of the NS subsystem essentially independent of η l and η a . Table 5: Lower bounds of the 1-norm condition numbers as computed by MATLAB's condest function for test case #1F in Table 1. Parameters of interest are liquid and ambient viscosities η l and η a .  Finally, we consider in test case 1G the effect of simultaneous variation of ε, m and L max , corresponding to a scenario in which the interface-thickness parameter ε is decreased to approximate the sharp-interface limit more closely, the mobility parameter m is reduced proportional to ε 3 to fix the diffusive time scale T diff (see Appendix A), and the maximum refinement level L max is raised to ensure adequate resolution of the diffuse interface at length scale ε. In particular, we regard with T diff = 2.686 × 10 −8 s. Other parameters are as listed in Table 1. Figure 5 plots the estimated condition number of the tangent matrix of the aggregated NSCH system and of the NS and CH subsystems versus the interface thickness, both for the original systems and their equilibrated counter parts. The condition numbers in Figure 5 pertain to the super-refined approximations. The largest value of ε in Figure 5 correspond to L max = 0. In this case, the super-refinement consists of a uniform mesh with mesh width h = h 0 /2 = 5 2 µm. Figure 5 conveys that the condition number of the aggregated NSCH system increases approximately as ε −2 as ε is reduced and m and L max are adapted accordingly. The condition number of the CH subsystem scales similarly with decreasing ε. The condition number of the NS subsystem is approximately proportional to ε −3/2 . Equilibration is effective, and renders the condition numbers of the NSCH system and the CH subsystem essentially independent of ε. The condition number of the equilibrated NS subsystem increases approximately proportional to ε −1 as ε decreases, and m and L max are adapted simultaneously. It is noteworthy that the results for the equilibrated condition numbers in Figure 5 are consistent with the results in Table 3 and Figure 3, in the sense that according to Table 3 the condition number of the equilibrated NSCH, NS and CH systems scale as ε −3 , ε −1 and ε −3 , respectively, under simultaneous ε and L max refinement, while according to Figure 3 the condition numbers of the equilibrated NSCH, NS and CH systems scale as m 1 , m 0 and m 1 , respectively, as m decreases.  Table 1. The interface thickness ε := ε k , mobility coefficient m := m k and the maximum refinement level L max := L max,k are adapted simultaneously according to (29).

preconditioned-GMRES convergence
To analyze the convergence behavior of GMRES with the P upp and P low preconditioners introduced in Section 4.3, we regard the residual reduction versus the number of GMRES iterations for the oscillating-droplet test case with parameter values as listed in Table 1, test case 2. Specifically, the parameters for the liquid and ambient fluid are selected such that the test case pertains to the oscillation of a water droplet suspended in air. We consider a relatively thin diffuse interface, corresponding to an interface-thickness parameter ε = 2 −10 ×100 µm. This implies that ε is approximately 150× smaller than the radius of the circular droplet with the same (2-dimensional) volume as the ellipse in (28). The maximum refinement level, L max = 7, is chosen such that the ratio of the interface thickness parameter, ε, to the minimal mesh width, h, is ε/h = 2 −2 ×10, and the interface is adequately resolved.
The convergence behavior of the preconditioned GMRES method displays some minor variability during the Newton iterations, during the refinement iterations, and in time. To illustrate the main aspects of the convergence behavior, Figure 6 plots the residual reduction r i / r 0 , versus the number of GMRES iterations, i, for the P upp and P low preconditioners, for the super-refinement of the maximum refinement level in the first time step and in time step 25, and for refinement level L = 5 in the first time step. For completeness, we mention that the results in Figure 6 pertain to the final step in the Newton procedure, but similar behavior is observed in other Newton iterations. One can observe that the upper-triangular preconditioner yields much more effective convergence behavior than the lower-triangular preconditioner (see also Remark 5), on account of the fact that P upp -preconditioned GMRES exhibits a significant decrease in the residual in the first iteration. In subsequent iterations, the residual reduction per iteration is very similar for P upp -preconditioned GMRES and P low -preconditioned GMRES, and both are essentially constant within each tangent problem, resulting in almost straight curves in Figure 6 for i > 1. Figure 6 conveys that, excepting the first iteration, the residual reduction per GMRES iteration lies between 0.5 and 0.6.  Table 1.
Remark 6. As a stopping criterion for the preconditioned GMRES procedure, we require that the residual reduction exceeds 10 8 . Such a stringent condition on the approximate solution of the linear tangent problems in (21) 1 is necessary to ensure proper (quadratic) convergence of the Newton process (21). If the stopping criterion for the GMRES procedure is too loose, the Newton method generally displays suboptimal convergence, or convergence can stagnate.
Remark 7. In the first time step of test case 2, time-step reduction (by a factor of 4) is invoked, to ensure adequate convergence of the Newton procedure; see Section 4. We conjecture that the poor convergence of the Newton method in the first time step with the original time step according to Table 1 is caused by the fact that the tanh profile of the initial data for the order parameter in (10) is not equilibrated for the elliptic initial configuration (28), which can lead to strong transient behavior. The results reported in Figure 6 for the first time step in fact pertain to the first time step which is completed at the original time step τ max in Tabel 1.

Comparison to sharp-interface model
A fundamental notion underlying diffuse-interface models for binary-fluid flows, is that in the sharp-interface limit ε → + 0 the solution approaches the solution of a corresponding sharpinterface model (provided that the limit exists), viz. a model comprising the incompressible Navier-Stokes equations on two complementary domains, separated by a manifold of codimension 1 which carries surface energy; see [1,3]. In actual computations of binary-fluid flows based on diffuse-interface models, the sharp-interface limit is inaccessible, i.e. one cannot pass to the limit ε → + 0, and one generally assumes that the diffuse-interface model forms an accurate representation of the sharp-interface model if the interface-thickness parameter ε is sufficiently small relative to other relevant length scales in the problem under consideration, e.g. the minimal radius of curvature of the fluid-fluid meniscus. It is noteworthy that analyses of sharp-interface limits of diffuse-interface models via matched-asymptotic expansions (see e.g. [3]) do not generally provide insight into the deviation between the two models, as the deviation pertains to higher-order terms.
In this section we compare results of the AGG NSCH diffuse-interface model with a relatively thin diffuse interface to corresponding sharp-interface results for an oscillating-droplet test case. The characteristic parameters are listed in Table 1 under test case 2. Note that this setup is the same as used for the investigation of the convergence properties of preconditioned GMRES in Section 5.2. For the ellipsoidal initial configuration of the droplet according to (28), the minimal radius of curvature occurs on the major axis of the ellipse and is 5 µm. The interfacethickness parameter ε is therefore approximately 70× smaller than this minimal radius of curvature, suggesting that the sharp-interface limit is adequately approximated. Let us note that the disparity in length scales is enabled by the adaptive refinement procedure described in Section 3. Rayleigh's classical theory [48] on the oscillation of inviscid droplets conveys that the period of oscillation of the considered droplet is approximately: with n = 2 for the considered ellipsoidal mode of oscillation and R = √ 200 µm the radius of the circular droplet with the same (2-dimensional) volume as the ellipse in (28). Equation (30) implies T osc ≈ 16 µs. The period of oscillation, T osc , is approximately 600× larger than the diffusive time scale T diff = ε 3 /σ la m.
We compare the results of the AGG NSCH model with those of a sharp-interface model based on the methodology described in [27]. In the sharp-interface model, the liquid and ambient fluids are represented by two systems of incompressible Navier-Stokes equations on complementary domains, Ω l and Ω a = Ω \ Ω l , respectively. The two NS systems are coupled at their mutual interface by dynamic and kinematic boundary conditions. The motion of the liquid and ambient domains is accommodated by means of an Arbitrary-Lagrangian-Eulerian formulation. The NS equations in ALE form are approximated by means of a finite-element method, using standard Taylor-Hood P2-P1 elements for velocity and pressure. The sharpinterface model employs a time step τ si = 2 −8 ×10 µs.
Intrinsically, diffuse-interface models do not provide an explicit representation of the interface. A comparison of the diffuse-interface results to the sharp-interface results in terms of the location of the interface is therefore non-trivial. We present a comparison in terms of geometric quantities that are well defined for both the diffuse-interface model and the sharp-interface model, viz. the second area moments m 20 and m 02 with respect to the axes of symmetry: m 02,nsch (t) = Note that these values differ by a factor 4 from those of a full droplet, because only a quarter of the droplet is considered. Figure 7 plots the second area moments of the NSCH and SI models versus time. The results in Figure 7 illustrate that the motion of the liquid-ambient interface, as expressed by the second area moments, is essentially indistinguishable for the diffuse-interface and sharp-interface models. The difference in oscillation period between both models, measured after two oscillations, is approximately 2 −8 × 10 µs, which is smaller than the time-step sizes used in both the NSCH and the sharp-interface simulations. Both the diffuse-interface and the sharp-interface binary-fluid model satisfy an energydissipation inequality. For the sharp-interface model, the energy-dissipation inequality is: ∇u : τ a dx + bnd (33) where 'bnd' denotes terms that are only supported on the external boundary Γ ext and E kin,si (t) = represent kinetic energy and interface energy, respectively. The viscous-stress tensors τ l and τ a in (33) are identical in form to (5) with the viscosity parameter according to η l and η a , respectively. One can infer that the integrals appearing on the right-hand side of (33) are non-negative, from which it follows that in the sharp-interface model, the energy decays by viscous dissipation, up to boundary terms. The energy-dissipation inequality corresponding to the NSCH model is with E kin,nsch (t) = where ρ is the mixture density according to (2). In the diffuse-interface model, the interface energy is represented by a Ginzburg-Landau-type functional, corresponding to the second functional in (36). Equation (35) conveys that in the diffuse-interface model, the energy decays by viscous dissipation and by gradients in the chemical potential, up to boundary terms.
To compare the energy-dissipation properties of the AGG NSCH model and the sharpinterface model, Figures 8 and 9 plot the evolution of the kinetic energy and the interface energy pertaining to both models, respectively. In addition, Figure 10 displays the evolution of the total energy for both models. From Figures 7 and 8, one can infer that the kinetic energy attains minima at the extrema of the second area moments, i.e. when the deformation of the droplet attains its (local in time) maximum. Figure 9 conveys that the interface energy is maximal at the minima of the kinetic energy, which illustrates the exchange of kinetic and interface energy during the oscillation of the droplet, up to dissipation. The minimal interface energy is approximately 1.617 µJ, which coincides with (one quarter of) the surface energy of a spherical droplet with the same volume as the ellipse in (28). From Figs. 8 and 10, one can observe that the energy dissipation is largest near the maxima of the kinetic energy, coinciding with the spherical configuration of the droplet. At these instants, the velocities in the liquid droplet are largest, inducing the largest viscous dissipation. Overall, the total energy as well as both its components exhibit very similar behavior in the NSCH model and in the sharp-interface model, except that the NSCH model exhibits slightly more dissipation.
From Figure 9 it can also be observed that the interface energy of the NSCH model deviates more pronouncedly from that of the sharp-interface model near the first minimum at approximately 4.4 µs, and a corresponding deviation occurs in the total energy in Figure 10. This deviation is caused by an initial transient, which is temporally under-resolved in the approximation of the NSCH model at the considered time step.  energy-dissipation relation (35). Despite the use of double-well splitting, neither the backward Euler scheme (16) nor the Crank-Nicolson scheme (17) do however preserve this property uniformly, i.e. independent of the time step. The results in Figure 10 convey that the discrete approximation provided by the Crank-Nicolson scheme with time step τ max according to Table 1 complies with the energy-dissipation relation. It is to be mentioned, however, that the energydissipation relation of the Crank-Nicolson scheme is conditional on the time step, and it is violated at a twice larger time step (results not displayed). Finally, we present a qualitative comparison of the the droplet configuration, velocity magnitude |u| and pressure field p. Figure 11 plots the computed |u| and p fields at t ≈ 33.8 µs, near the instant at which the second area moments reach their extrema after two full oscillations of the droplet; see Figure 7. In the NSCH results, the representation of the droplet geometry is indicated by a plot of the ϕ = 0 contour line. The droplet configurations pertaining to the NSCH and to the SI model are almost indistinguishable. The pressure fields of the SI and NSCH models are also virtually identical. The velocity-magnitude fields in the SI and NSCH models match closely: the local extrema are in nearly identical locations, both in the droplet and in the ambience. Minor deviations between the NSCH and SI results can be observed, which are accentuated by the contour lines. It is noteworthy that at the considered instant, the velocity in the droplet almost vanishes, while the ambient velocity is relatively large. Because the ambient density is 10 3 times smaller than the liquid density, however, the kinetic energy nearly vanishes; cf. Figure 8. Figure 11: Snapshot of the velocity magnitude (top) and pressure (bottom) profiles for the SI (left) and NSCH (right) models, after two oscillations at t = 3.38 · 10 −5 s. The colorbars are identical for both models; the contour lines for the velocity magnitude are divided into 9 intervals over the colorbar range. Γ la in (a) and a ϕ = 0 contour line in (b) have been added to highlight the interface locations.

Conclusion
In this work, we developed and analyzed an adaptive-approximation method for binary-fluidflow problems, based on the Abels-Garcke-Grün (AGG) Navier-Stokes-Cahn-Hilliard (NSCH) diffuse-interface model. To circumvent instabilities in the formulation for large density and viscosity contrasts due to range-violations of the phase field, we applied an Arrhenius interpolation for the mixture viscosity and a soft-clipped linear interpolation for the mixture density.
To resolve the spatial scale disparity between the diffuse-interface thickness and other representative length scales in the problem, we introduced an adaptive-refinement procedure directed by a two-level hierarchical a-posteriori error estimate. We regarded approximations based on truncated hierarchical B-splines, restricted to standard C 0 -continuous spline bases. The adaptive-approximation framework constructs a sequence of hierarchically refined approximation spaces within each time step by means of the Solve → Estimate → Mark → Refine (SEMR) procedure.
We introduced several enhancements of the basic SEMR process to improve the robustness and efficiency of the adaptive-refinement method and the corresponding solution procedure: • On coarse levels of the adaptive-refinement process, we apply a Backward Euler (BE) scheme with a second order splitting of the double-well potential. On the finest levels, we use a Crank-Nicolson (CN) scheme to recover second-order-in-time accuracy of the approximation. On the coarse levels, the computational complexity of the BE scheme is significantly lower than that of the CN scheme. In addition, the BE scheme is generally more stable on coarse meshes. The temporal accuracy of the BE scheme is too limited to apply it directly as a time integrator for diffuse-interface binary-fluid-flow problems with significant interface dynamics. However, the BE scheme is adequate to guide the initial refinement steps, and to provide an initial estimate for the Newton procedure on the next hierarchical level of refinement.
• We introduced an ε-continuation approach, in which both the interface-thickness parameter, ε, and the mobility, m := m ε , are adapted during the adaptive-refinement process, in such a manner that the diffuse-interface is adequately resolved at each refinement level, for the interface-thickness parameter at that level. The mobility is adapted in accordance with the interface thickness, such that the diffuse interface is suitably equilibrated within a time step. We established that to this purpose, the mobility needs to scale with the (level-dependent) diffuse-interface thickness as m ε ∝ ε 3 . The ε-continuation approach serves two purposes. First, it enhances the robustness of the solution procedure on coarse levels of approximation, where the diffuse interface corresponding to the original interfacethickness parameter would be severely under-resolved, leading to non-robustness of the Newton procedure. Second, it mitigates the conditions on the maximum time step in relation to robustness of the Newton procedure.
• We proposed a partitioned solution method for the linear tangent problems in Newton's method for the coupled NSCH system, to circumvent ill-conditioning of the monolithic system. The partitioned solution method corresponds to preconditioned GMRES, wherein the preconditioner is the original tangent matrix with either the CH-NS or the NS-CH submatrix suppressed. Suppression of either of these submatrices effectively decouples the NSCH system into its NS and CH constituents, which can then be solved separately. We conjectured that the preconditioner in which the NS-CH submatrix is retained (P upp ) is more effective than the preconditioner without the NS-CH submatrix (P low ), because the former provides an implicit treatment of the surface-tension terms. This conjecture was confirmed by the numerical experiments.
To determine the main properties of the AGG NSCH model, and of the enhanced adaptive approximation method, we conducted numerical experiments for an oscillating droplet in an ambient fluid. We empirically investigated the dependence of the condition number of the linear tangent matrices in Newton's method on the main problem parameters, both for the monolithic NSCH system and for the NS and CH submatrices. To account for the fact that contemporary direct solution methods generally apply equilibration, we also considered the effect of equilibration on the condition numbers of these (sub-)matrices. We found that equilibration is in general very effective, reducing the conditions numbers in many cases by orders of magnitude. Moreover, equilibration in many cases moderates or eliminates parameter dependence of the condition numbers. The effectiveness of equilibration however typically deteriorates as the system size increases, e.g. under mesh refinement. In the relevant range of the system-parameter values, we found that the condition numbers of the equilibrated NSCH system and of the equilibrated NS and CH subsystems depend on the minimal element size in the adaptively refined mesh, h min , the interface-thickness parameter, ε, the mobility, m, the time-step size, τ , the liquid and ambient densities, ρ l and ρ a , and the liquid and ambient viscosities, η l and η a , as: For simultaneous variations of the interface thickness ε, the minimal element size h min ∝ ε, and the mobility m ∝ ε 3 , corresponding to a practical scenario in which the interface-thickness parameter ε is decreased to approximate the sharp-interface limit more closely, the mobility parameter m is reduced proportional to ε 3 to fix the diffusive time scale, and the minimal element size is decreased proportional to ε to ensure adequate resolution of the diffuse interface, the condition numbers of the equilibrated NSCH, NS and CH (sub-)matrices are proportional to ε 0 , ε and ε 0 , respectively. The proposed partitioned solution method is generally very effective in solving the linear tangent problems in Newton's method. GMRES with the P upp preconditioner is more efficient than with the P low preconditioner. The P low preconditioner provides an essentially uniform convergence rate. The P upp preconditioner exhibits a similar convergence rate in all iterations except the first: in the first iteration, it yields a residual reduction of several orders of magnitude.
We presented a comparison of the numerical results of the AGG NSCH model with a thin diffuse interface with those of a corresponding sharp-interface model for the droplet-oscillation test case. The diffuse-interface results and the sharp-interface results are in excellent agreement, with regard to both the droplet-oscillation dynamics and the energy-dissipation characteristics. At the considered time step, the NSCH model with the CN time discretization was observed to be thermodynamically consistent in the sense that the pertinent free-energy-dissipation relation is satisfied, but this thermodynamic consistency is conditional on the time step and is violated at larger time steps.
The proposed enhancements of the SEMR method and the corresponding solution strategy enable the solution of the AGG NSCH diffuse-interface model for binary-fluid flow problems for physically relevant parameter values, i.e. with realistic values of the surface tension and of the liquid and ambient densities and viscosities. By virtue of the adaptive-refinement strategy, it is feasible to consider diffuse-interface thicknesses that are several orders of magnitude smaller than other characteristic length scales in the problem under consideration, thus enabling the investigation of sharp-interface limits. The proposed ε-continuation strategy is congruent with the adaptive-refinement procedure, because the continuation iterations can be merged with the adaptive-refinement iterations in the SEMR process. The fact that ε-continuation also effectively increases the admissible time step in the time-integration procedure, suggests that the continuation approach could also be beneficial in non-adaptive methods with implicit timeintegration schemes.
In the presented numerical simulations, the ε-continuation strategy and the partitioned solution method improve the robustness of the solution procedure to the extent that the accuracy and stability (in the sense of thermodynamic consistency) of the considered time-integration scheme becomes a limiting factor. This suggests that significant gains in efficiency can be achieved by means of high-order implicit time-integration schemes for NSCH systems, in conjunction with the proposed ε-continuation procedure. Nonetheless, we anticipate that scenarios exhibiting temporal multi-scale behavior, in particular, due to topological changes in the fluidfluid interface caused by break-up or coalescence, necessitate adaptive refinement of the time step.
with ε * = 2 −3 . Let us note that the initial dataφ 0 ,ū 0 corresponds for an equilibrium solution of the NSCH system (3) on Ω = R for ε * = 2 −3 , but not for the considered interface thickness ε = 1 in (B.1). The considered setting therefore essentially pertains to the expansion of the diffuse interface from length scale ε * = 2 −3 to length scale = 1. We assume that the effect of the domain truncation in combination with the boundary conditions (9) on the interface equilibration process is negligible and, accordingly, that the phase fieldφ approaches the equilibrium solutionφ ∞ (x) := − tanh(x/ √ 2) as t → ∞. To examine the interface-equilibration process, we consider a finite-element approximation on a uniform mesh with mesh width h = 10 −2 , without adaptive spatial refinement and εcontinuation. The time step size is restricted to τ max = 10 −1 , with an initial time-step-reduction factor m = 12, such that τ 0 = 2 −12 τ max . The small initial time-step size is selected to resolve the initial fast dynamics of the interface with sufficient accuracy.
The computed phase fieldφ(t, x) is presented in Figure B.12 for t = 10 {−1,0,1} . For reference, the equilibrium profileφ ∞ (x) := − tanh(x/ √ 2) is displayed as well. One can observe thatφ(t, x) indeed displays under-and overshoots, exceeding the bounds ±1. During the evolution, the velocity fieldū(t, x) and the corresponding kinetic energy E kin remain, up to numerical errors, equal to 0. This implies that the interface-expansion evolution is essentially independent of the fluid properties and, hence, that (φ,μ) corresponds to a solution of the generic CH equations (A.1). The results in Figure B.12 moreover convey that the diffuse interface is reasonably well equilibrated for t ≥ 1, and indistinguishable from the equilibrium profile for t ≥ 10.