Effective interfacial tension in flow-focusing of colloidal dispersions: 3D numerical simulations and experiments

TL;DR

This study combines 3D numerical simulations and experiments to investigate how an effective interfacial tension influences flow dynamics in microfluidic flow-focusing of colloidal dispersions, revealing its impact on thread formation and providing measurement methods.

Contribution

It introduces a combined numerical and experimental approach to quantify and analyze the effective interfacial tension in colloidal dispersion flows, a novel application in microfluidics.

Findings

Simulations closely match experimental flow and thread geometries.

Effective interfacial tension significantly affects threading dynamics.

Methodology enables measurement of ultra-low interfacial tensions.

Abstract

An interface between two miscible fluids is transient, existing as a non-equilibrium state before complete molecular mixing is reached. However, during the existence of such an interface, which typically be at short timescales, composition gradients at the boundary between the two liquids cause stresses effectively mimicking an interfacial tension. Here, we combine numerical modelling and experiments to study the influence of an effective interfacial tension between a colloidal fibre dispersion and its own solvent on the flow in a microfluidic system. In a flow-focusing channel, the dispersion is injected as core flow that is hydrodynamically focused by its solvent as sheath flows. This leads to the formation of a long fluid thread, which is characterised in 3D using Optical Coherence Tomography and simulated using a volume of fluid method. The simulated flow and thread geometries very closely reproduce the experimental results in terms of thread topology and velocity flow fields. By varying the effective interfacial tension numerically, we show that it clearly influences threading dynamics and that it can be described by an effective capillary number. Furthermore, we demonstrate that the applied methodology provide means to measure the ultra-low but dynamically highly significant effective interfacial tension.

Figures

Equations

$$Ca=\frac{\eta U}{\gamma },$$

$$Pe=\frac{T_d}{T_c}=\frac{hU}{D},$$

$$\gamma =K\frac{\Delta {\Phi }^2}{\delta },$$

$${\eta }_{\text{eff}}={\eta }_{\text{inf}}+({\eta }_0-{\eta }_{\text{inf}})[1+(\tau \dot{\gamma }{)}^2{]}^{\frac{n-1}{2}},$$

$$D=\frac{k_{B}T\ln (l/d)}{2\pi \eta l},$$

$$\frac{\partial{\rho_{b}\mbox{\boldmath$U$}}}{\partial t}+\mbox{\boldmath$\nabla$}\cdot(\rho_{b}\mbox{\boldmath$U$}{\mbox{\boldmath$U$}})=-{\mbox{\boldmath$\nabla$}}{p}+{\mbox{\boldmath$\nabla$}}\cdot{\mbox{\boldmath$T$}}+\rho_{b}{\mbox{\boldmath$g$}}+{\mbox{\boldmath$F_{s}$}},$$

$$\mbox{\boldmath$\nabla$}\cdot\,\mbox{\boldmath$U$}=0,$$

$$\frac{\partial\alpha}{\partial t}+\mbox{\boldmath$\nabla$}\cdot(\alpha\mbox{\boldmath$U$})=0.$$

$${\rho }_b={\rho }_1\alpha +{\rho }_2(1-\alpha ),$$

$${\mu }_b={\mu }_1\alpha +{\mu }_2(1-\alpha ),$$

$${\mbox{\boldmath$F_{s}$}=\gamma\kappa(\mbox{\boldmath$\nabla$}\alpha)},$$

$${\kappa=\mbox{\boldmath$\nabla$}\cdot\bigg{(}\frac{\mbox{\boldmath$\nabla$}\alpha}{|\mbox{\boldmath$\nabla$}\alpha|}\bigg{)}}.$$

$$\partial \alpha /\partial n=0,$$

$${\delta }_{{\epsilon }_z}(\gamma )=\left|\frac{{\epsilon }_z^{num}(\gamma )-{\epsilon }_z^{exp}}{{\epsilon }_z^{exp}}\right|\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}{\delta }_{L_w}(\gamma )=\left|\frac{L_w^{num}(\gamma )-L_w^{exp}}{L_w^{exp}}\right|.$$

$$\tau =\frac{{\eta }_1}{\delta P}\propto \frac{{\eta }_1}{\gamma },$$

$$x_c-L_w=U\tau \propto \frac{{\eta }_1{Q}_1}{\gamma h^2}.$$

$$U_{corr}(y,z)=U_{meas}(y,z)\frac{Q_1}{{\iint }_{\Sigma }{U}_{meas}(y,z)\text{d}y\text{d}z},$$

Full text

1 Introduction

Microfluidic techniques have facilitated the development of both fundamental and applied research in the field of chemistry, biology, medicine, material and physical sciences (Squires & Quake, 2005; Whitesides, 2006; Stone et al., 2004; Nunes et al., 2013). In particular, the flow-focusing configuration (illustrated in figure 1), where a core fluid is focused by two sheath flows has been successfully employed in a wide variety of applications such as bubble or droplet formation (Cubaud & Mason, 2008; Anna et al., 2003; Garstecki et al., 2004), micro and nano-particle production (Martín-Banderas et al., 2005), hydrodynamic assembly of nanoparticle dispersion into high-performance superstructures (Håkansson et al., 2014; Ekanem et al., 2015; Kamada et al., 2017; Mittal et al., 2018), cell patterning (Takayama et al., 1999), DNA stretching (Wong et al., 2003), and diffusion-mixers (Knight, 1998; Pollack et al., 1999). At the micrometre scale, phenomena such as interfacial tension and viscosity usually become dominant compared to the effects of gravity and inertia, which often are negligible. This accords unique characteristics to utilize multiphase flows in microfluidic devices and offers the possibility to thoroughly investigate the influence of fluid physiochemical properties on the flow.

For a pair of immiscible fluids, surface or interface properties are controlled by interfacial tension. The relevant dimensionless number controlling the ratio between viscous forces and interfacial tension is the capillary number

[TABLE]

