# Phase-Field Simulations for Dripping-to-Jetting Transitions: Effects of Low
Interfacial Tension and Bulk Diffusion

Fukeng Huang (黄富铿)  Weizhu Bao (包维柱)  Department of Mathematics, National
University of Singapore, Singapore, 119076  Tiezheng Qian (钱铁铮) Corresponding
author<EMAIL_ADDRESS>Department of Mathematics, The Hong Kong University of
Science and Technology
Clear Water Bay, Kowloon, Hong Kong, P. R. China

###### Abstract

The dripping-to-jetting transitions in coaxial flows have been experimentally
well studied for systems of high interfacial tension, where the capillary
number of the outer fluid and the Weber number of the inner fluid are in
control. Recent experiments have shown that in systems of low interfacial
tension, the transitions driven by the inner flow are no longer dominated by
the inertial force alone, and the viscous drag force due to the inner flow is
also quantitatively important. In the present work, we carry out numerical
simulations based on the Cahn-Hilliard-Navier-Stokes model, aiming for a more
complete and quantitative study to understand the effects of interfacial
tension when it becomes sufficiently low. The Cahn-Hilliard-Navier-Stokes
model is solved by using an accurate and efficient spectral method in a
cylindrical domain with axisymmetry. Plenty of numerical examples are
systematically presented to show the dripping-to-jetting transitions driven by
the outer flow and inner flow respectively. In particular, for transitions
dominated by inner flow, detailed results reveal how the magnitude of
interfacial tension quantitatively determines the relative importance of the
inertial and viscous forces due to the inner flow at the transition point. Our
numerical results are found to be consistent with the experimental
observation. Finally, the degree of bulk diffusion is varied to investigate
its quantitative effect on the condition for the occurrence of transition.
Such effect is expected for systems of ultralow interfacial tension where
interfacial motion is more likely to be driven by bulk diffusion.

††preprint: AIP/123-QED

##  I  Introduction

Dripping and jetting in coaxial flows of two immiscible fluids refer to the
phenomena in which two fluids are forced to flow through a cylindrical
conduit, with one fluid, namely the inner fluid, flowing at the center and the
other, namely the outer fluid, flowing around it in a coaxial manner. When the
flow rates of both fluids are low, dripping occurs due to the capillary
instability, with the inner fluid forming discrete drops close to the orifice.
On the other hand, when the flow rates are high enough, jetting occurs, with
the inner fluid forming a continuous jet that extends out of the orifice and
breaks into drops further downstream. The dripping and jetting of coaxial
flows of two immiscible fluids have many applications in ink jet printing,
biomedical engineering, and materials engineering whitesides2006origins;
stone2004engineering, and the transition from dripping to jetting is of
fundamental importance in these applications involving drop formation
marre2009dripping; kaufman2012structured. Extensive research efforts have been
made to study the dripping-to-jetting transition in coaxial geometry
utada2007dripping; guillot2007stability; castro2009scaling. Among these works,
Utada et al. utada2007dripping demonstrated that the transitions in coflowing
streams can be characterized by the capillary number of the outer fluid and
the Weber number of the inner fluid. The dripping-to-jetting transitions have
also been investigated in other geometries such as flow-focusing
ganan1998generation; cubaud2008capillary and T-junction thorsen2001dynamic;
abate2009impact. For a review summarizing the main observations and
understandings for common device geometries, we refer to nunes2013dripping and
the references therein.

To understand the hydrodynamics of the dripping-to-jetting transitions, most
of the previous studies have focused on systems, e.g., oil-water ones, of high
interfacial tension. For these systems, as demonstrated in utada2007dripping,
the dripping-to-jetting transitions can be described by a state diagram that
is controlled by the capillary number of the outer fluid and the Weber number
of the inner fluid. This means that for transitions driven by the inner fluid,
the effect of the viscous force due to the inner flow is negligible. However,
it has been shown experimentally in mak2017dripping that when the interfacial
tension is sufficiently low, the dripping-to-jetting transitions driven by the
inner fluid are no longer dominated by the inertial force alone, and the
viscous force due to the inner flow also plays a quantitatively important
role. Therefore, for a comprehensive understanding of the dripping-to-jetting
transitions, a more complete and quantitative study is needed to investigate
the effect of interfacial tension when it becomes sufficiently low. This will
help clarify the relative importance of the inertial and viscous forces due to
the inner flow in inducing the transitions. In addition, recent observations
in aqueous two-phase systems with ultra-low surface tension have revealed
novel and interesting pinch-off dynamics dominated by bulk diffusion Xu2019.
To the best of our knowledge, the effect of bulk diffusion on the dripping-to-
jetting transitions in systems of ultra-low interfacial tension has never been
investigated before.

The main purpose of the present work is to investigate the dripping-to-jetting
transitions in coaxial flows over a wide range of interfacial tension and with
variable bulk diffusion. Firstly, we aim to numerically observe and examine
the dripping-to-jetting transitions driven by the outer and inner fluids
respectively. This is to provide a crucial indicator to distinguish jetting
from dripping and hence locate the point of transition and the critical flow
rate. Secondly, regarding the contributions of inertial and viscous forces due
to the inner flow, we aim to establish a quantitative relationship between
them at the point of transition over a wide range of interfacial tension. Last
but not least, we will examine the quantitative effect of bulk diffusion on
the critical flow rate at the point of transition. This will also show if bulk
diffusion can change the relative importance of the inertial and viscous
forces due to the inner flow at the transition point.

While numerous experimental studies having been carried out on the dripping-
to-jetting transitions in coaxial flows, there have been very few works
focusing on numerical simulations. Guillaument et al. guillaument2013numerical
utilized the one-fluid model and the volume of fluid method to simulate
segmented micro coflows of CO2 and water in two dimensions. Lei et al.
lei2011dripping employed the phase-field model to investigate two types of
transitions driven by the outer flow and the inner flow in two dimensions.
Shahin et al. shahin2017three simulated dripping and jetting in a coflowing
system using a one-fluid model in three dimensions and developed a novel
algorithm to handle the topological change of the interface mesh.

To investigate the dripping-to-jetting transitions in immiscible two-phase
flows, we will employ the Cahn-Hilliard-Navier-Stokes (CHNS) model and carry
out numerical computation in a cylindrical domain with axisymmetry. The phase-
field methods have been widely used in the simulations of interfacial motion
in multiphase flows anderson1998diffuse; boettinger2002phase; yue2004diffuse;
Qian2003; Qian2004; yang2006JCP as they avoid the need of interface tracking
and can easily and efficiently accommodate topological changes such as pinch-
off, a key feature of the dripping and jetting phenomena.

