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

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 occurs at relatively short time scales, 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 characterized in three dimensions 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 interfacial tension numerically, we show that it controls the thread development, which can be described by an effective capillary number. Furthermore, we demonstrate that the applied methodology provide the means to measure the ultra-low but dynamically highly significant effective interfacial tension.

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 (Stone, Stroock & Ajdari 2004; Squires & Quake 2005; Whitesides 2006; 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 (Anna, Bontoux & Stone 2003; Garstecki et al. 2004; Cubaud & Mason 2008), microparticle and nanoparticle 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. 2017, 2018), cell patterning (Takayama et al. 1999), DNA stretching (Wong, Lee & Ho 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.

Figure 1. Schematic of the flow-focusing channel.

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

(1.1) $$\begin{eqnarray}Ca=\frac{\unicode[STIX]{x1D702}U}{\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where $U$ is the typical velocity of the flow, $\unicode[STIX]{x1D702}$ the dynamic viscosity and $\unicode[STIX]{x1D6FE}$ the interfacial tension between the two fluids. The capillary number acts as the most important dimensionless number for the 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 (Utada et al. 2007; Humphry et al. 2009). These flow patterns demonstrate the large variety of flow features typically observed in a microfluidic device. Moreover, much of the effort in the last two decades has been devoted to understanding 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 an interface whose thickness grows in time and 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

(1.2) $$\begin{eqnarray}Pe=\frac{T_{d}}{T_{c}}=\frac{hU}{D},\end{eqnarray}$$

which is the relevant dimensionless number to describe convective diffusive flows (Atencia & Beebe 2005; Cubaud & Mason 2006). 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 are short-lived and limited to short time scales before complete mixing, little is still known about the evolution 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. Nevertheless, the difference between the normal and tangential components of the stress tensor near an interface between fluids gives rise to an interfacial tension (Rowlinson & Widom 1982). Korteweg (1901) concluded that interfacial tension is present when two miscible fluids are brought in contact, arising from the composition gradients of a fluid property over the interface. Joseph & Renardy (1993) attributed the first observation to M. Bosscha, reported in a communication at the Academy of Sciences of Amsterdam in 1871. Thus, miscible fluids flowing under laminar flow conditions and high Péclet number may have a de facto interface (Joseph 1990; Joseph & Renardy 1993; Atencia & Beebe 2005). Such an interface exists only as long as the composition gradient of a fluid property is at hand before the interface smears out to an equilibrated homogeneous phase via diffusion or other mixing mechanisms. This de facto interface may also resemble a distinct interface present between two immiscible fluids (Garik et al. 1991; Anderson, McFadden & Wheeler 1998). As the interfacial stresses due to composition gradients over the interface region between two miscible liquids mimics the effect of interfacial tension, an effective interfacial tension, see e.g. Truzzolillo et al. (2014, 2016), can be introduced as

(1.3) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{e}=K\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}^{2}}{\unicode[STIX]{x1D6FF}},\end{eqnarray}$$

where $K$ is the Korteweg constant of the system, $\unicode[STIX]{x1D6FF}$ is the thickness of the interface and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$ is the change in the concentration or volume fraction $\unicode[STIX]{x1D6F7}$ . Truzzolillo et al. (2014, 2016) 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 effective interfacial tension is also a function of particle structure (basically size and shape) together with particle–particle and particle–solvent interactions (Truzzolillo et al. 2016). Typical values range from $10^{-4}$ to $10^{1}~\text{mN}~\text{m}^{-1}$ and correspond to ultra-low to low values of classical interfacial tensions.

Effective interfacial tension is expected to have relevance in many fields starting from material processing to multiphase complex fluid dynamic problems involving droplet and bubble 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 (when following the fluid) 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 effective interfacial tension 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 effective interfacial tension and its effects on the flow-focusing of a colloidal dispersion by its own solvent. Detailed three-dimensional (3-D) experimental and numerical investigations have been performed in a microfluidic flow-focusing channel of square cross-section using optical coherence tomography and the volume of fluid method, respectively. 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). For such a system, effective interfacial tension is expected to exist between the colloidal dispersion and its solvent (Truzzolillo & Cipelletti 2017). Here, we quantitatively investigate the characteristics of thread development such as 3-D topology and velocity fields as a function of interfacial tension and the ratio of the flow rates. 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 3-D shape of the thread. Indeed, there are only a few 3-D 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 shape. Thus, at large Péclet number, the role of effective interfacial tension potentially acting during the existence of an interface between two miscible fluids is not understood so far due to lack of 3-D experimental and numerical studies.

