Numerical investigation of the sharp-interface limit of the Navier–Stokes–Cahn–Hilliard equations

Abstract

In this article, we study the behaviour of the Abels–Garcke–Grün Navier–Stokes–Cahn– Hilliard diffuse-interface model for binary-fluid flows, as the diffuse-interface thickness passes to zero. For the diffuse-interface model to approach a classical sharp-interface model in the limit $\varepsilon \to +0$, the so-called mobility parameter $m$ in the diffuse-interface model must scale appropriately with the interface-thickness parameter $\varepsilon$. In the literature various scaling relations in the range $o(1)$ to $O(\varepsilon ^3)$ have been proposed, but the optimal order to pass to the limit has not been explored previously. Our primary objective is to elucidate this optimal order of the $m$–$\varepsilon$ scaling relation in terms of the rate of convergence of the diffuse-interface solution to the sharp-interface solution. Additionally, we examine how the convergence rate is affected by a sub-optimal parameter scaling. We centre our investigation around the case of an oscillating droplet. To provide reference limit solutions, we derive new analytical expressions for small-amplitude oscillations of a viscous droplet in a viscous ambient fluid in two dimensions. For two distinct modes of oscillation, we probe the sharp-interface limit of the Navier–Stokes–Cahn–Hilliard equations by means of an adaptive finite-element method. The adaptive-refinement procedure enables us to consider diffuse-interface thicknesses that are significantly smaller than other relevant length scales in the droplet-oscillation problem, allowing an exploration of the asymptotic regime.

1. Introduction

Binary-fluid flows in which the two fluid components are separated by a molecular transition layer are omnipresent in science and engineering. Examples are inkjet printing and additive manufacturing. Mathematical–physical models for binary-fluid flows generally fall under one of two categories, namely sharp-interface or diffuse-interface models. In sharp-interface models, the surface that separates the two fluid components is represented explicitly by a manifold of co-dimension one. This manifold carries kinematic and dynamic interface conditions, which act as boundary conditions on the initial boundary-value problems of the two contiguous fluid components and, in addition, determine the evolution of the manifold. Sharp-interface models are therefore of free-boundary type. In diffuse-interface models, the interface between the two fluid components is represented as a thin-but-finite transition layer, in which the two components are mixed in a proportion that varies continuously and monotonically between the two pure species across the layer. The strength of diffuse-interface models lies in their intrinsic ability to account for topological changes of the fluid–fluid interface due to coalescence or break-up of droplets or wetting, i.e. the propagation of the fluid–fluid front along a (possibly elastic) solid substrate (Seppecher 1996; Jacqmin 2000; Yue & Feng 2011; van Brummelen, Demont & van Zwieten 2021).

Diffuse-interface models for two immiscible incompressible fluid species are generally described by the Navier–Stokes–Cahn–Hilliard (NSCH) equations. The NSCH equations represent a class of models, of which various renditions have been proposed over the last half-century: by Hohenberg and Halperin in the late 1970s (Hohenberg & Halperin 1977), by Lowengrub and Truskinovsky in the late 1990s (Lowengrub & Truskinovsky 1998), by Shokrpour Roudbari et al. (2018) and by Abels, Garcke & Grün (2012). In this article, we focus on the latter model, in view of its thermodynamic consistency and its consistent reduction to the underlying single-fluid Navier–Stokes equations in the pure species setting.

NSCH models invariably contain three parameters related to the diffuse interface, viz. an interface-thickness parameter, $\varepsilon$, a mobility parameter, $m$, and a surface-tension parameter, $\sigma$. The interface-thickness parameter represents the transverse length scale of the transition layer between the two fluid components, and the transition layer collapses (specifically, is supposed to collapse) onto a manifold of co-dimension one in the so-called sharp-interface limit $\varepsilon \to +0$. The mobility parameter is responsible for the rate at which phase diffusion occurs in the vicinity of the diffuse interface. In the phase-separated regime in which the NSCH equations are typically applied as a binary-fluid model, the mobility parameter is responsible for the rate at which the interface equilibrates. In the mixture regime, it governs the dynamics of the Ostwald-ripening effect. The surface-tension parameter controls the excess free energy $\sigma _{\textrm {DA}}$ of the diffuse interface according to $2\sqrt {2}\sigma =3\sigma _{\textrm {DA}}$. It is to be noted that for the NSCH equations this proportionality holds independent of $\varepsilon$, as opposed to the Navier–Stokes–Korteweg equations.

Contemporary understanding of the sharp-interface limit of the NSCH equations is incomplete. An overview of known results and open questions is provided in Abels & Garcke (2018, § 4.3). One prominent open question pertains to the appropriate scaling of the mobility parameter in relation to the interface-thickness parameter, in the sharp-interface limit. The limit solution of the NSCH equations depends on the scaling $m:=m_{\varepsilon }$. Abels & Garcke (2018, § 4.1) establish that, if $m_{\varepsilon }={O}(1)$ as $\varepsilon \rightarrow +0$, their NSCH model converges to the non-classical sharp-interface Navier–Stokes/Mullins–Sekerka model; see also Jacqmin (1999, § 4). If, on the other hand, $m_{\varepsilon }$ vanishes suitably as $\varepsilon \to +0$, the classical sharp-interface binary-fluid model is obtained, where the interface is transported by the fluid velocity (Lowengrub & Truskinovsky 1998; Yue, Zhou & Feng 2010; Abels & Garcke 2018). However, the decay of the mobility cannot be too fast: if $m_{\varepsilon }=o(\varepsilon ^3)$ as $\varepsilon \to +0$, the resulting limit solution of the NSCH model generally violates the Young–Laplace condition on the pressure jump across the interface (Abels & Lengeler 2014). These results suggest that $m_{\varepsilon }\propto \varepsilon ^{a}$ with $0< a\leqslant 3$ as $\varepsilon \to +0$ represents a necessary and sufficient condition to achieve a classical sharp-interface solution. Still, the details of the approach of the diffuse-interface solution to the sharp-interface limit solution for these various admissible scalings of the mobility are not currently known, and different scalings have been proposed in the literature, in particular in the context of numerical simulation approaches. In Demont et al. (2022), the scaling $m_{\varepsilon } \propto \varepsilon ^3$ is considered, based on the argument that this proportionality fixes the diffusive time scale and, thus, the equilibration rate of the diffuse interface. This cubic scaling of the mobility with respect to the interface thickness (in terms of their usual dimensional forms) is also propounded in Khatavkar, Anderson & Meijer (2006), supported by numerical investigations. Based on partial matched-asymptotic-analysis arguments, Magaletti et al. (2013) concludes that $m \propto \varepsilon ^2$ is the appropriate scaling. However, because the matching procedure in this reference is incomplete, it is unclear whether this scaling relation in fact represents a necessary or sufficient condition. On the basis of a consideration of curvature-induced expansion/contraction modes at the diffuse interface, it is argued in Jacqmin (1999) that $m_{\varepsilon } \propto \varepsilon ^a$ with $1\leqslant a<2$. It is to be noted that the aforementioned scalings of the mobility pertain to situations without moving contact lines and topological changes; see, e.g. Yue et al. (2010).

In this article, we address the open question of the optimal scaling of the mobility parameter in terms of the convergence rate of the diffuse-interface solution to the sharp-interface solution. Specifically, if one considers a strictly decreasing sequence of interface-thickness parameters $\{\varepsilon _k\}$ (with only accumulation point zero) and for each $\varepsilon _k$ determines the mobility $m_k$ such that the deviation between the diffuse-interface solution corresponding to the parameters $\varepsilon _k,m_k$ and the sharp-interface solution (in some appropriate norm) is minimal, is there a proportionality $m_k\propto \varepsilon _k^{a_{opt}}$ as $k\to \infty$? Assuming the answer is affirmative, at which rate does the diffuse-interface solution approach the sharp-interface solution? Additionally, how does the convergence rate deteriorate for sub-optimal scaling relations? We consider these questions based on computational investigation of the Abels–Garcke–Grün NSCH model for different interface dynamics and different mobility parameters as the model limits toward a sharp-interface description of a two-dimensional oscillating droplet. To enable an exploration of the asymptotic regime, we apply an adaptive finite-element method, in which the adaptivity is guided by an a posteriori error estimate; see van Brummelen et al. (2021) and Demont et al. (2022) for details.

We conduct our analysis of the sharp-interface limit of the NSCH equations in the context of the prototypical oscillating-droplet problem, in two dimensions. Despite the fact that the oscillating-droplet problem is classical, it appears that the two-dimensional setting has not been extensively investigated, and that solutions of the two-dimensional problem have not been reported in the literature. The investigation of the oscillating-droplet problem dates back to Rayleigh, who presented the well-known frequency of oscillation of an inviscid droplet in vacuo in 1879 (Strutt 1879). This result was extended by Lamb in the 1930s to include the effect of an inviscid ambient fluid (Lamb 1932). In 1960, Reid generalized the theory of oscillating droplets in vacuo by including the effect of viscosity (Reid 1960). A complete theory, comprising solutions for small oscillations of a viscous droplet in a viscous ambient fluid, was then finally presented by Miller and Scriven in 1968 (Miller & Scriven 1968). The aforementioned references, however, exclusively consider the three-dimensional case and the results, especially those for the viscous solutions, do not trivially extend to the two-dimensional case. Clearly, the three-dimensional case is the practically relevant one, but the two-dimensional case has raison d’être independently as a means of verification for mathematical models and numerical methods. Our analysis of the sharp-interface limit of the NSCH equations requires access to closed-form solutions of the sharp-interface model, on the one hand to provide initial and boundary data for the NSCH equations, and on the other hand to systematically determine the deviation of the diffuse-interface solution relative to the sharp-interface solution. A secondary objective of this work is therefore to establish closed-form expressions for small-amplitude oscillations of a viscous droplet in a viscous ambient fluid. Our derivation follows that of Miller and Scriven, but we deviate from their derivation by a more complete elaboration of intermediate steps and assumptions and, in particular, an explicit accounting of the complex-valued nature of the different fields, and by presenting closed-form expressions of the final results.

The remainder of this article is structured as follows. In § 2, we lay out the Abels–Garcke–Grün NSCH model equations, and the coupled Navier–Stokes free-boundary problem that they should reduce to in the sharp-interface limit. In § 3, we derive a closed-form expression for the sharp-interface model corresponding to small-amplitude oscillations of a viscous droplet in a viscous ambient fluid in two dimensions. We make use of these expressions in § 4, where we study the approach of the NSCH solution to the sharp-interface solution in the limit $\varepsilon \to +0$, by means of systematic numerical experiments.

2. Governing equations