To the best of our knowledge, there has been no prior work that investigates
the dripping-to-jetting transitions in three dimensions using the phase-field
method. For the CHNS model employed here, a characteristic length scale has
been introduced in huang2022diffuse to measure the competition between
diffusion and viscous flow in interfacial motion. Parameters involved in
defining this length scale can be adjusted to tune the effect of bulk
diffusion in the simulated system. Numerically, we solve the CHNS model by
using the spectral method Shen97 for the spatial discretization and the
pressure-correction method shen1992error; lopez1998efficient for the temporal
discretization. These methods have been demonstrated to be accurate and
efficient in treating the phase-field models in cylindrical domains
yang2006JCP; lopez1998efficient; lopez2002efficient.

This paper is organized as follows. In Sec. II.1, the CHNS model is derived by
applying Onsager’s variational principle Qian2006; Onsager1931a; onsager1931b.
In Sec. II.2, the dimensionless equation system is presented with important
dimensionless parameters associated with the dripping-to-jetting transitions,
and the simulated systems are described in a cylindrical domain with necessary
boundary conditions for the inner and outer flows with adjustable flow rates.
In Sec. III, numerical results are presented to show the distinct between
dripping and jetting in the regime dominated by the outer flow and that by the
inner flow, respectively. Furthermore, in the regime dominated by the inner
flow, the relative importance of the inertial and viscous forces at the
transition point is investigated over a wide range of interfacial tension,
with numerical results showing agreement with recent experiments. Finally, the
quantitative effect of bulk diffusion on the critical flow rates at the
transition point is also measured. In Sec. IV, the paper is concluded with a
few remarks.

##  II  Modeling and simulation for immiscible two-phase flows

###  II.1  The Cahn-Hilliard-Navier-Stokes model

Consider a multi-component fluid with two co-existing immiscible phases. A
diffuse-interface model uses a Ginzberg-Landau-type free energy functional to
describe the thermodynamic properties of the fluid. Here we use the Cahn-
Hilliard (CH) free energy functional CH1958

$F_{\rm
CH}[\phi]=\int\left[\displaystyle\frac{K}{2}\left(\nabla\phi\right)^{2}+f(\phi)\right]d\mathbf{r},$
(1)

in which $\phi:=\phi(\mathbf{r})$ is the phase-field variable to measure the
local relative concentration, $f(\phi)$ is the Helmholtz free energy density
for a homogeneous phase, and $K$ is a positive material parameter. The free
energy density $f$ is given by
$f(\phi)=-\frac{\alpha}{2}\phi^{2}+\frac{\beta}{4}\phi^{4}$, which has a
double-well structure to stabilize the fluid-fluid interface between the two
co-existing phases around
$\phi_{\pm}=\pm\phi_{0}=\pm\sqrt{\frac{\alpha}{\beta}}$, where $\alpha$ and
$\beta$ are two positive parameters. Subject to appropriate boundary
conditions, $F_{\rm CH}[\phi]$ can be minimized to stabilize a flat interface
between the two equilibrium phases of $\phi=\pm\phi_{0}$. The interfacial
structure gives the interfacial tension
$\gamma=\frac{2\sqrt{2}\alpha^{2}\xi}{3\beta}$ and the characteristic length
scale $\xi=\sqrt{\frac{K}{\alpha}}$ for the interfacial thickness Qian2006.
Note that in many literatures, $\phi_{0}$ is made to equal $1$ through a
rescaling. Here $\phi_{0}=\sqrt{\frac{\alpha}{\beta}}$ is purposely retained
to measure the distance away from the critical point where $\phi_{0}$
vanishes.

For an incompressible fluid, the velocity field $\mathbf{v}$ is subject to the
incompressibility condition $\nabla\cdot\mathbf{v}=0$, and the phase field
$\phi$ satisfies the continuity equation

$\displaystyle\frac{\partial\phi}{\partial
t}=-\nabla\cdot\mathbf{J}=-\nabla\cdot\left(\phi\mathbf{v}+\mathbf{j}\right),$
(2)

where $\mathbf{J}=\mathbf{\phi\mathbf{v}+\mathbf{j}}$ is the total current
density, in which $\phi\mathbf{v}$ is contributed by the flow and $\mathbf{j}$
is the diffusive current density contributed by the bulk diffusion.

Hydrodynamic equations for immiscible two-phase flows can be derived by
applying Onsager’s variational principle (cf. appendix A in huang2022diffuse)
as follows. The Rayleighian ${\cal{R}}$ is given by ${\cal R}=\dot{F}_{\rm
CH}+\Phi$ in the bulk region. Here $\dot{F}_{\rm CH}$ is the rate of change of
$F_{\rm CH}[\phi]$, given by

$\dot{F}_{\rm CH}[\phi]=\int\mu\displaystyle\frac{\partial\phi}{\partial
t}d\mathbf{r}=\int\nabla\mu\cdot\left(\mathbf{\phi\mathbf{v}+\mathbf{j}}\right)d\mathbf{r},$
(3)

in which $\mu=\frac{\delta F_{\rm CH}}{\delta\phi}$ is the chemical potential,
given by $\mu=-K\nabla^{2}\phi+f^{\prime}(\phi)$, and the continuity equation
(2) has been used with the impermeability conditions for $\mathbf{v}$ and
$\mathbf{j}$ at the solid boundary. The other part in ${\cal R}$ is the
dissipation functional $\Phi$, which is half the rate of free energy
dissipation and given by

$\Phi=\int\displaystyle\frac{\eta}{4}\left[\nabla\mathbf{v}+\left(\nabla\mathbf{v}\right)^{T}\right]^{2}d\mathbf{r}+\int\displaystyle\frac{\mathbf{j}^{2}}{2M}d\mathbf{r},$
(4)

which is contributed by the viscous dissipation, with $\eta$ being the shear
viscosity, and the diffusive dissipation, with $M$ being the mobility
coefficient.

Subject to the incompressibility condition, the Rayleighian can be minimized
with respect to the rates $\mathbf{v}$ and $\mathbf{j}$. This gives the force
balance equation

$-\nabla p+\nabla\cdot\bm{\sigma}_{\rm visc}-\phi\nabla\mu=0$ (5)

for $\mathbf{v}$, and the constitutive equation

$\mathbf{j}=-M\nabla\mu$ (6)

for $\mathbf{j}$. Here $p$ is the pressure, which is the Lagrange multiplier
to locally impose $\nabla\cdot\mathbf{v}=0$, $\bm{\sigma}_{\rm visc}$ is the
Newtonian stress tensor given by $\bm{\sigma}_{\rm
visc}=\eta\left[\nabla\mathbf{v}+\left(\nabla\mathbf{v}\right)^{T}\right]$.
Equation (5) is the Stokes equation with the capillary force density, and it
can be readily generalized to the Navier-Stokes equation