where $U$ is the typical velocity of the flow, $\eta$ the kinematic viscosity, and $\gamma$ the interfacial tension between the two fluids. The capillary number acts as the most important dimensionless number for dynamics of bubbles or droplets in microfluidics (Thorsen et al., 2001; Anna et al., 2003). In the case of flow-focusing, the capillary number is the primary dimensionless number controlling the transition between different flow patterns denoted threading, tubing, dripping and jetting (Cubaud & Mason, 2008). At very low capillary numbers, in the dripping regime, elongated droplets of length larger than the channel width are formed. At large capillary numbers, the flow is in the threading regime where there is a formation of a continuous fluid thread. For intermediate capillary numbers, the thread can break up to form dripping and jetting flow patterns depending on the capillary instability mechanisms involved (Humphry et al., 2009; Utada et al., 2007). These flow patterns demonstrate the rich dynamics of flow features typically observed in a microfluidic device. Moreover, much of the effort in the last two decades has been devoted to understand the influence of system parameters like flow rates of two fluid phases, device geometry, surface chemistry, and material parameters, such as their viscosities, densities and interfacial tension, on the formation and motion of these flow patterns (Nunes et al., 2013; Anna, 2016).

When a pair of miscible fluids are in contact, there is a transient interface whose thickness is governed by the diffusion coefficient $D$. Indeed, as microfluidic flows are mainly laminar, the fluid streams remain parallel and the liquids mix only by diffusion. One can then define the diffusion time scale $T_{d}=h^{2}/D$ and the convective time scale $T_{c}=h/U$, where $h$ is the characteristic length scale of the system. These time scales can be compared using the dimensionless Péclet number

[TABLE]

which is the relevant dimensionless number to describe convective diffusive flows (Cubaud & Mason, 2006; Atencia & Beebe, 2004). When miscible fluids are subjected to flow-focusing, at high Péclet number, the main flow regime observed is a threading regime with common features similar to immiscible fluid threads while, at low Péclet number, the thread can be destabilized by diffusion instabilities (Cubaud & Notaro, 2014). However, as the existence of such thread structures is transient and limited to short timescales before complete mixing, little is still known about the dynamics of viscous threads in miscible environments.

From the above, it could be expected that the dynamics of immiscible and miscible fluids in microfluidic environments are controlled by the capillary and Péclet numbers, respectively. However, miscible fluids flowing under laminar flow conditions and high Péclet number have a boundary created by so called de facto dynamic interface, which gradually changes with time (Joseph, 1990). Such an interface may resemble a distinct interface present between two immiscible fluids (Garik et al., 1991; Anderson et al., 1998). The de facto dynamic interface arises due to non-equilibrium capillary forces generated by the interfacial stresses whenever gradients of a fluid property exist at the interface of two miscible fluids (Korteweg, 1901). Recently, Truzzolillo et al. (2014) coined the term effective interfacial tension (EIT) for this effect, and defined it as

[TABLE]

where $K$ is the Korteweg constant of the system, $\delta$ is the thickness of the interface and $\Delta\Phi$ is the change in the concentration or volume fraction $\Phi$. They noted that despite Korteweg’s theory being a century old, experimental investigations are limited to molecular fluids. They carried out experimental measurements of the Saffman-Taylor instability in a Hele-Shaw cell with miscible complex fluids including colloidal dispersions. Their measurements went beyond confirming Korteweg’s theory and showed that EIT is also a function of particle structure (basically size and shape) together with particle-particle and particle-solvent interactions (Truzzolillo et al., 2016).

Effective interfacial tension is expected to have relevance in many fields starting from material processing to multiphase complex fluid dynamic problems involving droplet and bubbles formation, jetting, coalescence and break-up of droplets (Truzzolillo & Cipelletti, 2017). Even though progress has been made on the understanding of the above aspects, transient mechanisms and their effect on the formation of the above flow patterns are still unexplored. To the best of our knowledge, there are no reports detailing the significance of EIT on the flow physics in microfluidic systems. Moreover, the recent review by Anna (2016) observed that many of the advancements in the field of microfluidics have been driven predominantly by experiments due to the complexity of geometries and fluid physics. Thus, there is a demand for numerical simulations that can provide insights on complex mechanisms, and to extract the variables that are difficult or impossible to measure in experiments.

We aim at addressing the influence of EIT and its effects on the flow-focusing of a colloidal dispersion by its own solvent. Detailed 3D experimental and numerical investigations have been performed in a microfluidic flow-focusing channel of square cross-section using Optical Coherence Tomography and Volume of Fluid method. The colloidal dispersion is formed by cellulose nanofibrils (CNF) dispersed in water. Such a colloidal fibre dispersion has been used to assemble very strong cellulose filaments through hydrodynamic focusing (Håkansson et al., 2014; Mittal et al., 2018). Furthermore,  Truzzolillo & Cipelletti (2017) noted that for such a system, EIT is expected to exist between the colloidal dispersion and its solvent. Here, we quantitatively investigate the characteristics of thread formation such as 3D topology and velocity fields as a function of EIT and flow-rate ratio. The previous experimental studies by Cubaud & Mason (2008, 2009) and Cubaud & Notaro (2014) for both immiscible and miscible high-viscosity threads in flow-focusing channels were limited to a top view and did not resolve the full 3D shape of the thread. Indeed, there are only a few 3D studies of the threading regime (Knight, 1998; Xu et al., 2016) but their focus are not on the effect of interfacial properties on the thread dynamics. Thus, at large Péclet number, the role of EIT potentially acting during the existence of a transient interface between two miscible fluids is not understood so far due to lack of 3D experimental and numerical studies.

This paper is organized as follows. In Section 2, experimental and numerical setups are presented together with the numerical validation. In Section 3, we compare the results of 3D computations with the experimental measurements in the threading regime. In particular, the focus is on the evolution of 3D thread shape, and morphology of the region wetted by the colloidal dispersion on the top and bottom walls. The effects of flow-rate ratio on variation in thread shape is investigated. We also examine numerically the effect of EIT on the dynamics of thread formation, allowing us to estimate the value of EIT for the present colloidal dispersion system. In Section 4, we look at the velocity fields both along the centreline as well as in cross-sectional planes of the channel. Finally, in Section 5, we present the conclusions.

2 Methodology

2.1 Experimental set-up

2.1.1 Flow cell