We consider a binary-fluid system, where both fluids are modelled as being incompressible, isothermal, immiscible and Newtonian with finite viscosity. In accordance with the later focus on a submerged droplet, we denote one of the fluids by $\textrm {D}$, for droplet, and the other by $\textrm {A}$, for ambient. Various modelling frameworks for describing the fluid motion exist. These make use of either a diffuse-interface representation or a sharp-interface representation. Figure 1 illustrates the different relevant domains and material parameters. We consider in this work the incompressible NSCH model – specifically, the model developed by Abels et al. (2012) – in order to describe the binary-fluid dynamics. Motivation for this choice lies in its thermodynamic consistency, incompressibility and consistent reduction to the underlying single-fluid Navier–Stokes equations in the pure species setting. Recently, the well posedness of the Abels–Garcke-Grün model in various settings has been shown (Abels, Depner & Garcke 2013a,b; Giorgini 2021).

Figure 1. Schematic of the physical setting of a submerged immiscible fluid, modelled with either a diffuse-interface representation or a sharp-interface representation.

2.1. Diffuse-interface representation

In diffuse-interface models, the two immiscible fluids are separated by a layer of finite thickness constituted by a mixture of both fluids, reflecting a gradual transition between fluid $\textrm {D}$ and fluid $\textrm {A}$. We consider an open time interval $(0, t_{fin}) \subseteq \mathbb{R}_{>0}$ and a spatial domain corresponding to a simply connected time-independent subset $\varOmega \subseteq \mathbb{R}^d$ ($d=2,3$). We make use of a NSCH type model that describes the evolution of a so-called order parameter $\varphi \in [-1,1]$ representing pure species $\textrm {D}$ and $\textrm {A}$ when $\varphi =1$ and $\varphi =-1$, respectively, and a mixture of both when $\varphi \in (-1,1)$, in addition to the velocity and pressure of the mixture. The NSCH model as presented by Abels, Garcke and Grün is given by (Abels et al. 2012)

(2.1a, 2.1b, 2.1c, and 2.1d)

where the volume-averaged velocity $\boldsymbol {u}$, the pressure $p$, the order parameter $\varphi$ and the chemical potential $\mu$ are the unknown fields. The closure relations for the relative mass flux $\boldsymbol {J}$, the viscous stress $\boldsymbol {\tau }$, the capillary stress $\boldsymbol {\zeta }$ and the mixture energy density $\varPsi$ are given as

(2.2a)$$\begin{gather} {} \boldsymbol{J} := m \frac{ \rho_{\textrm{A}} - \rho_{\textrm{D}} } { 2 } \boldsymbol{\nabla} \mu , \end{gather}$$

(2.2b)$$\begin{gather}\boldsymbol{\tau} := \eta ( \boldsymbol{\nabla} \boldsymbol{u} + ( \boldsymbol{\nabla} \boldsymbol{u} )^{\rm T} ) , \end{gather}$$

(2.2c)$$\begin{gather}\boldsymbol{\zeta} := - \sigma \varepsilon \boldsymbol{\nabla} \varphi \otimes \boldsymbol{\nabla} \varphi + \boldsymbol{I} \left( \frac{ \sigma \varepsilon } { 2 } | \boldsymbol{\nabla} \varphi |^2 + \frac{ \sigma } { \varepsilon } \varPsi \right) , \end{gather}$$

(2.2d)$$\begin{gather}\varPsi \left( \varphi \right) := \tfrac{1}{4} (\varphi^2 -1)^2 . \end{gather}$$

The remaining parameters are material and model parameters. The model parameters are the mobility parameter $m>0$ and the interface-thickness parameter $\varepsilon >0$, which affect the time and length scale of the diffuse interface, respectively. The material parameters are $\sigma$, a rescaling of the droplet-ambient surface tension $\sigma _\textrm {DA}$ according to $2\sqrt {2}\sigma =3\sigma _\textrm {DA}$, the mixture density $\rho$ and the mixture viscosity $\eta$. The mixture density and viscosity generally depend on $\varphi$. To ensure existence of a solution to the system of equations, we must allow $\varphi$ to take on values outside of $[-1,1]$ (Grün, Guillén-González & Metzger 2016). We include a density extension that ensures positive densities even for the non-physical scenario $\varphi \notin [-1,1]$ (Bonart, Kahle & Repke 2019)

(2.3) \begin{equation} \rho(\varphi) = \begin{cases} \frac{ 1 } { 4 } \rho_{\textrm{A}}, & \varphi \leqslant - 1 - 2 \lambda ,\\ \frac{ 1 } { 4 } \rho_{\textrm{A}} + \frac{ 1 } { 4 } \rho_{\textrm{A}} \lambda^{{-}2} ( 1 + 2 \lambda + \varphi )^2, & \varphi \in ( - 1 - 2 \lambda , - 1 - \lambda) ,\\ \dfrac{ 1 + \varphi } { 2 } \rho_{\textrm{D}} + \dfrac{ 1 - \varphi } { 2 } \rho_{\textrm{A}}, & \varphi \in [ - 1 - \lambda, 1 + \lambda ] ,\\ \rho_{\textrm{D}} + \frac{ 3 } { 4 } \rho_{\textrm{A}} - \frac{ 1 } { 4 } \rho_{\textrm{A}} \lambda^{{-}2} ( 1 + 2 \lambda - \varphi )^2, & \varphi \in ( 1 + \lambda, 1 + 2 \lambda ) ,\\ \rho_{\textrm{D}} + \frac{ 3 } { 4 } \rho_{\textrm{A}}, & \varphi \geqslant 1 + 2 \lambda , \end{cases} \end{equation}

where $\lambda = \rho _{\textrm {A}} / ( \rho _{\textrm {D}} - \rho _{\textrm {A}} )$. For the viscosity interpolation, we apply the Arrhenius mixture-viscosity model (Arrhenius 1887)

(2.4)\begin{equation} \log \eta( \varphi ) = \frac{ \left( 1 + \varphi \right) \varLambda \log \eta_{\textrm{D}} + \left( 1 - \varphi \right) \log \eta_{\textrm{A}} } { \left( 1 + \varphi \right) \varLambda + \left( 1 - \varphi \right) } \,, \end{equation}

where $\varLambda = { \rho _{\textrm {D}} M_{\textrm {A}} }/{\rho _{\textrm {A}} M_{\textrm {D}} }$ is the intrinsic volume ratio between the two fluids (with $M_{\textrm {A}}$ and $M_{\textrm {D}}$ their respective molar masses).

Remark 2.1 To eliminate the Ostwald-ripening effect in the pure species, a degenerate dependence of the mobility on the phase field can be introduced, according to $m(\varphi ) \geqslant 0$ with inequality if and only if $|\varphi |<1$. However, as a degenerate mobility introduces complications with regard to numerical-approximation procedures (Barrett, Blowey & Garcke 1999), we opt for a constant mobility parameter.

Remark 2.2 In this work, we use a volume-fraction-based Arrhenius relation, i.e. $\varLambda =1$. Because the denominator in (2.4) then reduces to a non-zero constant, this choice eliminates singularities, and the mixture viscosity is bounded away from zero in a finite interval including $[-1,1]$. See Remark 2 in van Brummelen et al. (2021) for further details.

2.2. Sharp-interface limit

As $\varepsilon \rightarrow +0$, the width of the diffuse interface in the NSCH model reduces to zero. As pointed out in the introduction, the particular model that arises in this limit depends on the scaling relation of the mobility $m$. If the mobility also tends to zero appropriately, then the following classical sharp-interface model is obtained:

(2.5a)$$\begin{gather} \rho_{i}\partial_t \boldsymbol{u}_i + \rho_{i} \left( \boldsymbol{u}_i \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u}_i - \eta_{i} \Delta \boldsymbol{u}_i + \boldsymbol{\nabla} p_i = 0\quad \textrm{in}\ \varOmega_{i} = \varOmega_{i}(t) , \end{gather}$$

(2.5b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_i = 0 \quad \textrm{in } \varOmega_{i} , \end{gather}$$

(2.5c)$$\begin{gather}\boldsymbol{u}_{i}\boldsymbol{\cdot}\boldsymbol{n} = \mathcal{V}\quad \textrm{on } \varGamma = \varGamma(t) , \end{gather}$$

(2.5d)$$\begin{gather} {[}\![ \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{t}_j ]\!] = 0 \quad \textrm{on } \varGamma \text{, for }j=1,\ldots,d-1 , \end{gather}$$

(2.5e)$$\begin{gather} {[}\![ -( \boldsymbol{\nabla} \boldsymbol{u} + ( \boldsymbol{\nabla} \boldsymbol{u} )^{\rm T} ) \boldsymbol{n} + p \boldsymbol{n} ]\!] = \sigma_{\textrm{D}\textrm{A}} \kappa \boldsymbol{n} \quad \textrm{on } \varGamma , \end{gather}$$

for $i \in \{\textrm {D},\textrm {A}\}$, and with $\boldsymbol {n}$ the unit normal vector on $\varGamma$ external to $\varOmega _{\textrm {D}}$, $\mathcal {V}$ the interface normal velocity, $\kappa$ the (additive) curvature of the interface and $[\![ {\cdot } ]\!]$ the interface jump operator $[\![ g ]\!] = g|_{\textrm {D}} - g|_{\textrm {A}}$. We adhere to the convention that curvature is negative if the centre of the osculating circle in the normal plane is located in the droplet domain.

As opposed to the NSCH model, the sharp-interface model (2.5) represents a set of equations for each fluid species separately, complemented by appropriate coupling conditions at the interface. The sharp-interface model represents a free-boundary problem. The domains on which the various fields are defined evolve in time, as reflected by the time dependence of $\varOmega _i(t)$ and $\varGamma (t)$. There is an intrinsic coupling between the velocity field and the evolution of $\varOmega _i$ and $\varGamma$, according to (2.5c). Because this equation holds on both sides of the interface and $\mathcal {V}$ is single valued, (2.5c) implies $[\![ \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n} ]\!] = 0$.

3. Response of an oscillating droplet

To provide a reference solution for the sharp-interface limit of the NSCH equations, we now consider solutions of the free-boundary problem (2.5) corresponding to small perturbations of a circular droplet set in an ambient fluid, in two dimensions. Readers mainly interested in the asymptotic results of the NSCH model may wish to proceed directly to § 4.

Denoting by $R_0$ the radius of the droplet, one can verify that

(3.1a)$$\begin{gather} \varOmega_{\textrm{D},0}=\{\boldsymbol{x}\in\mathbb{R}^2:|\boldsymbol{x}|< R_0\}, \varOmega_{\textrm{A},0}=\mathbb{R}^2\setminus\overline{\varOmega_{\textrm{D},0}}, \end{gather}$$

(3.1b)$$\begin{gather}\boldsymbol{u}_{\textrm{D},0}=0, \boldsymbol{u}_{\textrm{A},0}=0, \end{gather}$$

