Generalized-α methods for the Navier–Stokes–Cahn–Hilliard equations

Diffuse-interface models offer a powerful modeling framework for intricate binary-fluid flows in science and engineering, in view of their ability to accurately capture complex interfacial dynamics. A key advantage of diffuse-interface models lies in their inherent capability to handle topological changes, such as merging and breakup of interfaces (Khanwale et al., 2022), without the need for ad hoc treatments. In addition, they provide a natural framework for simulating contact line dynamics (Xu et al., (2018; Demont et al., 2026) and can robustly accommodate geometrically complex features (Stoter et al., 2023).

Diffuse-interface models for binary fluids find their origin in the seminal work of van der Waals (1893, 1979), who was the first to describe an interface carrying surface energy as a continuous-transition layer, and of Cahn and Hilliard (1958), who provided the first thermodynamically consistent model for the composition evolution of non-uniform mixtures. The first systematic extension to binary fluids, referred to as model-H, was presented by Hohenberg and Halperin (1977) in the late 1970s. Model-H is an incompressible binary-fluid model, i.e. the mixture velocity is solenoidal. However, it has been specifically derived for binary fluids with uniform mass density, and it is only thermodynamically consistent in this restricted setting. Various extensions of model-H have been presented to account for non-matching component densities, e.g. the quasi-incompressible models by Lowengrub and Truskinovsky (1998) and Shokrpour Roudbari et al. (2018), and the incompressible models by Ding et al. (2007) and Abels et al. (2012). The aforementioned models are collectively referred to as Navier–Stokes–Cahn–Hilliard (NSCH or CHNS) models. An overview and classification of the main characteristics of existing models can be found in Eikelder et al. (2023). In addition, various reformulations (e.g. using different state variables) of these models have been presented, often with the aim of deriving tailored numerical approximation methods that retain certain fundamental properties of the binary-fluid system, such as species conservation and energy dissipation.

Numerical simulation of diffuse-interface models carries several fundamental challenges (Demont et al., 2022; Han and Wang, 2015): (i) spatial multiscale behavior arises due to the extremely thin fluid–fluid interface compared to macroscopic length scales, demanding high resolution; (ii) temporal multiscale behavior stems from the disparity between the fast diffusion dynamics at the interface, slower macroscopic processes and time singularities in topological changes; (iii) $ε$-conditional stability imposes severe time-step restrictions tied to a CFL-type condition where the interface displacement per time step is limited to its thickness (⁠$ε$⁠), leading to prohibitive computational costs as the sharp-interface limit is approached, i.e. as $ε→+0$⁠; (iv) heterogeneous parameters such as large density and viscosity contrasts can cause instability and ill-conditioning; (v) incompatibility and spurious velocities corresponding to discretization errors that manifest as parasitic flows, with a velocity magnitude that is inversely proportional to the Ohnesorge number and which are, therefore, particularly problematic in low-viscosity regimes; (vi) Ill-conditioning of linear systems complicates the solution of the coupled nonlinear equations due to asymmetry, parameter disparity and subsystem aggregation. Together, these issues present significant barriers to practical and predictive simulations of multiphase systems based on diffuse-interface models, particularly in engineering and scientific applications involving significant parameter variations.

Various methodologies have been proposed to address the aforementioned challenges. Adaptive mesh refinement has proven to be an effective approach to resolving the spatial multiscale behavior of diffuse-interface models (Demont et al., 2022; Garcke et al., 2016; Van Brummelen et al., 2021). To address non-robustness emerging from the temporal multiscale behavior, specific energy-stable time-integration methods have been proposed (Shen and Yang, 2010, 2013; Shokrpour Roudbari et al., 2018; Khanwale et al., 2022; El Haddad and Tierra, 2022; Liu et al., 2023; Garcke et al., 2016). By construction, the energy stability is unconditional in some schemes (Han and Wang, 2015; Garcke et al., 2016; Shokrpour Roudbari et al., 2018; Shen and Yang, 2010, 2013) and conditional in others (El Haddad and Tierra, 2022; Khanwale et al., 2022; Liu et al., 2023). Robustness of the solution process, in the sense of unique solvability of the algebraic system within each time step, is ensured for some of these approaches by monotonicity properties (Han and Wang, 2015) or linearity (Shokrpour Roudbari et al., 2018; Liu et al., 2023; El Haddad and Tierra, 2022; Garcke et al, 2016). In Demont et al. (2022), it was noted that the standard Newton solution procedure fails if the diffuse-interface displacement within one time step exceeds the interface thickness, owing to the fact that the error in the initial estimate provided by the solution in the previous time step then becomes excessively large. This non-robustness can be effectively resolved by means of $ε$-continuation, i.e. a continuation procedure in the diffuse-interface thickness, and this $ε$-continuation approach is synergetic with adaptive-refinement procedures. The issue of ill-conditioning in NSCH systems can generally be bypassed using a partitioned solution procedure that decomposes the NSCH system into its NS and CH parts (Demont et al., 2022; El Haddad and Tierra, 2022; Garcke et al., 2016; Liu et al., 2023; Shen and Yang, 2013).

By virtue of the preceding advances in computational methodologies, the robustness of time-integration methods for diffuse-interface binary-fluid models has improved significantly. However, the accuracy of energy-stable time-integration schemes is generally limited to first (El Haddad and Tierra, 2022; Shen and Yang, 2010, 2013; Shokrpour Roudbari et al., 2018; Liu et al., 2023) or second (Khanwale et al., 2022; Han and Wang, 2015) order and the artificial numerical dissipation in these schemes in many cases surpasses physical dissipation, unless time steps are applied that are exceedingly small compared to relevant problem-specific time scales, e.g. the period of oscillation of a droplet (Demont et al., 2022). In Demont et al. (2022) the standard second-order Crank-Nicolson scheme, acclaimed for its favorable low numerical dissipation, was applied to the problem of an oscillating droplet, and it was found that without $ε$-continuation the robustness of the Newton procedure limits the admissible time step, while with $ε$-continuation the accuracy of the Crank–Nicolson scheme is restrictive. This reveals the need for high-order time-integration methods for diffuse-interface binary-fluid models.

In this work, we consider a higher-order generalized-$α$ time-integration method for the NSCH equations, with the objective of assessing the properties of such an approach for diffuse-interface binary-fluid flow models. The class of generalized-$α$ methods offers higher-order alternatives to conventional two-step schemes such as the $θ$-method, which can in principle be extended up to arbitrary order (Behnoudfar et al., 2023, 2021). Generalized-$α$ methods hold several advantageous properties that render them promising in the context of the NSCH equations. First, their two-step, single-stage form, as opposed to multi-stage methods such as Runge–Kutta schemes or multi-step methods such as BDFs or Adams–Moulton schemes, facilitates temporally varying spatial adaptivity. Second, in generalized-$α$ schemes, the higher-order accuracy is accomplished by means of correction terms based on auxiliary variables corresponding to higher-order derivatives. These correction terms provide an intrinsic a-posteriori error estimate, which can be leveraged in adaptive time-step refinement procedures. Third, the correction stages are linearly implicit and unidirectionally coupled, i.e. each higher-order correction depends on the previous, but not vice-versa. As a result, generalized-$α$ schemes generally admit a very efficient implementation. Finally, generalized-$α$ methods have tunable numerical dissipation properties.

The remainder of this article is arranged as follows. In Section 2 we introduce the NSCH problem and review relevant numerical solution methods to motivate consideration of the generalized-$α$ scheme. Section 3 presents the generalized-$α$ method and the resulting weak forms. In Section 4, we present numerical experiments for a prototypical one-dimensional diffuse-interface flow problem pertaining to a uniform translation of the interface, to establish the properties of the third-order generalized-$α$ scheme. In addition, we presented corresponding results for two $θ$-methods to enable a comparison. Section 5 presents concluding remarks.

This section presents the NSCH model as derived by Abels, Garcke and Grün (AGG) (Abels et al., 2012) to describe the dynamics of a mixture of two immiscible incompressible Newtonian fluids separated by a thin-but-finite transition layer. Next, we briefly review the $ε$-continuation method used in combination with spatial adaptivity to resolve the spatial multiscale behavior of the NSCH equations and bypass the $ε$-dependent time-step restriction in the time-integration process for the NSCH equations.

In diffuse-interface models, the interface between two immiscible fluids is represented by a finite-thickness transition layer composed of a mixture of both constituents. The fluid species are regarded as being omnipresent and the local composition of the binary fluid is described by an order parameter (or phase-field variable), $φ∈[−1,1]$⁠, where $φ=1$ denotes pure species 1 and $φ=−1$ denotes pure species 2. The AGG NSCH model describes the evolution of the multiphase-flow system in terms of the volume-averaged velocity field, $u$⁠, the pressure, p, the order parameter, $φ$ and the chemical potential, $μ$⁠.

On an open time interval $(0,tfin)⊆R>0$ and a spatial domain $Ω⊆Rd$ (⁠$d=2,3$⁠), the AGG NSCH equations are given by the following:

The constitutive relations for the relative mass flux $J$⁠, specific to the AGG NSCH formulation, the viscous stress $τ$⁠, the capillary stress $ζ$ and the chemical potential $μ$ are defined as follows:

with $∇s$ the symmetric gradient for vector fields, i.e. $∇su=12(∇u+(∇u)T)$⁠. This NSCH system generally admits non-matching densities and viscosities and, accordingly, $ρ=ρ(φ)$ and $η=η(φ)$⁠. However, to facilitate the presentation in this precursory investigation of the properties of generalized-$α$ methods for the NSCH system, we restrict ourselves to matched-density and matched-viscosity scenarios, alluding to the effect of generalizations if appropriate. In the matched-density case, the relative mass flux $J$ vanishes and the AGG NSCH equations reduce to the standard model-H (Hohenberg and Halperin, 1977).

The model parameters in equations (1) and (2) involve the interface thickness parameter, $ε>0$ and the mobility parameter, $m>0$⁠. Respectively, these control the diffuse-interface length-scale and the diffusive time scale. The material parameters in equations (1) and (2) are the fluid-fluid surface tension, $σ12=232<private-char specific-use="loader-generated">&#xef07;</private-char>σ$⁠, mass density, $ρ=ρ1=ρ2>0$ and dynamic viscosity, $η=η1=η2>0$⁠.

The NSCH system can generally be equipped with various types of boundary conditions, e.g. no-slip conditions (Jacqmin, 2000), (generalized-α) Navier slip conditions (Xu et al., 2018; Demont et al., 2026), dynamic and static contact angle conditions (Jacqmin, 2000; Xu et al., 2018; Demont et al., 2026), etc. For transparency, in the present work, we restrict ourselves to settings in which both the velocity and phase field are subject to stationary Dirichlet conditions on some parts of the boundary, $ΓDu$ and $ΓDφ$⁠, respectively, homogeneous traction conditions hold on the complementary part $ΓNu=∂Ω∖ΓDu$⁠, the phase field satisfies homogeneous Neumann conditions on $ΓNφ=∂Ω∖ΓDφ$ and $μ$ satisfies homogeneous Neumann conditions throughout. Denoting by $g={uD,φD}$ the collected Dirichlet data for $u$ and $φ$⁠, we define the trial space for the weak formulation of (1) subject to the aforementioned boundary conditions as:

and the test space $V0$ as its counterpart with homogeneous traces of $u$ and $φ$ on $ΓDu$ and $ΓDφ$⁠, respectively. The weak formulation of equation (1) subject to the aforementioned boundary conditions can then be condensed into:

where the bilinear form $[·,·]$ corresponds to the usual $L2(Ω)$ inner-product of its arguments, i.e. $[u,v]=∫Ω(uv)$ for scalar-valued arguments and $[u,v]=∫Ω(u·v)$ for vector-valued arguments, and:

Remark 1. For matched densities, the volume term in the trilinear form$T$inequation (5a)is skew symmetric in$u$and$w$⁠. Owing to its skew-symmetric structure, this formulation of the convection operator mitigates spurious energy generation (or dissipation) when the transport field$w$is not strictly solenoidal, for instance, when the velocity and pressure are approximated using Taylor–Hood elements; seeLayton (2008).

Numerical approximation of equation (4) by a standard Rothe-type approach generally leads to a sequence of nonlinear algebraic equations associated with each time step. Numerical solution of each of these nonlinear algebraic systems by means of a Newton method requires an initial estimate. Conventionally, this initial estimate is obtained from the previous time step. However, this standard procedure leads to $ε$-conditional stability: if the diffuse-interface displacement within one time step is large relative to the interface thickness, then the error in the initial estimate is too large, causing breakdown of the Newton procedure; see the illustration in Figure 1. The error in the phase field obtained from a previous time step is in fact approximately $uτ/2ε$ for $0≤uτ/ε≤1$⁠, with u as the transversal velocity of the interface and $τ$ as the time step. In the limit $uτ/ε→+0$⁠, the error is $O(uτ/ε)$⁠; see Figure 1. Hence, to control the error in the initial estimate in the Newton procedure, the time step must be reduced if the interface thickness is reduced, to ensure that $uτ/ε$ remains sufficiently small.

In Demont et al. (2022), a procedure has been presented to bypass the $ε$-dependent time step restriction imposed by the robustness of the Newton procedure. This $ε$-continuation method applies a continuation procedure in the interface-thickness parameter, in which the initial estimate for the largest interface thickness parameter, $ε0$⁠, is obtained from the previous time step, but the initial estimate for each of the subsequent reduced interface thicknesses, $εl=εl−1/2$ (⁠$l=1,2,…L$ with $εL=ε$⁠) corresponds to the solution in the current time step at $εl−1$⁠, instead. Accordingly, the time step restriction imposed by the Newton procedure takes the form $τ ≲ ε0/u$ instead of $τ ≲ ε/u$⁠. The $ε$-continuation procedure is synergetic with adaptive spatial refinement based on an a-posteriori error estimate, in that the continuation steps can be combined with the refinement steps in the adaptive process. Such adaptive spatial refinement is instrumental in numerical methods for the NSCH equations, in view of the large gradients that occur in the vicinity of the diffuse interface.

Equipped with the $ε$-continuation method, the most stringent time-step restriction derives from the temporal discretization error, related to the temporal-multiscale behavior of the diffuse-interface problem. For instance, for the second-order Crank–Nicolson scheme, the discretization error visibly impacts the results or the robustness of the solution procedure, well before $uτ/ε=1$⁠; see Demont et al. (2022) and also Section 4.2.

Higher-order (⁠$k>2$⁠) time-integration methods can potentially avoid the stringent time-step restrictions of second (or first) order methods. In principle, higher-order methods offer higher convergence rates in the limit $τ→+0$⁠. Conversely, for suitably small time steps, this implies that higher-order methods also provide a significant reduction in the discretization error relative to a lower-order scheme. However, when adaptive refinement is used in space to resolve the diffuse interface, the feasibility of various higher-order time-integration schemes can be significantly affected. This is because the diffuse interface must be resolved at every time step and/or stage of the time-integration process, requiring fine spatial resolution over the entire integration sequence. This causes spatial refinement in the length scale $ε$ to resolve the diffuse interface, in multiple areas near the stages or steps, leading to excessive computational complexity.

To illustrate the effect of local spatial refinement near the diffuse interface on the computational complexity of higher-order time-integration schemes, Figure 2 illustrates the necessary common refinements for a multi-stage method, a multi-step method and a two-step method such as the generalized-$α$ scheme. In multi-stage methods such as Runge–Kutta schemes, the diffuse interface must be resolved in each stage. However, because the approximation space in all stages within one time step must be identical, a common refinement is required that unites the approximation spaces (meshes) for all stages; see the left panel in Figure 2. In addition, a projection step is required at the end of each time step, to avoid the auxiliary refinements for the intermediate stages carrying over to the next time step and, recursively, to all subsequent time steps. For multi-step methods such as BDF, the approximation at time step n involves the computed solutions in time steps $n−i$⁠, $i=0,1,…,s$ with s the number of steps; see Figure 2 (center). This implies that a common refinement is required that unites the approximation spaces for time step n and all s previous time steps, or the bi-linear and semi-linear forms in (4) must be evaluated for test functions at time level n and trial functions at all time levels $n−i$ (⁠$i=0,…,s$⁠). Both approaches carry a significant computational cost. In addition, auxiliary projection steps may be required to avoid adaptive refinements from carrying over to subsequent time steps. Two-step methods only require a common refinement that unites the approximation spaces for two time levels, namely, n and $n−1$⁠. Therefore, two-step methods have the potential of being computationally more efficient in comparison to multi-stage and multi-step methods for diffuse-interface models of binary-fluid flows.

Generalized-$α$ methods represent a class of higher-order (⁠$≥2$⁠) time-integration methods, alternative to (at most) second-order accurate $θ$-methods, designed to provide stable and accurate solutions for time-dependent differential equations. In this work, we focus specifically on the third-order generalized-$α$ scheme. This section considers the application of the third-order generalized-$α$ scheme to the NSCH equations and elaborates on the different steps in the solution procedure.

To provide a setting for a Rothe-type discretization of equation (4), with a generalized-$α$ discretization in the temporal dependence, we consider a partition ${tn}$ of the time interval $[0,tfin]$⁠. For convenience, we assume the partition to be uniform and denote by $τ=tn−tn−1$ the time step. To discretize (4) in the spatial dependence, in turn, we regard a conforming approximation space $Vgh,n⊂Vg$ for time step $tn$⁠, subordinate to a mesh-resolution parameter $h>0$⁠.

In generalized-$α$ methods, the state variables of the problem under consideration, $U=(u,p,φ,μ)$⁠, are treated as the primary field, while auxiliary variables are introduced associated with higher-order time derivatives of U. These auxiliary variables are defined by update equations based on Taylor-series expansion and state equations obtained by differentiation of the original weak formulation of the problem under consideration with respect to the temporal dependence. The conforming approximation of the primary variables for the NSCH equations in time step $tn$ then corresponds to $Un=(un,pn,φn,μn)∈Vgh,n$⁠. For the third-order generalized-$α$ scheme, the auxiliary variables correspond to the first, second and third-order time derivatives of U. Their conforming FEM approximation at time-step $tn$ is denoted by $U˙n∈V0h,n$⁠, $U¨n∈V0h,n$ and $U⃛n∈V0h,n$⁠.

In each time step, $n=1,2,…$⁠, the primary variable $Un$ and auxiliary variables $U˙n$⁠, $U¨n$ and $U⃛n$ in the generalized-$α$ scheme are obtained in four sequential steps, which will be elaborated on below.

Step 1

In the first step, the time derivative $U˙n$ is extracted from the weak formulation of the NSCH problem (4) with $∂tu$ and $∂tφ$ replaced by suitable approximations:

where we have introduced the condensed notation $Ψn=Ψ(φn)$ and $Ψn′=Ψ′(φn)$⁠. The occurrences of $Un=(un,pn,φn,μn)∈Vgh,n$ in (6) are implicitly defined in terms of $U˙n$ and $Un−1,U˙n−1,U¨n−1$ and $U⃛n−1$ via the first update equation:

Hence, $U˙n$ is the only unknown in (6). The parameters $α1∈[1,2]$ and $γ1∈[12,32]$ in (6) and (7), respectively, are certain user-defined parameters; see Remark 6.

Remark 2. The Dirichlet data for$U˙n$is homogeneous under the standing assumption that the Dirichlet data$g$for U is stationary. If the Dirichlet data is time-dependent, then the trial space for$U˙n$is subject to non-homogeneous Dirichlet data, which are determined by the prescribed Dirichlet data for$Un$ via the update equation (7).

Remark 3. The update equation (7) requires data$Un−1,U˙n−1,U¨n−1,U⃛n−1$from the previous time step. This implies that for$n=1$, consistent initial data must be constructed for these items. Noting that the NSCH system (1) contains first-order time derivatives for$u$and$φ$, the specification of the initial-boundary-value problem requires initial data$u0=u(t=0)$and$φ0=φ(t=0)$⁠. Given$φ0$, the initial chemical potential,$μ0$, can be determined from (1d). In turn,$φ˙0$can be derived fromequation (1c). Initial data for pressure and acceleration,$p0$and$u.0$, can be extracted from (1a) subject toequation (1b)differentiated with respect to time, i.e. ∇· $u.0$ = 0. Consistent initial data for higher-order time derivatives can be obtained similarly, based on conditions obtained by differentiating the NSCH equations with respect to time.

Step 2

In the second step, $Un$ is extracted from the update equation (7), based on data from the previous time step, with index $n−1$ and $U˙n$ from step 1. If $Vh,n⊇Vh,n−1$ is a refinement of the approximation space in the previous time step, then the addition in (7) can be executed directly. Otherwise, the identity in equation (7) must be interpreted as a projection. Specifically, in the case of an $L2$-projection, $Un$ follows from:

One may note that the approximation space $Vh,n$ is not constrained by the Dirichlet data $g$ for $Un$⁠. The essential boundary conditions on $Un$ are implicitly accounted for in the construction of $U˙n$⁠; see Remark 2.

Step 3

The third step produces $U⃛n$⁠, based on the weak form equation (4), differentiated in time twice. The differentiation in time leads to a new weak formulation for $U⃛n$⁠, which is however significantly more convoluted due to the nonlinearity of equation (4):

where $Ψ̈αf$ and $Ψ̈αf′$ are abbreviations according to:

The occurrences of $U¨n=(u¨n,p¨n,φ¨n,μ¨n)$ in (9) are implicitly defined in terms of $U⃛n$⁠, $U¨n−1$ and $U⃛n−1$ via the second update equation:

so that $U⃛n$ is the only unknown in (9). The parameters $α2∈[12,2]$⁠, $αf∈[12,1]$ and $γ2∈[12,32]$ are user-defined; see Remark 6. Similarly to step 1, if the boundary data $g$ is time dependent, then the Dirichlet data for $U⃛n$ is constructed via the update equation (11) and the second-order time derivative of the data, $g¨(tn)$⁠; cf. Remark 2.

Remark 4. For non-matching mass density and viscosity of the fluid components, the weak formulation (9) for$U⃛n$becomes significantly more complicated, on account of the associated additional nonlinearities. Firstly, the additional term$∇·(u⊗J)$related to the relative mass flux$J$according to (2a) must be incorporated in the weak formulation (4), which in turn results in various additional terms inequation (9). Second, the mass density assumes a dependence on the phase field. To account for the fact that the range of the phase field is not generally restricted to$[−1,1]$, the dependence$φ↦ρ(φ)$is commonly defined by a so-called soft-clipped linear interpolation (Bonart et al., 2019):

with$δ=ρ2/(ρ1−ρ2)$⁠. The mass density appears in the weak formulation as a product$ρ(φ)u$with$u$, which in turn generates a significant number of additional terms inequation (9), requiring higher-order derivatives of the piece-wise expression (12). Third, for non-matching species viscosities, the mixture viscosity$φ↦η(φ)$is commonly modeled by the Arrhenius mixture-viscosity relation (Arrhenius, 1887;Van Brummelen et al., 2021;Demont et al, 2022):

where$Λ=ρ1M2ρ2M1$with$M1$and$M2$as the molar masses. The nonlinear dependence between$φ,u$and the viscous stress tensor,$(φ,u)↦η(φ)∇su$, generates several additional terms in (9). The multitudinous additional terms that occur in (9) as a result of the higher-order derivatives of the additional nonlinear relations associated with the extension to non-matching densities and viscosities, render this extension exceedingly complicated or even unfeasible. It is to be remarked that this complexity is further compounded for higher-order generalized-$α$schemes.

Remark 5. In this manuscript we restrict ourselves to scenarios in which the velocity$u$is sufficiently small to obviate convection stabilization, and we apply the weak formulations in equation (6)and, correspondingly,equation (9)directly. It is to be noted, however, that for larger$u$, the weak formulation (6) must be endowed with additional stabilization terms, which in turn lead to a proliferation of terms inequation (9).

Step 4

The fourth step is analogous to the second step, comprising the construction of $U¨n$ via the second update equation (11). This can be accomplished by the $L2$-projection:

Analogously to step 2, the approximation spaces do not explicitly incorporate Dirichlet boundary conditions, and these conditions are implicitly accounted for via the construction of boundary conditions on $U⃛n$ in equation (9) in step three via the second update equation.

Remark 6. The parameters$α1,<private-char specific-use="loader-generated">&#xef07;</private-char>α2,<private-char specific-use="loader-generated">&#xef07;</private-char>αf$and$γ1,<private-char specific-use="loader-generated">&#xef07;</private-char>γ2$in the third-order generalized-$α$method are implicitly defined via two user-defined parameters, namely,$ρ1∞∈[0,1]$and$ρ2∞∈[0,1]$, which determine the numerical dissipation and thereby the high-frequency damping (Behnoudfar et al., 2020). The$α$parameters:

determine the implicit versus explicit character of the scheme. The$α$parameters in turn determine the two$γ$parameters, which are used in the update equations in solution steps two and four, and which ensure the third-order accuracy of the scheme:

This choice of parameters guarantees an optimal convergence rate and unconditional stability for linear systems (Behnoudfar et al., 2020). For linear problems, the parameters$ρ1∞$and$ρ2∞$govern the spectral radius of the amplification matrix in the high-frequency limit, controlling the numerical dissipation. For nonlinear systems such as the NSCH equations, they act as tunable filters for the high-frequency components of the approximate solution that typically arise from temporal and spatial discretization errors. Smaller values of$ρi∞$introduce stronger damping, which helps stabilize stiff or convection-dominated flows, while values closer to one preserve more of the transient behavior with minimal artificial dissipation. Hence, the choice of$ρi∞$provides a flexible means to balance stability, accuracy and temporal smoothness in unsteady flow computations. Jansen et al. (2000) and Liu et al. (2021)proclaim that a spectral radius of$12$in high-frequency regions generally provides a suitable balance for incompressible Navier–Stokes simulations. In our method, this is obtained by setting$ρ1∞=ρ2∞=12$⁠. However, for accuracy-sensitive cases,$ρ1∞=ρ2∞=1$can preserve transient details, at the expense of reduced high-frequency damping. It is to be noted, however, that for sufficiently small time steps, the algorithm’s spectral radius does not bifurcate and the order of accuracy is independent of the selected parameters. To illustrate this aspect,Table A3in  Appendix A reports discretization errors for various settings of the stabilization parameters, for the test case presented in Section 4. These results convey that for the considered spatial and temporal resolution, the discretization error is insensitive to$ρ1∞$and$ρ2∞$⁠.

Remark 7. The correction terms in the generalized-$α$scheme can be leveraged to obtain an a-posteriori error estimate. Specifically, the expression in the right member of the update equation (8) can be conceived of as a Taylor series expansion of$U(tn)$around$tn−1$⁠. The Taylor series expansion has remainder of$O(τ5)$as$τ→+0$for$γ1=γ1*:=14.$However, for admissible values of$ρi∞$, it holds that$γ1∈[12,32]$; see Remark 6. Denoting by$|||·|||$the norm associated with the approximation space$V0h,n$, the expression:

therefore provides an estimate of the local truncation error. More precisely, under the assumption that the data at$tn−1$is exact, i.e. $(U,U˙,Ü,U⃛)n−1=(U,U˙,Ü,U⃛)(tn−1)$and U is sufficiently smooth on the interval$[tn−1,tn]$, it holds that:

while$|||U(tn)−Un|||=O(τ4)$⁠. The error indicator$ιn$according to (17) is therefore asymptotically exact, i.e. the effectivity index$Inloc:=ιn/|||U(tn)−Un|||=1$in the limit$τ→+0$⁠. It is to be noted that this asymptotic exactness only holds for the local truncation error, as accumulated errors in the data at$tn−1$will propagate into the error estimate$ιn$⁠. We refer to  Appendix B for further elaboration and for an analysis of the performance of the error indicator (17) for the numerical experiment in Section 4.

In this section, we assess the properties of the third-order generalized-$α$ scheme for a prototypical model problem, and we compare its performance to two $θ$-methods, namely, the implicit-Euler and Crank–Nicolson methods. The model problem pertains to a one-dimensional periodic setting of the NSCH system, with a solution corresponding to a translating (double, to enforce periodicity) interface. For this model, we can construct an analytical reference solution to verify the accuracy of the various time-integration methods. In addition, as the reference solution is characterized by a constant velocity, it is well-suited to assess the impact of the translation velocity of the diffuse interface on the robustness and accuracy of the time-integration schemes.

The considered numerical experiment concerns the convective transport of an equilibrium solution of the NSCH system on a one-dimensional periodic domain $Ω=(−ℓ/2,ℓ/2)$⁠. Such equilibrium solutions are hyperbolic tangent functions. On a periodic domain, we require a pair of hyperbolic tangent profiles to enforce periodicity in the translating solution. We define initial conditions for $φ(t=0)$ and $u(t=0)$ in accordance with the reference solution:

and $uref(t,x)=u0$ (uniformly constant), where $t˜$ is the remainder time, used to enforce periodicity of the reference solution:

The repeating pairs of hyperbolic tangent functions in equation (19) mimic the influence of neighboring profiles of virtually adjacent interfaces. They mitigate continuity issues that would arise due to the mismatch between the hyperbolic tangent profile within the finite domain size and its asymptotic values of $−1$ and 1.

In addition to the canonical initial conditions for $φ$ and $u$⁠, the third order generalized-$α$ scheme requires initial data for p and $μ$⁠, to complete $U(t=0)$ and for $U˙(t=0),U¨(t=0)$ and $U⃛(t=0)$⁠. The initial data for the time derivatives of $φ$ are obtained via direct time differentiation of the reference solution of the phase field in equation (19). The reference solution of the velocity field is uniform in $(t,x)$⁠. Accordingly, the initial conditions for the time derivatives of $u$ are homogeneous. Moreover, for $u$ uniform in $(t,x)$⁠, it follows from equation (1a) that p is uniform in $(t,x)$ (set to zero) and, hence, the initial conditions for p and its time derivatives are homogeneous. Finally, because $φref$ corresponds to an equilibrium solution, $μ$ vanishes uniformly in $(t,x)$⁠. Accordingly, homogeneous initial conditions are prescribed for $μ$ and its time derivatives. It is to be noted that the derivation of the consistent auxiliary initial data that are required by the third-order generalized-$α$ scheme is straightforward for the considered reference test case, but not generally; see Remark 3.

The physical parameters and time-stepping parameters used in the numerical experiment are shown in Table 1. They resemble the properties of water at room temperature with a length-scale corresponding to a picoliter droplet. Regarding spatial discretization, we use $C1$-continuous cubic splines for the velocity field, phase field and chemical potential, and $C1$-continuous quadratic splines for the pressure field. To ensure that the spatial discretization error is at least two orders of magnitude smaller than the temporal discretization error, we use a highly refined mesh consisting of 12,000–48,000 elements. The numerical implementation is performed in Nutils (Van Zwieten et al., 2022).

To assess the convergence behavior of the third-order generalized-$α$ method in comparison to the Crank-Nicolson and implicit-Euler methods, Figure 3 displays the $L2$-norm of the discretization error of the three schemes, $||φn−φ(tn,·)||L2(Ω)$⁠, at $tn=48ε/u0$⁠, versus the normalized time-step size, $u0τ/ε$⁠. One can observe that for sufficiently small time steps, the convergence rate of all three methods is in accordance with theory, i.e. $O(τ)$⁠, $O(τ2)$ and $O(τ3)$ for the implicit-Euler, Crank-Nicolson and generalized-$α$ method, respectively. For both $θ$-methods, the asymptotic convergence behavior occurs for time steps $τ≲ε/u0$⁠, which corresponds to the time scale at which the diffuse interface moves by approximately its thickness per time step. For the generalized-$α$ scheme, however, the asymptotic optimal convergence rate only manifests for time-step sizes $u0τ/ε≲10−1$⁠. At larger time steps, the error increases rapidly as $τ$ increases, until it roughly matches that of the Crank-Nicolson and implicit-Euler schemes at $u0τ/ε≃1$⁠. Furthermore, the Crank–Nicolson and implicit-Euler methods appear more robust than the generalized-$α$ method for large $τ$⁠: for $u0τ/ε≥2$⁠, Newton’s method fails to converge for the generalized-$α$ scheme, while it converges effectively for the Crank–Nicolson and implicit-Euler schemes. This nonrobustness of the (linearly) unconditionally stable generalized-$α$ scheme is further elaborated in Remark 8 below.

Remark 8. High-order generalized-$α$methods are unconditionally stable for linear problems. This guarantee does however not carry over to nonlinear settings. In practice, we observe that for nonlinear problems and large time steps, Newton’s method (with safeguards) fails to converge. We conjecture that this is primarily due to the breakdown of the high-order Taylor expansions used in the formulation. The auxiliary variables in the generalized-$α$method, corresponding to higher-order time derivatives of the state variables, depend nonlinearly on the state variables and can grow rapidly or lead to severe ill-conditioning of the linear tangent problems in Newton’s method. For instance, the time derivatives of the reference solution in (19) scale as$||∂t(m)φref(t,·)||L∞(Ω)∝ε−m$as$ε→+0$⁠. Consequently, the Taylor series can diverge, and the residuals in the Newton solver can grow or oscillate, preventing convergence.

To elucidate the discretization error of the generalized-$α$ scheme in comparison to that of the $θ$-methods, Figure 4 plots the error in the approximation of the phase field at $tn=48ε/u0$ for time-step sizes $u0τ/ε=1$ and $u0τ/ε=1/4$⁠. The time-step size $τ=ε/u0$ in Figure 4 (left) corresponds to the largest step size for which the Newton procedure converges for the generalized-$α$ method, and where the three methods yield similar accuracy; cf. Figure 3. It is notable that the implicit-Euler method shows a significant error throughout the entire domain, while for the Crank–Nicolson method, the error remains confined to a narrow region around the phase transition. For the generalized-$α$ method, the error is largest in the vicinity of the diffuse interface, but the error also extends into the single-phase regions. At the smaller time-step size $u0τ/ε=14$ in Figure 4 (right), the benefit of higher-order convergence in the generalized-$α$ scheme becomes apparent: its error is substantially lower than that of both $θ$-methods.

The computational cost of the generalized-$α$ scheme generally exceeds that of the two $θ$-methods on account of the auxiliary problems that must be solved. However, as these auxiliary problems are linear, the associated cost is limited. To provide insight into the differences in computational cost of the three schemes, we monitor the wall-clock time and the number of Newton iterations spent per time step. For the generalized-$α$ scheme, we moreover record the time spent in solving the auxiliary linear systems. We acknowledge that wall-clock time is only an indirect indicator of computational cost, sensitive to numerous hidden factors. However, the code base for the considered schemes as well as the underlying numerical methods such as integration procedures and linear solvers, are essentially identical. Wall-clock time can therefore serve as a relative indicator for computational cost, enabling a comparison between the different schemes. The results for a single time-step of size $u0τ/ε=14$ are listed in Table 2. We note that the variability of this wall-clock time over the time steps is negligible. Solving the nonlinear NSCH system using a Newton method takes essentially an equal amount of time for the generalized-$α$ method and both $θ$-methods. The implicit Euler scheme is less expensive per Newton iteration, by virtue of its simple structure, but in this case requires one more Newton iteration per time step than the Crank-Nicolson and generalized-$α$ schemes. For the $θ$-methods, the nonlinear NSCH system is the only system solved at each time step. Consequently, the total runtime of these schemes corresponds directly to the cost of solving the NSCH system itself. In contrast, the generalized-$α$ method requires three additional linear solves. The two update equations associated with equations (8) and (14), incur only minor computational cost by virtue of their simple structure. The time-differentiated auxiliary problem (9), however, requires substantially more time, reflecting its significantly more complex form. Overall, the computational cost of the generalized-$α$ method is approximately twice that of the $θ$-methods.

Based on the results of this experiment, we conclude that by virtue of its higher-order asymptotic convergence rate, the generalized-$α$ method outperforms the Crank–Nicolson method for sufficiently small time-step sizes, in particular, for the considered test case, for $u0τ/ε≤1/4$⁠. However, due to a steep deterioration of the accuracy of the generalized-$α$ scheme in the pre-asymptotic regime, the scheme fails to achieve comparable performance to the Crank–Nicolson method in terms of balancing accuracy and computational expense for $u0τ/ε≥12$⁠. For such larger time-step sizes, both methods yield comparable errors, but the generalized-$α$ scheme is more computationally demanding on account of its auxiliary variables. This suggests that the generalized-$α$ scheme is only appropriate for applications that require high accuracy. For such applications, its tunable numerical dissipation and built-in error estimation further enhance its potential.

The computational approximation of diffuse-interface binary-fluid models presents several fundamental challenges. In this work, we focus on the challenge of efficient and accurate time integration in an effort to combat its inherent temporal multiscale behavior. In earlier work, we introduced the concept of $ε$-continuation as a remedy for potential non-robustness of iterative solvers. When this approach is adopted, lower-order time-integration methods transition from being robustness-bound to accuracy-bound in terms of the admissible time-step size. Accordingly, in this work we investigated the feasibility of higher-order time-stepping methods to achieve accuracy when using moderate time-step sizes.

We explored the use of the higher-order generalized-$α$ method. In the context of diffuse-interface models, generalized-$α$ methods are promising candidates as they are two-step methods, only requiring data from the previous time step to compute the solution in the new time step. This is a crucial property, since diffuse-interface models are commonly combined with spatial adaptivity. A two-step time-integration method requires storage and cross-communication of only two adaptively refined meshes – unlike multi-step or multi-stage methods, where each additional step or stage introduces another refined mesh. Additional properties that render higher-order generalized-$α$ schemes a compelling option are their tunable parameters, which provide control over high-frequency numerical dissipation, and the inherent temporal-error estimator, which arises from manipulating the higher-order time derivatives solved for in the scheme.

From an implementation perspective, however, application of higher-order generalized-$α$ methods to the NSCH equations proves challenging. Since generalized-$α$ methods involve auxiliary (linear) problems based on repeated time differentiation of the governing PDEs with respect to the temporal dependence, the potentially large number of nonlinearities in NSCH models – exacerbated when material properties depend on the order parameter or additional stabilization terms are required – leads to a proliferation of terms. Moreover, the auxiliary state variables resulting from the time differentiation require consistent initial conditions, the derivation of which is unwieldy for nontrivial initial data, in particular for higher-order generalized-$α$ schemes. Similarly, for transient data for Dirichlet boundary conditions, the derivation of consistent data for the auxiliary variables is tedious.

To assess the performance of the generalized-$α$ scheme for the NSCH equations, we compared the third-order generalized-$α$ method with the implicit-Euler and Crank–Nicolson methods for a prototypical one-dimensional NSCH problem corresponding to uniform translation of the diffuse interface. A convergence study confirms their expected asymptotic convergence rates as the time step vanishes. The implicit-Euler and Crank-Nicolson methods (first- and second-order, respectively) each require one solve of a nonlinear system of equations. The third-order generalized-$α$ method additionally requires the solution of a linear system of equations, thus resulting in a higher computational cost. As expected, the Crank–Nicolson and generalized-$α$ methods consistently outperform the implicit-Euler scheme for sufficiently small time steps. By virtue of its higher-order asymptotic convergence rate, the generalized-$α$ scheme delivers superior accuracy compared to both $θ$-methods for sufficiently small time-step sizes. However, this superior accuracy of the generalized-$α$ scheme does not carry over to larger time steps, due to a steep deterioration of its accuracy in the pre-asymptotic regime. The approximation errors of the three schemes are comparable for time-step sizes for which the diffuse interface moves by approximately its thickness within one time step. Moreover, the generalized-$α$ method is less robust, in the sense that Newton’s method fails to converge at larger time-step sizes.

For applications that require very high accuracy, the superior accuracy of the generalized-$α$ method, deriving from its higher-order asymptotic convergence rate, outweighs its methodological and computational complexity, compared to the $θ$-methods. However, the generalized-$α$ method seems non-competitive if moderate accuracy is required and, accordingly, larger time-step sizes are applied. Its additional features – user-controlled dissipation and the inherent error estimation for adaptive timestepping – are potential advantages, the performance and impact of which warrant further investigation.

T.B. van Sluijs gratefully acknowledges the financial support through the Industrial Partnership Program Fundamental Fluid Dynamics Challenges in Inkjet Printing (FIP), a joint research program of Canon Production Printing, Eindhoven University of Technology, University of Twente, and the Netherlands Organization for Scientific Research (NWO).

The Netherlands Organization for Scientific Research (NWO).