In the present work, we use a flow-focusing channel system with square cross-sections having $h=1$ mm sides as shown schematically in figure 1. The channel system comprises four channels, with three inlet channels intersecting to form a cross-junction, and a main outlet channel for the fluids to flow out. A Cartesian coordinate system ($x,y,z$) is introduced, with the $x$ and $y$ axes in directions of the core and sheath flow inlets, respectively. The $z$ axis is perpendicular to the $(x,y)$-plane representing the third direction. Here, a colloidal dispersion in the core flow is focused by two water sheath flows. The flow-focusing channel is fabricated out of a $1$ mm thick stainless-steel plate (Håkansson et al., 2014). The steel plate is covered on both sides with a $140\leavevmode\nobreak\ \mu$m Cyclic Olefin Copolymer (COC) film providing optical access. The COC - channel plate - COC sandwich is held together by two aluminium plates with screws. Two syringe pumps (WPI, AI-$4000$) were used to drive the core and sheath flows at constant volumetric flow rates. In most of the experiments, flow rates are $Q_{1}$ = 6.5 mm${}^{3}\cdot$s*-1* for the colloidal dispersion and $Q_{2}$ =  $7.5$ mm${}^{3}\cdot$s*-1* for the sheath flows, respectively. The ratio of these flow rates is defined as $\phi=Q_{1}/Q_{2}$ and is set to $\phi=0.8667$ in the standard case, except when explicitly varied in Section 3.2. The indices 1 and 2 will be used to denote the core and sheath fluid/flow respectively, throughout the paper.

2.1.2 Colloidal dispersion

The colloidal dispersion used as the core fluid is cellulose nanofibrils (CNF) suspended in water. Cellulose nanofibrils were prepared by liberating fibrils from bleached softwood fibres (Domsjö Fabriker AB, Sweden). Before liberation, the fibrils were carboxymethylated (Wågberg et al., 1987) to a degree of substitution of $0.1$. The fibrils were then liberated from the fibre wall following a protocol described in Fall et al. (2011). Post-liberation unfibrillated fibre fragments were removed by centrifuging the dispersions at 4500 rpm. A transparent colloidal dispersion of concentration $3$ g$\cdot$dm*-3* was obtained with typical fibril length $l$ of $700$ nm and fibril diameter $d$ of $2$ to $3$ nm.

The rheological characterization of the colloidal dispersion was carried out using a Kinexus pro+ rheometer (Malvern) with a modified Couette geometry, suitable for accurate measurement of viscosity at low shear rates. The shear viscosity behaviour of the colloidal dispersion is indicated by red dots as shown in figure 2. It exhibits a classical shear thinning behaviour at moderate and high shear rates, and constant viscosity at low shear rates (Martoïa et al., 2016; Nechyporchuk et al., 2016). The rheological data can be fitted well by a Carreau model

[TABLE]

where $\eta_{\textrm{{eff}}}$ is the shear viscosity, $\dot{\gamma}$ the shear rate, $\eta_{0}$ the zero shear viscosity, $\eta_{\textrm{{inf}}}$ the infinite shear viscosity, $\tau$ the relaxation time, and $\mathit{n}$ the power index. The best fit is shown by the black solid line in figure 2 and gives $\eta_{\textrm{{inf}}}=5$ mPa$\cdot$s, $\eta_{0}=1756$ mPa$\cdot$s, $\tau=16.16$ s, and $n=0.56$. The viscosity ratio between the core and sheath flow defined as $\chi=\eta_{1}/\eta_{2}=1756$ is arrived by considering the viscosity of colloidal dispersion at low shear rate $\eta_{1}$ = $\eta_{0}=1756$ mPa$\cdot$s and viscosity of sheath fluid $\eta_{2}=1$ mPa$\cdot$s, respectively.

2.1.3 System time scales

Due to their dimensions, the fibrils are Brownian and diffuse in the solvent. The diffusion coefficient $D$ is given by (Doi & Edwards, 1986)

[TABLE]

with the temperature $T=300$ K, the Boltzmann constant $k_{B}=1.38\times 10^{-23}$ J$\cdot$K*-1*, and the solvent viscosity $\eta=1$ mPa$\cdot$s. Plugging in the all the values, this leads to $D\approx 5\times 10^{-12}$ m${}^{2}\cdot$s*-1* and to a typical diffusion time scale in the system $T_{d}\approx 2\times 10^{5}$ s. With a velocity of the flow $U\approx 10$ mm$\cdot$s*-1*, the convective time scale can be estimated to be $T_{c}\approx 10^{-1}$ s. Using these values and equation (2), we get $Pe\approx 2\times 10^{6}$. The high Péclet number implies that the system is predominantly convective driven and diffusion is very weak. Therefore, it is possible for a de facto dynamic interface to exist between the two fluids (Joseph, 1990; Korteweg, 1901). The properties of a de facto dynamic interface may resemble a distinct interface between two immiscible fluids and thus, the two fluids, the colloidal dispersion and water, exhibit a transient immiscible behaviour. In the following two sections, the experimental method and numerical setup will be presented. Experimental method involves more about the description of optical coherence tomography (OCT) employed for topology and velocity measurements.

2.1.4 Optical Coherence Tomography

Huang et al. (1991) first demonstrated the three-dimensional non-invasive imaging technique known as Optical Coherence Tomography (OCT). Based on low-coherence interferometry, OCT is used to produce an image of optical scattering within transparent or turbid millimetric samples with a few micrometre resolutions. Usually designed for biomedical applications, the characteristics of the OCT, namely high-resolution, depth-resolved information, contact-free imaging and velocity acquisition, make it a pertinent measurement technique for a large range of other applications (Stifter, 2007), e.g. material characterisation and microfluidics (Xi et al., 2004; Ahn et al., 2008).