(3.1c)$$\begin{gather}p_{\textrm{D},0}=\sigma_{\textrm{D}\textrm{A}}/R_0, p_{\textrm{A},0}=0, \end{gather}$$

represents a stationary solution to (2.5). We will use (3.1) as a generating solution. We consider perturbations of the solution (3.1) that are suitably bounded and vanish toward infinity. The latter condition can be expressed as

(3.2)\begin{equation} \lim_{|\boldsymbol{x}|\rightarrow\infty} \left(\boldsymbol{u},p\right)(\boldsymbol{x},t) = 0 . \end{equation}

Our derivation of the natural response of such a droplet follows that of Miller & Scriven (1968), except that we provide a more complete elaboration of intermediate steps and assumptions and, in particular, an explicit accounting of the complex-valued nature of the different fields. The approach essentially comprises four steps. First, we linearize the governing equations around the generating solution (3.1), perturbed by a small deformation of the interface. Second, we derive the general solutions corresponding to the natural response in both domains separately. In the third step, we incorporate the interface coupling conditions by constraining the free parameters in the general solutions. Finally, the characteristic temporal response (frequency of oscillation and rate of damping) of the assumed interface displacement, as well as the corresponding shapes, follow from a solution-existence condition.

3.1. Formal linearization

We consider small-amplitude perturbations of the interface that are sinusoidal along the circumference of the droplet. To facilitate the presentation, we introduce polar coordinates $r\in \mathbb {R}_{\geqslant 0}$ and $\theta \in [0,2{\rm \pi} )$ and the coordinate transformation $\boldsymbol {x}=(x_1,x_2)=r(\cos \theta,\sin \theta )$. We regard perturbations of the interface $\varGamma _0=\partial \varOmega _{\textrm {D},0}$ corresponding to the following parametrization:

(3.3)\begin{equation} \varGamma_{\delta}(t)=\{\boldsymbol{x}\in\mathbb{R}^2: \boldsymbol{x}=R_{\delta}(\theta,t)\,(\cos(\theta),\sin(\theta)),\theta\in[0,2{\rm \pi})\} , \end{equation}

where

(3.4) \begin{align} R_{\delta} ( \theta , t ) & = R_0 + R_0\, \delta( \beta \cos ( k \theta ) + \sqrt{1-\beta^2 } \, \sin ( k \theta ) ) \cos(\nu t) {\rm e}^{- \alpha t} \nonumber\\ & = {\rm Re} ( R_0 + R_0\, \delta ( \beta \cos ( k \theta ) + \sqrt{1-\beta^2 } \, \sin ( k \theta ) ) {\rm e}^{- \gamma t} ) . \end{align}

The interface configuration (3.3)–(3.4) represents a damped oscillation of the droplet with a mode shape described by the mode number $k \in \mathbb {N}$, with angular orientation dependent on $0\leqslant \beta \leqslant 1$ and with an amplitude described by $\delta \ll {}1$ as the fraction of the droplet radius $R_0$. Our interest is restricted to droplet configurations for which $\operatorname {meas}(\varOmega _{\textrm {D}})=\operatorname {meas}(\varOmega _{\textrm {D},0})+{O}(\delta ^2)$ and the barycentre of $\varOmega _{\textrm {D}}$ is located at the origin. This implies that $k \in \mathbb {N}_{\geqslant 2}$. The damping rate $\alpha \geqslant 0$ and the frequency of oscillation $\nu > 0$ implicitly depend on the mode number and will follow from the subsequent analysis. The second expression in (3.4) provides a representation of $R_{\delta }$ as the real part of a complex-valued function, with $\gamma := \alpha - {\rm i} \nu$. This form enables us to condense some of the expressions that appear in the sequel.

In conjunction with the interface configuration (3.3)–(3.4), we consider linear asymptotic solutions of the sharp-interface problem of the form

(3.5)\begin{equation} \left(\boldsymbol{u}_{i},p_i,\boldsymbol{n}, \boldsymbol{t},\kappa,\mathcal{V}\right)= \left(\boldsymbol{u}_{i},p_i,\boldsymbol{n}, \boldsymbol{t},\kappa,\mathcal{V}\right)_0 +\delta\left(\boldsymbol{u}_{i},p_i,\boldsymbol{n}, \boldsymbol{t},\kappa,\mathcal{V}\right)_1, \end{equation}

i.e. functions conforming to (3.5) that satisfy (2.5) modulo terms of $o(\delta )$ as $\delta \to {}0$. Substituting (3.5) into the sharp-interface equations (2.5), collecting terms of distinct orders in $\delta$ and noting that all terms of ${O}(1)$ vanish on account of the fact that the first term in (3.5) represents a solution to (2.5), we obtain the following infinitesimal conditions on the second term in (3.5):

(3.6a)$$\begin{gather} \rho\partial_t \boldsymbol{u}_{i,1} - \eta \Delta \boldsymbol{u}_{i,1} + \boldsymbol{\nabla} p_{i,1} = 0 \quad \textrm{in } \varOmega_{i,0} , \end{gather}$$