This paper is organized as follows. In § 2, experimental and numerical set-ups are presented together with the numerical validation. In § 3, we compare the results of 3-D computations with the experimental measurements in the threading regime. In particular, the focus is on the 3-D thread shape, and morphology of the region wetted by the colloidal dispersion on the top and bottom walls. The effects of flow rate ratio variations on the thread shape is investigated experimentally and numerically as well as the influence of interfacial tension through numerical simulations. Lastly, we discuss the physical significance of interfacial tension, allowing us to estimate a value of the effective interfacial tension for the present system. In § 4, we look at the velocity fields both along the centre line as well as in the cross-sectional planes of the channel. Finally, in § 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~\text{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 the 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~\unicode[STIX]{x03BC}\text{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 (standard case), flow rates are $Q_{1}=6.5~\text{mm}^{3}~\text{s}^{-1}$ for the colloidal dispersion and $Q_{2}=7.5~\text{mm}^{3}~\text{s}^{-1}$ for the sheath flows, respectively. The ratio of these flow rates is defined as $\unicode[STIX]{x1D719}=Q_{1}/Q_{2}$ and $\unicode[STIX]{x1D719}=0.8667$ in the standard case, except when explicitly varied in § 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 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 revolutions per minute. A transparent colloidal dispersion of concentration $3~\text{g}~\text{dm}^{-3}$ was obtained with typical fibril length  $l$ of 700 nm and fibril diameter  $d$ of 2–3 nm.

Figure 2. Shear viscosity measurements of the colloidal dispersion represented by red dots, and the non-Newtonian Carreau model fit (2.1) indicated by the black solid line.

The rheological characterisation 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 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, Belgacem & Pignon 2016). The rheological data can be fitted well by a Carreau model

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{eff}=\unicode[STIX]{x1D702}_{inf}+(\unicode[STIX]{x1D702}_{0}-\unicode[STIX]{x1D702}_{inf})[1+(\unicode[STIX]{x1D70F}\dot{\unicode[STIX]{x1D6FE}})^{2}]^{n-1/2},\end{eqnarray}$$

where $\unicode[STIX]{x1D702}_{eff}$ is the shear viscosity, $\dot{\unicode[STIX]{x1D6FE}}$ the shear rate, $\unicode[STIX]{x1D702}_{0}$ the zero shear viscosity, $\unicode[STIX]{x1D702}_{inf}$ the infinite shear viscosity, $\unicode[STIX]{x1D70F}$ the relaxation time and $\mathit{n}$ the power index. The best fit is shown by the black solid line in figure 2 and gives $\unicode[STIX]{x1D702}_{inf}=5~\text{mPa}~\text{s}$ , $\unicode[STIX]{x1D702}_{0}=1756~\text{mPa}~\text{s}$ , $\unicode[STIX]{x1D70F}=16.16~\text{s}$ and $n=0.56$ . The viscosity ratio between the core and sheath flow defined as $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D702}_{1}/\unicode[STIX]{x1D702}_{2}=1756$ is arrived at by considering the viscosity of colloidal dispersion at low shear rate $\unicode[STIX]{x1D702}_{1}=\unicode[STIX]{x1D702}_{0}=1756~\text{mPa}~\text{s}$ and the viscosity of sheath fluid $\unicode[STIX]{x1D702}_{2}=1~\text{mPa}~\text{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)

(2.2) $$\begin{eqnarray}D=\frac{k_{B}T\ln (l/d)}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}l},\end{eqnarray}$$

with the temperature $T=300~\text{K}$ , the Boltzmann constant $k_{B}=1.38\times 10^{-23}~\text{J}~\text{K}^{-1}$ and the solvent viscosity $\unicode[STIX]{x1D702}=1~\text{mPa}~\text{s}$ . This gives $D\approx 5\times 10^{-12}~\text{m}^{2}~\text{s}^{-1}$ , and thus a typical diffusion time scale in the system $T_{d}=h^{2}/D\approx 2\times 10^{5}~\text{s}$ . For a flow velocity $U\approx 10~\text{mm}~\text{s}^{-1}$ , the convective time scale can be estimated to be $T_{c}=h/U\approx 10^{-1}~\text{s}$ . Using the above values and (1.2), we get $Pe\approx 2\times 10^{6}$ . The high Péclet number implies that the system is predominantly convective driven and the diffusion is very weak. Therefore, since the time scale for diffusion of the nanofibrils from the core flow into the sheath flows being substantially larger than the convective time scale of the two fluids in the channel, it is possible for a de facto interface to exist between the two fluids (Joseph 1990; Joseph & Renardy 1993; Atencia & Beebe 2005). Explicitly, in the context of the present experiments, the de facto interface exists due to the presence of a region of composition gradient at the boundary between the colloidal dispersion and its solvent. Here, interdiffusion between colloidal dispersion and its solvent is almost negligible compared to convection thereby keeping the interface sharply defined and preventing it from smearing (Atencia & Beebe 2005). As mentioned in the introduction, the properties of such an interface may resemble a distinct interface between two immiscible fluids and thus, the two fluids, the colloidal dispersion and water, exhibit a near-immiscible behaviour. In the following two sections, the experimental method and numerical set-up will be presented. The experimental method section involves more about the description of optical coherence tomography employed for topology and velocity measurements.