Since optical coherence tomography might be new to the reader, a brief description is in place. The interferometric technique used in OCT is based on the Michelson interferometer. Its principle is sketched in figure 3. A broadband light source is divided into two parts by a beam splitter. One part is sent in a reference arm while the other goes through the sample. The reflected light from the sample is then recombined in the beam splitter with the light from the reference arm, providing an interference signal. Interferences are only observed when the reference and sample arm optical path lengths are matched to within the coherence length of the light (Tomlins & Wang, 2005). Therefore, the depth resolution of the system is determined by the coherence of the light source. A full interference pattern, i.e. interferences as a function of depth, can be obtained by analysing the output spectrum, thanks to the spectrometer. In the interference pattern, refractive index variations or differences of optical scattering between different layers in the sample correspond to intensity peaks or to different magnitudes of intensity. There, the interference pattern contains information about the depth and the optical scattering within the sample, with a micrometric resolution. Moreover, the depth and lateral resolutions are decoupled.

This imaging technique can be combined with a Doppler acquisition when there is a flow in the sample (Chen et al., 1997). The backscattered light from moving particles experiences a double Doppler shift (one from the source, one from the particle back to the probe). Thereby, the velocity of moving particles can be measured by using the Doppler shift and the relative angle between the optical beam and the direction of the flow. Again, the depth resolution depends on the coherence length of the light source and is decoupled from the lateral resolution. The velocity resolution depends on the acquisition time and the scan angle. The range of velocities that can be measured is wide, varying from several micrometres per second to tens of millimetre per second.

The OCT apparatus used in this study is a commercial spectral domain OCT (SD-OCT) Telesto II from Thorlabs. The light source has a central wavelength of $1310$ nm and a bandwidth of $270$ nm, giving a resolution in both the axial and transverse directions of 3 $\mu$m. Here, OCT together with Doppler acquisition is employed for 3D measurements of thread topology and velocity fields. More details on the OCT acquisition and data processing can be found in Appendix A.

2.2 Numerical set-up

2.2.1 Numerical method

The numerical simulation of multiphase flows with fluid-fluid interfaces faces two major challenges. First, the accurate estimation and advection of the complex interface between two incompressible fluids. Second, suppression of spurious currents/velocities arising from inaccurate calculations of interface curvature. Thus, the demand for highly resolved simulations in terms of spatial resolution at the liquid-liquid interface has become pertinent to capture a sharp interface and analyse the physics of such flows.

Volume of Fluid (VOF) (Hirt & Nichols, 1981) and level set methods(Osher & Sethian, 1988), have received significant attention and are the most popular state-of-the-art tools in simulating multiphase flow problems involving extensive interface movement and topological changes. Here, we employ the VOF methodology, where the interface is implicitly represented using an indicator function, a fluid fraction parameter $\alpha$ via the volume fractions of the two fluids in the computational cells; $\alpha$ = 0 correspond to fluid 1, $\alpha$ = 1 to fluid 2, with 0 $<$ $\alpha$ $<$ 1 in the interface region. The advection of the interface is achieved through the redistribution of $\alpha$ between neighbouring cells. Various VOF methods have been developed since its inception, and are broadly categorised into algebraic and geometric methods. Algebraic VOF methods are developed for general mesh types, easier to implement and faster (Deshpande et al., 2012). However, they are less accurate and often require higher order schemes to retain the sharpness and boundedness of fluid fraction $\alpha$. Geometric VOF methods, on the other hand, use geometric operations to reconstruct the interface inside a computational cell, and are more accurate, but also complex to implement and the execution is relatively slow. In the present work, we have chosen a new geometric VOF approach called IsoAdvector (Roenby et al., 2016) which combines high accuracy with advection of interface.

The full system of equations being solved consists of the Navier-Stokes equation

[TABLE]

the continuity equation

[TABLE]

and the equation for the advection of fluid fraction $\alpha$

[TABLE]

Here U is the velocity vector field, $\rho_{b}$ the density, $p$ the pressure field, T the deviatoric stress tensor (\mbox{\boldmathT}=2\mu_{b}\mbox{\boldmathS}-2\mu_{b}(\mbox{\boldmath\nabla}\cdot\mbox{\boldmathU})\textit{{I}}/3) , where S= 0.5(\mbox{\boldmath\nabla}\mbox{\boldmathU}+\mbox{\boldmath\nabla}\mbox{\boldmathU}^{T}), I is the identity matrix and g is the gravitational force. The bulk parameters density $\rho_{b}$ and viscosity $\mu_{b}$ are computed based on the weighted average distribution of the fluid fraction $\alpha$

[TABLE]

[TABLE]

where $\rho_{1}$, $\rho_{2}$, $\mu_{1}$, ${\mu_{2}}$ are the densities and the viscosities of the two phases. $F_{s}$ is the body force per unit mass, which accounts for the surface tension forces present only at the interface between two fluids. The surface tension force $F_{s}$ is modelled as a volumetric force using Continuum Surface Force (CSF) method (Brackbill et al., 1992), and is defined as

[TABLE]

where $\gamma$ is the interfacial tension and $\kappa$ is the curvature of the interface

[TABLE]

The numerical simulations were performed with the finite-volume-based open-source code OpenFOAM (Weller et al., 1998). The PISO (pressure-implicit with splitting of operators) scheme (Rusche, 2003) is applied for solving the momentum balance equation (6) in conjunction with the continuity equation (7). More detailed information about the implementation of the above equations in the OpenFOAM framework can be found elsewhere (Deshpande et al., 2012; Nekouei & Vanapalli, 2017). A first order implicit Euler scheme was used for transient terms, and the time step was controlled by setting the maximum Courant number to 0.3. A second order Total Variation Diminishing (TVD) scheme with van Leer limiter was used for the spatial discretization. Boundedness of the fluid phase fraction $\alpha$, was controlled by the IsoAdvector (Roenby et al., 2016) scheme implemented in the interIsoFoam solver of OpenFOAM-v1706 (OpenCFD, 2017).

Our system has three inlets, two for the sheath flows and one for the core flow, and one outlet as depicted in figure 1. The length of three inlets and outlet channel correspond to $10$ mm and $30$ mm with square cross-section $1\times 1$ mm2, respectively. The non-Newtonian rheology of the core fluid has been implemented numerically using the Carreau model depicted in figure 2. We assume that the viscosity $\eta_{\textrm{{eff}}}$ is sufficient to describe the relevant material properties of the colloidal dispersion. To initialize the simulation, we set the fluid phase fraction, $\alpha=0$ in the core flow inlet channel, $\alpha=1$ in the sheath flow inlet channel and for the outlet channel. A uniform velocity was set at the inlets calculated as per the flow rates mentioned in Section 2.1.1. At the walls, we apply the no-slip boundary condition and use an equilibrium contact angle $\theta=0^{\circ}$ considering full wetting by the sheath fluid. The contact angle boundary condition was employed to correct the surface normal vector, and thereby adjusts the curvature of the interface in the vicinity of the wall. At the outlet, we prescribe a reference pressure (atmospheric) and zero gradient of the volume fraction