$\rho\left[\displaystyle\frac{\partial\mathbf{v}}{\partial
t}+\left(\mathbf{v}\cdot\nabla\right)\mathbf{v}\right]=-\nabla
p+\nabla\cdot\bm{\sigma}_{\rm visc}-\phi\nabla\mu.$ (7)

Combining equations (2) and (6) gives the advection-diffusion equation for the
phase field $\phi$:

$\displaystyle\frac{\partial\phi}{\partial
t}+\mathbf{v}\cdot\nabla\phi=-\nabla\cdot\mathbf{j}=M\nabla^{2}\mu,$ (8)

which is the CH equation for a constant mobility $M$. Equations (7) and (8)
govern the hydrodynamics of immiscible two-phase flows.  In the present work,
the simplest situation is treated with the two fluids having equal density,
equal viscosity and equal mobility.

###  II.2  Dimensionless equations and simulated systems

Numerical simulations are carried out by solving the CHNS system:

$\displaystyle\frac{\partial\phi}{\partial
t}+\mathbf{v}\cdot\nabla\phi=M\nabla^{2}\mu,$ (9a)
$\displaystyle\mu=-K\nabla^{2}\phi-\alpha\phi+\beta\phi^{3},$ (9b)
$\displaystyle\rho\big{(}\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{v}\big{)}=-\nabla
p+\eta\nabla^{2}\mathbf{v}+\mu\nabla\phi,$ (9c)
$\displaystyle\nabla\cdot\mathbf{v}=0,$ (9d)

in a cylindrical domain
$\Omega=\\{\mathbf{r}=(x,y,z):x^{2}+y^{2}<L^{2},z\in(0,H)\\}$. Here $M$, $K$,
$\alpha$, $\beta$, $\rho$, and $\eta$ are material parameters introduced in
Sec. II.1. Note that the pressure $p$ in equation (9c) is different from that
in equation (7), with $-\phi\nabla\mu$ there being replaced by $\mu\nabla\phi$
here. The boundary conditions on $x^{2}+y^{2}=L^{2}$ are

$\frac{\partial\phi}{\partial\mathbf{n}}=0,\quad\frac{\partial\mu}{\partial\mathbf{n}}=0,\quad\mathbf{v}=0.$
(10)

In our simulations, two immiscible phases flow into the cylinder on the
boundary $z=0$ and out of the cylinder on the boundary $z=H$. The boundary
conditions there for $\phi$ and $\mu$ are given by

$\phi=\tanh\Big{(}\frac{r-R}{\sqrt{2}\xi}\Big{)},\quad\mu=0,$ (11)

on $z=0$, with $r=\sqrt{x^{2}+y^{2}}$, $R$ being the radius of the inner tube,
and

$\frac{\partial\phi}{\partial\mathbf{n}}=0,\quad\frac{\partial\mu}{\partial\mathbf{n}}=0,$
(12)

on $z=H$. Finally, the boundary conditions for
$\mathbf{v}:=(v_{x},v_{y},v_{z})$ on $z=0$ and $z=H$ are given by