(3.6b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_{i,1} = 0 \quad \textrm{in } \varOmega_{i,0} , \end{gather}$$

(3.6c)$$\begin{gather}\boldsymbol{u}_{i,1} \boldsymbol{\cdot} \boldsymbol{n}_0 = \mathcal{V}_1 \quad \textrm{on } \varGamma_0 , \end{gather}$$

(3.6d)$$\begin{gather} {[}\![ \boldsymbol{u}_{1} \boldsymbol{\cdot} \boldsymbol{t}_0 ]\!] = 0 \quad \textrm{on } \varGamma_0 , \end{gather}$$

(3.6e)$$\begin{gather} {[}\![ - \eta \left(\boldsymbol{\nabla} \boldsymbol{u}_{1} + \left(\boldsymbol{\nabla} \boldsymbol{u}_{1} \right)^{\rm T} \right) \boldsymbol{n}_0 + p_{1} \boldsymbol{n}_0 ]\!] = \sigma_{\textrm{D}\textrm{A}} \kappa_1 \boldsymbol{n}_0 \quad \textrm{on } \varGamma_0 , \end{gather}$$

for $i \in \{\textrm {D},\textrm {A}\}$. The nonlinear advective term has dropped since the only nonlinear first-order perturbation terms are cross-terms between $\boldsymbol {u}_{i,1}$ and $\boldsymbol {u}_{i,0} = 0$. Similarly, only $\boldsymbol {n}_0$ appears in the first-order conditions (3.6), because its multiplication with $\boldsymbol {u}_{i,1}$ is a second-order term, its dot product with $\boldsymbol {u}_{i,0}$ vanishes and $[\![ p_0 \boldsymbol {n}_1 ]\!] = \sigma _{\textrm {D}\textrm {A}} R^{-1}_0 \boldsymbol {n}_1$ cancels with the right-hand side $\sigma _{\textrm {D}\textrm {A}} \kappa _0 \boldsymbol {n}_1$.

Remark 3.1 The first-order conditions are set on the stationary generating domains, $\varOmega _{i,0}$, and the generating interface, $\varGamma _0$. This is a universal characteristic of linearizations of free-boundary problems. The perturbation of the interface according to (3.4) appears implicitly in (3.6) in the interface-velocity perturbation, $\mathcal {V}_1$, and the curvature perturbation, $\kappa _1$.

The remainder of this section is devoted to finding general solutions to (3.6) for different perturbation wavenumbers $k$. For purposes of readability, henceforth we suppress the subscripts corresponding to the order of perturbation.

3.2. First-order solutions in the droplet and ambient domain

Next, we derive the general first-order solutions $(\boldsymbol {u},p)_{i,1}$ in accordance with the differential equations (3.6a)–(3.6b) for both the droplet $\textrm {D}$ and ambient $\textrm {A}$ domains. We proceed by deriving the vorticity equation corresponding to (3.6a), which we solve by means of separation of variables. The first-order velocity field is subsequently retrieved from the vorticity solutions. The corresponding first-order pressure fields are derived as those that yield balance of linear momentum.

3.2.1. Vorticity solution

The pressure may be eliminated from the governing equations by taking the curl of (3.6a) and by using the identity $\boldsymbol {\nabla } \times \boldsymbol {\nabla }({\cdot }) = 0$. Introducing the vorticity $\omega =\boldsymbol {\nabla } \times \boldsymbol {u}$, we infer from (3.6a) that

(3.7)\begin{equation} \rho \partial_t \left( \boldsymbol{\nabla} \times \boldsymbol{u} \right) - \eta \varDelta \left( \boldsymbol{\nabla} \times \boldsymbol{u} \right) + \boldsymbol{\nabla} \times\boldsymbol{\nabla}p= \rho \partial_t \omega - \eta \Delta \omega = 0. \end{equation}

Let us note that in a two-dimensional setting, vorticity can be represented as a scalar-valued field. The evolution equation for this scalar vorticity field can be recognized as a diffusion equation.

To determine the general solution to (3.7), we assume the following separation of variables form:

(3.8)\begin{equation} \omega ( r,\theta , t ) := P(r) \varTheta(\theta) T(t). \end{equation}

Substitution in the polar coordinate representation of the diffusion equation gives

(3.9)\begin{equation} P(r) \varTheta(\theta) T'(t) = \frac{\eta}{\rho r^2}P(r) \varTheta''(\theta) T(t)+ \frac{\eta}{\rho r}P'(r) \varTheta(\theta) T(t) + \frac{\eta}{\rho}P''(r) \varTheta(\theta) T(t) , \end{equation}

where primes denote differentiation. The usual separation of variables argument leads to

(3.10a)$$\begin{gather} T'(t) ={-}\frac{\eta}{ \rho} m^2 T(t) \quad m \in \mathbb{C}, \end{gather}$$

(3.10b)$$\begin{gather}\varTheta''(\theta) ={-}n^2 \varTheta(\theta) \quad n \in \mathbb{C} , \end{gather}$$

(3.10c)$$\begin{gather}r^2 P''(r) + r P'(r) + (m^2 r^2 - n^2) P(r) = 0 , \end{gather}$$

with $\mathbb {C}$ the set of complex numbers.

The general solutions of (3.10a) and (3.10b) consist of complex-valued exponential functions, according to

(3.11a)$$\begin{gather} T (t) = c\, {\rm e}^{- ({\eta}/{ \rho}) m^2 t}, \end{gather}$$

(3.11b)$$\begin{gather}\varTheta (\theta) = c_1\, {\rm e}^{ i n t} + c_2 {\rm e}^{ {-}i n t}. \end{gather}$$

From the periodicity of the droplet perturbations in the angular dependence, conforming to (3.3), we infer that $n\in \mathbb {Z}_{\geqslant {}0}$. The arbitrary constants $c_1$ and $c_2$ can then be selected such that (3.11b) reduces to the sum of two real-valued trigonometric functions

(3.12)\begin{equation} \varTheta (\theta) = C\cos( n \theta) + D\sin( n \theta)\quad (n \in \mathbb{Z}_{\geqslant0}) , \end{equation}

where $C,D \in \mathbb {R}$ are coefficients that determine the angular orientation of the solution. Regarding (3.10c), we note that this equation corresponds to Bessel's equation with a complex-valued scaling $m\in \mathbb {C}$. Solutions of (3.10c) therefore consist of extensions of Bessel functions to the complex plane. Such extensions of Bessel functions are well defined, by virtue of the fact that Bessel functions are analytic functions on $\mathbb {R}$ and can hence be extended to analytic functions on $\mathbb {C}$ via their power-series expansion. The general solution of (3.10c) consists of a linear combination of two Bessel functions of order $n\in \mathbb {Z}_{\geqslant 0}$. For reasons that will become clear once we consider the boundary conditions, we choose to work with a Bessel function of the first kind, ${\rm J}_n$, and a Hankel function of the second kind, $H^{(2)}_n$,

(3.13)\begin{equation} P(r) = A m^2 {\rm J}_n ( m r ) + B m^2 H^{(2)}_n ( m r ) , \end{equation}

with $A,B \in \mathbb {C}$. Substitution of (3.11a), (3.12) and (3.13) into (3.8) gives the general rotationally periodic solution of the vorticity equation (3.7)

(3.14)\begin{equation} \omega(r,\theta,t) = [A m^2 {\rm J}_n ( m r ) + B m^2 H^{(2)}_n ( m r )] \left[ C\cos( n \theta) + D\sin( n \theta) \right] {\rm e}^{- ({\eta}/{ \rho}) m^2 t} , \end{equation}

for arbitrary $A,B\in \mathbb {C}$, $C,D\in \mathbb {R}$ and $m\in \mathbb {C}$.

3.2.2. Velocity solutions

To obtain the velocity fields from the general vorticity solution (3.14), we introduce a streamfunction $\chi$ according to

(3.15)\begin{equation} \Delta \chi ={-} \omega . \end{equation}

The velocity can be retrieved from this streamfunction as

(3.16)\begin{equation} \boldsymbol{u} = \frac{1}{r} \frac{\partial}{\partial\theta} \chi\, \boldsymbol{e}_r - \frac{\partial}{\partial r} \chi \, \boldsymbol{e}_\theta. \end{equation}

Based on the expression for $\omega$ in (3.14), we anticipate that a particular solution to (3.15) is of the form

(3.17)\begin{equation} \chi (r,\theta,t) = \varUpsilon ( r ) \left[ C\cos( n \theta) + D\sin( n \theta) \right] {\rm e}^{- ({\eta}/{ \rho}) m^2 t} . \end{equation}

Substitution of (3.17) into (3.15) leads to the following ordinary differential equation for $\varUpsilon$:

(3.18)\begin{equation} \varUpsilon'' ( r ) + \frac{1}{r} \varUpsilon' ( r ) - \frac{n^2}{r^2} \varUpsilon ( r ) ={-} A m^2 {\rm J}_n ( m r ) - B m^2 H^{(2)}_n ( m r ) . \end{equation}

The general solution to the non-homogeneous ordinary differential equation (3.18) is given by

(3.19)\begin{equation} \varUpsilon(r) = E r^{n} + F r^{{-}n} + A {\rm J}_n ( m r) + B H^{(2)}_{n} ( m r ) . \end{equation}

The first and second terms in (3.19) constitute the homogeneous part of the solution. The third and fourth terms represent the particular part. From (3.16) it then follows that the $r$ and $\theta$ components of the velocity field are given by

(3.20a)$$\begin{gather} u_r = \frac{1}{r} \frac{\partial}{\partial\theta} \chi = \frac{n}{r} \varUpsilon ( r ) \left[ D \cos( n \theta) - C \sin( n \theta) \right] {\rm e}^{- ({\eta}/{\rho}) m^2 t} \,, \end{gather}$$

(3.20b)$$\begin{gather}u_\theta ={-} \frac{\partial}{\partial r} \chi ={-} \varUpsilon' ( r ) \left[ C \cos( n \theta) + D \sin( n \theta) \right] {\rm e}^{- ({\eta}/{ \rho}) m^2 t} . \end{gather}$$

The velocity solutions (3.20) are not generally bounded in the limits $r\to {}0$, $r\to \infty$ and $t\to \infty$. Auxiliary conditions must be imposed on the coefficients in (3.20) to extract general solutions that are bounded in the droplet domain $\varOmega _{\textrm {D},0}$ (respectively ambient domain $\mathbb {R}^2\setminus \overline {\varOmega _{\textrm {D},0}}$) as $r\to {}0$ (respectively $r\to \infty$) and as $t\to \infty$. To ensure boundedness of the solutions in the limit $t\to \infty$, we insist that ${\rm Re}(m^2)\geqslant {}0$. To assess the boundedness of the solutions (3.20) in the spatial dependence, we note that the Bessel function ${\rm J}_n(mr)$ is singular at $r\to \infty$ for all $m\in \mathbb {C}$ such that $|{\rm Im}(m)|>0$, and the Hankel function $H^{(2)}_{n}(mr)$ is singular at the origin and at $r\to \infty$ for all $m\in \mathbb {C}$ with ${\rm Im}(m)>0$. In addition, in relation to (3.2), we note that $H^{(2)}_{n}(mr)$ vanishes in the limit $r\to \infty$ if ${\rm Im}(m)\leqslant {}0$. Moreover, $r^n$ (respectively $r^{-n}$) is singular in the limit $r\to \infty$ (respectively $r\to {}0$). On account of their singularity at the origin, $H^{(2)}_{n}$ and $r^{-n}$ are inadmissible in the droplet domain $\varOmega _{\textrm {D},0}\supset \{0\}$. Hence, in the droplet domain, it must hold that $B=0$ and $F=0$. Conversely, ${\rm J}_n(mr)$ and $r^n$ are inadmissible in the ambient domain, in view of their singularity at $r\to \infty$. Hence, in the ambient domain, it must hold that $A=0$ and $E=0$.

Summarizing, we obtain the following general bounded complex-valued velocity solutions in the droplet and ambient domains:

(3.21a)\begin{align} u_{\textrm{D},r} &= \exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) \left[ D_\textrm{D} \cos( n_\textrm{D} \theta) - C_\textrm{D} \sin( n_\textrm{D} \theta) \right]\nonumber\\ &\quad\times\frac{n_\textrm{D}}{r} [ A {\rm J}_{n_\textrm{D}} ( m_\textrm{D} r ) + E r^{n_\textrm{D}}], \end{align}

(3.21b)\begin{align} u_{\textrm{D},\theta}&= \exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) \left[ C_\textrm{D} \cos( n_\textrm{D} \theta) + D_\textrm{D} \sin( n_\textrm{D} \theta) \right]\nonumber\\ &\quad \times \left[- A \left( m_\textrm{D} {\rm J}_{n_\textrm{D}-1} ( m_\textrm{D} r ) - \frac{n_\textrm{D}}{r} {\rm J}_{n_\textrm{D}} ( m_\textrm{D} r ) \right) - E n_\textrm{D} r^{n_\textrm{D}-1} \right], \end{align}

(3.21c)\begin{align} u_{\textrm{A},r}&= \exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right) \left[ D_\textrm{A} \cos( n_\textrm{A} \theta) - C_\textrm{A} \sin( n_\textrm{A} \theta) \right]\nonumber\\ &\quad \times\frac{n_\textrm{A}}{r}[B H^{(2)}_{n_\textrm{A}} ( m_\textrm{A} r ) + F r^{{-}n_\textrm{A}}], \end{align}

(3.21d)\begin{align} u_{\textrm{A},\theta}&= \exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right) \left[ C_\textrm{A} \cos( n_\textrm{A} \theta) + D_\textrm{A} \sin( n_\textrm{A} \theta) \right]\nonumber\\ &\quad \times \left[ - B \left( m_\textrm{A} H^{(2)}_{n_\textrm{A}-1} ( m_\textrm{A} r ) - \frac{n_\textrm{A}}{r} H^{(2)}_{n_\textrm{A}} ( m_\textrm{A} r ) \right) + F n_\textrm{A} r^{{-}n_\textrm{A}-1} \right] , \end{align}

subject to ${\rm Re}(m_i^2)\geqslant {}0$ $(i\in \{\textrm {D},\textrm {A}\})$ and ${\rm Im}(m_{\textrm {A}})<0$.

3.2.3. Pressure solutions

To facilitate the derivation of the infinitesimal pressure solutions associated with (3.21), we first note that, by virtue of (3.6a) and (3.6b), the pressure solutions are harmonic functions. Considering functions that are sinusoidal and periodic in the angular dependence, that conform to (3.21) in the temporal dependence, and that are appropriately bounded, we find the following general expression for $p_{\textrm {D}}$:

(3.22)\begin{equation} p_{\textrm{D}}(r,\theta,t)= \exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) r^{\tilde{n}_{\textrm{D}}}(\tilde{D}_{\textrm{D}}\cos(\tilde{n}_{\textrm{D}}\theta)+\tilde{C}_{\textrm{D}}\sin(\tilde{n}_{\textrm{D}}\theta)), \end{equation}

with $\tilde {n}_{\textrm {D}}\in \mathbb {Z}_{\geqslant 0}$ and $\tilde {C}_{\textrm {D}},\tilde {D}_{\textrm {D}}\in \mathbb {C}$ arbitrary constants. By substituting (3.21) and (3.22) into (3.6a), we deduce that compatibility between (3.22) and (3.21) imposes that the index $\tilde {n}_{\textrm {D}}\in \mathbb {Z}_{\geqslant 0}$ and coefficients $\tilde {C}_{\textrm {D}},\tilde {D}_{\textrm {D}}\in \mathbb {C}$ satisfy

(3.23a–c)\begin{equation} \tilde{D}_{\textrm{D}}=D_{\textrm{D}}E\eta_{\textrm{D}}m_{\textrm{D}}^2,\quad \tilde{C}_{\textrm{D}}={-}C_{\textrm{D}}E\eta_{\textrm{D}}m_{\textrm{D}}^2,\quad \tilde{n}_{\textrm{D}}=n_{\textrm{D}}. \end{equation}

The infinitesimal pressure solution in the ambient domain can be determined similarly. In summary, we obtain

(3.24a)$$\begin{gather} p_{\textrm{D}} = E \eta_\textrm{D} m_\textrm{D}^2\exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) r^{n_\textrm{D}} \left[ D_\textrm{D} \cos(n_\textrm{D} \theta) - C_\textrm{D} \sin(n_\textrm{D} \theta) \right] , \end{gather}$$