[TABLE]

where $n$ is the unit vector normal to the wall.

2.2.2 Validation and convergence

The validation of the numerical scheme is performed by simulating a reference case, which also demonstrates the rich dynamics of flow in the flow-focusing geometry. The Cubaud & Mason (2008) experiments involve two immiscible Newtonian fluids subjected to hydrodynamic focusing achieved through a square microfluidic flow-focusing channel of side, with $h=100\leavevmode\nobreak\ \mu$m. They have shown that the different regimes can be represented in a capillary number based flow map. Their experiments have been reproduced numerically here using glycerol and PDMS oil as core ($Q$1) and sheath ($Q$2) flows. These fluids have viscosities of $\eta_{1}=1214$ mPa$\cdot$s, and $\eta_{2}=4.59$ mPa$\cdot$s, respectively. The interfacial tension between the two fluids is $\gamma_{12}=27.0$ mN$\cdot$m*-1*. The capillary numbers of core and sheath flows are defined as $Ca_{1}$ = ${\eta_{1}Q_{1}}/{(\gamma_{12}h^{2})}$ and $Ca_{2}$ = ${\eta_{2}Q_{2}}/{(\gamma_{12}h^{2})}$, respectively. Figure 4 (a) shows the typical flow patterns viz., (i) threading, (ii) jetting, (iii) dripping, and (iv) tubing obtained through our 3D simulations. The state space in figure 4 (b) shows that the agreement between flow regimes observed by Cubaud & Mason (2008) and in our simulations is very good. In particular, the transition between threading and tubing (from circles to diamonds) is well captured. In addition, the agreement with the experiments is further confirmed quantitatively by measuring the width $\varepsilon_{y}/h$ and height $\varepsilon_{z}/h$ of the thread in the (i) threading regime as a function of the flow rate-ratio $\phi$ as depicted in figure 4 (c). Indeed, the thread cross-section is circular and scales as $\varepsilon_{y}/h=\varepsilon_{z}/h=\leavevmode\nobreak\ (\phi/2)^{1/2}$, observed experimentally and expected theoretically by Cubaud & Mason (2008, 2009). The pink solid line represents the fit of $\varepsilon_{y}/h=\varepsilon_{z}/h=(\phi/2)^{1/2}$.

All the validation studies and upcoming results of the present study in Sections 3 and 4 were obtained at a resolution of $h/\Delta=50$, where $\Delta$ is the grid size, i.e., $50\times 50$ cells in a square channel cross-section. Figure 5 shows the sensitivity of the grid size $\Delta$ on width $\varepsilon_{y}/h$, height $\varepsilon_{z}/h$ of the thread, and on the flow field velocity $\mathit{U}$. We observe that for $h/\Delta=50$, the variation with respect to the finest resolution $h/\Delta=80$ in the width ($\delta_{\varepsilon_{y}}$), height ($\delta_{\varepsilon_{z}}$) of the thread, and on the flow field velocity ($\delta_{U}$) is smaller than $2\%$, while for $h/\Delta>60$ it is smaller than $1\%$. It will be seen that $h/\Delta=50$ is sufficient to capture the physics of flow with desired quality both in quantitative and qualitative forms. Typical clock time for the simulations range from $24$ to $48$ hours using $64$ to $256$ processors going from coarse to fine grid.

3 Thread shape

In this section, we focus on the shape of the thread through a detailed comparison between experimental measurements and numerical simulations. In particular, the influence of flow-rate ratio $\phi$ and effective interfacial tension $\gamma$ on thread formation is investigated. The flow-rate ratio and the effective interfacial tension (EIT) for the standard case are respectively set to $\phi=0.8667$ and $\gamma\leavevmode\nobreak\ =\leavevmode\nobreak\ 0.054$ mN$\cdot$m*-1* in all simulations, except when explicitly varied in Sections 3.2 and  3.3.

3.1 Shape of the thread in the standard case

The 3D shapes of the thread region occupied by the colloidal dispersion measured using OCT and obtained by numerical simulations are shown in figure 6. The qualitative agreement is very good (quantitative comparisons will follow later). For $x/h<0$, the colloidal dispersion attach to the channel walls. The beginning of focusing region where the entry of sheath flows commence corresponds to $x/h=0$. By inspecting figure 6 carefully at the position where the sheath flows are injected (0 $\leq$ $x/h$ $\leq$ 1), it can be noted that the colloidal dispersion does not lose contact with the walls located at $z/h=\pm 0.5$ immediately. Instead, there are regions wetted by the colloidal dispersion before detaching completely from the top and bottom walls as highlighted by cyan colour in figure 6. Further downstream at around $x/h=2$, the cross-section of the colloidal dispersion thread forms a nearly elliptically shaped thread. The very good agreement between numerical and experimental results is again demonstrated in figure 7 where the experimentally and numerically obtained morphologies of the region wetted by the colloidal dispersion on the wall at $z/h=0.5$ are compared. The distance from $x/h=0$ to the point where the colloidal dispersion detaches from the wall is named wetted length and is denoted by $L_{w}/h$. All features such as wetted length $L_{w}/h$ and the morphology of the wetted region observed in experiments are very well captured and reproduced by the numerical simulations.

