Non-local model for surface tension in fluid-fluid simulations

TL;DR

This paper introduces a non-local integral model for surface tension in fluid-fluid simulations, capturing nanoscale effects and improving the accuracy of droplet behavior predictions compared to traditional models.

Contribution

The paper presents a novel non-local surface tension model that accounts for size-dependent effects and demonstrates its effectiveness through analytical and numerical comparisons.

Findings

The model agrees with the Young-Laplace law for large droplets.

Surface tension decreases with smaller droplet radii, matching experimental observations.

The non-local model accurately predicts droplet dynamics and coalescence behaviors.

Abstract

We propose a non-local model for surface tension obtained in the form of an integral of a molecular-force-like function with support $3.5\varepsilon$ added to the Navier-Stokes momentum conservation equation. We demonstrate analytically and numerically that with the non-local model interfaces with a radius of curvature larger than the support length behave macroscopically and microscopically, otherwise. For static droplets, the pressure difference $P_{\varepsilon, in} - P_{\varepsilon, out}$ satisfies the Young-Laplace law for droplet radius greater than $3.5\varepsilon$ and otherwise deviates from the Young-Laplace law. The latter indicates that the surface tension in the proposed model decreases with decreasing radius of curvature, which agrees with molecular dynamics and experimental studies of nanodroplets. Using the non-local model we perform numerical simulations of droplets under dynamic conditions, including a rising droplet, a droplet in shear flow, and two colliding droplets in shear flow, and compare results with a standard Navier-Stokes model subject to the Young-Laplace boundary condition at the fluid-fluid interface implemented via the Conservative Level Set (CLS) method. We find good agreement with existing numerical methods and analytical results for a rising macroscopic droplet and a droplet in a shear flow. For colliding droplets in shear flow, the non-local model converges (with respect to the grid size) to the correct behavior, including sliding, coalescence, and merging and breaking of two droplets depending on the capillary number. In contrast, we find that the results of the CLS model are highly grid-size dependent.

Tables

Table 1. Table 1: Parameters for benchmark test cases 1 and 2 from [ 50 ] .

Test case	${\rho }_1$	${\rho }_2$	${\mu }_1$	${\mu }_2$	$g$	$\sigma$
1	1000	100	10	1	0.98	24.5
2	1000	1	10	0.1	0.98	1.96

Table 2. Table 2: Comparison between benchmark values from the local models in [ 50 ] and the non-local model with Δ ​ x = 1 / 80 , 1 / 160 , Δ 𝑥 1 80 1 160 \Delta x=1/80,1/160, and 1 / 320 1 320 1/320 . The definitions for the center of mass y c subscript 𝑦 𝑐 y_{c} , rise velocity v c subscript 𝑣 𝑐 {v}_{c} , and circularity c ^ ^ 𝑐 \hat{c} are given in eqs. 22 – 24 .

	$\Delta x=\frac{1}{80}$	$\Delta x=\frac{1}{160}$	$\Delta x=\frac{1}{320}$	[50]
$y_c\left(t=3.0\right)$	1.0863	1.0824	1.0807	1.081 $\pm$ 0.001
$v_{c,max}$	0.2438	0.2426	0.2418	0.2419 $\pm$ 0.0002
$t_{v_{c,max}}$	0.9199	0.9030	0.9174	0.9263 $\pm$ 0.005
${\hat{c}}_{min}$	0.8961	0.8967	0.8968	0.9012 $\pm$ 0.0001
$t_{{\hat{c}}_{min}}$	1.9206	1.9314	1.9125	1.89 $\pm$ 0.015

Table 3. Table 3: Comparison between values from [ 50 ] and the non-local model with Δ ​ x = 1 / 80 Δ 𝑥 1 80 \Delta x=1/80 and 1 / 160 1 160 1/160 . The definitions for the center of mass y c subscript 𝑦 𝑐 y_{c} and the rise velocity v c subscript 𝑣 𝑐 {v}_{c} are given in eqs. 22 and 23 .

	$\Delta x=\frac{1}{80}$	$\Delta x=\frac{1}{160}$	[50]
$y_c\left(t=3.0\right)$	1.102	1.119	1.134 $\pm$ 0.009
$v_{c,max}$	0.250	0.250	0.252 $\pm$ 0.002
$t_{v_{c,max}}$	0.674	0.750	0.731$\pm$ 0.003

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} and L ∞ subscript 𝐿 L_{\infty} error in the velocity field generated by the parasitic currents calculated with the non-local model and the CLS method.

$\Delta x$	Non-local ${\Vert 𝐮\Vert }_2$	Non-local ${\Vert 𝐮\Vert }_{\infty }$	CLS ${\Vert 𝐮\Vert }_2$	CLS ${\Vert 𝐮\Vert }_{\infty }$
1/16	$2.35\times {10}^{-4}$	$3.74\times {10}^{-2}$	$2.31\times {10}^{-4}$	$3.93\times {10}^{-2}$
1/32	$1.98\times {10}^{-5}$	$1.17\times {10}^{-2}$	$7.18\times {10}^{-5}$	$3.40\times {10}^{-2}$
1/64	$3.94\times {10}^{-6}$	$5.93\times {10}^{-3}$	$3.49\times {10}^{-5}$	$5.39\times {10}^{-2}$
1/128	$3.14\times {10}^{-7}$	$1.26\times {10}^{-3}$	$1.78\times {10}^{-5}$	$9.37\times {10}^{-2}$
1/256	$5.13\times {10}^{-9}$	$3.82\times {10}^{-5}$	$3.08\times {10}^{-6}$	$4.83\times {10}^{-2}$

Equations

$$\nabla \cdot {\mathbf{u}}_i=0,\mathbf{x}\in {\Omega }_i(t),$$

$${\rho }_i\left(\frac{\partial {\mathbf{u}}_i}{\partial t}+{\mathbf{u}}_i\cdot \nabla {\mathbf{u}}_i\right)=-\nabla P_i+\nabla \cdot {\mathbf{\tau }}_i+{\rho }_i\mathbf{g},\mathbf{x}\in {\Omega }_i(t),$$