2.1.4 Optical coherence tomography

Figure 3. Schematic of the working principle of optical coherence tomography, adapted from Tomlins & Wang (2005).

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 resolution. 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, Jung & Chen 2008).

Since OCT might be new to the reader, a brief description is given. 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 and 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 millimetres per second.

The OCT apparatus used in this study is a commercial spectral domain 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~\unicode[STIX]{x03BC}\text{m}$ . Here, OCT together with Doppler acquisition is employed for 3-D 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 $\unicode[STIX]{x1D6FC}$ via the volume fractions of the two fluids in the computational cells; $\unicode[STIX]{x1D6FC}=0$ corresponds to fluid 1, $\unicode[STIX]{x1D6FC}=1$ to fluid 2, with $0<\unicode[STIX]{x1D6FC}<1$ in the interface region. The advection of the interface is achieved through the redistribution of $\unicode[STIX]{x1D6FC}$ between neighbouring cells. Various VoF methods have been developed since its inception, and are broadly categorized into algebraic and geometric methods. Algebraic VoF methods are developed for general mesh types, easier to implement and faster (Deshpande, Anumolu & Trujillo 2012). However, they are less accurate and often require higher-order schemes to retain the sharpness and boundedness of fluid fraction  $\unicode[STIX]{x1D6FC}$ . 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, Bredmose & Jasak 2016) that combines high accuracy with advection of the interface.

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

(2.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{b}\boldsymbol{U}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{b}\boldsymbol{U}\boldsymbol{U})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}+\unicode[STIX]{x1D70C}_{b}\boldsymbol{g}+\boldsymbol{F}_{\boldsymbol{s}},\end{eqnarray}$$

the continuity equation

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0,\end{eqnarray}$$

and the equation for the advection of fluid fraction $\unicode[STIX]{x1D6FC}$

(2.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}\boldsymbol{U})=0.\end{eqnarray}$$

Here $\boldsymbol{U}$ is the velocity vector field, $\unicode[STIX]{x1D70C}_{b}$ the density, $p$ the pressure field, $\unicode[STIX]{x1D64F}$ the deviatoric stress tensor $(\unicode[STIX]{x1D64F}=2\unicode[STIX]{x1D707}_{b}\boldsymbol{S}-2\unicode[STIX]{x1D707}_{b}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U})\unicode[STIX]{x1D644}/3)$ , where $\boldsymbol{S}=0.5(\unicode[STIX]{x1D735}\boldsymbol{U}+\unicode[STIX]{x1D735}\boldsymbol{U}^{\text{T}})$ , $\unicode[STIX]{x1D644}$ is the identity matrix and $\boldsymbol{g}$ is the gravitational acceleration. The density $\unicode[STIX]{x1D70C}_{b}$ and viscosity $\unicode[STIX]{x1D707}_{b}$ are computed based on the weighted average distribution of the fluid fraction $\unicode[STIX]{x1D6FC}$

(2.6) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70C}_{b}=\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D70C}_{2}(1-\unicode[STIX]{x1D6FC}), & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D707}_{b}=\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D707}_{2}(1-\unicode[STIX]{x1D6FC}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{1}$ , $\unicode[STIX]{x1D70C}_{2}$ , $\unicode[STIX]{x1D707}_{1}$ , $\unicode[STIX]{x1D707}_{2}$ are the densities and the viscosities of the two phases. Here $\boldsymbol{F}_{\boldsymbol{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 $\boldsymbol{F}_{\boldsymbol{s}}$ is modelled as a volumetric force using the continuum surface force (CSF) method (Brackbill, Kothe & Zemach 1992), and is defined as

(2.8) $$\begin{eqnarray}\boldsymbol{F}_{\boldsymbol{s}}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the interfacial tension and $\unicode[STIX]{x1D705}$ is the curvature of the interface

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D705}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}|}\right).\end{eqnarray}$$