The cross-sectional views of the thread shape at various downstream positions $x/h=1.5$, $x/h=3$ and $x/h=10$ are illustrated in figure 8, and figure 9 depicts the length of the ellipsoid axes $\varepsilon_{y}/h$, $\varepsilon_{z}/h$ and aspect ratio $\varepsilon_{z}/\varepsilon_{y}$ of the cross-section as a function of downstream positions $x/h$. The agreement between the numerical and experimental measurements is superb. One could also notice that at $x/h=1.5$, the colloidal dispersion is still being in contact with the top and bottom walls (figure 8 (a),(d)) and then after the detachment, the thread becomes thinner and stays nearly stable (figure 8 (b),(e) and (c),(f)). This is reflected in the streamwise development of the height $\varepsilon_{z}/h$ and the width $\varepsilon_{y}/h$ of the thread (see figure 9). First, the width decreases much faster than the height and also reaches a constant value at around $x/h\approx 5$ whereas there is a very slow change in the height of thread. The values of width and height reached far downstream ($\varepsilon_{y}/h\approx 0.5$ and $\varepsilon_{z}/h\approx 0.8$ at $x/h=20$) are different, confirming that the cross-section of the thread is nearly elliptical as visible in figure 6. These clear differences indicate that the thread formation in the $y$ and $z$ directions could be governed by different physical mechanisms.

3.2 Shape variations with flow rate ratio

In order to investigate the role played by the sheath flows on the thread formation, a parametric variation of the flow-rate ratio $\phi$ over two orders of magnitude has been performed. This is an aspect that is of immense interest when the flow-focusing is used for hydrodynamic assembly of nanofibrils (Håkansson et al., 2014), since the shape determines the cross-section of an assembled material and could also provide means for synthesizing materials of complex shapes such as rods, discs, ellipsoids (Xu et al., 2005). The dimensionless width $\varepsilon_{y}/h$ and height $\varepsilon_{z}/h$ of the thread measured far downstream at  $x/h=20$ are plotted in log-log scale in figure 10 as a function of flow-rate ratio $\phi$. The agreement between experimental and numerical points is very good over the full range of the flow-rate ratio variation. Both curves show two different regimes: for low flow-rate ratio, there is a power law regime $\alpha\phi^{k}$ corresponding to the threading regime, while for high flow-rate ratio, there is saturation at $1$ corresponding to the tubing regime. The fit of the threading regime as a function of flow-rate ratio $\phi$, represented by pink solid lines in figure 10, gives $\varepsilon_{y}/h=0.58\phi^{0.5}$ for the width and $\varepsilon_{z}/h=0.89\phi^{0.28}$ for the height. A power law in $\phi^{0.5}$ is predicted for a circular thread only controlled by the ratio of the flow rates (Cubaud & Mason, 2009) and has been verified experimentally only by measuring the width of the thread assuming a circular cross section (Cubaud & Mason, 2008; Cubaud & Notaro, 2014). The coefficient $0.58$ in front of the power law is also in good agreement with previous studies (Cubaud & Mason, 2008, 2009). This implies that the width of the thread is purely controlled by the flow-rate ratio $\phi$. However, the height of the thread exhibits a different power law, in $\phi^{0.28}$. This indicates that the thread cross section is not circular and that the height of the thread is driven only partly by the flow-rate ratio $\phi$ but could also be dependent on other parameters such as viscosity ratio $\chi$ and interfacial tension $\gamma$. Later, in the upcoming Section 3.3, the effect of EIT ($\gamma$) on the height of thread $\varepsilon_{z}/h$ is discussed.

One can also notice that the transitions between the power law regimes happen at different flow-rate ratio values for the width and the height of the thread. There are therefore three different characteristic regimes that can are observed within the range of the flow-rate ratio sampled in this study. They are illustrated in figure 10 (c). For $\phi<1.5$, the colloidal dispersion forms a thread with an elliptical cross-section (i). In this regime, both width and height of the thread enlarge when increasing $\phi$. In the range $1.5<\phi<3$, the height of the thread saturates and reaches the size of the channel $h$ in $z$ direction while the width keeps increasing. This is the first regime of tubing (ii). Above $\phi>3$, both height and width have saturated and reaches the size of the channel $h$ in the $y$ and $z$ directions. This defines a second tubing regime (iii) where the sheath flows are only localized in the corners of the channel.

3.3 Effect of interfacial tension

The transient immiscible behaviour between the two fluids, the colloidal dispersion and its own solvent, provides an opportunity to study the effect and measure the value of so-called effective interfacial tension (EIT) (Truzzolillo et al., 2014) for a dynamic interface. Since EIT is a ultra-low and transient phenomena, it is difficult to measure experimentally. We utilize the numerical simulations to quantify values within the range $0.024$ mN$\cdot$m${}^{-1}<\leavevmode\nobreak\ \gamma<\leavevmode\nobreak\ 3.000$ mN$\cdot$m*-1*. The height of the thread $\varepsilon_{z}/h$ at $x/h=20$ together with the wetted length $L_{w}/h$ are now used to illustrate how EIT affects the shape in the calculations by comparing with the experimental values. This has been done using the following definitions for the differences,

[TABLE]

The differences $\delta_{\varepsilon_{z}}(\gamma)$ and $\delta_{L_{w}}(\gamma)$ are plotted in figure 11. Both curves show significant variations i.e., $\varepsilon_{z}/h$ and $L_{w}/h$ are controlled by $\gamma$ for a fixed viscosity ratio $\chi$ and flow-rate ratio $\phi$. The minima of $\delta_{\varepsilon_{z}}(\gamma)$ and $\delta_{L_{w}}(\gamma)$, are indicated by red filled symbols in figures 11 (a) and (b) and occur at $\gamma=0.054$ mN$\cdot$m*-1*. This value is used for the comparison of numerical results with the experimental measurements in this study. Moreover, the order of magnitude of $\gamma\sim 10^{-2}$ mN$\cdot$m*-1*, is in line with the recent experimental observations of EIT made by Truzzolillo et al. (2016) for miscible complex fluids like colloidal and polymer dispersions where there is a presence of dynamic interface due to weak diffusion. Thus, the presence of EIT in this study is a consequence of the high Péclet number leading to the existence of a de facto dynamic interface.