$$(P_{\alpha }-P_{\beta })\mathbf{n}=({\mathbf{\tau }}_{\alpha }-{\mathbf{\tau }}_{\beta })\cdot \mathbf{n}+\kappa {\sigma }_0\mathbf{n},\mathbf{x}\in {\Omega }_{\alpha }\cap {\Omega }_{\beta },$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla P+\nabla \cdot \mathbf{\tau }+\rho \mathbf{g}+\mathbf{F},\mathbf{x}\in \Omega ,$$

$$\mathbf{F}={\sigma }_0\kappa \hat{\mathbf{n}}$$

$$\phi(\mathbf{x})=\left\{\begin{array}[]{ll}0,&\mathbf{x}\in\Omega_{\alpha},\\ 1,&\mathbf{x}\in\Omega_{\beta}.\end{array}\right.$$

$$\mathbf{F}={\sigma }_0\left(-\nabla \cdot \frac{\nabla \varphi }{|\nabla \varphi |}\right)\nabla \varphi$$

$$\mathbf{F}=-{\int }_{\Omega }s(\mathbf{x},\mathbf{y})f_{\epsilon }(|\mathbf{x}-\mathbf{y}|)\frac{\mathbf{x}-\mathbf{y}}{|\mathbf{x}-\mathbf{y}|}d\mathbf{y},\mathbf{x}\in \Omega ,$$

$$s(\mathbf{x},\mathbf{y})=\left\{\begin{array}[]{ll}s_{\alpha\alpha},&\mathbf{x}\in\Omega_{\alpha},\;\mathbf{y}\in\Omega_{\alpha},\\ s_{\alpha\beta},&\mathbf{x}\in\Omega_{\alpha},\;\mathbf{y}\in\Omega_{\beta},\\ s_{\beta\beta},&\mathbf{x}\in\Omega_{\beta},\;\mathbf{y}\in\Omega_{\beta}.\\ \end{array}\right.$$

$$s_{\alpha \alpha }=s_{\beta \beta }=\frac{1}{2(1-10^{-k})}\frac{{\sigma }_0}{\lambda },$$

$$f_{\epsilon }(|\mathbf{x}-\mathbf{y}|)=|\mathbf{x}-\mathbf{y}|\left[-A{e}^{-\frac{|\mathbf{x}-\mathbf{y}{|}^2}{2{\epsilon }_0^2}}+e^{-\frac{|\mathbf{x}-\mathbf{y}{|}^2}{2{\epsilon }^2}}\right]$$

$$P(r)-P(\infty )=-4\pi (s_{aa}-s_{ab})\left({\epsilon }^{4}G(r,\epsilon )-A{\epsilon }_0^{4}G(r,{\epsilon }_0)\right)$$

$$G(r,\epsilon )$$

$$\frac{{\mathbf{u}}^{\ast }-{\mathbf{u}}^n}{\Delta t}={\nabla }_h\cdot {\mathbf{u}}^n{\mathbf{u}}^n+\mathbf{g}+\frac{1}{{\rho }^n}{\mathbf{F}}^n+\frac{1}{{\rho }^n}{\nabla }_h^2{\mathbf{u}}^n$$

$${\nabla }_h\cdot \left(\frac{1}{{\rho }^n}{\nabla }_h{p}^n\right)=\frac{1}{\Delta t}{\nabla }_h\cdot {\mathbf{u}}^{\ast }$$

$$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^{\ast }}{\Delta t}=\frac{-{\nabla }_h{p}^n}{{\rho }^n}$$

$${\mathbf{F}}^n(\mathbf{x})$$

$$=\sum _j{f}_{\epsilon }(|\mathbf{x}-\mathbf{y}|)\frac{\mathbf{x}-\mathbf{y}}{|\mathbf{x}-\mathbf{y}|}{V}_j\left(s_{\alpha \alpha }[{\varphi }^n(\mathbf{x}){\varphi }^n(\mathbf{y})+(1-{\varphi }^n(\mathbf{x}))(1-{\varphi }^n(\mathbf{y}))]+s_{\alpha \beta }[{\varphi }^n(\mathbf{x})(1-{\varphi }^n(\mathbf{y}))+(1-{\varphi }^n(\mathbf{x})){\varphi }^n(\mathbf{y})]\right),$$

$${\varphi }_t+\text{u}\cdot \nabla \varphi =0$$

$${\varphi }_{\tau }+\nabla \cdot (\varphi (1-\varphi )\mathbf{n})=\kappa \Delta \varphi$$

$${\mu }^{n+1}={\mu }_1+({\mu }_2-{\mu }_1){\varphi }^{n+1},$$

$${\rho }^{n+1}={\rho }_1+({\rho }_2-{\rho }_1){\varphi }^{n+1}.$$

$${\mathbf{x}}_c=(x_c(t),y_c(t))=\frac{{\int }_{{\Omega }_2(t)}\mathbf{x}d\mathbf{x}}{{\int }_{{\Omega }_2(t)}1d\mathbf{x}},$$

$${\mathbf{v}}_c=(u_c(t),v_c(t))=\frac{{\int }_{{\Omega }_2(t)}\mathbf{v}(\mathbf{x},t)d\mathbf{x}}{{\int }_{{\Omega }_2(t)}1d\mathbf{x}},$$

$$\hat{c}(t)=\frac{2\pi a_{circle}(t)}{P_b(t)},$$

$$D=\text{Ca}\frac{(19\lambda +16)}{(16\lambda +16)}$$

$$\theta =\frac{\pi }{4}-\text{Ca}\frac{(19\lambda +16)(2\lambda +3)}{80(\lambda +1)}.$$

Full text

††journal: Journal of Computational Physics

\floatsetup

[figure]style=plain,subcapbesideposition=top

1 Introduction

Since the work of Young, Laplace, and Gauss, multiphase flow at the continuum scale has been modeled almost exclusively by the Young-Laplace (YL) law, imposed as a boundary condition for the Navier-Stokes (NS) equations at the fluid-fluid interface. In this work, we propose an alternative continuum description of the multiphase flow in the form of an integral of a molecular-force-like function with support $\varepsilon$ added to the Navier-Stokes momentum conservation equation.

At the nanoscale, multiphase flow is traditionally described by the equations of molecular dynamics (MD). MD simulations show a sharp density and pressure drop across the fluid-fluid interface in a region of around ten nanometers [1, 2]. Outside of this region the pressure satisfies the YL law, i.e., the pressure difference across the interface is linearly proportional to the interface curvature.

Another molecular-scale feature of the fluid-fluid interface is the surface tension dependence on the curvature radius for curvature radii smaller than 10 nanometers [3]. As the curvature radius increases, surface tension asymptotically approaches its “macroscopic” value $\sigma_{0}$. Therefore, one can conclude that the NS-YL model is suitable for describing interfacial dynamics on scales larger than 10 nanometers given that the thermal fluctuations are properly accounted for [4]. Comparisons between continuum and MD simulations have shown disagreement below a droplet diameter of 36 nanometers, with drastic differences at 10 nanometers [5].

One example of a physical system that cannot be described by the YL law are nanobubbles. Bulk nanobubbles, also known as ultrafine bubbles, show extraordinarily long-term stability compared to that which would be expected from the Young-Laplace Law with macroscopic surface tension [6], up to three months in a laboratory setting [7]. This long-term stability leads to a number of applications, for example long term stability of bulk nanobubbles is an advantage when used as an ultrasound contrast agents over ultrasound contrast agents with larger bubbles, which have a limited half-life [8]. Additionally, nanobubbles have a higher surface area than macro-scale bubbles for the same volume fraction of bubbles, and change the physiochemical properties of the fluid-bubble system [9, 10]. Bulk nanobubbles show strong promise for biomedical applications in ultrasound [11, 12], radiofrequency ablation of tumors [13], and for mediating drug delivery [14], as well as water and waste treatment [15, 16, 17, 18, 19, 20], prevention of membrane and surface fouling and cleaning [21, 22, 23], froth flotation [24], and plant and animal growth [25, 26]. Numerical simulations of a large number of nanobubbles over large timescales are essential for these and other important applications.

On the molecular scale, surface tension results from the broken symmetry in the molecular interactions near the interface, i.e., the molecular forces acting between like molecules differ from forces acting between unlike molecules. In a similar manner, the molecular-like-forces generate surface tension in the non-local model.

Many interfacial phenomena are multiscale in nature. For example, modeling colliding macroscale droplets requires resolving a large range of relevant length scales and topology changes that occur at the droplet interfaces and a thin fluid film forming between two colliding droplets, and is another example where the YL law fails under certain capillary numbers [27]. The interfaces of droplets become locally flat as they approach each other, so the surface tension force due to the YL law becomes zero. The film of surrounding fluid forms between droplets and reaches $100-1000\AA$ before rupturing. This film drains as the droplets come closer to each other, eventually rupturing due to the inter-molecular van der Waals forces, which are absent in the macroscale models [27].

In traditional local front capturing Level Set and volume-of-fluid methods, droplet coalescence will occur if the droplets are separated by less than about one grid point, a phenomenon known as numerical (or artificial) coalescence [28]. Simulations of binary head-on collisions by [29] using the Level Set method were unable to capture the small deformation regime of bubble merging because they could not resolve the film drainage and rupturing. The Coupled Level-Set/Volume-of-Fluid (CLSVOF) method was proposed to overcome numerical coalescence, at the expense of prohibiting even physically-feasible coalescence [30]. In the simulations by [31], coalescence occurred at timing determined by experiments or by a van der Waals force with augmented range such that the length scale was large enough to be resolved by the simulations, implemented by adding a surface force at the interface. This study, and others, showed that the behavior of the droplets after coalescence is sensitive to the timing of the front rupture, and experimental data may not be available in all cases to determine the coalescence time [31, 32]. To resolve this challenge, sub-grid-scale (SGS) models have been proposed, where a semi-analytical model for the thin film dynamics is coupled with the local model at the droplet lengthscale to determine if and when the front will rupture [32, 33]. This results in a predictive method, removing the necessity of prescribing the time for the fronts to merge. SGS models were considered in the CLSVOF method by [28] and [34] and in the front tracking method by [35]. These models require coupling knowledge of the dynamics of the nano-scale gap with the macro-scale bubble dynamics as well as computationally expensive periodic searching through the domain to find interfaces close to collision [36]. Another approach by [27] included van der Waals forces in the momentum conservation equation in addition to the YL surface tension force. This approach was limited to the head-on collisions with equal sized droplets.

To demonstrate the multiscale nature of the non-local model, we present a semi-analytical steady-state solution for the fluid pressure across a fluid-fluid interface. This solution shows nanoscale behavior for the radius of curvature smaller than $3.5\varepsilon$ and macroscopic behavior (i.e., the solution follows the YL law) for the radius of curvature larger than $3.5\varepsilon$.

Using the non-local model, we perform numerical simulations of droplets under dynamic conditions, including a rising droplet, a droplet in a shear flow, and two colliding droplets in a shear flow, and compare results with standard NS model subject to the YL boundary condition at the fluid-fluid interface implemented via the Conservative Level Set (CLS) method. We find good agreement with the CLS method for a rising macroscopic droplet and a droplet in a shear flow. For a small capillary number (Ca = 0.24), we find that the non-local model predicts a microscopic droplet in shear flow to form “ears”. The CLS method predicts the same (macroscopic) shape of the droplet regardless of the droplet size.

Finally, for colliding droplets in shear flow we find that the non-local model converges (with respect to the grid size) to the correct behavior, including (depending on capillary number) sliding, coalescing, and temporary bridging (coalescing and then separating) of two droplets. On the other hand, in our simulations the CLS method results are highly grid-size dependent.

2 Non-local surface tension model

We consider flow of two incompressible Newtonian fluids, denoted $\alpha$ and $\beta$ in a fixed domain $\Omega=\Omega_{\alpha}(t)\cup\Omega_{\beta}(t)$. The macroscopic (hydrodynamic) model for two-phase flow includes the continuity equation for phase $i=\alpha$ or $\beta$

[TABLE]

and the momentum conservation equation

[TABLE]

subject to the YL boundary condition at the fluid-fluid interface

[TABLE]

where the subscript $i$ denotes phase $i$, $\rho_{i}$ is the density, $\mathbf{u}_{i}$ is the velocity, $P_{i}$ is the pressure, $\mathbf{g}$ is the gravitational acceleration, $\tau_{i}=[\mu_{i}(\nabla\mathbf{u}_{i}+\nabla\mathbf{u}_{i}^{T})]$ is the viscous stress tensor with the dynamic viscosity $\mu_{i}$, $\sigma_{0}$ is the macroscopic surface tension, and $\mathbf{n}$ is the normal vector.

To simplify numerical treatment of these equations, it is common to replace eqs. 2 and 3 with [37]

[TABLE]

where the surface force $\mathbf{F}$ is given by the YL law as

[TABLE]

where $\kappa$ is the interface curvature. We introduce the color function $\phi$:

[TABLE]

The color function is advected with the velocity $\mathbf{u}$ as $\frac{\partial\phi}{\partial t}+\mathbf{u}\cdot\nabla\phi=0.$ Then, from [37], the force

[TABLE]

gives the same total force as eq. 5, but spread over the interface width.

Here, we propose to replace the YL definition of the surface tension force from eq. 5 with the non-local model

[TABLE]

where $f_{\varepsilon}(|\mathbf{x}-\mathbf{y}|)$ is the force shape function and $s(\mathbf{x},\mathbf{y})$ is the force strength. We obtain eq. 8 as the continuous limit of the so-called pairwise surface tension force that is used in multiphase Smoothed Particle Hydrodynamics [38, 39], a fully Lagrangian particle method. The force strength is given by

[TABLE]

To ensure that $\sigma_{0}$ is positive, the coefficients must satisfy $s_{\alpha\alpha}+s_{\beta\beta}>2s_{\alpha\beta}$. For convenience, we take $s_{\alpha\alpha}=s_{\beta\beta}=10^{k}s_{\alpha\beta}$ with $k=3$. Then, $s_{\alpha\alpha}$ and $s_{\beta\beta}$ can be found as a function of $\sigma_{0}$:

[TABLE]

where $\lambda=\frac{1}{8}\pi\int\limits_{0}^{\infty}z^{4}f_{\varepsilon}(z)dz$ and $\lambda=\frac{1}{3}\int\limits_{0}^{\infty}z^{3}f_{\varepsilon}(z)dz$ in three and two spatial dimensions, respectively.

As in molecular dynamics, $f_{\varepsilon}(|\mathbf{x}-\mathbf{y}|)$ must be repulsive for small $|\mathbf{x}-\mathbf{y}|$, attractive for large $|\mathbf{x}-\mathbf{y}|$, and, for computational efficiency, should have compact support $h=O(\varepsilon)$ or become negligibly small for $|\mathbf{x}-\mathbf{y}|\geq h$. Several forms of $f_{\varepsilon}$ have been proposed [39]. Here, we use

[TABLE]

and $\lambda=\frac{1}{2}\sqrt{2\pi}(-A\varepsilon_{0}^{5}+\varepsilon^{5}).$

A fundamental difference between the non-local model in eq. 8 and the YL law is that the former has an internal length scale $\varepsilon$, while the latter does not have any internal length scale. Because of this, the YL law predicts the same behavior (for the same dimensionless numbers) regardless of the problem’s length scale. In the following, we obtain an analytical solution for pressure that demonstrates that the non-local model behaves “macroscopically” (follows the YL law) on the scale larger than $\varepsilon$ and “microscopically” (deviates from the YL law in a way consistent with molecular dynamics simulations of droplets), otherwise.

Under static conditions, eqs. 4 and 8 can be solved analytically for a circular interface separating two fluids in two dimensions. The pressure as a function of the distance $r$ from the droplet center can be found as [40]:

[TABLE]

where $G(r,\varepsilon)$ is given by eq. 13:

[TABLE]

The pressure in eq. 12 as a function of $r$ is plotted in Fig. 1$A$ for a droplet with the radius $a$. The pressure profile qualitatively agrees with behavior observed in molecular dynamics simulations [41, 1, 2, 42]. The pressure is zero far from the droplet, then at the interface the pressure drops and becomes negative, as seen in MD simulations. At the center of the droplet the pressure reaches a constant value. The pressure jump from inside to outside the interface is given by the YL law (eq. 5).

Fig. 1$B$ demonstrates that the pressure difference $P_{\varepsilon,in}-P_{\varepsilon,out}$ as a function of curvature agrees with the YL law (eq. 3) for droplet radius $a>3.5\varepsilon$ ($P_{\varepsilon,in}$ and $P_{\varepsilon,out}$ are pressures inside and outside of the droplet at the distance greater than $3.5\varepsilon$ from the interface). For smaller droplets with $a<3.5\varepsilon$, the pressure jump begins to deviate from the YL law. For “tiny” droplets with $a<0.3\varepsilon$, the pressure difference begins to decrease (the triangle symbols in Fig. 1$B$). This indicates the limit of incompressible treatment of small (nano) droplets.

3 Numerical implementation of the non-local model

Eq. 2 is discretized with a prediction-correction scheme based on the finite volume method presented in [43]. A temporary velocity is first calculated without the pressure,

[TABLE]

where $\mathbf{u}^{n}$ is the velocity field at time step $n$, $\mathbf{u}^{*}$ is the temporary velocity, $\mathbf{g}$ represents external forces such as gravity, and $\rho^{n}$ is the density field. The subscript $h$ denotes discrete finite volume operators. The pressure $p^{n}$ is then calculated so that $\mathbf{u}^{n+1}$ satisfies the incompressibility condition, $\nabla_{h}\cdot\mathbf{u}^{n+1}=0$, which gives

[TABLE]

for the pressure and

[TABLE]

for the updated velocity at time step $n+1$. We use a staggered mesh, with the velocity calculated on the mesh edges and pressure, density, viscosity, and the color function updated on the cell centers. The pressure $p^{n}$ is solved using a successive over relaxation (SOR) scheme [43]. Following [43], we set the time step to $\Delta t=\frac{\Delta x^{2}\min(\rho_{\alpha},\rho_{\beta})}{4.25\max(\mu_{\alpha},\mu_{\beta})}$.

We use a linear interpolation within elements and the color function at time step $n$, $\phi^{n}$, to discretize eq. 8:

[TABLE]

where $V_{j}$ is the size of element $j$. In the numerical simulations, $\varepsilon/\varepsilon_{0}=2$, $A=(\varepsilon/\varepsilon_{0})^{3}$, and $s_{\alpha\alpha}={s}_{\beta\beta}=10^{5}{s}_{\alpha\beta}$. Unless otherwise noted we take $\varepsilon=4\Delta x/3.5$ where $\Delta x$ is the grid resolution and $V_{j}=\Delta x^{2}$.

When advecting the color function, $\phi$, it is important to use an accurate method that is conservative and also preserves the sharpness of the front. As noted in [44], implementation of high order upwind schemes is difficult on a staggered mesh. We use the method presented in [45] and updated in [46] for a conservative advection scheme of the color function. To update the color function, we first solve

[TABLE]

using the conservative Total Variation Diminishing (TVD) method with Superbee limiter from [45] with a second-order Runge-Kutta scheme in time. As noted in [45], with eq. 18 alone the interface width and profile will not remain constant over the course of a simulation. To remedy this, we solve the compression-diffusion equation from [45] after eq. 18. The compression-diffusion equation is given by

[TABLE]

where $\mathbf{n}$ is the normal of the interface and $\kappa=\frac{1}{2}(\Delta x)^{1-d}$. $d$ is a parameter chosen so that $d\ll 1$ [45]. In the paper we take $d=0.$ The time parameter is denoted by $\tau$ to distinguish it from the simulation time step, and is chosen as $\Delta\tau=\frac{1}{2}(\Delta x)^{1+d}$. Eq. 19 is solved until the solution reaches steady state, defined by $\int|\phi^{n+1}-\phi^{n}|\leq\text{TOL}\cdot\Delta\tau$ for a given tolerance TOL, where $\phi^{n}$ denotes the solution after $n$ iterations. To update $\phi^{n+1}$ we use a second-order Runge-Kutta scheme. The intermediate compression step provides an artificial compression normal to the front interface to maintain a sharp front, with constant thickness proportional to $\kappa$, preventing diffusion as the front is advected. The added viscosity term, $\kappa\Delta\phi$, prevents discontinuities at the interface. The compression flux $(\phi(1-\phi)\mathbf{n})$ acts in regions where $0<\phi<1$ and provides compression normal to the interface, sharpening the interface profile.

Because the area each phase is defined as the area inside and outside of the $\phi=0.5$ contour, conservation of $\phi$ does not necessarily imply conservation of each phase. However, our results in sec. 4.2 indicate that the area loss is minimal. Once the color function is calculated, the density and viscosity are found directly by

[TABLE]

[TABLE]

The numerical algorithm described in this section is implemented in a parallel C++ code with OpenMP. In this paper, simulations with the non-local methods are compared with simulations calculated with the CLS method. To provide a direct comparison, the simulations using the CLS method are completed with the same code, with the surface tension force in eq. 3 replaced by eq. 5.

4 Numerical results

In this section, we validate the numerical implementation of the non-local model against analytical solutions, solutions obtained with the CLS [47, 48, 45, 49, 46] finite-volume discretization of the Navier-Stokes equations subject to the YL boundary condition, and published numerical simulations. In the following, we refer to CLS simulations of the Navier-Stokes equations subject to the YL boundary condition as a local model. Comparisons with the results of the local model will be used to illustrate advantages of the proposed here non-local model.

4.1 Static pressure in a droplet

We first test the code against analytical non-local and YL solutions for a static droplet with radius $a$. In Fig. 2$A$, the pressure across the center of the droplet is plotted for varying resolution $\Delta x$. As the resolution across the domain increases, $P_{\varepsilon,in}-P_{\varepsilon,out}$ converges to the value expected by the YL law, with an error of less than $1\%$ for $a/\Delta x=64$. In Fig. 2$B$, the resolution remains constant while the parameter $\varepsilon$ is varied relative to $a$. Here, $P_{\varepsilon,in}-P_{\varepsilon,out}$ is in excellent agreement with the YL law when $a/\varepsilon>3.5$.

4.2 Rising bubble

A set of quantitative benchmarks for a rising bubble was proposed in [50] for validation and comparison between numerical methods for interfacial flow. [50] considered a bubble rising due to buoyancy effects in a two-dimensional domain with $L_{x}=1$, $L_{y}=2$, no-slip boundary conditions at the top and bottom walls, and free slip conditions on the vertical walls. The droplet has an initial radius $a=0.25$ and is centered at $(x_{c},y_{c})=(0.5,0.5)$. Values of the viscosity, density, surface tension, and gravitational acceleration are given in table 1 for the two test cases, denoted Case 1 and Case 2. In [50], three numerical discretizations of local models were compared, two Eulerian Level Set finite-element codes (TP2D, denoted by Group 1, and FreeLIFE, denoted by Group 2) and an arbitrary Lagrangian-Eulerian moving grid method (MooNMD, denoted by Group 3.) The results of [50] have been compared with other numerical local models by many subsequent authors, including a conservative Level Set method with front sharpening [51], the volume of fluid method implemented in OpenFOAM® [52], a diffuse interface model [53], and a finite element based Level Set method [54]. In [51] and [52], the results from [50] were matched to within $2\%$ and $4\%$, respectively, while in [54], it was found that higher resolution in their model is needed to achieve the accuracy in [50].

We run the non-local model for Case 1 and Case 2 until the dimensionless time $t=3$ with the parameters given in table 1 and $\varepsilon=\frac{8}{3.5}\Delta x$, where $\Delta x=1/80$, $1/160$, or $1/320$ is the grid spacing. The characteristic length scale is $L=2a$, where $a$ is the bubble radius, and time scale $L/\sqrt{2ga}$ where $g$ is the gravitational velocity. In [50], data are provided for the bubble area $A(t)$, the center of mass $\mathbf{x}_{c}$

[TABLE]

where $\Omega_{2}(t)$ is the domain occupied by the bubble, the rise velocity $\mathbf{v}_{c}$

[TABLE]

where $\mathbf{v}=(u,v)$ is the fluid velocity, and the circularity $\hat{c}(t)$

[TABLE]

where $a_{circle}(t)=\sqrt{A(t)/\pi}$ and $P_{b}(t)$ is the bubble perimeter.

A comparison of the front locations given by [50] and from the non-local model at time $t=3.0$ for Case 1 and Case 2 is given in Figs. 3a and 4a. For Case 1 there is excellent agreement between the non-local model and the benchmark data from [50], both in terms of the final bubble location and shape. Case 2 is computationally more challenging due to the topological changes that occur as the bubble breaks up to form satellite droplets. The results from [50] disagree on the point of breakup and the final bubble shape, as shown in Fig. 3$c$.

Tables 2 and 3 give values for Case 1 and Case 2 of $y_{c}(t=3)$, the maximum rise velocity $v_{c,max}$ and the time $t_{v_{c,max}}$ at which the maximum occurs, and the minimum circularity $\hat{c}_{min}$ and the time $t_{\hat{c}_{min}}$ at which the minimum occurs for the range of values presented by the three groups in [50] and the non-local model. Our results show that the non-local model converges to the benchmark values as its resolution increases. For Case I is an excellent agreement between the non-local models and the results from [50], with a relative error for $\Delta x=1/320$ of $0.02\%$ for $y_{c}(t=3.0)$, $0.04\%$ for ${v}_{c,max}$, $0.96\%$ for $t_{{v}_{c,max}}$, $0.488\%$ for $\hat{c}_{min}$, and $1.19\%$ for $t_{\hat{c}_{min}}$. In each case the error is less than or on the same order of magnitude as the errors found by [51] and [52]. The inset of Fig. 3$d$ shows that the final center of mass in the vertical direction $y_{c}$ for the non-local model with $\Delta x=1/320$ is bounded by the benchmark values from [50]. For Case 2 the final bubble shape after breakup is inconclusive, however Table 3 shows good agreement between the non-local model and [50]. Due to computational limitations, we are unable to run higher resolution simulations of Case 2, which may improve the agreement between [50] and the non-local model.

While the numerical method used for advecting the color function $\phi$ is conservative, the bubble area is calculated as the contour corresponding to $\phi=0.5$, and the area inside this contour is not necessarily conserved. Nevertheless, we see good conservation of the bubble area to within $0.226\%$ for Case 1 with $\Delta x=1/80$, $0.077\%$ for $\Delta x=1/160$, and $0.027\%$ for $\Delta x=1/320$. The change in area for Case 2 is $0.283\%$ for $\Delta x=1/80$ and $0.116\%$ for $\Delta x=1/160$. The change in the bubble area is plotted in Figs. 3$b$ and 4$b$ for Case 1 and Case 2, respectively.

Figs. 3$e$ and 4$e$ show that the initial circularity is slightly under the expected value of one for a circular droplet at the beginning of the non-local simulations. This occurs because the droplet originally is fit to the square mesh. Over the first several time steps the compression algorithm for the advection of the color function rounds the edges of the droplet increasing the circularity to one.

4.3 Circular bubble in shear flow

A “macroscopic” spherical bubble in a linear shear flow $\mathbf{u}=(\dot{\gamma}y,0)$ will form a steady elongated ellipsoidal shape. Denoting the major and minor axes of the ellipse as $l$ and $b$, the deformation ratio $D=\frac{l-b}{l+b}$ of an initially spherical bubble with radius $a$ was found by [55, 56] using the local model as

[TABLE]

where $\lambda=\mu_{d}/\mu_{f}$ is the viscosity ratio between the drop and the surrounding fluid and $\text{Ca}=\mu_{f}a\dot{\gamma}/\sigma_{0}$ is the capillary number, the ratio of the magnitude of the viscous forces to the capillary forces. This formula has been validated for periodic suspensions of two-dimensional macroscopic droplets in a channel [57], and corrections have been made for macroscopic bubbles in confinement [58]. The three-dimensional results of [56] have been shown to agree well with two-dimensional simulations for small and moderate capillary number $\text{Ca}<1$ [57, 59]. The results of [55, 56] have been extended to second order [60, 61] to describe the angle the major axis of the deformed elliptical bubble makes with the horizontal axis, called the orientation angle and denoted by $\theta$, to give

[TABLE]

We simulate a macroscale bubble ($a/\varepsilon=5$) in shear flow with the non-local and local (CLS) models for $\text{Ca}=0.24$ and $\lambda=1.4$ to provide a comparison to the experiments of [59]. Results for the bubble deformation, $D$, and the location of the fronts are given in Fig. 5. At times $0\leq t\leq 4.0$, the non-local and local models agree well with the analytical results from [55, 56] and the experimental results [59]. The orientation angle $\theta$ agrees as well with the second-order analytical results from [60]. The relative error between the exact deformation, $D_{exact}=0.26625$, and the non-local model with $a/\varepsilon=5.0$ is $0.74\%$, and the relative error for the local model is $0.39\%$. The error between the analytical value for the orientation angle $\theta_{exact}$ and the non-local model is $4.13\%$, compared with $7.15\%$ for the CLS method.

Next, we model a “microscopic” bubble with the initial radius $a/\varepsilon=2.5$, $\text{Ca}=0.24$, and $\lambda=1.4$, the same Ca and $\lambda$ values as in the simulations of the macroscopic bubble. We observe that the deformation of the microscopic bubble is significantly larger than that of the macroscopic bubble and the resulting bubble has a pronounced sigmoid shape, see Fig. 5$C$. In this case, the curvature of the front is significantly lower at the rounded corners of the droplet, which in turn decreases the force due to the surface tension at those points. The droplet then continues to elongate, coming into equilibrium with a higher deformation than for macroscale droplets. Mesoscale Dissipative Particle Dynamics (DPD) simulations of a spherical microscale droplet in shear flow also predicted a larger droplet deformation for capillary number $\text{Ca}>0.2$ than that predicted by the Taylor theory [62, 63]. The deformation values found in these DPD studies correspond well with the deformation in the non-local model for microscale bubble. Note that the local model, the behavior depends only on Ca and $\lambda$ and not on the size of the bubble. We examined the droplet circulation for each case and found very similar results for all simulations in Fig. 5.

4.4 Droplet collisions in shear flow

In this section we use the non-local model and the local CLS model to simulate the collision of two equal size droplets in shear flow. This system has been studied previously, both experimentally and, more recently, numerically. Results have shown that the droplet behavior falls into one of three regimes: at low capillary number the droplets coalesce to form one large droplet, and at higher capillary number the droplets slide past each other [64, 21, 65]. In some cases, particularly with small droplets, a third regime is possible at moderate capillary number where the droplets temporarily coalesce before breaking apart [65]. These three regimes are illustrated in Fig. 6, simulated with the non-local model for $a=0.15$ at three capillary numbers. At the lowest capillary number, the droplets coalesce, forming one large droplet. At long times, this large coalesced droplet will find a final shape as described for a single droplet in shear flow in sec. 4.3. As the capillary number increases there is a transitional regime where the droplets coalesce and then separate, known as temporary bridging. Finally, at high enough capillary number the droplets slide past each other without coalescence. In [65], all three regimes were simulated using a three-dimensional free-energy lattice Boltzmann method (LBM). The LBM offers an advantage over many front capturing models, including the Volume of Fluid (VoF) method and CLS method, where coalescence of droplets will depend on the grid resolution. As the grid is further refined, the film between the droplets will be better resolved, which increases the time needed for the film to fully drain and delays coalescence. [65] shows the LBM does offer grid independent results, but the behavior depends on the thickness of the diffuse droplet interface relative to the droplet size.

Here, we choose parameters to match the results in [65], including the Reynolds number, $Re=1$. In comparison, the Reynolds number in the experiments in [21] is $Re<10^{-7}.$ Further studies are needed to provide a detailed understanding of how the coalescence behavior may depend on the Reynolds number. In our simulations, the droplets have initial separation in both the $x$ and $y$ directions given by $|x_{1}-x_{2}|/a=2.52$ and $|y_{1}-y_{2}|/a=0.78$, where $(x_{1},y_{1})$ and $(x_{2},y_{2})$ are the centers of the two droplets. The relative channel width is $H/a=5.1282$. The droplets have the same density and viscosity as the surrounding fluid, with $\rho=6\mu$. We find that the transitions between the three regimes occur at capillary numbers of the same order of magnitude as presented by [65], although we do not expect exact agreement due to the differences between two-dimensional and three-dimensional simulations.

A key advantage of the non-local model for this problem is that droplet coalescence is controlled by the parameter $\varepsilon$, not the grid resolution. In [65], the critical capillary number for the transition from coalescing to sliding droplets was found to decrease with the droplet radius, and follow a power-law with exponent that increases with increasing Péclet number. In our CLS simulations, the droplets coalesce at different times depending on the resolution and do not coalesce at all at the highest resolution considered with interface smoothing length $\kappa=\Delta x/2$, as shown in Fig. 7$A$ at $\text{Ca}=0.1$ and Fig. 8$A$ at $\text{Ca}=0.15$. This is because in the local model, the droplets will coalesce when they are within one grid point and there are no forces interacting between the droplets to assist coalesce when the droplets are separated by one more than one grid point. When the smoothing length is kept fixed, $\kappa=0.0125$, the droplets continue to coalesce at slightly different times, although the final behavior is the same, as shown in Fig. 7$B$ at $\text{Ca}=0.1$. In contrast, when the resolution is sufficient ($\Delta x\leq 1/80$), the behavior of the non-local method is similar regardless of resolution as long as $\varepsilon$ is kept constant with respect to the droplet radius as the mesh is refined. This is demonstrated in Fig. 7$C$ at $\text{Ca}=0.1$ and Fig. 8$B$ at $\text{Ca}=0.15$. At $\text{Ca}=0.1$, small satellite droplets form as two droplets coalesce. Similar phenomenon can be seen in [65].

The effect of $a/\varepsilon$ (the relative droplet size) is studied in Fig. 9, which shows the droplet trajectories for five values of $\varepsilon$ at resolution $\Delta x=1/160$ and droplet radius $a=0.15$. The resulting droplets range from micro-droplets ($a/\varepsilon=2.5$) to mesoscopic droplets ($a/\varepsilon=6,10$, and $12$) and macroscopic droplets ($a/\varepsilon=15$, and $20$). We find that the relatively smaller droplets (i.e., droplets with smaller $a/\varepsilon$) coalesce sooner. The droplets with $a/\varepsilon\geq 15$ do not coalesce during the simulation and instead slide past each other. The behavior of macroscopic bubbles does not depend on the $a/\varepsilon$ as long as $a/\varepsilon\geq 15$. The micro-scale droplets with $a/\varepsilon=2.5$ display distinct characteristics, although the overall behavior is similar to that of the mesoscopic droplets. The fronts between these droplets flatten significantly at $\dot{\gamma}t=1.2$, to a much larger degree than is seen in the other droplets. This results in many satellite droplets forming, which then migrate to the edge of the coalesced droplet over time. A zoomed-in view of the front locations and the corresponding color functions for macroscopic droplets can be seen in Fig. 10 for $a/\varepsilon=10$ and $6$, showing that the front flattening occurs to a lesser degree than the micro-scale droplet in Fig. 9 (see Fig. 10$B$.) The effect of increasing $\varepsilon$ relative to the droplet radius, for the same resolution and time, is to cause the two droplets to be significantly closer together leading to earlier coalescence.

An example of the coalescence process with the non-local simulation is shown in detail in Fig. 11, which plots the droplet fronts as well as the color function. In this example $\varepsilon/\Delta x=2.4$, and $a/\varepsilon=10$. When the droplets are sufficiently close an asymmetry forms in the front which becomes the initial location of coalescence. As the simulation continues, this neck widens until the droplets form a single large droplet.

Conclusion

In this paper we propose a non-local PDE model for multiphase flow that replaces the Young-Laplace law and is valid at both nano and macroscopic scales. Our model allows for simulations of multiscale systems with curvature radii ranging from nano to macro (micron and larger) scales, without the need to couple MD simulations with continuum NS models. Nanoscale multiphase flows, including nanobubbles, are important in biomedical applications, water and waste treatment, and membrane and surface defouling and cleaning. These applications typically include suspensions of many densely packed nanobubbles over macroscopic length scales, so the ability to use a continuum NS model will allow for efficient simulation of such systems in a way not possible with MD.

The non-local model limits numerical error from interpolating a sharp interface across a grid. Calculating the surface tension by integrating over a local neighborhood at each point on the surface reduces the spurious, or parasitic, currents that occur in other methods. It also handles merging interfaces without the behavior depending on the grid resolution. Instead, the dynamics are controlled by the relative curvature radius $a/\varepsilon$, where $\varepsilon$ is a parameter in the non-local surface tension force. The parameter $\varepsilon$ defines the relative width of the interface. Physically, the interface width is of the order of the support of the molecular forces. Therefore, when modeling nanodroplets, $\varepsilon$ should be on the order of the support of molecular forces, i.e., on the order of several $r{A}$. The parameter $\varepsilon$ also define the resolution of numerical simulations as the grid size should be smaller than $\varepsilon$. Therefore, using $\varepsilon$ on the order of several ångströms for modeling macroscopic droplets is computationally infeasible. Our results show that interfaces behave macroscopically when the radius of curvature larger than $3.5\varepsilon$ and the interface merging dynamics becomes independent of $\varepsilon$ for $a/\varepsilon>15$. Thus, it is sufficient to set $\varepsilon=a/3.5$ for modeling individual macroscopic droplets and $\varepsilon=a/15$ for modeling multiple interacting droplets.

When compared with existing benchmarks for macroscale droplets, the proposed non-local method matches or exceeds the accuracy of local numerical methods for multiphase flow. However, the non-local method can also capture mesoscale features that cannot be resolved by local methods, as demonstrated by a droplet in shear flow. Calculating the surface tension with the non-local model does require taking an integral at each point, which does make the surface tension calculation more computationally expensive, however each integral is independent and therefore lends itself well to massive parallelization. While not considered here, the non-local model can be easily implemented in existing codes by changing only the surface tension calculation. Additionally, the number of points included in each integral can be truncated to those within a distance of 3.5$\varepsilon/\Delta x$ from the target point, reducing the additional amount of work per time step to be $\mathcal{O}\left(N\varepsilon^{2}/\Delta x^{2}\right)$. For a fixed mesh, the force shape function $f_{\varepsilon}(|{\mathbf{x}}-{\mathbf{y}}|)$ can be precomputed before the simulation and used each time step. Our future work will include full three-dimensional simulations. Because the force due to surface tension only needs to be calculated near the interface, fast Level Set methods such as the one proposed by [66] could be adapted to the non-local model to reduce the computational time. In addition, simulations of multiscale systems are well fitted for adaptive mesh refinement, particularly when considering the interface between merging bubbles.

5 Acknowledgements

This work was supported by the U.S. Department of Energy (DOE) Office of Science, Office of Advanced Scientific Computing Research as part of the New Dimension Reduction Methods and Scalable Algorithms for Nonlinear Phenomena project. Pacific Northwest National Laboratory is operated by Battelle for the DOE under Contract DE-AC05-76RL01830.

The authors wish to thank Professor S. Guido for the data provided in Fig. 5.

Appendix A Spurious currents

Many numerical schemes for surface tension in two-phase flows struggle to capture equilibrium solutions exactly, resulting in what are known as parasitic or spurious currents, non-zero velocity fields when the system is in a static equilibrium and thus the velocity field should be identically zero. The origin of these currents comes from errors in approximating continuous quantities by discrete operators [67]. The parasitic currents may not converge with spatial resolution [68] and scale with the surface tension and viscosity [69]. Some numerical methods have been proposed to limit or eliminate parasitic currents, see for example [68, 70].

The simplest system for studying spurious currents is that of a circular droplet with $a>3.5\varepsilon$ suspended in another fluid in the absence of gravity. In this case, the velocity field is zero and the pressure jump satisfies eq. 3 exactly. We model a droplet at the equilibrium using the finite-volume (local) CLS and non-local models. The magnitude of the velocity fields are given in Table 4 and Fig. 12$A$. Because the expected velocity field is zero, any non-zero velocity is taken to be a spurious current. The magnitude and directions of the velocity field for both the non-local and CLS model are shown in Figs. 12$B$ and 12$C$. The non-local model has spurious currents that are at least two orders of magnitude smaller than the CLS method for all quantities considered.