(3.24b)$$\begin{gather}p_{\textrm{A}} ={-} F \eta_\textrm{A} m_\textrm{A}^2\exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right) r^{{-}n_\textrm{A}} \left[ D_\textrm{A} \cos(n_\textrm{A} \theta) - C_\textrm{A} \sin(n_\textrm{A} \theta) \right] . \end{gather}$$

Remark 3.2 It is noteworthy that the Bessel and Hankel functions that appear in the velocity solutions (3.21) are absent from the pressure solutions (3.24). This can be rationalized by noting that these special functions originated directly from the (pressure free) vorticity equation and thus satisfy the momentum equation for a uniform pressure field.

3.3. Interface conditions

The general solutions for the pressure and velocity fields in the droplet and ambient domains according to (3.21) and (3.24), involve twelve unknown coefficients: $A$, $B$, $C_\textrm {D}$, $C_\textrm {A}$, $D_\textrm {D}$, $D_\textrm {A}$, $E$, $F$, $m_\textrm {D}$, $m_\textrm {A}$, $n_\textrm {D}$ and $n_\textrm {A}$. We next extract from the general solutions the subspace that complies with the kinematic interface conditions (3.6c) and (3.6d), and the dynamic condition (3.6e), by introducing auxiliary conditions on the coefficients.

3.3.1. Kinematic compatibility conditions

The kinematic conditions (3.6c)–(3.6d) can be equivalently reformulated as

(3.25a)$$\begin{gather} \boldsymbol{u}_\textrm{D} \boldsymbol{\cdot} \boldsymbol{n}_0 = \boldsymbol{u}_\textrm{A} \boldsymbol{\cdot} \boldsymbol{n}_0= \mathcal{V}_1 \quad \textrm{on } \varGamma_0 , \end{gather}$$

(3.25b)$$\begin{gather}\boldsymbol{u}_\textrm{D} \boldsymbol{\cdot} \boldsymbol{t}_0 = \boldsymbol{u}_\textrm{A} \boldsymbol{\cdot} \boldsymbol{t}_0 \quad \textrm{on } \varGamma_0 . \end{gather}$$

The generating solutions of the interface normal and tangent vectors corresponding to the circular droplet, $\boldsymbol {n}_0$ and $\boldsymbol {t}_0$, simply coincide with the radial and angular basis vectors, respectively. The kinematic condition (3.25a) (respectively (3.25b)) thus pertains to the radial (respectively angular) components of $\boldsymbol {u}$ in (3.21a) and (3.21c) (respectively (3.21b) and (3.21d)). To impose (3.25a), we require the infinitesimal interface velocity $\mathcal {V}_1$ corresponding to (3.4)

(3.26)\begin{equation} \mathcal{V}(\theta,t) = \partial_t{}R_{\delta}(\theta,t) ={-}\delta R_0 \gamma\, {\rm e}^{- \gamma t}( \beta \cos ( k \theta ) + \sqrt{1-\beta^2 } \sin ( k \theta ) ) , \end{equation}

where we retain the entire complex form, to facilitate the exposition. Noting that $\mathcal {V}_0$ vanishes and, hence, $\mathcal {V}=\delta \mathcal {V}_1$, we infer from (3.25a) that

(3.27) \begin{align} &\exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right)\left[ D_\textrm{D} \cos( n_\textrm{D} \theta) - C_\textrm{D} \sin( n_\textrm{D} \theta) \right] [ R^{{-}1}_0 A n_\textrm{D} {\rm J}_{n_\textrm{D}} ( m_\textrm{D} R_0 ) + E n_\textrm{D} R_0^{n_\textrm{D}-1} ]\nonumber\\ &\quad= \exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right)\left[ D_\textrm{A} \cos( n_\textrm{A} \theta) - C_\textrm{A} \sin( n_\textrm{A} \theta) \right]\nonumber\\ &\qquad \times [ R^{{-}1}_0 B n_\textrm{A} H^{(2)}_{n_\textrm{A}} ( m_\textrm{A} R_0 ) + F n_\textrm{A} R_0^{{-}n_\textrm{A}-1}]\nonumber\\ &\quad={-}R_0 \gamma\, {\rm e}^{- \gamma t} (\beta \cos ( k \theta ) + \sqrt{1-\beta^2 }\sin ( k \theta )) . \end{align}

The equalities in (3.27) must hold for all $\theta \in [0,2{\rm \pi} )$ and all $t\in \mathbb {R}_{>0}$. Keeping $\theta$ fixed and varying $t$, one can infer that the temporal exponents must coincide. Subsequently, by fixing $t$ and varying $\theta$, it follows that all parameters that characterize the trigonometric terms must be the same. Hence,

(3.28a)$$\begin{gather} \frac{\eta_\textrm{D}}{\rho_\textrm{D}} m_\textrm{D}^2 = \frac{\eta_\textrm{A}}{\rho_\textrm{A}} m_\textrm{A}^2 = \gamma , \end{gather}$$

(3.28b)$$\begin{gather}n_\textrm{D} = n_\textrm{A} = k , \end{gather}$$

(3.28c)$$\begin{gather}D_\textrm{D} = D_\textrm{A} = \beta , \end{gather}$$

(3.28d)$$\begin{gather}C_\textrm{D} = C_\textrm{A} ={-}\sqrt{1-\beta^2} . \end{gather}$$

In the sequel, we continue to use $m_\textrm {D}$ and $m_\textrm {A}$ in the arguments of the Bessel and Hankel functions, but we tacitly suppose the relation to $\gamma$ per (3.28a). Substitution of (3.28) in (3.25) yields the following three conditions on $A,B,E$ and $F$:

(3.29a)$$\begin{gather} A k R^{{-}1}_0 {\rm J}_k (m_\textrm{D} R_0) + E k R_0^{k-1} ={-} \gamma R_0 , \end{gather}$$

(3.29b)$$\begin{gather}B k R_0^{{-}1} H^{(2)}_{k} ( m_\textrm{A} R_0 ) + F k R_0^{{-}k-1} ={-} \gamma R_0 , \end{gather}$$

(3.29c) \begin{align} &- A ( m_\textrm{D} {\rm J}_{k-1} ( m_\textrm{D} R_0) - k R_0^{{-}1} {\rm J}_k ( m_\textrm{D} R_0 ) ) -E k R_0^{k-1}\nonumber\\ &\quad + B ( m_\textrm{A} H_{k-1}^{(2)}(m_\textrm{A} R_0) - k R_0^{{-}1} H_k ^{(2)}(m_\textrm{A} R_0 )) - F k R_0^{{-}k-1} = 0 . \end{align}

3.3.2. Dynamic compatibility conditions

The dynamic condition (3.6e) can be separated into radial and angular components to form the following two conditions:

(3.30a)\begin{align} &- \eta_{\textrm{D}} \boldsymbol{n}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} )^{\rm T} ) \boldsymbol{n}_0 + p_{\textrm{D}} \nonumber\\ &\quad + \eta_{\textrm{A}} \boldsymbol{n}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} )^{\rm T} ) \boldsymbol{n}_0 - p_{\textrm{A}} = \sigma_{\textrm{D}\textrm{A}} \kappa_1 \quad \textrm{on } \varGamma , \end{align}

(3.30b)\begin{align} &- \eta_{\textrm{D}} \boldsymbol{t}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} )^{\rm T} ) \boldsymbol{n}_0 + \eta_{\textrm{A}} \boldsymbol{t}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} )^{\rm T} ) \boldsymbol{n}_0 = 0 \quad \textrm{on } \varGamma. \end{align}

To elaborate on the dynamic interface condition (3.30a), we require the first-order perturbation of the interface curvature. From the postulated interface displacement (3.4), the complex-valued form of the curvature can be derived up to second-order terms

(3.31)\begin{align} \kappa(\theta,t) & = \kappa_0+\delta\kappa_1(\theta,t)+{O}(\delta^2)\nonumber\\ & =R_0^{{-}1} + \delta R^{{-}1}_0 e^{-\gamma t} ( k^2 - 1 ) ( \beta \cos ( k \theta ) - \sqrt{1-\beta^2} \sin ( k \theta)) + {O} (\delta^2 ) . \end{align}

From the expressions for the velocity (3.21) and pressure (3.24), the relation between the coefficients in (3.28), and the dynamic conditions (3.30), it then follows that

(3.32a) \begin{align} &- A \eta_\textrm{D} [ 2 m_\textrm{D} k {R_0}^{{-}1} {\rm J}_{k-1}(m_\textrm{D} {R_0}) - 2 (k+1) k {R_0}^{{-}2} {\rm J}_k (m_\textrm{D} {R_0})]\nonumber\\ &\qquad+ B \eta_\textrm{A} [ 2 m_\textrm{A} k {R_0}^{{-}1} H_{k-1}^{(2)}(m_\textrm{A} {R_0}) - 2 (k+1) k {R_0}^{{-}2} H_k ^{(2)}(m_\textrm{A} {R_0}) ]\nonumber\\ &\qquad+ E\eta_\textrm{D} [ m_\textrm{D}^2 {R_0}^k - 2 (k-1) k {R_0}^{k-2} ]+ F \eta_\textrm{A} [ m_\textrm{A}^2 {R_0}^{{-}k} - 2 (k+1) k {R_0}^{{-}k-2} ]\nonumber\\ &\quad= \sigma_{\textrm{D}\textrm{A}} {R_0}^{{-}1} (k^2 - 1) , \end{align}

(3.32b) \begin{align} &-A \eta_\textrm{D} [ 2 m_\textrm{D} {R_0}^{{-}1} {\rm J}_{k-1}(m_\textrm{D} {R_0}) + ( m_{\textrm{D}}^2 - 2 (k+1) k {R_0}^{{-}2} ) {\rm J}_k (m_\textrm{D} {R_0}) ]\nonumber\\ &\qquad+ B \eta_\textrm{A} [ 2 m_\textrm{A} {R_0}^{{-}1} H_{k-1}^{(2)}(m_\textrm{A} {R_0}) + ( m_\textrm{A}^2 - 2 (k+1) k {R_0}^{{-}2} ) H_k ^{(2)}(m_\textrm{A} {R_0}) ]\nonumber\\ &\qquad+ 2 E \eta_\textrm{D} (k-1) k {R_0}^{k-2} + 2 F \eta_\textrm{A} (k+1) k {R_0}^{{-}k-2}\nonumber\\ &\quad= 0 . \end{align}

3.4. Dispersion relation