$v_{x}=v_{y}=0,\quad v_{z}=\left\\{\begin{array}[]{lr}a(R^{2}-r^{2}),\quad
0<r<R,\\\
-b(r^{2}-R^{2})+\frac{b(L^{2}-R^{2})}{\ln\frac{L}{R}}\ln\frac{r}{R},\quad
R\leq r<L,\end{array}\right.$ (13)

on $z=0$, where $a$ and $b$ are the parameters determining the mean velocities
(i.e., flow rates) of the inner and outer phases, respectively, and

$v_{x}=v_{y}=0,\quad v_{z}=c(L^{2}-r^{2}),$ (14)

on $z=H$, where $c$ is the parameter determining the mean velocity of the flow
out of the cylinder. Note that these flow profiles are based on the Poiseuille
profile, and the parameters $a$, $b$ and $c$ satisfy

$aR^{4}+b(L^{4}-R^{4})-b\frac{(L^{2}-R^{2})^{2}}{\ln\frac{L}{R}}=cL^{4},$ (15)

for the volume conservation.

To nondimensionalize the above system, we use the radius $L$ of the
computational domain $\Omega$ as the length unit, $u=\frac{\gamma}{\eta}$ as
the velocity unit, $\tau=\frac{L}{u}$ as the time unit, and
$p_{0}=\frac{\eta}{\tau}$ as the pressure unit. We also define the following
quantities:

  * •

$\bar{H}=\frac{H}{L}$ as the dimensionless length of the computational domain,

  * •

$\bar{R}=\frac{R}{L}$ as the dimensionless radius of the inner tube,

  * •

$\phi_{0}=\sqrt{\frac{\alpha}{\beta}}$, with the two equilibrium phases
separated by a flat interface being of $\phi=\pm\phi_{0}$,

  * •

$\varepsilon=\frac{\xi}{L}=\frac{1}{L}\sqrt{\frac{K}{\alpha}}$ as the
dimensionless interfacial thickness of the diffuse interface,

  * •

$D=2M\alpha$ as the diffusion coefficient for $\phi$ close to $\pm\phi_{0}$
far away from the interface,

  * •

$l_{c}=\frac{\sqrt{M\eta}}{\phi_{0}}$ as the characteristic length scale,
determined from the competition between diffusion and viscous flow
huang2022diffuse,

  * •

$\gamma=\frac{2\sqrt{2}}{3}\alpha\phi_{0}^{2}{\xi}$ as the interfacial
tension,

  * •

$Re_{\gamma}=\frac{\rho uL}{\eta}$ as the Reynolds number defined from the
velocity unit $u=\frac{\gamma}{\eta}$ and the length unit $L$,

  * •

$B=\frac{\eta D}{\alpha\phi_{0}^{2}L^{2}}=\frac{2l_{c}^{2}}{L^{2}}$ as the
dimensionless parameter measuring the characteristic length scale $l_{c}$ with
respect to $L$.

Dimensionless variables, denoted by using overbar, are defined as follows:

  * •

$\bar{\phi}=\frac{\phi}{\phi_{0}}$, $\bar{\mathbf{v}}=\frac{\mathbf{v}}{u}$,
$\bar{\mu}=\frac{\mu}{\alpha\phi_{0}\varepsilon}$, $\bar{p}=\frac{p}{p_{0}}$,

and the dimensionless operators:

  * •

$\frac{\partial}{\partial\bar{t}}=\tau\frac{\partial}{\partial{t}}$,
$\bar{\nabla}=L{\nabla}$.

Using the above definitions, we obtain the dimensionless CHNS system in the
cylindrical domain
$\bar{\Omega}=\\{(\bar{x},\bar{y},\bar{z}):\bar{x}^{2}+\bar{y}^{2}<1,\bar{z}\in(0,\bar{H})\\}$
as

$\displaystyle\frac{\partial\bar{\phi}}{\partial\bar{t}}+\bar{\mathbf{v}}\cdot\bar{\nabla}\bar{\phi}=\frac{3}{4\sqrt{2}}B\bar{\nabla}^{2}\bar{\mu},$
(16a)
$\displaystyle\bar{\mu}=-\varepsilon\bar{\nabla}^{2}\bar{\phi}+\frac{1}{\varepsilon}(-\bar{\phi}+\bar{\phi}^{3}),$
(16b) $\displaystyle
Re_{\gamma}\left(\frac{\partial\bar{\mathbf{v}}}{\partial\bar{t}}+\bar{\mathbf{v}}\cdot\bar{\nabla}\bar{\mathbf{v}}\right)=-\bar{\nabla}\bar{p}+\bar{\nabla}^{2}\bar{\mathbf{v}}+\frac{3}{2\sqrt{2}}\bar{\mu}\bar{\nabla}\bar{\phi},$
(16c) $\displaystyle\bar{\nabla}\cdot\bar{\mathbf{v}}=0.$ (16d)

The boundary conditions are

$\frac{\partial\bar{\phi}}{\partial\mathbf{n}}=0,\quad\frac{\partial\bar{\mu}}{\partial\mathbf{n}}=0,\quad\mathbf{\bar{v}}=0,$
(17)

on $\bar{x}^{2}+\bar{y}^{2}=1$,

$\begin{array}[]{l}\bar{\mu}=0,\quad\bar{\phi}=\tanh(\frac{\bar{r}-\bar{R}}{\sqrt{2}\varepsilon}),\\\
\bar{v}_{x}=\bar{v}_{y}=0,\quad\bar{v}_{z}=\left\\{\begin{array}[]{lr}\bar{a}(\bar{R}^{2}-\bar{r}^{2}),\quad
0<\bar{r}<\bar{R},\\\
-\bar{b}(\bar{r}^{2}-\bar{R}^{2})+\frac{\bar{b}(1-\bar{R}^{2})}{\ln\frac{1}{\bar{R}}}\ln\frac{\bar{r}}{\bar{R}},\quad\bar{R}\leq\bar{r}<1,\end{array}\right.\end{array}$
(18)

on $\bar{z}=0$ with $\bar{r}=\sqrt{\bar{x}^{2}+\bar{y}^{2}}$, and

$\begin{array}[]{l}\frac{\partial\bar{\phi}}{\partial\mathbf{n}}=0,\quad\frac{\partial\bar{\mu}}{\partial\mathbf{n}}=0,\\\
\bar{v}_{x}=\bar{v}_{y}=0,\quad\bar{v}_{z}=\bar{c}(1-\bar{r}^{2}),\end{array}$
(19)

on $\bar{z}=\bar{H}$, with the dimensionless parameters $\bar{a}$, $\bar{b}$
and $\bar{c}$ satisfying

$\bar{a}\bar{R}^{4}+\bar{b}(1-\bar{R}^{4})-\bar{b}\frac{(1-\bar{R}^{2})^{2}}{\ln\frac{1}{\bar{R}}}=\bar{c}.$
(20)

Here the dimensionless $\bar{a}$, $\bar{b}$ and $\bar{c}$ are obtained by
multiplying the dimensional ones by $\frac{L^{2}}{u}$. The above dimensionless
CHNS system involves the dimensionless parameters $\varepsilon$,
$Re_{\gamma}$, $B$, $\bar{R}$, $\bar{H}$, $\bar{a}$, and $\bar{b}$. Here
$\varepsilon$ is the dimensionless interfacial thickness, which is the
smallest length to be resolved, $Re_{\gamma}$ is the Reynolds number defined
from the velocity unit $u=\frac{\gamma}{\eta}$, $B$ controls the competition
between bulk diffusion and viscous flow, $\bar{R}$ measures the size of the
orifice (i.e., the radius of the inner tube), $\bar{H}$ measures the length of
the computational domain, and $\bar{a}$ and $\bar{b}$ control the flow rates
of the inner and outer fluids.

From our three-dimensional (3D) simulations, it is verified that given an
axisymmetric initial condition in the cylindrical domain, the axisymmetry can
be accurately preserved during the whole dynamic process. Therefore, in the
absence of any evidence for non-axisymmetric modes, we treat the axisymmetric
3D problem as a reduced two-dimensional (2D) problem by making use of the
cylindrical coordinates to improve the computational efficiency
huang2022diffuse.  Technically, we first transform the 3D problem into a 2D
problem using the cylindrical coordinates Shen97. We then adopt the usual
semi-implicit scheme to solve the phase-field variable and the spectral-
projection method to solve the velocity and pressure fields for the Navier-
Stokes equation in cylindrical geometry lopez1998efficient. At each time step,
we can efficiently solve a series of Poisson-type equations with constant
coefficients.

##  III  Results and discussion

With the dimensionless parameters introduced in the previous section, the
average velocity of the inner flow $\bar{v}_{\rm in}$ and that of the outer
flow $\bar{v}_{\rm out}$ are given by

$\bar{v}_{\rm in}=\frac{1}{2}\bar{a}\bar{R}^{2},\quad\bar{v}_{\rm
out}=\frac{\bar{b}}{2}(1+\bar{R}^{2})+\frac{\bar{b}}{2}\frac{1-\bar{R}^{2}}{\ln\bar{R}}.$
(21)

Using $\bar{v}_{\rm in}$ and $\bar{v}_{\rm out}$, the capillary number of the
outer flow $\mathcal{C}_{\rm out}$ and the Weber number of the inner flow
$\mathcal{W}_{\rm in}$ can be expressed as

$\mathcal{C}_{\rm out}=\bar{v}_{\rm out},\quad\mathcal{W}_{\rm
in}=\bar{v}_{\rm in}^{2}Re_{\gamma}\bar{R}.$ (22)

Here $\mathcal{C}_{\rm out}$ is defined by $\mathcal{C}_{\rm
out}=\frac{\eta(\bar{v}_{\rm out}u)}{\gamma}$, and $\mathcal{W}_{\rm in}$ is
defined by $\mathcal{W}_{\rm in}=\frac{\rho(\bar{v}_{\rm in}u)^{2}R}{\gamma}$,
where $\bar{v}_{\rm in}u$ and $\bar{v}_{\rm out}u$ are the dimensional average
velocities with $u$ being the velocity unit. Physically, the capillary number
measures the viscous drag force, and the Weber number measures the inertial
force relative to the interfacial tension force.

It has been well established that there are two classes of dripping-to-jetting
transitions in coflowing streams utada2007dripping. The first one is driven by
strong outer flows and will be numerically investigated in Section III.1 by
fixing a small $\bar{v}_{\rm in}$ and varying $\bar{v}_{\rm out}$. The second
one is driven by strong inner flows, and will be numerically investigated in
Section III.2 by fixing a small $\bar{v}_{\rm out}$ and varying $\bar{v}_{\rm
in}$. In this regime, our numerical results show that in addition to the
inertial force measured by $\mathcal{W}_{\rm in}$, the viscous force due to
the inner flow, measured by the capillary number $\mathcal{C}_{\rm
in}=\bar{v}_{\rm in}$, also contributes to the occurrence of dripping-to-
jetting transition when the interfacial tension is sufficiently low. This
numerical observation is in agreement with recent experiments mak2017dripping.
Finally, Section III.3 demonstrates the quantitative effect of bulk diffusion
on the critical flow rates at the transition point. Such effect is expected
for systems of ultralow interfacial tension where interfacial motion is more
likely to be driven by bulk diffusion Xu2019.

###  III.1  Transitions dominated by outer flows

In this subsection, we investigate the first class of dripping-to-jetting
transitions driven by strong outer flows. For this purpose, the value of
$\bar{v}_{\rm in}$ is fixed to be small, and the value of $\bar{v}_{\rm out}$
is increased to induce the transition.

We start by demonstrating the dripping-to-jetting transitions in the regime
dominated by strong outer flows. Let $Z_{p}$ denote the distance between the
pinch-off position and the boundary of $\bar{z}=0$ (the orifice). For a slow
inner flow with $\bar{a}=15$ being fixed, $Z_{p}$ is expected to increase with
the increasing outer flow rate, i.e., the increasing $\bar{b}$. Figure 1 shows
two different pinch-off positions for two different outer flow rates. It is
clearly observed from figure 2 that $Z_{p}$ exhibits a sharp increase from
$\bar{b}=0.575$ (figure 1(a)) to $\bar{b}=0.6$ (figure 1(b)), indicating a
transition from a dripping state to a jetting state as $\mathcal{C}_{\rm out}$
is increased from $0.1450$ to $0.2030$. This critical magnitude of
$\mathcal{C}_{\rm out}$ is in agreement with the experimental results in
utada2007dripping.  According to the state diagram reported in
utada2007dripping, for the dripping-to-jetting transitions dominated by outer
flows, the critical values of $\mathcal{C}_{\rm out}$ are typically
distributed between $0.2$ and $0.4$.

Figure 1:  Two different pinch-off positions for two different outer flow
rates. (a) A dripping state for $\bar{b}=0.575$. (b) A jetting state for
$\bar{b}=0.6$. Other parameter values used in simulations are
$\varepsilon=0.01$, $Re_{\gamma}=500$, $B=0.0002$, $\bar{R}=0.1$, $\bar{H}=6$,
and $\bar{a}=15$. Figure 2:  Variation of the pinch-off position $Z_{p}$ with
the parameter $\bar{b}$ which controls the outer flow rate. A transition is
noted to occur between $\bar{b}=0.575$ (dripping in figure 1(a)) and
$\bar{b}=0.6$ (jetting in figure 1(b)). Other parameter values used in
simulations are $\varepsilon=0.01$, $Re_{\gamma}=500$, $B=0.0002$,
$\bar{R}=0.1$, $\bar{H}=6$, and $\bar{a}=15$.

A jetting state maintained by a strong outer flow is characterized by a long,
narrow jet and small drops utada2007dripping. In fact, a stronger outer flow
results in a narrower jet and smaller drops. Here we present some quantitative
results on the relationship between the outer flow rate and the size of the
corresponding jet, i.e., the radius of the jet. To obtain reliable data, we
have ensured that the jets are long and wide enough by using values of
$\bar{b}$ and $\bar{R}$ that are sufficiently large. Figure 3(a) presents a
jetting state obtained from our simulations, and figure 3(b) shows the
dependence of the jet radius $r_{j}$ on the outer flow rate ($\bar{v}_{\rm
out}\propto\bar{b}$), with the jet radius $r_{j}$ being measured at the plane
of $\bar{z}=\frac{\bar{H}}{2}$. When $\bar{a}$ and $\bar{R}$ are both fixed,
the total flux of inner fluid is given, and a faster outer flow (with a larger
$\bar{b}$) leads to a thinner jet in which the inner fluid flows with a larger
average velocity ($\propto\bar{b}$). According to mass conservation,
$\bar{b}r_{j}^{2}$ must be a constant in order to maintain the total flux of
inner fluid, as shown in figure 3(b).

We can measure both the jet diameter $d_{\rm jet}$ and the drop diameter
$d_{\rm drop}$ in the jetting regime. From these two diameters, we obtain
$\lambda_{f}$, the wavelength of the fastest growing mode of the Rayleigh-
Plateau instability, through the relation $\frac{\pi}{4}d_{\rm
jet}^{2}\lambda_{f}=\frac{\pi}{6}d_{\rm drop}^{3}$ for drop volume. Using the
simulation results shown in figure 3(a), we obtain $d_{\rm drop}\approx
2d_{\rm jet}$ and hence $\lambda\approx 5.3d_{\rm jet}$, which is in the
physically reasonable range. It is noted that for dripping-to-jetting
transitions dominated by outer flows, $d_{\rm drop}\approx 2d_{\rm jet}$ has
been experimentally observed utada2007dripping; mak2017dripping.

Figure 3:  (a) A jetting state obtained from our simulations. Note that the
radius of the jet $r_{j}$ is measured at the plane of
$\bar{z}=\frac{\bar{H}}{2}$. (b) Log-log plot of the jet radius $r_{j}$ versus
$\bar{b}$ which controls the outer flow rate. Here the mass conservation of
the inner fluid is ensured by $\bar{b}r_{j}^{2}\approx 0.011$. The jetting
state in (a) is obtained for $\bar{b}=0.5$. Other parameter values used in
simulations are $\varepsilon=0.01$, $Re_{\gamma}=10$, $B=0.0002$,
$\bar{R}=0.2$, $\bar{H}=8$, and $\bar{a}=6$.

###  III.2  Transitions dominated by inner flows

In this subsection, we investigate the second class of dripping-to-jetting
transitions driven by strong inner flows. For this purpose, the value of
$\bar{v}_{\rm out}$ is fixed to be small, and the value of $\bar{v}_{\rm in}$
is increased to induce the transition.  For $\bar{R}=0.1$, we have
$\bar{v}_{\rm out}=0.29\bar{b}$ and $\mathcal{C}_{\rm out}=0.029$ for the
typical value $0.1$ used for $\bar{b}$.

We start from the drop size in the dripping regime. When $\bar{v}_{\rm in}$ is
not large enough, the system is in a dripping state in which drops of the same
size are periodically generated at the same pinch-off position. From the
periodic dynamics and mass conservation, we obtain
$\frac{4}{3}\pi(\frac{d_{e}}{2})^{3}=\pi\bar{R}^{2}\bar{v}_{\rm in}t_{p}$,
where $\bar{R}$ is the radius of the orifice, $\bar{v}_{\rm in}$ is the
average velocity of the inner fluid, $t_{p}$ is the time period of the
periodic generation of drops, and $d_{e}$ is the diameter of the drops
expected from mass conservation. Figure 4(a) shows a comparison between the
expected diameter $d_{e}$ and the diameter $d_{m}$ which is measured in our
numerical simulations. It is noted that in each simulation, $d_{m}$ is
slightly smaller than $d_{e}$ expected from mass conservation. This is
attributed to the bulk diffusion which continuously reduce the size of drops.

To understand how the drop size is controlled by the inner and outer flows, we
show that the time period $t_{p}$ can be related to the drop diameter $d$ as
follows:

$t_{p}\approx\frac{\kappa d}{\bar{v}_{\rm
out}+\nu\frac{\bar{R}^{2}}{d^{2}}\bar{v}_{\rm in}}$ (23)

where $\mu$ and $\nu$ are two adjustable parameters of the order of magnitude
of $1$, and $\bar{v}_{\rm in}$ and $\bar{v}_{\rm out}$ have been defined in
(21). For $\bar{R}=0.1$, we have

$t_{p}\approx\frac{\kappa d}{(0.29\bar{b}+\nu\frac{10^{-4}}{2d^{2}}\bar{a})},$
(24)

which has been numerically verified by figure 4(b) in which the measured
diameter $d_{m}$ is used for the drop diameter $d$.  Physically, equation (23)
describes the advection of a growing drop, with the advected distance being
$\sim d$ and the velocity being $\sim\bar{v}_{\rm
out}+\nu\frac{\bar{R}^{2}}{d^{2}}\bar{v}_{\rm in}$, in which the contribution
of $\bar{v}_{\rm in}$ is rescaled by a factor $\sim\frac{\bar{R}^{2}}{d^{2}}$.
Equation (24) is then obtained by using equation (21) to express $\bar{v}_{\rm
in}$ and $\bar{v}_{\rm out}$ for $\bar{R}=0.1$. From our simulation results,
the data points in figure 4(b) are produced by using optimal values for the
adjustable parameters $\kappa$ and $\nu$ to best fit the solid line
representing equation (24).  Furthermore, it is seen from the inset to figure
4(b) that the contribution of $\bar{v}_{\rm out}$ is much larger than that of
$\nu\frac{\bar{R}^{2}}{d^{2}}\bar{v}_{\rm in}$ in equation (23), i.e., the
contribution of $0.29\bar{b}$ is much larger than that of
$\nu\frac{10^{-4}}{2d^{2}}\bar{a}$ in equation (24) for $\bar{R}=0.1$. This
means that for the advection of a growing drop, the distance is typically
$\sim d$, and the velocity is predominantly $\sim\bar{v}_{\rm out}$. It
follows that the time period $t_{p}$ of drop generation is
$\sim\frac{d}{\bar{v}_{\rm out}}$. Combining $t_{p}\sim\frac{d}{\bar{v}_{\rm
out}}$ and $\frac{4}{3}\pi(\frac{d}{2})^{3}=\pi\bar{R}^{2}\bar{v}_{\rm
in}t_{p}$ from mass conservation, we have
$d\sim\bar{R}\sqrt{\frac{\bar{v}_{\rm in}}{\bar{v}_{\rm out}}}$,  which has
been experimentally verified utada2007dripping.

Figure 4:  (a) A comparison between the diameter $d_{e}$ expected from mass
conservation and the diameter $d_{m}$ measured in our simulations. Note that
$d_{m}$ is always slightly smaller than $d_{e}$ due to the bulk diffusion. (b)
The relation between the drop diameter $d_{m}$ and the time period $t_{p}$ of
drop generation. Here the parameters $\bar{a}$ and $\bar{b}$, which control
$\bar{v}_{\rm in}$ and $\bar{v}_{\rm out}$, are also involved according to
equations (23) and (24), with $\kappa=3.1$ and $\nu=3$. The inset shows that
the contribution of $0.29\bar{b}$ is much larger than that of
$\nu\frac{10^{-4}}{2d^{2}}\bar{a}$ in equation (24), indicating that the
growing drop is mainly advected by the outer flow. The data are obtained by
using $\varepsilon=0.01$, $B=0.0002$, $\bar{R}=0.1$, $\bar{H}=4$, and
different combinations of $Re_{\gamma}$, $\bar{a}$ and $\bar{b}$.

Now we focus on the dripping-to-jetting transitions dominated by inner flows.
Same as done in the previous subsection, we use $Z_{p}$ to denote the distance
between the pinch-off position and the boundary of $\bar{z}=0$ (the orifice).
For a slow outer flow fixed at $\bar{b}=0.1$, $Z_{p}$ is expected to increase
with the increasing inner flow rate, i.e., the increasing $\bar{a}$. Figure
5(a) shows the pinch-off position in a dripping state for $\bar{a}=24$ just
before the transition, and figure 5(b) shows the pinch-off position in a
jetting state for $\bar{a}=25$ just after the transition. It is clearly
observed that there is a sharp increase of $Z_{p}$ from figure 5(a) to 5(b),
indicating the occurrence of a dripping-to-jetting transition. Here the value
of $Re_{\gamma}$ is $500$, and we have $\mathcal{W}_{\rm in}=0.72$ and
$\mathcal{C}_{\rm in}=0.12$ for $\bar{a}=24$, and $\mathcal{W}_{\rm in}=0.781$
and $\mathcal{C}_{\rm in}=0.125$ for $\bar{a}=25$. It is noted that the value
of $Re_{\gamma}$ used here is large enough to let $\mathcal{W}_{\rm in}$ be in
control, with $\mathcal{C}_{\rm in}$ being less important. It is also noted
that the critical magnitude of $\mathcal{W}_{\rm in}$ is in agreement with the
experimental results in utada2007dripping.  According to the state diagram
reported in utada2007dripping, for the dripping-to-jetting transitions
dominated by inner flows, the critical values of $\mathcal{W}_{\rm in}$ are
typically distributed around $1$.

To understand the underlying mechanism of the dripping-to-jetting transitions
dominated by inner flows, we use figures 5(c) and 5(d) to show the variation
of the neck radius $r_{n}$ with the neck position $z_{n}$ as time goes on.
Note that as the neck radius approaches $0$, i.e., $r_{n}\to 0$, pinch-off
occurs with the neck position approaching the pinch-off position, i.e.,
$z_{n}\to Z_{p}$. In the dripping regime, it is observed that $r_{n}$
decreases monotonically to $0$ (as shown in figure 5(c)), while in the jetting
regime, $r_{n}$ exhibits a transient increase before it eventually decreases
to $0$ (as shown in figure 5(d)). It is this transient increase of $r_{n}$
that leads to a visible jump in the value of $Z_{p}$ that marks the transition
from dripping to jetting.

Figure 5: (a)-(b) Two different pinch-off positions for two different inner
flow rates, with $\bar{a}=24$ for a dripping state in (a) and $\bar{a}=25$ for
a jetting state in (b). (c)-(d) Variation of the neck radius $r_{n}$ with the
neck position $z_{n}$ as time goes on. Here the third inset to (c) corresponds
to (a) for dripping, and the third inset to (d) corresponds to (b) for
jetting. It is noted in (d) that before $r_{n}$ eventually decreases to $0$,
it exhibits a transient increase that leads to a visible jump in the value of
$Z_{p}$. The data are obtained by using $\varepsilon=0.01$, $Re_{\gamma}=500$,
$B=0.0002$, $\bar{R}=0.1$, $\bar{H}=4$, and $\bar{b}=0.1$.

Figure 6 shows the variation of the pinch-off position $Z_{p}$ with the
parameter $a$ which controls the inner flow rate. For each value of
$Re_{\gamma}$, a transition is noted around a critical value of $a$.
Furthermore, this critical value of $a$ increases with the decreasing
$Re_{\gamma}$.

Figure 6:  Variation of the pinch-off position $Z_{p}$ with the parameter
$\bar{a}$ which controls the inner flow rate. For each value of $Re_{\gamma}$,
a transition is noted around a critical value of $\bar{a}$. The data are
obtained by using $\varepsilon=0.01$, $B=0.0002$, $\bar{R}=0.1$, $\bar{H}=4$,
and $\bar{b}=0.1$.

The jump in the pinch-off position is a clear indicator that can be used to
locate the dripping-to-jetting transition. In the following, we focus on the
critical velocity of the inner flow that is needed to induce the transition,
with the interfacial tension $\gamma$ being varied for nearly two orders of
magnitude. The dripping-to-jetting transitions in systems of high interfacial
tension have been extensively studied utada2007dripping. In particular, when
the transition is dominated by the inner flow (with the outer flow rate
measured by $\mathcal{C}_{\rm out}$ being negligible), the inertial force due
to the inner flow, measured by $\mathcal{W}_{\rm in}$, plays a dominant role
in systems of high interfacial tension. However, when the interfacial tension
is continuously lowered, the viscous force due to the inner flow, measured by
$\mathcal{C}_{\rm in}$, becomes more and more important in driving the
transition. This trend has been reported experimentally mak2017dripping, and a
theoretical understanding can be described as follows. The Weber number of the
inner flow is given by $\mathcal{W}_{\rm in}=\bar{v}_{\rm
in}^{2}Re_{\gamma}\bar{R}$, where the interfacial tension $\gamma$ is involved
in the Reynolds number $Re_{\gamma}$ defined by $Re_{\gamma}=\frac{\rho
uL}{\eta}$, with $u=\frac{\gamma}{\eta}$ being the velocity unit. Let’s
suppose that the transition occurs at $\mathcal{W}_{\rm in}\approx 1$, with
the interfacial tension force being balanced by the inertial force due to the
inner flow. If $Re_{\gamma}$ is made sufficiently small by a sufficiently low
interfacial tension, then the value of $\bar{v}_{\rm in}$ corresponding to
$\mathcal{W}_{\rm in}=\bar{v}_{\rm in}^{2}Re_{\gamma}\bar{R}\approx 1$ can be
made large enough to be comparable to $\mathcal{W}_{\rm in}$. Note that the
capillary number of the inner flow is given by $\mathcal{C}_{\rm
in}=\bar{v}_{\rm in}$. With $\mathcal{C}_{\rm in}$ being comparable to
$\mathcal{W}_{\rm in}\approx 1$ for sufficiently low interfacial tension, it
is deduced that the viscous force due to the inner flow is no longer
negligible compared to the inertial force in driving the transition in systems
of low interfacial tension.

Let $\bar{V}_{\rm in}$ denote the critical velocity of the inner flow. In
systems of high interfacial tension, the inertial force due to the inner flow
is dominant, and hence the transition occurs at $\bar{V}_{\rm
in}^{2}Re_{\gamma}\bar{R}\approx 1$ for the critical Weber number
$\mathcal{W}_{\rm in}\approx 1$. As a result, $\bar{V}_{\rm
in}^{2}Re_{\gamma}={\rm const.}$ is expected for large $Re_{\gamma}$. This is
indeed observed in figure 7(a). When $Re_{\gamma}$ is no longer large enough,
deviation from $\bar{V}_{\rm in}^{2}Re_{\gamma}={\rm const.}$ does show up.
From figure 7(a), it is seen that toward the low end of the range of
$Re_{\gamma}$, the critical $\bar{V}_{\rm in}$ is actually below that
predicted by $\bar{V}_{\rm in}^{2}Re_{\gamma}={\rm const.}$, which only
considers the inertial force due to the inner flow.  As explained above, when
the interfacial tension is low and hence $Re_{\gamma}$ is small, the value of
$\bar{V}_{\rm in}$ predicted by $\bar{V}_{\rm in}^{2}Re_{\gamma}={\rm const.}$
is large. This means a large viscous force due to the inner flow. As a result,
the viscous force and inertial force due to the inner flow are added up to
jointly balance the interfacial tension force. Consequently, the critical
$\bar{V}_{\rm in}$ becomes smaller than that predicted by $\bar{V}_{\rm
in}^{2}Re_{\gamma}={\rm const.}$, which only considers the inertial force due
to the inner flow.

For $Re_{\gamma}$ being varied between $100$ and $4000$, numerical simulations
have been carried out to determine the critical velocity of the inner flow
$\bar{V}_{\rm in}$ at which the transition occurs. The data obtained for
$\bar{V}_{\rm in}$ are used to produce a formula that describes the
contributions of the Weber number of the inner flow $\mathcal{W}_{\rm in}$ and
the capillary number of the inner flow $\mathcal{C}_{\rm in}$ at the
transition. Figure 7(b) shows that $\mathcal{W}_{\rm in}$ and
$\mathcal{C}_{\rm in}$ at the transition satisfy a linear relation given by
$\mathcal{W}_{\rm in}+2.6\mathcal{C}_{\rm in}=1.13$ approximately. It is worth
emphasizing that this equation holds for the interfacial tension $\gamma$
being varied for nearly two orders of magnitude. Note that from the upper left
to the lower right, the value of $Re_{\gamma}$ decreases and consequently the
relative importance of $\mathcal{C}_{\rm in}$ increases. Therefore, it is
numerically verified that the viscous force due to the inner flow plays a
quantitatively important role in driving the dripping-to-jetting transitions
in systems of low interfacial tension.

Figure 7: (a) Log-log plot for the critical velocity of the inner flow
$\bar{V}_{\rm in}$ versus the Reynolds number $Re_{\gamma}$. For large
$Re_{\gamma}$, $\bar{V}_{\rm in}^{2}Re_{\gamma}={\rm const.}$, while for small
$Re_{\gamma}$, the critical $\bar{V}_{\rm in}$ becomes smaller than that
predicted by a constant $\bar{V}_{\rm in}^{2}Re_{\gamma}$. (b) The Weber
number of the inner flow $\mathcal{W}_{\rm in}$ and the capillary number of
the inner flow $\mathcal{C}_{\rm in}$ at the transition satisfy
$\mathcal{W}_{\rm in}+2.6\mathcal{C}_{\rm in}=1.13$ approximately. The data
are obtained by using $\varepsilon=0.01$, $B=0.0002$, $\bar{R}=0.1$,
$\bar{H}=4$, and $\bar{b}=0.1$, with $Re_{\gamma}$ being varied between $100$
and $4000$.  A thin dotted line with arrow is used to indicate the direction
of change of $Re_{\gamma}$.  Note that $\bar{b}=0.1$ used here gives
$\mathcal{C}_{\rm out}=0.029$, which is much smaller than the typical values
of $\mathcal{W}_{\rm in}$ and $\mathcal{C}_{\rm in}$ at the transition.

###  III.3  Effect of bulk diffusion

In this subsection, we investigate the effect of bulk diffusion on the
condition for the occurrence of transition. Physically, bulk diffusion is a
dissipative process that can lower the interfacial energy and lead to the
breakup of a liquid thread Xu2019; huang2022diffuse. Therefore, adding bulk
diffusion to the system will facilitate the pinch-off dynamics and hence
hinder the development of jetting state. As a result, a larger critical
velocity $\bar{V}_{\rm in}$ is needed to induce the dripping-to-jetting
transition. In an earlier work huang2022diffuse, we demonstrate that the
effect of bulk diffusion can be enhanced by increasing the characteristic
length scale $l_{c}$, which enters into the dimensionless system through the
parameter $B=\frac{2l_{c}^{2}}{L^{2}}$. Figure 8 shows that at different
levels of bulk diffusion controlled by $B$, $\mathcal{W}_{\rm in}$ and
$\mathcal{C}_{\rm in}$ at the transition always satisfy a linear relation for
transitions dominated by inner flow. Two important observations are made from
Figure 8. (i) Stronger diffusion indeed necessitates a larger critical
velocity $\bar{V}_{\rm in}$ to induce the transition. (ii) The three fitting
lines are parallel, showing that the relative contributions of the inertial
force and viscous force due to the inner flow remain the same regardless of
the variation of bulk diffusion.

Figure 8:  The Weber number of the inner flow $\mathcal{W}_{\rm in}$ and the
capillary number of the inner flow $\mathcal{C}_{\rm in}$ at the transition
always satisfy a linear relation for transitions dominated by inner flows. In
addition, the three fitting lines for three different values of $B$ are
parallel. The data are obtained by using $\varepsilon=0.01$, $\bar{R}=0.1$,
$\bar{H}=4$, $\bar{b}=0.1$, with $Re_{\gamma}$ being varied between $100$ and
$2000$, and $B=0.0002$, $0.0005$, and $0.00075$. The value of $Re_{\gamma}$
decreases from the upper left to the lower right along each line. Note that
$\bar{b}=0.1$ used here gives $\mathcal{C}_{\rm out}=0.029$, which is much
smaller than the typical values of $\mathcal{W}_{\rm in}$ and
$\mathcal{C}_{\rm in}$ at the transition.

##  IV  Concluding remarks

The CHNS model has been solved in a cylindrical domain with axisymmetry to
investigate the dripping-to-jetting transitions in coaxial flows of two
immiscible fluids. Numerous numerical examples are presented to demonstrate
that the distance between the orifice and pinch-off position increases when
either the outer or the inner flow rate is enhanced. It is observed that there
is an apparent jump in this distance when the outer or the inner flow rate
reaches the critical value for the dripping-to-jetting transition to occur.
The critical flow rates numerically obtained for both the outer and inner
flows are consistent with the corresponding experimental results in order of
magnitude. For transitions dominated by outer flows, a thin and long jet is
generated when jetting occurs, and our numerical results for the jet radius
are validated by its dependence on the outer flow rate according to the mass
conservation. For transitions dominated by inner flows, the interfacial
tension is varied for nearly two orders of magnitude, and a quantitative
relation is established between the contributions of the inertial and viscous
forces due to the inner flow at the transition point. Finally, the degree of
bulk diffusion is varied to show its quantitative effect on the critical flow
rate at the transition point.

To the best of our knowledge, there has been no prior work that employs a
phase-field model to investigate the dripping-to-jetting transitions in three
dimensions, with a focus on the effects of low interfacial tension and bulk
diffusion. In the present work, we have considered the simplest situation in
which the two fluids have equal density, equal viscosity and equal diffusion
coefficient. Actually, these restrictions can be lifted in both experiments
utada2007dripping; mak2017dripping; Xu2019 and numerical simulations dong2012.
Although the dripping-to-jetting transitions for high interfacial tension have
been extensively studied in the past two decades, low interfacial tension and
bulk diffusion may inject new ingredients into this classical problem. In this
regard, quantitative effects of density ratio, viscosity ratio and diffusivity
ratio largely remain to be explored in both experiments and numerical
simulations.

###### Acknowledgements.

The work of F. Huang and W. Bao was supported by the Ministry of Education of
Singapore under its AcRF Tier 2 funding MOE-T2EP20122-0002(A-8000962-00-00),
and the work of T. Qian was supported by the Hong Kong RGC grants CRF No.
C1006-20WF and GRF No. 16306121. T. Qian was also supported by the Key Project
of the National Natural Science Foundation of China (No. 12131010). Part of
this work was done when the first two authors were visiting the Institute for
Mathematics Sciences at the National University of Singapore in February 2023.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References