Figure 12 (a) shows the height of the thread $\varepsilon_{z}/h$ as a function of downstream position $x/h$ for various $\gamma$. First, one can note that the wetted length $L_{w}/h$ is slightly influenced by $\gamma$, as the detachment of the various threads from the top and bottom walls of the channel takes place in the range $x/h=1$ - $3$. Once the thread is detached, the larger the $\gamma$, the faster the thread converges towards a value of $0.63$, corresponding to the width of the thread $\varepsilon_{y}/h$ for the simulation with the highest $\gamma$. Far downstream, for $x/h>20$, the height of the threads with low $\gamma$ is still decreasing. This strongly suggests that the near-elliptical shape of the thread cross-section is only a transient state, before flattering towards circular shape. Indeed, the thread will eventually become circular as the effective interfacial tension tends to minimize the contact surface. However, the time scale for the thread to converge towards a circular cross-section appears to be dependent on the EIT $\gamma$.

For instance, let us consider the two-dimensional problem of the thread cross-section, in the $y$-$z$ plane as shown in Figure 8. Just after the detachment, at $x=L_{w}$, the cross-section has an near-elliptical shape, with $\varepsilon_{z}/h\approx 1$ and $\varepsilon_{y}/h<1$. The thread reorganizes itself by converging towards a circular cross-section, with $\varepsilon_{z}/h=\varepsilon_{y}/h<1$. This dynamics is driven by  $\gamma$ at the interface. Indeed, as the interface between the core and sheath fluids is curved, there is a Laplace pressure difference between the two fluids, proportional to the $\gamma$ and to the local curvature of the interface. As $\varepsilon_{z}/h>\varepsilon_{y}/h$ and assuming a constant pressure in the sheath flow, the pressure on the top and bottom of the thread is thus larger than the pressure on the sides. This in turn leads to a pressure gradient $\delta P$ within the thread, proportional to EIT ($\gamma$) and dependent on the geometry of the thread, driving the thread cross-section from an near-elliptical to a circular shape.

As the viscosity ratio $\chi$ in the system is very high, one can neglect the influence of the sheath flows on this dynamics. Consequently, it is mainly dependent on the viscosity $\eta_{1}$ of the core fluid. The typical time scale $\tau$ for this dynamics can be estimated via the Stokes equation as

[TABLE]

where the thread geometry dependence of $\delta P$ is hidden in the proportional symbol. During this time scale, the thread cross-section is advected downstream with a velocity $U$. As the flow rate is conserved, this velocity is dependent on the shape of the thread $U=4Q_{1}/(\pi\varepsilon_{z}\varepsilon_{y})$, but we neglect these variations for simplification and consider that $U\approx Q_{1}/h^{2}$. Accordingly, the typical distance $x_{c}$ needed for the thread to converge from near-elliptical to circular cross-section can be estimated by

[TABLE]

Similar to immiscible fluids, one can therefore define an effective capillary number for the thread, $Ca_{1}=\eta_{1}Q_{1}/(h^{2}\gamma)$, and normalize the modified downstream location $(x-L_{w})/h$ with it. This is shown in Figure 12 (b). All the numerical data collapse well on a master curve, showing that the dynamics of the height of the thread is governed by this effective capillary number. To measure the EIT acting between two miscible fluids with high Péclet and high viscosity ratio, it is therefore only necessary to measure the evolution of the height of the thread $\varepsilon_{z}/h$ formed by these two fluids and to compare it with the master curve as depicted in Figure 12 (b).

It is worthwhile to note that Cubaud & Mason (2009) made a circular prediction for a square-cross section channel of width $h$ assuming a circular annular flow, which means that the way to pour liquids in the channel was not taken into account. In a flow-focusing geometry, there is no rotational symmetry of $\frac{\pi}{2}$ as for a square duct. As a consequence, the thread only respects the two planes of symmetry of the channel at $y/h=0$ and $z/h=0$. The thread therefore detaches from the channel wall with an elliptical cross-section, thinner in the direction of the sheath inlets, i.e. the $y$ direction because the sheath has higher momentum in that direction.

4 Velocity fields

In addition to the comparison between numerical simulations and experiments with respect to shape and topology of the thread, we will also compare the velocity fields. These are responsible for the alignment of the nanofibrils in the dispersion (Jeffery, 1922) and are therefore of importance to material fabrication.

In figure 13, the centerline velocity $U_{c}$ is shown as a function of the downstream location $x/h$, obtained by both experimental measurements and numerical simulations. The velocity is first constant before the focusing region ($x/h<0$) and right after the injection of the sheath flows at $x/h=0$, there is a minor deceleration followed by rapid acceleration, before almost stabilizing at around $x/h=10$. The increasing velocity implies that there is an acceleration of the core flow, along the direction ($x$-direction) of the channel. Thus the flow-focusing geometry through sheath flows creates an extensional flow with a velocity gradient in the direction of the core flow, leading to the alignment of the nanofibrils in the colloidal dispersion (Håkansson et al., 2016).

Figure 14 shows experimental and numerical maps and profiles of the $x$-component of the velocity in different cross-sectional $x/h$ planes: before the focusing region ($x/h=\leavevmode\nobreak\ -3$), just after the focusing region ($x/h=1$) and at far downstream ($x/h=10$). Both the velocity maps and profiles show again an excellent quantitative agreement. Before the focusing point, the velocity profile is near parabolic. Just after the focusing, at the beginning of the acceleration, one can see that the colloidal dispersion is still attached to the walls at $z/h=\pm 0.5$, in the wetted region, but is surrounded by the sheath flows in the $y$-direction. The sheath flows are visible on the velocity maps as the region with the highest velocities on the sides. The thread exhibits a uniform velocity profile in the centre. Far downstream, the thread has attained its nearly elliptical shape and steady velocity flow-field. The velocity profile within the thread is completely uniform indicating plug flow. Surrounding the colloidal dispersion, the sheath flow velocity profiles are parabolic, connecting the thread to the walls of the channel. This has been predicted by Cubaud & Mason (2009) for a circular thread in a square channel.

5 Conclusions