For each wavenumber $k\in \mathbb {N}_{\geqslant {}2}$, the corresponding characteristic temporal response coefficient of the solution, $\gamma \in \mathbb {C}$, as well as the mode shapes, encoded in the remaining free parameters, follow from a solution-existence condition. To elucidate this condition, we first recall that the general bounded complex-valued velocity and pressure solutions of the partial-differential equations (3.6a)–(3.6b), subject to the limit condition (3.2), are given by (3.21) and (3.24). These general solutions contain twelve coefficients. Eight of these coefficients are determined by the kinematic interface condition (3.6c), in accordance with (3.28). The kinematic conditions (3.6c) and (3.6d) imply that the remaining four coefficients, $A,B,E$ and $F$ must satisfy the three identities in (3.29). The dynamic condition (3.6e) demands that, in addition, these four coefficients satisfy the two identities in (3.32). The remaining five conditions on the coefficients can be cast in the form

(3.33)\begin{equation} \underbrace{ \begin{pmatrix} {a}_{11} & 0 & {a}_{13} & 0 \\ 0 & {a}_{22} & 0 & {a}_{24} \\ {a}_{31} & {a}_{32} & {a}_{33} & {a}_{34} \\ {a}_{41} & {a}_{42} & {a}_{43} & {a}_{44} \\ {a}_{51} & {a}_{52} & {a}_{53} & {a}_{54} \end{pmatrix}}_{\boldsymbol{\mathsf{A}}(k,\gamma)} \begin{pmatrix} A \\ B \\ E \\ F \end{pmatrix} = \underbrace{ \begin{pmatrix} -\gamma R_0 \\ -\gamma R_0 \\ 0 \\ \sigma_{\textrm{D}\textrm{A}} R^{{-}1}_0 (k^2 - 1) \\ 0 \end{pmatrix}}_{{\boldsymbol{b}}(k,\gamma)} \,, \end{equation}

in such a manner that the first three equations in (3.33) represent (3.29) and the latter two represent (3.32). Noting the dependence of (3.29) and (3.32) on the wavenumber $k$, and recalling the dependence of $m_{\textrm {D}}$ and $m_{\textrm {A}}$ in these equation on the temporal response coefficient $\gamma$ via (3.28a), we infer that the entries of ${\boldsymbol{\mathsf{A}}}$ depend on $k$ and $\gamma$. With five constraints and four unknowns, the system of equations (3.33) is formally over-constrained, and a solution is non-existent unless the right-hand side vector ${\boldsymbol {b}}(k,\gamma )$ is in the column space of ${\boldsymbol{\mathsf{A}}}(k,\gamma )$. The relation between the existence of a solution and the condition

(3.34)\begin{equation} {\boldsymbol{b}}(k,\gamma)\in\operatorname{span}({\rm{col}}({\boldsymbol{\mathsf{A}}}(k,\gamma))) , \end{equation}

is indicative of the fact that only specific combinations of the wavenumber $k\in \mathbb {N}_{\geqslant {}2}$ and the temporal response coefficient $\gamma \in \mathbb {C}$ in the postulated interface configuration (3.3) correspond to a natural response of the droplet.

To determine the combinations $(k,\gamma )$ for which the existence condition (3.34) is fulfilled, we note that (3.34) is equivalent to

(3.35)\begin{equation} \det\left((\boldsymbol{\mathsf{A}}\mid\boldsymbol{b})(k,\gamma)\right)=0 , \end{equation}

where $(\boldsymbol{\mathsf{A}}\mid \boldsymbol {b})$ corresponds to $\boldsymbol{\mathsf{A}}$ augmented by $\boldsymbol {b}$. The equivalence between (3.34) and (3.35) follows from the fact that the column vectors of $\boldsymbol{\mathsf{A}}$ are linearly independent for all $(k,\gamma )$ and, hence, the augmented matrix is singular if and only if the vector $\boldsymbol {b}$ resides in the column space of $\boldsymbol{\mathsf{A}}$. Moreover, by virtue of the linear independence of the columns of $\boldsymbol{\mathsf{A}}$, if (3.35) holds, then (3.33) has a unique solution. This solution corresponds to the coefficients $(A,B,E,F)$ that, in combination with (3.28), define the droplet and ambient velocity-pressure pairs corresponding to $(k,\gamma )$ according to (3.21) and (3.24).

To facilitate and generalize the root finding of the determinant in (3.35), we non-dimensionalize the matrix entries of the augmented matrix based on the droplet density, $\rho _{\textrm {D}}$, droplet viscosity, $\eta _{\textrm {D}}$, and droplet radius $R_0$. The non-dimensionalized parameters are indicated with a tilde diacritic. Additionally, we introduce the following condensed notation:

(3.36)\begin{equation} \left.\begin{gathered} \mathcal{J} = {\rm J}_k ( \tilde{m}_\textrm{D} (\tilde{\gamma}) ),\quad \mathcal{H} = H^{(2)}_k ( \tilde{m}_\textrm{A} (\tilde{\gamma}) ) ,\quad \zeta = 2(k-1)k , \\ \hat{\mathcal{J}} = \tilde{m}_\textrm{D} (\tilde{\gamma}) {\rm J}_{k-1} ( \tilde{m}_\textrm{D} (\tilde{\gamma}) ) ,\quad \hat{\mathcal{H}} = \tilde{m}_\textrm{A} (\tilde{\gamma}) H^{(2)}_{k-1} ( \tilde{m}_\textrm{A} (\tilde{\gamma}) ) ,\quad \xi = 2 (k+1) k . \end{gathered}\right\} \end{equation}

The non-dimensionalized augmented matrix can then be expressed as

(3.37) \begin{align} & (\tilde{\boldsymbol{\mathsf{A}}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma})\nonumber\\ &\quad = \begin{pmatrix} k \mathcal{J} & 0 & k & 0 & - \tilde{\gamma} \\ 0 & k\mathcal{H} & 0 & k & - \tilde{\gamma} \\ - \hat{\mathcal{J}} + k \mathcal{J} & \hat{\mathcal{H}} - k \mathcal{H} & -k & -k & 0 \\ -2 k \hat{\mathcal{J}} + \xi \mathcal{J} & \tilde{\eta}_\textrm{A} \big[ 2 k \hat{\mathcal{H}} - \xi \mathcal{H} \big] & \tilde{m}^2_\textrm{D} - \zeta & \tilde{\eta}_\textrm{A} \big [ \tilde{m}^2_\textrm{A} - \xi \big] & \tilde{\sigma}_{\textrm{D}\textrm{A}}\left(k^2 - 1\right) \\ -2\hat{\mathcal{J}} - \left(\tilde{m}^2_\textrm{D} - \xi \right) \mathcal{J} & \tilde{\eta}_\textrm{A} \big[ 2 \hat{\mathcal{H}} + \left(\tilde{m}^2_\textrm{A} - \xi \right) \mathcal{H}\big] & \zeta & -\tilde{\eta}_\textrm{A} \xi & 0 \end{pmatrix}. \end{align}