The numerical simulations were performed with the finite-volume-based open-source code OpenFOAM (Weller et al. 1998). The pressure-implicit with splitting of operators scheme (Rusche 2003) is applied for solving the momentum balance equation (2.3) in conjunction with the continuity equation (2.4). 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 scheme with van Leer limiter was used for the spatial discretization. Boundedness of the fluid phase fraction $\unicode[STIX]{x1D6FC}$ , was controlled by the IsoAdvector (Roenby et al. 2016) scheme implemented in the interIsoFoam, a two-phase incompressible immiscible flow solver in 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~\text{mm}^{2}$ , 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 $\unicode[STIX]{x1D702}_{eff}$ is sufficient to describe the relevant material properties of the colloidal dispersion. To initialize the simulation, we set the fluid phase fraction, $\unicode[STIX]{x1D6FC}=0$ in the core flow inlet channel, $\unicode[STIX]{x1D6FC}=1$ in the sheath flow inlet channel and for the outlet channel. A uniform velocity flow profile was set as boundary condition at the entrance of the core and sheath flows channel inlets. The velocities were calculated based on the flow rates $Q_{1}=6.5~\text{mm}^{3}~\text{s}^{-1}$ for the core flow and $Q_{2}=7.5~\text{mm}^{3}~\text{s}^{-1}$ for the sheath flows, respectively. At the walls, we apply the no-slip boundary condition and use an equilibrium contact angle $\unicode[STIX]{x1D703}=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

(2.10) $$\begin{eqnarray}\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}/\unicode[STIX]{x2202}n=0,\end{eqnarray}$$

where $n$ is the unit vector normal to the outflow boundary.

2.2.2 Validation and convergence

The validation of the numerical scheme is performed by simulating a reference case, which also demonstrates the large range of flow patterns 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 with sidelength $h=100~\unicode[STIX]{x03BC}\text{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 $\unicode[STIX]{x1D702}_{1}=1214~\text{mPa}~\text{s}$ , and $\unicode[STIX]{x1D702}_{2}=4.59~\text{mPa}~\text{s}$ , respectively. The interfacial tension between the two fluids is $\unicode[STIX]{x1D6FE}_{12}=27.0~\text{mN}~\text{m}^{-1}$ . The capillary numbers of core and sheath flows are defined as $Ca_{1}=\unicode[STIX]{x1D702}_{1}Q_{1}/(\unicode[STIX]{x1D6FE}_{12}h^{2})$ and $Ca_{2}=\unicode[STIX]{x1D702}_{2}Q_{2}/(\unicode[STIX]{x1D6FE}_{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 3-D 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 $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{z}/h$ of the thread in the threading regime (i) as a function of the flow rate ratio  $\unicode[STIX]{x1D719}$ as depicted in figure 4(c). Indeed, the thread cross-section is circular and scales as $\unicode[STIX]{x1D700}_{y}/h=\unicode[STIX]{x1D700}_{z}/h=(\unicode[STIX]{x1D719}/2)^{1/2}$ , as observed experimentally and expected theoretically by Cubaud & Mason (2008, 2009). The pink solid line represents the fit of $\unicode[STIX]{x1D700}_{y}/h=\unicode[STIX]{x1D700}_{z}/h=(\unicode[STIX]{x1D719}/2)^{1/2}$ .

Figure 4. Validation of three-dimensional simulations in a microfluidic flow-focusing channel comparing with the observations of Cubaud & Mason (2008) experiments. (a) Flow-patterns viz., (i) threading, (ii) jetting, (iii) dripping and (iv) tubing; (b) state space of flow regimes mapped based on capillary numbers of the core $Ca_{1}$ and sheath $Ca_{2}$ fluid phases, adapted from Cubaud & Mason (2008). The markers show simulated cases: ○ threading, ▹ jetting, ▫ dripping and ♢ tubing. (c) Dimensionless width $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{z}/h$ of the thread versus flow rate ratio  $\unicode[STIX]{x1D719}$ . Inset: top and cross-sectional view of flow-focusing channel depicting the flow-pattern of threading regime (i), channel width  $\mathit{h}$ , width of thread $\unicode[STIX]{x1D700}_{y}/h$ and height of thread $\unicode[STIX]{x1D700}_{z}/h$ .

Figure 5. Effect of grid size $\unicode[STIX]{x1D6E5}$ on the (a) width $\unicode[STIX]{x1D700}_{y}/h$ , height $\unicode[STIX]{x1D700}_{z}/h$ of the thread, and (b) velocity in a square channel cross-section of width $h$ . The variation with respect to the finest resolution $h/\unicode[STIX]{x1D6E5}=80$ in $\unicode[STIX]{x1D700}_{y}/h$ , $\unicode[STIX]{x1D700}_{z}/h$ and velocity $\mathit{U}$ due to grid size $\unicode[STIX]{x1D6E5}$ is indicated by $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{y}}$ , $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{z}}$ and $\unicode[STIX]{x1D6FF}_{U}$ . Red markers denote that a grid size $h/\unicode[STIX]{x1D6E5}=50$ which is chosen for the main results in the present study.

All the validation studies and upcoming results of the present study in §§ 3 and 4 were obtained at a resolution of $h/\unicode[STIX]{x1D6E5}=50$ , where $\unicode[STIX]{x1D6E5}$ 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 $\unicode[STIX]{x1D6E5}$ on width $\unicode[STIX]{x1D700}_{y}/h$ , height $\unicode[STIX]{x1D700}_{z}/h$ of the thread, and on the flow field velocity $\mathit{U}$ . We observe that for $h/\unicode[STIX]{x1D6E5}=50$ , the variation with respect to the finest resolution $h/\unicode[STIX]{x1D6E5}=80$ in the width ( $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{y}}$ ), height ( $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{z}}$ ) of the thread, and on the flow field velocity ( $\unicode[STIX]{x1D6FF}_{U}$ ) is smaller than 2 %, while for $h/\unicode[STIX]{x1D6E5}>60$ it is smaller than 1 %. It will be seen that $h/\unicode[STIX]{x1D6E5}=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 h 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 $\unicode[STIX]{x1D719}$ and interfacial tension $\unicode[STIX]{x1D6FE}$ on thread formation is investigated. The flow rate ratio for the standard case is respectively set to $\unicode[STIX]{x1D719}=0.8667$ in the experiments and simulation. The interfacial tension is set to $\unicode[STIX]{x1D6FE}=0.054~\text{mN}~\text{m}^{-1}$ in the simulations, except when explicitly varied in §§ 3.2 and 3.3.