In this paper, we present a comprehensive comparison of experimental and numerical results on 3D thread formation under the influence of EIT with a miscible fluid pair of colloidal dispersion and its own solvent as a model system. With a high Péclet number, this system is a weakly diffusive i.e., the timescale for diffusion of the nanofibrils from the core flow into the sheath flows being substantially larger than the residence time of the two fluids in the channel. Therefore, the system exhibits a near-immiscible behaviour with a de facto dynamic interface where EIT is acting. Three-dimensional versions of the finite volume based open-source code OpenFOAM (OpenCFD, 2017) with volume of fluid (VOF) method have been used for the numerical simulations and experimental measurements were carried out with Optical Coherence Tomography (OCT). The validation of the numerical scheme was performed by choosing and reproducing a reference case for immiscible fluids reported in Cubaud & Mason (2008) that demonstrates the rich dynamics of the flow features typically seen in a microfluidic flow-focusing configuration. Our simulations capture all the features of various flow patterns: threading, jetting, dripping and tubing regimes, in very good agreement with the experiments of Cubaud & Mason (2008) both in terms of flow visualisations and in the state space mapped based on capillary numbers of two fluids.

In the threading regime, the 3D numerical simulations of the present study correctly reproduce the experimental observations in terms of thread topology and morphology of the wetted region. The measurements and simulations of the cross-sectional width $\varepsilon_{y}/h$ and height $\varepsilon_{z}/h$ show the same behaviour where the thread width reaches a constant value faster than the height, and where the cross-sectional shape of the thread far downstream is nearly elliptic. The shape variation of width $\varepsilon_{y}/h$ and height $\varepsilon_{z}/h$ of the thread as a function of flow-rate ratio $\phi$ show good agreement between experiments and numerical simulations revealing interesting features: the width and height of thread both exhibit power-law behaviour with flow-rate ratio $\phi$. The width of the thread scales purely with the flow-rate ratio as $\phi^{0.5}$ in good agreement with previous studies (Cubaud & Mason, 2009, 2008). However, the height of thread scales as $\phi^{0.28}$ confirming the nearly elliptical cross-sectional shape, thus indicating the dependence of height on other parameters like EIT along with the flow-rate ratio $\phi$.

The three-dimensional numerical simulations gave access to probe for a detailed investigation of the influence of EIT by comparing with the elliptical thread shape given by experiments. Thus, it provided a framework to deduce the values of EIT by minimizing the differences between simulations and experiments. Effective interfacial tension was determined by varying the surface tension in the simulations and comparing the numerically generated data with the experimental wetted length $L_{w}/h$ and height $\varepsilon_{z}/h$ of the thread. The minimum differences are found to be very small, typically about less than 1%. Furthermore, by estimating the time and length scales, we show that the transition of the thread shape from elliptical to circular is governed by so called effective capillary number. This proves the significance of EIT on the thread dynamics in the microfluidic flow system.

Also, we show that the numerically simulated velocity flow fields along the centreline and across the cross-section of the channel, show an excellent agreement with the experimental measurements obtained through Doppler OCT. The transversal profile of the streamwise velocity divulges the parabolic nature of core flow before the focusing. Once the dispersion thread detaches from the walls and attains a near-elliptical shape, its velocity profile is flat. The water sheath flow close to the walls has a parabolic profile.

Summing up, our three-dimensional experiments and numerical simulations in a microfluidic system have shed new light on the significance and impact of EIT present in miscible complex fluids such as colloidal and polymer fluids, unveiling a distinct framework for assembly, transport and control of material properties through the optimal usage of the surface area of thread shape and wetted region along with other flow parameters such as flow-rate and viscosity ratios. Two-dimensional observations of microfluidic experiments cannot capture this fascinating flow physics. Furthermore, the developed methodology combining experiments and numerical simulations can be used to measure effective interfacial tension acting between two miscible fluids.

Acknowledgements

Financial support by the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (FORMAS), & Knut and Alice Wallenberg foundation through WWSC is gratefully acknowledged. Computer time was provided by the Swedish National Infrastructure for Computing (SNIC). The authors thank Dr. Mathias Kvick, Calvin Brett and Dr. Jingmei Li for experimental assistance.

Appendix A

Optical coherence tomography measurements and analysis

This appendix gives some details about the data acquisition and processing using OCT. The simplest acquisition is a scan of the sample along the depth, i.e. the direction of the light in the sample (vertical in figure 3). Such a scan is called an A-scan and the frequencies given below corresponds to the acquisition frequencies of this type of scans. To obtain a two-dimensional (B-scan) or a three-dimensional image of the sample, it is necessary to perform multiple A-scans by moving laterally the beam in either one or two orthogonal directions. Three-dimensional scans have been used to measure the shape of the thread while the velocity acquisition is only limited to two-dimensional scans.

The OCT apparatus is placed vertically and the channel is inclined with an angle of $7^{\circ}$ with respect to the horizontal in order to limit the apparition of multiple reflections of light rays created by the sample and reverberated back to OCT. As the colloidal dispersion and water exhibit differences in optical scattering, the different fluid phases are separated after image binarization using a contrast threshold. The frequency of acquisition is $f=5.5$ kHz for 3D acquisitions to determine the shape of the thread, leading to a total scanning time of a few minutes. The scanning time is longer with respect to the typical convective time scale $T_{c}$ but it is acceptable since the system is in a stationary state. The Doppler acquisitions have been done at $f=5.5$ kHz or at $28$ kHz, depending on the maximum velocity in the measured region of the thread. Typical scanning time is less than a second and the velocity measurements have been averaged over $100$ repetitions.

In order to measure the velocity fields using Doppler velocimetry, it is necessary to seed the fluids with scattering particles. For the water sheath flows, this is done by adding commercial milk to reach $10\%$ concentration, only when measuring velocity profiles. The properties of the sheath flows are not changed by this small modification. For the colloidal dispersion, light scattering already exists without any addition of particles but is relatively low. However, we have chosen not to add particles in the dispersion in order to preserve its physiochemical properties. This low scattering in the dispersion leads to incorrect estimations of the absolute velocity magnitude but to correct relative velocity variations, i.e. the measured velocity profiles are fully valid but need to be rescaled. This rescaling is done by integrating the measured velocity fields $U_{meas}(y,z)$ in the thread and normalizing them by the flow rate $Q_{1}$ prescribed by the syringe-pump:

[TABLE]

where $U_{corr}(y,z)$ are the corrected velocity fields and $\Sigma$ represents the thread cross-section.