It is not generally feasible to determine the roots of $\det ((\tilde {\boldsymbol{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$ with respect to $\tilde {\gamma }$ in closed form and, in practice, it is necessary to revert to a numerical root-finding algorithm. Once a root has been determined, one can extract the kernel of the augmented matrix (3.37) and scale the corresponding vector such that its fifth entry is minus one, to obtain the coefficients $\tilde {A},\tilde {B},\tilde {E},\tilde {F}$.

Remark 3.3 The roots of $\det ((\tilde {\boldsymbol{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$ are not unique: one can infer that

(3.38) \begin{equation} [(\tilde{\boldsymbol{\mathsf{A}}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}{}^*)] = [(\tilde{\boldsymbol{\mathsf{A}}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma})]^*, \end{equation}

where $({\cdot })^*$ denotes complex conjugation. Because the eigenvalues of the complex conjugate of a matrix are the complex conjugates of the original eigenvalues, it follows that if $\tilde {\gamma }$ is a root of $\det ((\tilde {\boldsymbol{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$, then so is $\tilde {\gamma }{}^*$. Since, in addition, it must hold that $\textrm {Re}(\tilde {\gamma })>0$, it suffices to consider roots in the fourth quadrant of the complex plane. Noting that the entries of the augmented matrix are analytic functions, one can infer that so is its determinant. This implies that the roots of $\det ((\tilde {\boldsymbol{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$ form a totally disconnected set and, accordingly, for each root there exists a neighbourhood in which that root is unique. A detailed investigation of the uniqueness of the roots of $\det ((\tilde {\boldsymbol{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$ in the fourth quadrant is beyond the scope of this work. In our numerical root-finding procedure, we have verified that there are no other roots in a region around the found root.

By virtue of the complex representation of the interface parametrization (3.4), we obtain the real-valued velocity and pressure fields by taking the real parts of (3.21) and (3.24) after substitution of the relations (3.28), and the temporal response coefficient $\gamma$ and the corresponding coefficients $A$, $B$, $E$ and $F$. Table 1 provides computed parameter values for the physical setting outlined in table 2, representing a water-in-air picolitre-sized droplet. For completeness, we mention that we have applied Mathematica's root-finder to determine $\tilde {\gamma }:=\tilde {\gamma }_k$ in the fourth quadrant of the complex plane such that $\det ((\tilde {\boldsymbol{\mathsf{A}}}| \tilde {\boldsymbol {b}} )(k,\tilde {\gamma }_k))=0$. The dimensions of parameters $E$ and $F$ depend on the mode number $k$. As a result, as $k$ increases, the values of $E$ and $F$ grow rapidly, conveying that these parameters are ill conditioned in terms of $k$. Furthermore, as the mode number $k$ increases, the frequency and damping rate of the corresponding oscillation, both encoded in $\gamma$, increase. This implies that if a droplet sustains an initial perturbation that is characterized by multiple modes, the higher wavenumber modes decay quickly, and low-order modes dominate the long-term dynamics of viscous-in-viscous oscillating droplets.

Table 1. Modal solution parameter values for a water droplet of radius $R= \sqrt {2}\times 10^{1}\ \mathrm {\mu } \textrm {m}$ suspended in air.

Table 2. Physical and numerical parameter values of the considered numerical experiments. Entries marked with the symbol $\ast$ indicate a range of values, which will be specified in the text.

4. Numerical experiments

The free-boundary problem (2.5) formally represents the sharp-interface limit of the Abels–Garcke–Grün NSCH model (2.1), provided that the mobility is appropriately scaled in the limit $\varepsilon \to {}+0$. For sufficiently small $\delta$, the oscillating-droplet solutions derived in § 3 can therefore serve to investigate the approach of the diffuse-interface solution to the sharp-interface limit solution. In this section, we investigate this sharp-interface limit numerically, by means of an adaptive finite-element method. Specifically, we focus on the scaling of the mobility parameter $m:=m_{\varepsilon }$ in the limit $\varepsilon \to {}+0$, and investigate the deviation of the diffuse-interface solution from the sharp-interface solution in relation to $m$. As reference solutions, we consider the lowest mode of oscillation ($k=2$), as well as the next higher doubly symmetric mode ($k=4$), for a viscous droplet in a viscous ambient with parameter values according to table 2. The corresponding coefficients of the velocity solution (3.21) and pressure solution (3.24) fields are presented in table 1.

4.1. Set-up and discretization

The oscillating-droplet test cases that we consider pertain to doubly symmetric modes. The set-up of the test cases is similar to that in Demont et al. (2022, § 5). To reduce computational expense, we exploit the symmetry of the configurations and consider only one quarter of the droplet-ambient domain. We regard a domain $\varOmega =(0,50)^2\ \mathrm {\mu }\textrm {m}^2$ and prescribe symmetry conditions on $\varGamma _{{sym}}:=\{(x_1,x_2)\in \partial \varOmega :\{x_1=0\}\cup \{x_2=0\}\}$. Because the linear sharp-interface solution is in fact defined on the generating circular droplet domain and the corresponding ambient domain according to (3.1a), while the diffuse-interface model exhibits a moving interface, we prescribe auxiliary conditions in accordance with an initially circular droplet. Specifically, with reference to (3.4), we select $t_0$ such that $-\textrm {Im}(\gamma )t_0={\rm \pi} /2$ and, hence, $R_{\delta }(\theta,t_0)=R_0$, and prescribe initial data corresponding to the reference solution at $t_0$ and boundary data corresponding to $t+t_0$. The complementary part of the boundary, $\varGamma _{{ext}}:=\partial \varOmega \setminus \varGamma _{{sym}}$, is furnished with Dirichlet conditions for velocity and homogeneous Neumann conditions for the order parameter and the chemical potential

(4.1) \begin{equation} \left.\begin{array}{c@{}} \boldsymbol{u}({\cdot},t)=\delta\boldsymbol{u}_{\textrm{A}}({\cdot},t_0+t)\\ \partial_n\varphi=0\\ \partial_n\mu=0 \end{array}\right\}\quad\text{on}\ \varGamma_{{ext}}, \text{ for }t\in[0,T), \end{equation}

where $\delta \boldsymbol {u}_{\textrm {A}}$ corresponds to the ambient velocity solution (3.21) with appropriate coefficients and scaling $\delta$, and $[0,T)$ denotes the time interval under consideration. For the aforementioned combination of boundary conditions, the pressure variable $p$ is only determined up to a constant. We impose the auxiliary condition that $p$ vanishes on average.

We impose an initial condition for the order parameter corresponding to a circular interface, in accordance with the initial configuration of the sharp-interface reference solution, viz.

(4.2)\begin{equation} \varphi(\boldsymbol{x},0)=\varphi_0(\boldsymbol{x}):=\tanh\left(\frac{d_{{\pm}}(\boldsymbol{x},\varGamma_0)}{\sqrt{2}\varepsilon}\right), \end{equation}

where $d_{\pm }(\boldsymbol {x},\varGamma _0)$ represents the signed distance from $\boldsymbol {x}$ to $\varGamma _0$. The function $s\mapsto \tanh (s/\sqrt {2}\varepsilon )$ corresponds to an equilibrium solution of the Cahn–Hilliard equations for the phase field in one spatial dimension and, accordingly, the phase field (4.2) is meta-stable if $\varepsilon$ is sufficiently small compared with the radius of curvature of $\varGamma _0$. In conjunction with (4.2), we impose the following initial condition for velocity:

(4.3)\begin{equation} \boldsymbol{u}(\boldsymbol{x},0)= \begin{cases} \delta\boldsymbol{u}_{\textrm{D}}(\boldsymbol{x},t_0) & \text{if }\boldsymbol{x}\in\varOmega_{\textrm{D},0},\\ \delta\boldsymbol{u}_{\textrm{A}}(\boldsymbol{x},t_0) & \text{if }\boldsymbol{x}\in\varOmega_{\textrm{A},0}, \end{cases} \end{equation}

where the data in the right member of (4.3) correspond to the velocity solutions according to (3.21) in the droplet and ambient domains. We select the perturbation magnitude $\delta =10^{-2}$, after verifying that this choice renders the linearization error negligible in comparison with the deviation between the diffuse-interface and the sharp-interface solutions, for the $\varepsilon$ considered below. Hence, the selected value of $\delta$ is suitable for our investigation of the sharp-interface limit. The characteristic parameters pertaining to the droplet and ambient fluids, and to the interface are reported in table 2.

To perform the numerical simulations, we make use of the adaptive finite-element approximation method presented in Demont et al. (2022). For coherence, we present a concise overview of the numerical methodology. The weak form of the NSCH equations (2.1) is discretized with respect to the spatial dependence with $P3-P2$ (Taylor–Hood) $\mathcal {C}^0$ truncated hierarchical B-splines (see Hughes, Cottrell & Bazilevs 2005; Cottrell, Hughes & Bazilevs 2009; Giannelli, Jüttler & Speleers 2012) for the velocity and pressure fields, and $P3$ $\mathcal {C}^0$ truncated hierarchical B-splines for the order parameter and chemical potential; see van Brummelen et al. (2021, § 3.1) for further details. The adaptive-refinement procedure is guided by a two-level hierarchical a posteriori error estimate, and follows the standard SEMR (Solve $\rightarrow$ Estimate $\rightarrow$ Mark $\rightarrow$ Refine) process (Dörfler 1996; Bertoluzza et al. 2012). To improve the robustness of the solution procedure on the coarse meshes that occur in the sequence of adaptive refinements within each time step, an $\varepsilon$-continuation process is introduced, in which the thickness parameter $\varepsilon$ (and, in conjunction, the mobility $m \propto \varepsilon ^3$) is enlarged for the first $K$ iterations of the adaptive-refinement process; see Demont et al. (2022) and van Brummelen et al. (2021) for details. In each time step, the fluid domain is initially covered with a uniform mesh comprising $10\times 10$ elements, corresponding to an initial mesh width $h_0=5\ \mathrm {\mu }\textrm {m}$, and we perform $L_{max}$ refinement steps. Refinement steps $L=0, 1, \ldots, K-1$ make use of the $\varepsilon$- and $m$-continuation process, while in refinement steps $L=K, \ldots, L_{max}$ the original parameter values for $\varepsilon$ and $m$ are used. A skew-symmetric formulation according to Layton (2008) is used for the convective term in the Navier–Stokes equations, enhancing the stability of the discrete approximation by eliminating potential artificial energy production due to deviations from solenoidality in pure-species regions. On the coarse meshes, a first-order backward Euler scheme with second-order contractive–expansive splitting of the double-well potential with stabilization (Wu, van Zwieten & van der Zee 2014) is employed. On the finest mesh, a second-order Crank–Nicolson scheme is applied with implicit treatment of the double-well potential. The second-order Crank–Nicolson scheme provides significant better accuracy than the backward Euler scheme; cf. e.g. John, Matthies & Rang (2006). For the temporal discretization, we employ a time-step size $\tau = 2^{-7} \times 10\ \mathrm {\mu }\textrm {s}$. The parameter setting of the numerical procedure is also summarized in table 2. The nonlinear algebraic systems corresponding to the discretized NSCH equations, are solved with a Newton procedure, in which the linear tangent problems are solved with the generalized minimal residual method (GMRES) with a preconditioner based on a partition of the NSCH system into NS and CH subsystems; see Demont et al. (2022) for further details.

To illustrate the set-up of the numerical experiments, and the resemblance between the analytic sharp-interface solution and the numerical approximation of the diffuse-interface solution for sufficiently small $\varepsilon$, we conduct numerical experiments with interface-thickness parameter $\varepsilon =2^{-10}\times 10^2\ \mathrm {\mu } \textrm {m}$, mobility $m=9.5272\times 10^{-13} \textrm {m}^2\textrm {s}\ \textrm {kg}^{-1}$, maximum number of refinement levels $L_{max}=7$ and number of continuation levels $K=5$. Figures 2 and 3 display snapshots of the velocity field and pressure field at six time instants. The top (respectively bottom) half of each panel displays the velocity (respectively pressure) field. The right (respectively left) half of each panel depicts the diffuse-interface simulation (respectively sharp-interface solution). The figures convey that the sharp-interface solutions and the diffuse-interface solutions are visually indistinguishable.

Figure 2. Snapshots of the magnitude of the velocity field $|\boldsymbol {u}|$ (top) and pressure field $p$ (bottom) throughout half a droplet oscillation of mode $k=2$. The left half of each panel displays the analytical sharp-interface solution, while the right half displays the numerical diffuse-interface solution.

Figure 3. Snapshots of the magnitude of the velocity field $|\boldsymbol {u}|$ (top) and pressure field $p$ (bottom) throughout half a droplet oscillation of mode $k=4$. The left half of each panel displays the analytical sharp-interface solution, while the right half displays the numerical diffuse-interface solution.

4.2. Optimal mobility scaling

To elucidate the dependence of the diffuse-interface solution in the sharp-interface limit $\varepsilon \to +0$ on the scaling of the mobility $m:=m_{\varepsilon }$, we conduct numerical experiments for a range of combinations of $\varepsilon$ and $m$. For each combination of $\varepsilon$ and $m$, we determine the deviation relative to the sharp-interface solution according to

(4.4)

where the considered length of the time interval, $T$, corresponds to half a period of oscillation. We regard a set of decreasing interface-thickness parameters $\varepsilon \in \mathscr {E}:=\{2^0,\ldots,2^{-3}\}\varepsilon _{max}$ relative to the baseline interface thickness $\varepsilon _{max} = 2^{-7}\times 10^{2}\ \mathrm {\mu } \textrm {m} = 7.8125\times 10^{-1}\ \mathrm {\mu } \textrm {m}$. The baseline interface-thickness parameter corresponds to approximately 5 % of the droplet radius. For each $\varepsilon$, we consider mobility parameters in (a relevant subset of) the set $m\in \mathscr {M}:=\{2^{0},2^{-1},\ldots,2^{-12}\}m_{max}$ with $m_{{max}} = 2.4389632\times 10^{-10}\ \textrm {m}^d \textrm {s}\ \textrm {kg}^{-1}$. The range of mobility parameters has been determined empirically such that $[2^{-12},1]m_{{max}}$ includes the optimal mobility, i.e. the one for which $\operatorname {dev}(\varepsilon,m)$ is minimal, for all $\varepsilon \in \mathscr {E}$. It is to be noted that $\mathcal {E}\times \mathcal {M}$ contains various monomial scalings of the mobility with respect to the interface thickness, viz. $m\propto \varepsilon ^l$ with $l\in \{0,1,2,3\}$.

Tables 3 and 4 present the deviations $\operatorname {dev}(\varepsilon,m)$ for the two modes of oscillation, $k=2$ and $k=4$, respectively. For each $\varepsilon \in \mathcal {E}$, the entry corresponding to the mobility $m\in \mathcal {M}$ that yields the smallest deviation, is highlighted. One can observe that, indeed, the mobility corresponding to the minimal deviation decreases as $\varepsilon$ decreases. More precisely, for both modes, the optimal scaling of the mobility parameter with the interface-thickness parameter appears to lie between $m\propto \varepsilon$ and $m\propto \varepsilon ^2$. One may moreover note that the entries corresponding to the optimal mobility decrease by a factor of approximately two if $\varepsilon$ is halved, which indicates that for the considered droplet-oscillation case, the diffuse-interface solution approaches the sharp-interface solution at rate ${O}(\varepsilon )$, provided that the mobility in the diffuse-interface model is appropriately scaled.

Table 3. Deviation between the diffuse-interface solution and the sharp-interface solution according to (4.4), for $k=2$ for a half-period of oscillation.

Table 4. Deviation between the diffuse-interface solution and the sharp-interface solution according to (4.4), for $k=4$ for a half-period of oscillation.

To provide a more precise assessment of the optimal scaling relation $m:=m_{\varepsilon }$, we determine for each $\varepsilon$ the optimal value of $m$ based on a quadratic log–log interpolation around the minimal values in tables 3 and 4. Figure 4 plots the optimal value of $m$ vs $\varepsilon$. For both wavenumbers, we observe an optimal scaling $m\propto \varepsilon ^{a_{opt}}$ with $a_{opt}\approx 1.7$. It is noteworthy that the constant of proportionality in the scaling relation is different for the two modes, and that the graphs are offset in the $\varepsilon$-dependence by a factor of approximately two, i.e. the optimal mobility for $k=4$ is approximately $2^{a_{opt}}$ larger than the optimal mobility for $k=2$. This suggests that the optimal mobility in fact scales with $(\varepsilon /\ell )^{a_{opt}}$, where $\ell$ represents another characteristic length scale of the interface which, for the considered droplet-oscillation test case, is proportional to the wavelength of the perturbation. This observation calls for further investigation, but we consider a detailed analysis of this aspect beyond the scope of the present work.

Figure 4. Optimal mobility $m$ obtained from quadratic interpolation around the minima in tables 3 and 4.

4.3. Sensitivity to the proportionality constant

In the previous section, we established optimal values of the mobility parameter and inferred an optimal scaling $m\propto \varepsilon ^{a_{opt}}$ in the sharp-interface limit. The results in § 4.2 also convey that the constant of proportionality in the scaling relation $m_{\varepsilon }=\mathscr {C}\varepsilon ^{a_{opt}}$ depends on the configuration and dynamics of the interface. This raises the question how sensitive the solution is to suboptimality of the proportionality constant in the scaling relation.

To elucidate the sensitivity of the deviation of the diffuse-interface solution to the sharp-interface solution with respect to the mobility in the limit $\varepsilon \to {}+0$, figure 5 (respectively figure 6) plots for each $\varepsilon \in \mathcal {E}$ the ratio of the deviation $\operatorname {dev}(\varepsilon,m)$ in the columns of able 3 (respectively table 4) to the minimal deviation $\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })$ vs the ratio $m/m_{opt,\varepsilon }$. Noting that the curves in figures 5 and 6 exhibit a vanishing slope near $m/m_{opt,\varepsilon }=1$, one can conclude that in the vicinity of the optimal mobility, the relative deviation is essentially independent of the mobility. However, for larger departures from the optimal mobility and sufficiently small $\varepsilon$, the relative deviation increases almost linearly in $\max (m/m_{opt,\varepsilon },m_{opt,\varepsilon }/m)$. For the largest $\varepsilon \in \mathcal {E}$, the relative deviation appears to be less sensitive to underestimation than to overestimation of the mobility. However, the results plotted with solid markers in figures 5 and 6 indicate that in the sharp-interface limit, the relative deviation is equally sensitive to under- and overestimation of the mobility.

Figure 5. Normalized deviation $\operatorname {dev}(\varepsilon,m)/\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })$ vs normalized mobility $m/m_{opt,\varepsilon }$ for $k=2$.

Figure 6. Normalized deviation $\operatorname {dev}(\varepsilon,m)/\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })$ vs normalized mobility $m/m_{opt,\varepsilon }$ for $k=4$.

4.4. Suboptimal mobility scaling and convergence in the sharp-interface limit

To illustrate the effect of the scaling of the mobility on the approach to the sharp-interface limit solution, figures 7 and 8 plot the deviation $\operatorname {dev}(\varepsilon,(\varepsilon /\varepsilon _{max})^a \, m_0)$ vs $\varepsilon$ for $a\in \{0,\ldots,3\}$. Herein, $m_0$ corresponds to the optimal sampled mobility for $\varepsilon _{max}$, viz. $m_0=2^{-3}m_{max}$ for $k=2$ and $m_0=2^{-2}m_{max}$ for $k=4$; cf. tables 3 and 4. For reference, the figures also contain the estimated minimal deviation obtained by minimization of the quadratic interpolation, $\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })$. The figures convey that, for the optimal scaling of the mobility, the diffuse-interface solution approaches the sharp-interface solution essentially at order $\varepsilon$, i.e. $\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })={O}(\varepsilon )$ as $\varepsilon \to {}+0$. For the linear and quadratic scaling of the mobility with the interface thickness, $m\propto \varepsilon ^a$ with $a\in \{1,2\}$, i.e. the integer scalings of the mobility closest to the optimum, we observe convergence to the sharp-interface solution, but at a suboptimal (sublinear) rate. For the constant and cubic scalings of the mobility, $m \propto \varepsilon ^a$ with $a\in \{0,3\}$, the deviation $\operatorname {dev}(\varepsilon,m_0(\varepsilon /\varepsilon _{max})^a)$ does not vanish as $\varepsilon \to {}+0$, i.e. the diffuse-interface solution does not convergence to the sharp-interface solution. For $m \propto \varepsilon ^0$, this confirms the known result that for constant mobility, the Abels–Garcke–Grün NSCH model converges to the non-classical sharp-interface Navier–Stokes/Mullins–Sekerka model; see Abels & Garcke (2018).

Remark 4.1 It is noteworthy that the subcubic scaling of the mobility, $m\propto \varepsilon ^a$ with $0<{}a<3$, which is necessary to converge to the classical sharp-interface solution in the limit $\varepsilon \to {}+0$, implies that the characteristic diffusive time scale $T_{{diff}}:=\varepsilon ^3/\sigma {}m$ associated with the diffuse interface, approaches zero in the sharp-interface limit. Consequently, for any $\varepsilon$-independent characteristic time scale $T_*$ in the problem under consideration, e.g. the period of oscillation of a droplet, it holds that the ratio $T_{{diff}}/T_*\to {}+0$ as $\varepsilon \to {}+0$. For numerical time-integration methods for the NSCH equations, it is therefore essential that such methods are robust in the limit $T_{{diff}}/\tau \to +0$ (with $\tau$ denoting the time-step size), to avoid excessive computational complexity in the sharp-interface limit.

Figure 7. Error convergence rates for optimal and suboptimal scalings of mobility $m$, for mode of oscillation $k=2$.

Figure 8. Error convergence rates for optimal and suboptimal scalings of mobility $m$, for mode of oscillation $k=4$.

5. Conclusions

Diffuse-interface binary-fluid models bear significant potential for describing complex phenomena in fluid mechanics, such as topological changes of the fluid–fluid interface and dynamic wetting, by virtue of their implicit representation of the interface. In the absence of topological changes of the interface, diffuse-interface models should reduce to corresponding classical sharp-interface models in the so-called sharp-interface limit, viz. if the interface-thickness parameter, $\varepsilon$, passes to zero. Contemporary understanding of the sharp-interface limit is, however, incomplete and, in particular, the scaling of the mobility parameter, $m$, with $\varepsilon$ as $\varepsilon \to {}+0$ is incompletely understood. In this article, we investigated the limit behaviour of the Abels–Garcke–Grün NSCH model for the classical case of an oscillating droplet in two dimensions, by means of an adaptive finite-element methodology.

To provide reference sharp-interface solutions, we derived new two-dimensional analytical expressions for the velocity and pressure fields for small-amplitude oscillations of a viscous droplet in a viscous ambient fluid with different densities and viscosities.

For mode numbers $k=2,4$ of the droplet oscillation, we compared the solutions of the NSCH model with the corresponding analytical solutions for a decreasing sequence of interface-thickness parameters and a suitably chosen sequence of mobility parameters. Based on an analysis of the deviation between the diffuse-interface solution and the sharp-interface solution, we inferred that $m_{opt,\varepsilon }\propto \varepsilon ^{a_{opt}}$ with ${a_{opt}}\approx 1.7$ corresponds to the optimal scaling of the mobility in the sharp-interface limit. We found that this optimal scaling is universal for $k=2$ and $k=4$. However, we also observed that for $k=4$ the factor of proportionality is approximately $2^{a_{opt}}$ larger than for $k=2$, suggesting that the optimal mobility in fact scales with the ratio of $\varepsilon$ to another characteristic length scale, proportional to the wavelength of the perturbation. However, as the wavelength is specific to the droplet-oscillation case, the question is what the underlying generic length scale is. The intrinsic dependence of the optimal mobility on the configuration and motion of the interface is non-trivial and warrants further investigation.

For the optimal scaling of the mobility parameter, we observed that the deviation between the diffuse-interface solution and its sharp-interface limit decreases according to ${O}(\varepsilon )$ in the limit $\varepsilon \to {}+0$. Our investigation of suboptimal integer scalings of the mobility conveyed that the sharp-interface limit is also attained for a linear and quadratic scaling of the mobility, but not for a constant or cubic scaling. For the linear and quadratic scaling, the approach of the diffuse-interface solution to the sharp-interface solution occurs at a suboptimal (sublinear) rate. The fact that the scaling of the mobility with $\varepsilon$ must be subcubic, implies that the characteristic diffusive time scale $\varepsilon ^3/\sigma {}m$ (with $\sigma$ as surface tension) passes to zero in the sharp-interface limit.

An important question pertains to generalization of the aforementioned results beyond the oscillating-droplet scenario for which they have been determined. A main conclusion that we expect to hold universally, is that the sharp-interface limit is attained for different scalings $m_{\varepsilon }\propto \varepsilon ^a$ of the mobility, and that the value of $a$ determines the rate of convergence to the sharp-interface solution. Based on the presented results for the oscillating-droplet case and consistent with other results in the literature, we conjecture that $0< a<3$ is necessary and sufficient. Whether the observed optimal scaling rate $a_{opt}\approx 1.7$ for the oscillating-droplet case extends to other scenarios including, specifically, three-dimensional settings and problems involving contact lines, cannot be convincingly asserted without further investigation.