Figure 6. 3-D view of the experimental and numerical threads. The black dashed lines defines the beginning of the focusing region at $x/h=0$ , and cyan colour at the top plane $z/h=0.5$ denotes the wetted region.

3.1 Shape of the thread in the standard case

The 3-D 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\leqslant x/h\leqslant 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 the experiments are very well captured and reproduced by the numerical simulations.

Figure 7. Comparison of experimental and numerical wetted regions of the core fluid in the plane $z/h=0.5$ . The length of this region in the downstream direction $x/h$ is defined as $L_{w}/h$ .

Figure 8. Experimental (a–c) and numerical (d–f) cross-sections of the thread. The Black colour represents the core fluid while the white colour indicates the sheath fluid. The cross-sectional views are taken at $x/h=1.5$ (a,d), $x/h=3$ (b,e) and $x/h=10$ (c,f).

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 $\unicode[STIX]{x1D700}_{y}/h$ , $\unicode[STIX]{x1D700}_{z}/h$ and aspect ratio $\unicode[STIX]{x1D700}_{z}/\unicode[STIX]{x1D700}_{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 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 $\unicode[STIX]{x1D700}_{z}/h$ and the width $\unicode[STIX]{x1D700}_{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 ( $\unicode[STIX]{x1D700}_{y}/h\approx 0.5$ and $\unicode[STIX]{x1D700}_{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.

Figure 9. Thread width $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{z}/h$ in (a), and aspect ratio $\unicode[STIX]{x1D700}_{z}/\unicode[STIX]{x1D700}_{y}$ (b) as a function of downstream locations $x/h$ . The experimental data are shown in red while the numerical data are represented by solid black lines. The vertical dashed-dotted lines at $x/h=2$ indicate the end of the wetted region.

3.2 Shape variations with flow rate ratio

Figure 10. Width $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{z}/h$ of the thread as a function of flow rate ratio $\unicode[STIX]{x1D719}$ . Red square symbols represent experiments while diamond symbols denote numerical simulations. The fit of the first part of the curves by the relation $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D719}^{k}$ are shown by the pink lines and the fitted values of $k$ are written in each panel. (c) Illustration using numerical data of the different regimes separated in (a,b) by vertical dashed-dotted lines: threading regime (i) and tubing regime (ii, iii).

In order to investigate the role played by the sheath flows on the thread formation, a parametric variation of the flow rate ratio $\unicode[STIX]{x1D719}$ over two orders of magnitude is performed with the same rheology of colloidal dispersion as given by the Carreau model in figure 2. 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 the means to synthesize materials of complex shapes such as rods, discs or ellipsoids (Xu et al. 2005). The dimensionless width $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{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 $\unicode[STIX]{x1D719}$ . 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 $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D719}^{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 $\unicode[STIX]{x1D719}$ , represented by pink solid lines in figure 10, gives $\unicode[STIX]{x1D700}_{y}/h=0.58\unicode[STIX]{x1D719}^{0.5}$ for the width and $\unicode[STIX]{x1D700}_{z}/h=0.89\unicode[STIX]{x1D719}^{0.28}$ for the height. A power law in $\unicode[STIX]{x1D719}^{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 $\unicode[STIX]{x1D719}$ . However, the height of the thread exhibits a different power law, in $\unicode[STIX]{x1D719}^{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 $\unicode[STIX]{x1D719}$ but could also be dependent on other parameters such as viscosity ratio $\unicode[STIX]{x1D712}$ and interfacial tension $\unicode[STIX]{x1D6FE}$ . Later, in the upcoming § 3.3, the effect of interfacial tension on the height of thread $\unicode[STIX]{x1D700}_{z}/h$ is discussed using numerical simulations.

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 are observed within the range of the flow rate ratio sampled in this study. They are illustrated in figure 10(c). For $\unicode[STIX]{x1D719}<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 $\unicode[STIX]{x1D719}$ . In the range $1.5<\unicode[STIX]{x1D719}<3$ , the height of the thread saturates reaching the size of the channel $h$ in $z$ direction while the width keeps increasing. This is the first regime of tubing (ii). Above $\unicode[STIX]{x1D719}>3$ , both height and width have saturated reaching 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

Figure 11. Variations in per cent of the error $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{z}}$ on the height of the thread (a) and of the error $\unicode[STIX]{x1D6FF}_{L_{w}}$ on the wetted length (b) as a function of interfacial tension $\unicode[STIX]{x1D6FE}$ obtained from numerical simulations. The red points in (a) and (b) indicate the interfacial tension that gives the best match with the experimental data.

Figure 12. Variations in the height of the thread $\unicode[STIX]{x1D700}_{z}/h$ (a) with downstream locations $x/h$ for various interfacial tensions $\unicode[STIX]{x1D6FE}$ (b) as a function of scaled downstream location for different interfacial tensions $\unicode[STIX]{x1D6FE}$ . $Ca_{1}$ corresponds to the capillary number of core fluid thread.

We utilize the numerical simulations to study the influence of interfacial tension on the thread shape within the range $0.024~\text{mN}~\text{m}^{-1}<\unicode[STIX]{x1D6FE}<3.000~\text{mN}~\text{m}^{-1}$ . The height of the thread $\unicode[STIX]{x1D700}_{z}/h$ at $x/h=20$ together with the wetted length $L_{w}/h$ are now used to illustrate how interfacial tension affects the shape in the calculations by comparing with the experimental values. This has been done using the following definitions for the differences

(3.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{z}}(\unicode[STIX]{x1D6FE})=\left|\frac{\unicode[STIX]{x1D700}_{z}^{num}(\unicode[STIX]{x1D6FE})-\unicode[STIX]{x1D700}_{z}^{exp}}{\unicode[STIX]{x1D700}_{z}^{exp}}\right|,\quad \text{and}\quad \unicode[STIX]{x1D6FF}_{L_{w}}(\unicode[STIX]{x1D6FE})=\left|\frac{L_{w}^{num}(\unicode[STIX]{x1D6FE})-L_{w}^{exp}}{L_{w}^{exp}}\right|.\end{eqnarray}$$

The differences $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{z}}(\unicode[STIX]{x1D6FE})$ and $\unicode[STIX]{x1D6FF}_{L_{w}}(\unicode[STIX]{x1D6FE})$ are plotted in figure 11. Both curves show significant variations i.e. $\unicode[STIX]{x1D700}_{z}/h$ and $L_{w}/h$ are controlled by $\unicode[STIX]{x1D6FE}$ for a fixed viscosity ratio $\unicode[STIX]{x1D712}$ and flow rate ratio $\unicode[STIX]{x1D719}$ . The minima of $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}_{z}}(\unicode[STIX]{x1D6FE})$ and $\unicode[STIX]{x1D6FF}_{L_{w}}(\unicode[STIX]{x1D6FE})$ , are indicated by red filled symbols in figure 11(a,b) and occur at $\unicode[STIX]{x1D6FE}=0.054~\text{mN}~\text{m}^{-1}$ . This value is used in the current study for the comparison of numerical results with the experimental measurements presented in §§ 3.1, 3.2 and 4.

Figure 12(a) shows the height of the thread $\unicode[STIX]{x1D700}_{z}/h$ as a function of downstream position $x/h$ for various $\unicode[STIX]{x1D6FE}$ . First, one can again note that the wetted length $L_{w}/h$ is influenced by $\unicode[STIX]{x1D6FE}$ , as the detachment of the various threads from the top and bottom walls of the channel takes place in the range $1\leqslant x/h\leqslant 3$ (see also figure 11). Once the thread is detached, the larger the $\unicode[STIX]{x1D6FE}$ , the faster the thread converges towards a value of 0.63, corresponding to the width of the thread $\unicode[STIX]{x1D700}_{y}/h$ for the simulation with the highest $\unicode[STIX]{x1D6FE}$ . Far downstream, for $x/h>20$ , the height of the thread with low $\unicode[STIX]{x1D6FE}$ is still decreasing. This strongly suggests that the near-elliptical shape of the thread cross-section is only an intermediate state, before flattening towards a circular shape far downstream. 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 with three inlets in the $(x,y)$ -plane, there is no rotational symmetry of $\unicode[STIX]{x03C0}/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.

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 a near-elliptical shape, with $\unicode[STIX]{x1D700}_{z}/h\approx 1$ and $\unicode[STIX]{x1D700}_{y}/h<1$ . Figure 12(a) suggests that the thread reorganizes itself by converging towards a circular cross-section, with $\unicode[STIX]{x1D700}_{z}/h=\unicode[STIX]{x1D700}_{y}/h<1$ . This spatial evolution depends on the interfacial tension and can therefore be assumed to be mainly driven by $\unicode[STIX]{x1D6FE}$ at the interface, as the interfacial tension tends to minimize the contact surface. Indeed, as the interface between the core and sheath fluids is curved, there is a Laplace pressure difference between the two fluids, proportional to $\unicode[STIX]{x1D6FE}$ and to the local curvature of the interface. As $\unicode[STIX]{x1D700}_{z}/h>\unicode[STIX]{x1D700}_{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 $\unicode[STIX]{x1D6FF}P$ within the thread, proportional to the interfacial tension $\unicode[STIX]{x1D6FE}$ and dependent on the geometry of the thread, driving the thread cross-section from a near-elliptical to a circular shape.

As the viscosity ratio $\unicode[STIX]{x1D712}$ in the system is very high, one can neglect the influence of the sheath flow viscosity on this evolution. Consequently, it is mainly dependent on the viscosity $\unicode[STIX]{x1D702}_{1}$ of the core fluid. The typical time scale $\unicode[STIX]{x1D70F}$ for this dynamics can be estimated via the Stokes equation as

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=\frac{\unicode[STIX]{x1D702}_{1}}{\unicode[STIX]{x1D6FF}P}\propto \frac{\unicode[STIX]{x1D702}_{1}}{\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where the thread geometry dependence of $\unicode[STIX]{x1D6FF}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}/(\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{z}\unicode[STIX]{x1D700}_{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

(3.3) $$\begin{eqnarray}x_{c}-L_{w}=U\unicode[STIX]{x1D70F}\propto \frac{\unicode[STIX]{x1D702}_{1}Q_{1}}{\unicode[STIX]{x1D6FE}h^{2}}.\end{eqnarray}$$

One can therefore normalize the modified downstream location $(x-L_{w})/h$ by the capillary number $Ca_{1}=\unicode[STIX]{x1D702}_{1}Q_{1}/(h^{2}\unicode[STIX]{x1D6FE})$ . This is shown in figure 12(b). All the numerical data collapse well on a master curve, showing that the shape of the thread is governed by the capillary number.

Nevertheless, given the presence of high viscosity contrast between colloidal dispersion and its solvent, there could be a possibility for the shape of the thread to be driven through other hydrodynamic phenomena, independent of interfacial tension. Indeed, one can argue that the shape of the thread naturally evolves in the parabolic velocity profile of the sheath flow in order to minimize the viscous dissipation. Thus, as the isovelocity contours of the sheath flow are almost circular near the centre line of a square duct, one should expect the thread to progressively adopt a circular shape further downstream to reduce high-viscosity shear. However, the master curve (see figure 12 b) obtained through the simple model shows only the interfacial tension ((3.2) and (3.3)) describing very well the development of the thread. Therefore, there is no evidence of major effects of other hydrodynamic phenomena on the thread shape, even though viscous dissipation and lubrication at the interface might play a role in the system.

3.4 Physical significance of interfacial tension

The numerical results in § 3.3 affirm the strong effects of interfacial tension $\unicode[STIX]{x1D6FE}$ on the shape of the thread, particularly at ultra-low values ( $0.024~\text{mN}~\text{m}^{-1}<\unicode[STIX]{x1D6FE}<0.5~\text{mN}~\text{m}^{-1}$ ) where the thread cross-section remains elliptical far downstream. Indeed, the interfacial tension $\unicode[STIX]{x1D6FE}$ controls the spatial evolution of thread height along the channel length and this can be well understood by the simple model ((3.2) and (3.3)) leading to a master curve depicted in figure 12(b). Moreover, the detailed discussion of numerical and experimental results in §§ 3.1 and 3.2 shows that when we choose the appropriate interfacial tension in the numerics ( $\unicode[STIX]{x1D6FE}=0.054~\text{mN}~\text{m}^{-1}$ , which is three orders of magnitude lower than classical ones), all the features such as the morphology of the wetted region and topology of the experimental thread for different flow rates are very well captured and reproduced. Furthermore, it will be shown in § 4 that the agreement between the experimental and numerically obtained velocity fields is also excellent for the same value of interfacial tension, i.e. $\unicode[STIX]{x1D6FE}=0.054~\text{mN}~\text{m}^{-1}$ .

As the Péclet number is high in the system, the presence of a de facto interface between the colloidal dispersion and its solvent is expected in the experiments. Such a de facto interface can be associated with a near-immiscible behaviour of the two fluids and can lead to a weak but effective interfacial tension $\unicode[STIX]{x1D6E4}_{e}$ between the core and sheath fluid. Thus, the results presented in § 3, using an immiscible fluid solver with ultra-low values of interfacial tension $\unicode[STIX]{x1D6FE}$ , allow us to claim that an effective interfacial tension exists in the system and that in turn controls the thread morphology. From a physical interpretation perspective, the interfacial tension  $\unicode[STIX]{x1D6FE}$ in the numerical simulations models the ‘effective’ interfacial tension in the experiments. Hence, we refer to $\unicode[STIX]{x1D6E4}_{e}=\unicode[STIX]{x1D6FE}=0.054~\text{mN}~\text{m}^{-1}$ . Note that the value of the effective interfacial tension found in the present system is of the same order as experimental observations of effective interfacial tensions ( $10^{-2}~\text{mN}~\text{m}^{-1}$ ) made by Truzzolillo et al. (2016) between complex fluids like colloidal or polymer dispersions and its own solvent.

Since, the interface between miscible fluids is in a non-equilibrium state, it is often difficult to measure the effective interfacial tension experimentally (Truzzolillo et al. 2014, 2016; Truzzolillo & Cipelletti 2017). As an alternative, we show that numerical simulations can be used to extract this variable, by reproducing the thread shape obtained in the experiments. Indeed, as the experimental shape is controlled by an ‘effective’ capillary number $Ca_{e}=\unicode[STIX]{x1D702}_{1}Q_{1}/(h^{2}\unicode[STIX]{x1D6E4}_{e})$ , it is only necessary to measure the evolution of the height of the thread $\unicode[STIX]{x1D700}_{z}/h$ formed by these two fluids and to compare it with the master curve depicted in figure 12(b).

4 Velocity fields

Figure 13. Experimental (red points) and numerical (solid black line) centre line velocity $U_{c}$ as a function of the downstream position $x/h$ . Error bars for experimental data are shown and are around 2 %.

In addition to the comparison between numerical simulations and experiments with respect to the 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 this Section, the flow rate ratio is fixed to $\unicode[STIX]{x1D719}=0.8667$ in the experiments and simulation while the interfacial tension is set to $\unicode[STIX]{x1D6FE}=0.054~\text{mN}~\text{m}^{-1}$ in the simulation, as in the standard case in § 3.1.

In figure 13, the centre line 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. Experimental and numerical velocity maps (a,d,g and b,e,h) and profiles (c,f,i) at three different downstream locations: $x/h=-3$ (a–c), $x/h=1$ (d–f) or $x/h=10$ (g–i). The velocity profiles in (c,f,i) show the $x$ -component of the velocity at $z/h=0$ . The experimental data are represented by the red line while the numerical data are shown with the black line.

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=-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 3-D thread formation under the influence of effective interfacial tension 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 time scale 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 interface where effective interfacial tension is acting. Three-dimensional versions of the finite-volume-based open-source code OpenFOAM (OpenCFD 2017) with VOF method have been used for the numerical simulations and experimental measurements were carried out with optical coherence tomography. The validation of the numerical scheme was performed by choosing and reproducing a reference case for the immiscible fluids reported in Cubaud & Mason (2008) that demonstrates the diversity 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 3-D 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 $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{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 $\unicode[STIX]{x1D700}_{y}/h$ and height $\unicode[STIX]{x1D700}_{z}/h$ of the thread as a function of flow rate ratio $\unicode[STIX]{x1D719}$ 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 $\unicode[STIX]{x1D719}$ . The width of the thread scales purely with the flow rate ratio as $\unicode[STIX]{x1D719}^{0.5}$ in good agreement with previous studies (Cubaud & Mason 2008, 2009). However, the height of thread scales as $\unicode[STIX]{x1D719}^{0.28}$ confirming the nearly elliptical cross-sectional shape, thus indicating the dependence of height on other parameters like interfacial tension along with the flow rate ratio $\unicode[STIX]{x1D719}$ .

The three-dimensional numerical simulations gave access to probe for a detailed investigation of the influence of interfacial tension by comparing with the elliptical thread shape given by experiments. Thus, it provided a framework to deduce the values of interfacial tension by minimizing the differences between simulations and experiments. The experimental effective interfacial tension $\unicode[STIX]{x1D6E4}_{e}$ was determined by varying the interfacial tension $\unicode[STIX]{x1D6FE}$ in the simulations and comparing the numerically generated data with the experimental wetted length $L_{w}/h$ and height $\unicode[STIX]{x1D700}_{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 mainly governed by so called effective capillary number. This proved the significance of interfacial tension on the thread morphology in the microfluidic flow system.

Finally, we show that the numerically simulated velocity flow fields along the centre line 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 effective interfacial tension 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.