Formation of colloidal threads in geometrically varying flow-focusing channels

When two miscible fluids are brought into contact with each other, the concentration gradients induce stresses. These are referred to as Korteweg stresses and are analogous to interfacial tension between two immiscible fluids, thereby acting as an effective interfacial tension (EIT) in inhomogeneous miscible systems. EIT governs the formation of a viscous thread in flow-focusing of two miscible fluids. To further investigate its significance, we have studied thread formation of a colloidal dispersion focused by its own solvent. Experiments are combined with three-dimensional numerical models to systematically expand previous knowledge utilising different flow-focusing channel setups. In the reference setup, the sheath flows impinge the core flow orthogonally while in four other channel setups, the sheath flows impinge the core flow at an oblique angle that is both positive and negative with respect to the reference sheath direction. As an initial estimate of the EIT, we fit the experimentally determined thread shape in the reference setup to a master curve that depends on EIT through an effective capillary number. By numerically reproducing these experimental results, it is concluded that the estimated EIT is within 25% of the optimal EIT value that can be deduced by iteratively fitting the numerical results to the experimental measurements. Regardless of channel setups, further numerical calculations performed using the optimal EIT evaluated from the reference setup show good agreement with the experimental findings in terms of thread shapes, wetted region morphologies, and velocity flow fields. The one-to-one comparison of numerical and experimental findings unveil the crucial role of EIT on the thread detachment from the top and bottom walls of the channel, bringing useful insights to understand the physical phenomenons involved in miscible systems with a high-viscosity contrast.


I. Introduction
A boundary between two fluid phases is identifiable, irrespective of whether the fluids are immiscible or miscible. In the case of immiscible fluids, the boundary zone is clearly distinguished through a distinct interface created by the equilibrium interfacial tension acting between the two fluids. On the other hand, in the context of miscible fluids, no such distinct interface exists as the two fluids are fully mixed at equilibrium. However, when two miscible fluids come into contact, a boundary between them can persist till the time the two fluids eventually form an equilibrated homogeneous mixture due to diffusive mixing. This boundary zone can act as a de-facto interface with some of its properties resembling that of a distinct interface observed between immiscible fluids [1][2][3][4][5]. In 1901, Korteweg [6] first proposed that when two miscible fluids are brought into contact, the composition inhomogeneties/gradients of the fluid property at the zone of contact gives rise to additional stresses (so-called Korteweg stresses). These stresses effectively mimic capillary-like stress effects across the boundary zone that can be seen as a sharp de-facto interface. Accordingly, analogous to the interfacial tension γ in the immiscible fluids, an effective interfacial tension (EIT) for miscible fluids in non-equilibrium state, commonly denoted as Γ e , can be written as: where K is the Korteweg factor accounting for the relevant interaction effects (e.g. particle-solvent interactions in miscible complex fluids) between the two fluids [7,8], δ is the interface thickness and ∆Φ is the variation in composition or volume fraction Φ. As can be seen from Eq. (1), Γ e exists as long as the composition gradients persist at the de-facto interface, and as the interface smears out due to diffusion over time, the EIT goes to zero.
The role of these transient stresses or EIT in non-equilibrium miscible fluids is evident in numerous multi-fluid dynamical processes occuring at short times. Examples are microgravity experiments [9], evolution of miscible droplets [10,11], modeling of hydrodynamic instabilities like viscous fingering or in Hele-Shaw flows [12,13], stabilization of Rayleigh-Taylor instabilities induced by evaporation between a polymer solution and its own solvent [14] and so on. Some of the recent experimental techniques explored to measure the EIT between miscible fluids are through the evolution of drop shape [15,16], examination of hydrodynamic instabilities [7,17], and probing of capillary waves by light scattering [18,19]. In spite of these attempts, measuring the EIT between miscible fluids is intrinsically difficult due to ultra-low values (Γ e ∼ 10 −4 −10 1 mN m −1 ) and absence of a distinct interface [7,8,20].
Lately, in the framework of Korteweg's theory for miscible fluids, employing a microfluidic flow-focusing setup, we have proposed a methodology with a possibility to determine the EIT between two miscible fluids. In short, Γ e can be estimated by measuring the spatial evolution of the thread shape formed by these two fluids, and equating it with a master curve obtained based on a simple scaling model [21]. However, this method was not yet fully exploited.
Flow-focusing essentially comprises of a central core fluid focused by two side (sheath) fluids as shown schematically in Fig. 1(a). For miscible systems, the Péclet number P e = U h/D, which measures the relative importance of convective and diffusive effects, is the relevant non-dimensional quantity to describe the flow with an average velocity U in a microchannel of width h, and D being the diffusion coefficient between the two fluids [5,22]. When a pair of miscible fluids is passed through the flow-focusing channel, the Péclet number dictates the flow-patterns. At high Péclet number under laminar flow conditions, where diffusion is almost negligible, viscous threads are formed, while at low Péclet number, the threads undergo diffusive instabilities leading to a wide diversity of flow-patterns [23][24][25]. It is even observed that the miscible threads at high Péclet number is equivalent to the multi-fluid viscous flow problem, and could be employed to examine and realise the role of diffusion [26,27].
For a high Péclet number system, with a pair of miscible fluids comprising of a colloidal dispersion and its own solvent, we have detected and reported that the characteristics of spatial thread evolution could be accurately captured and modelled as an immiscible fluid problem with a very weak interfacial tension γ for a set of flow rates and specified rheology of the two fluids [21]. Ideally, at high Péclet number, the time scale for interdiffusion between the colloidal dispersion and its own solvent is almost negligible compared to the convective time scale of the fluids in the channel. In such a scenario, the presence of a sharp de-facto interface due to composition gradients is expected in the experiments [5,7,8]. However, it is extremely difficult to access such a sharp de-facto interface between the two fluids experimentally. The miscible viscous thread structures formed at high Péclet number are in an out-of-equilibrium state and occur before the two fluids are fully mixed.
In such occurences, it is possible to reckon the de-facto interface by appyling a weak interfacial tension γ in numerical modelling. In reality, the weak interfacial tension γ accounts for the Korteweg stresses induced by the composition gradients in the experiments. The experimental observations were reproduced numerically with a very weak interfacial tension γ (O∼10 −2 mN m −1 ). Moreover, the magnitude of γ corroborated with the experimental measurements of previous studies conducted for miscible complex fluids involving colloidal dispersions or polymer solutions and its own solvent [8,20]. Thus, it was established that the very weak interfacial tension γ acts as an effective interfacial tension in the experiments with γ ≡ Γ e . Furthermore, in microfluidic systems, the interfacial tension and viscous dissipation of energy dominate over inertial effects [28][29][30] due to the small length scales. Pertaining to miscible systems, and in the context of weakly diffusive threads, the phenomenon of viscous dissipation of energy is pertinent as the viscosity contrast between the fluid pairs is fairly large [23,26]. In addition to these effects, we noted that the effect of EIT due to concentration gradients also play a significant role in such systems at small scales [21]. For a finite length of the channel, the streamwise spatial evolution of the high-viscosity thread cross-section from non-circular (near-ellipsoidal) to a circular shape was found to be dependent on EIT. Indeed, it is plausible for a high-viscosity non-circular thread to evolve naturally and adapt a circular shape through minimization of energy via viscous dissipation. But, such a natural process was found to be inconsequential compared to the effect of EIT. This was evident from the estimation of time and length scales derived from a simple scaling model [21], wherein, the spatial evolution of the high-viscosity thread from near-ellipsoidal crosssection to circular shape could be quantified by an effective capillary number. Analogous to the capillary number in immiscible systems, the effective capillary number based on Γ e for the non-equilibirum miscible system could be defined as, where η is the dynamic viscosity and Q is the flow rate of the core fluid thread, and h is the width of channel. Thus, when the streamwise spatial coordinate is normalised with the above Ca e , all the spatial evolution of the thread heights at different γ (Γ e ) collapse on to a master curve.
In the present work, we widen the studies of viscous thread formation at high Péclet number in miscible environments into a systematic investigation employing geometrically varying flow-focusing setups with a confluence angle β as x z y Figure 1: (a) -(e) Schematic illustration of the top view of geometrically varying flow-focusing channel configurations with various confluence angles β. The red color represents the core fluid entering the central inlet channel with a volumetric flow rate Q 1 . The light cyan color denotes the sheath fluid entering from the side inlet channels at an angle β with a flow rate of Q 2 /2 in each channel. The central channel has a square cross-section with sidelength h. illustrated in Fig. 1. The Confluence angle β refers to the angle made by the side (sheath) flow channel inlets with the central (core) flow channel inlet. We characterise the flow in three dimensions (3-D) both experimentally and numerically employing optical coherence tomography (OCT) and volume of fluid (VoF) method [31] implemented in the open-source code OpenFOAM [32].
There are three profound objectives for this study. The first is to investigate the effect of confluence angles on the 3-D shape of the thread structures. Here, we aim to understand the influence of sheath fluid impingement and geometrical channel effects on the thread formation and wetted region morphologies. Through this, we address an important aspect: whether the process of thread detachment from the top and bottom walls of the channel in miscible systems [25,26] occurs (i) naturally through the self-lubrication principle [23], a phenomenon associated with the effect of viscous dissipation of energy and originally observed in core-annular flow involving high-viscosity contrast immiscible fluid pairs [3,33], or (ii) through some other mechanisms with respect to high-viscosity contrast miscible fluid pairs. The effects of confluence angle β could also be potentially useful to understand how the sheath flow momentum affects the system, and identify the means to achieve efficient extensional flows. The second objective is, to exploit the experimentally measured 3-D spatial evolution of the thread shape together with an effective capillary number dependent master curve to estimate an experimentally intractable variable, namely the Γ e between two miscible fluids. The final goal is, to bring out the implications of EIT in microfluidic flow systems. In addition to these, a valuable feature of microfluidic flows is demonstrated, where the 3-D flow characteristics of different confluence angle geometries can be replicated with a geometry of β = 90°by altering the sheath flow channel widths. The present study also illustrates how additional insights into the physical mechanisms acting experimentally in the complex microfluidic systems can be gained from a diligent analysis of numerical calculations. These insights are difficult to ascertain in experiments due to inherent composition-dependent fluid properties.
As the experimental fluids, similar to our previous work [21], we use a colloidal dispersion for the core and its own solvent as the sheath fluid. The ingredients of the colloidal dispersion consists of cellulose nanofibrils (CNF) dispersed in water. Such nanofibrils have been assembled into high-performance structural cellulose filaments via hydrodynamic focusing [34,35]. Understanding the underlying flow behaviour among various geometrical configurations is critical for controlling the hydrodynamic assembly [36,37]. Moreover, these colloidal systems exhibit variant compositiondependent fluid properties based on particle structure and interparticle interactions, and the Γ e for such colloidal dispersion and its solvent system is expected to vary considerably [8,38].
The organisation of the paper is as follows. In Sec. II, we give an overview of the experimental setup and a brief outline of the numerical method. In Sec. III A, we compare and discuss the numerical and experimental 3-D thread topologies and wetted region morphologies on the top and bottom walls for various flow-focusing configurations emphasizing the role of confluence angle β. In Sec. III B, we briefly recall the scaling model, and the master curve, and use it to estimate the Γ e for the present experimental fluids. The estimated Γ e , in turn, is verified with the value of Γ e obtained through an optimization procedure. In Sec. III C, numerical and experimental results of the centreline velocity is compared. Further, a numerical comparison of strain rates along the centreline in different geometrical configurations is undertaken to dissect the effects of confluence angle β on the flow-field. In Sec. IV, we present numerical results of cross-sectional velocity distributions for all the confluence angle geometries. In Sec. V, replicability of the 3-D flow features of different confluence angles β in modified flow-focusing channels with β = 90°i s presented. In Sec. VI, a remark highlighting the significance of EIT in microfluidic channels is elucidated, and finally a brief summary of the conclusions is provided in Sec. VII.

II. Experimental and numerical setups
A. Experimental setup

Flow-focusing geometries
In this work, we employ five distinct types of flow-focusing channel setups. The setups are planar with square cross-sections of sidelength h = 1 mm and have the geometrical configurations as shown in Fig. 1. All the geometries have one main central inlet channel for the core flow and two sheath flow inlets inclined with a confluence angle β. The confluence angle is varied between 30°and 150°with an interval of 30°, so as to systematically investigate the impact of sheath flow impingement and channel effect on the core fluid 3-D thread topology and flow-field.
All the channel geometries are built using a stainless-steel plate of 1 mm thickness similar to our previous works involving flow-focusing configuration with confluence angle β = 90° [21,34,36,37]. The steel channel plate is enclosed on both sides with layers of aluminum plates and a cyclic olefin copolymer (COC) film forming an assembly of 'aluminum plate−COC−steel channel−COC−aluminum plate' sandwich. The fluids are injected at constant volumetric flow rates into the core and sheath flow channel inlets by two syringe pumps (WPI, AI-4000). Furthermore, all the geometrical configurations (Figs. 1(b), 1(c) and 1(d), 1(e)) will be discussed in relative to the reference configuration ( Fig. 1(a)), and this choice will be clarified in Sec. III A.
In all the configurations, the experimental measurements are performed at constant volumetric flow rates (see Fig. 1) with Q 1 = 6.5 mm 3 s −1 and Q 2 = 7.5 mm 3 s −1 , respectively. The core fluid is a colloidal dispersion exhibiting a non-Newtonian viscosity behaviour as shown in Fig. 2. The sheath fluid is a de-ionized (DI) water of viscosity η 2 = 1 mPa s.

Dispersion material and its rheological properties
The colloidal dispersion is composed of cellulose nanofibrils (CNF) suspended in water. Cellulose nanofibrils were obtained by liberating fibrils from never dried sulfite softwood pulp (Domsjö Fabriker AB, Sweden) supplied by RISE (Research Institute of Sweden). Before defibrillation, never dried sulfite softwood pulp was subjected to TEMPOmediated oxidation following the protocol described elsewhere [39,40]. Thereafter, defibrillation of the oxidized pulp was carried out by passing through a high pressure (1600 bars) Microfluidizer (M-110EH, microfluidics) with 400/200 µm (1 pass) and 200/100 µm (4 passes) wide chambers connected in series. The resulting output is a cellulose nanofibril dispersion of 1 wt%. Finally, a transparent colloidal dispersion of concentration 3 g dm −3 was obtained by further dilution from 1 wt% to 0.3 wt% and homogenization by an Ultra-Turrax dispersing tool (IKA, Sweden) for 10 min at 12000 revolutions per minute. The typical fibril lengths L vary from 100−1500 nm, and the average fibril diameter d is 2.3 ± 0.7 nm as determined by Atomic Force Microscopy (AFM, Dimension 3100 SPM, Veeco, USA) measurements.
The rheological characterization of the colloidal dispersion was performed using a bob and cup Kinexus pro+ rheometer (Malvern). This rheometer is well suited for accurate measurements of viscosity (indicated by red dots) over a range of shear rates, as observed in Fig. 2. The dispersion displays a non-Newtonian shear thinning behaviour [41][42][43].
The rheological data can be described well by a Carreau model where η eff is the shear viscosity, η inf the infinite shear viscosity, η 0 the zero shear viscosity, τ the relaxation time, γ the shear rate, and n the power index. The parameters of the Carreau model fit denoted by the solid black line in Fig. 2 are as follows: η inf = 12 mPa s, η 0 = 4500 mPa s, τ = 1.306 s, and n = 0.16. The shear thinning behaviour is the consequence of micro-structural re-arrangements observed in the colloidal dispersion resulting from fibril/solvent molecular interactions, electro-static interactions or due to effects of Brownian motion and fibre entanglement [44]. As the dimensions of the fibrils are very small, Brownian motion along with the entanglement effects is dominant over other interactions. At low shear rates, the micro-structure of the fibrils are in disordered or isotropic arrangement. Once the shear rates are high enough to overcome Brownian effects, fibrils tend to reorganize by aligning and orienting towards the flow direction. As a result, the viscosity starts to decrease at a particular critical shear rateγ crit indicating the onset of shear thinning. The critical shear rate can be obtained by taking the inverse of the orientational relaxation time of fibrils,γ crit = τ −1 r where τ r =1/(6D r ) [45]. The orientational relaxation time of a Brownian system can be estimated by calculating the rotational diffusion coefficient D r , which is a measure of the rate at which a anisotropic system relaxes towards isotropy.
The rotational diffusional coefficient for a fibril of length L in a polydisperse system close to isotropy, could be appoximated as [36,37,46,47] whereβ is a numerical factor 1000 [48][49][50], the Boltzmann constant k B = 1.38 × 10 −23 J K −1 , the temperature T = 300 K, the solvent viscosity η s = 1 mPa s. The entanglement length L * is defined as [46,47] wherec is the concentration distribution dependent on fibril length L. The entanglement length L * indicates if the fibrils are in a dilute (not entangled, cL 3 1, c being concentration of fibrils , L < L * ) or semi-dilute regime (1 cL 3 L/d, entangled, L > L * ). More details related to L * can be found in Refs. [36,37]. For the present colloidal dispersion, L * ≈ 45 nm, and all the fibril lengths L > L * , fall in the semi-dilute regime. Substituting all the values in Eq. (4), for an average fibril length L ∼ 750 nm, the rotational diffusional coefficient close to isotropy becomes D r ≈ 0.126 rad 2 s −1 . Thus,γ crit is approximated to be aroundγ crit = 6D r 0.80 s −1 as depicted in Fig. 2.
Thus, the low shear viscosity i.e. zero shear viscosity η 0 of the colloidal dispersion is considered near to the critical shear rateγ crit estimated as per the Eq. (4) (see Fig. 2). The viscosity ratio between the core and sheath flow defined as χ = η 1 /η 2 = 4500 is arrived based on the zero shear viscosity η 0 of the core fluid i.e. η 1 = η 0 = 4500 mPa s, and η 2 = 1 mPa s, being the viscosity of sheath fluid.

Péclet number estimation
In order to verify whether a de-facto interface persists in the experimental system, time scale analysis is carried out by estimating the Péclet number based on the characteristic length scale h of the channel system [5,22]. The translational diffusion coefficient D for Brownian fibrils of average length L ∼ 750 nm and diameter d ∼ 2.3 nm diffused in a solvent of water is given as [45] where the temperature T = 300 K, the Boltzmann constant k B = 1.38 × 10 −23 J K −1 , and the solvent viscosity η = 1 mPa s. Substituting all the values, D becomes approximately 5 × 10 −12 m 2 s −1 . The Péclet number is estimated as P e = U h/D ≈ 2 × 10 6 for an average flow velocity of U ≈ 10 mm s −1 . Thus, as the Péclet number is very large, the time scale for the interdiffusion between the colloidal dispersion and its solvent is almost negligible compared to the convection time scale of the two fluids in the channel. Therefore, in such an event, a sharp de-facto interface between the two fluids is likely to exist in the experiments [3,5]. In the context of the present experimental fluids, at such a de-facto interface, EIT is expected to be present [8,38]. In addition, the rheological behaviour of the colloidal dispersion is also controlled by the fact that the nanofibrils form a percolating volume spanning arrested state at very low concentrations [51]. This, in turn, makes the diffusion of nanofibrils into the surrounding sheath flows very slow compared to the time scales of the dynamics in the flowfocusing channel, making the Korteweg stresses long-lived.

Data acquisition method
Three-dimensional core fluid thread topologies and the velocity field measurements are carried out by employing an light-based spectral domain optical coherence tomography. Optical coherence tomography (OCT) is a non-invasive volumetric imaging technique that uses a broadband light source, and operates based on the principle of low-coherence interferometry [52]. Utilising the Doppler principle, OCT can simultaneously capture the structural properties as well as the motion of opaque and turbid media sample with micron-level resolution [53,54]. The wavelength of the light source of the spectral domain OCT used in the present work is 1310 nm with a bandwidth of 270 nm, and a resolution of ∼3 µm in both the axial and transverse directions. More details related to the working principle of OCT, subsequent data acquisition and processing employed for the present study is described in Ref. [21].
In the present work, the 3-D experimental measurements are performed for all the flow-focusing configurations illustrated in Fig. 1 with the above mentioned fluid properties. These experimental measurements will be used in Secs. III A and III C for cross-comparison with the numerical observations, and in Sec. III B to measure the experimentally acting Γ e between the colloidal dispersion-solvent system.

B. Numerical setup
The numerical computations have been performed by utilising a recently developed finite volume based geometric volume of fluid (VoF) method, an interface advection algorithm called isoAdvector [31,55]. The implementation of the algorithm is incorporated in interIsoFoam, a two-phase incompressible immiscible open-source flow solver which is a part of the OpenFOAM ® community [32]. The algorithm accurately captures and advects the sharp interface, a key aspect in the numerical computation of multiphase flows. The rationale behind choosing the immiscible fluid solver was elucidated in the introduction, and will be again clarified in the upcoming Sec. III B. The set of equations being solved for an immiscible system of two fluids are, the continuity equation the Navier-Stokes equation together with the continuum representation of an interfacial tension force F s [56] where γ is the interfacial tension and κ is the curvature of the interface and the equation for the advection of phase/volume fraction α ∂α ∂t + ∇ · (αU ) = 0.
Here T represents the deviatoric stress tensor, U is the velocity vector field and p is the pressure field. The density ρ b and viscosity µ b are computed as based on the weighted average distribution of the volume fraction α of fluid where ρ 1 , ρ 2 , µ 1 , µ 2 are the densities and the viscosities of the two fluids, respectively.
In the present study, 3-D numerical computations are performed for the geometrical configurations illustrated in Fig. 1. At the channel walls, a no-slip velocity boundary condition and a contact angle of θ = 0 • is imposed for the phase field. A uniform velocity flow profile is prescribed at the core and sheath flow channel inlets based on the flowrates of Q 1 = 6.5 mm 3 s −1 and Q 2 = 7.5 mm 3 s −1 . At the channel outlet, pressure is set to atmospheric, and zero gradient for the volume fraction. The non-Newtonian Carreau model depicted in Fig. 2 is implemented for the rheology of the core fluid while the viscosity of water η 2 = 1 mPa s is set for the sheath fluid.
The numerical computations are performed on Cartesian meshes. The size of the computational domain comprises of three inlet channels of length 5h and an outlet channel of length 30h with a square cross-section of width h = 1 mm. The inlet channel domains are discretized with a unidirectional grid along the channel length and an equidistant grid spacing across the cross-section. The outlet channel domain has an equidistant grid spacing of ∆ = 2×10 −5 m (∆x = ∆y = ∆z) both along the channel length and across the square cross-section. An adaptive time step method is used to achieve the stability and convergence of the computations. The computations were run on 256 processors, and the simulations took 48 h. For additional details concerning to numerics and grid convergence, we guide the interested reader to look upon Ref. [21].
In application to two-phase microfluidic flows via hydrodynamic focusing, the solver has been thoroughly tested and validated through an extensive comparison with the experimental data. Excellent qualitative and quantitative agreement between numerical and experimental results [21] is demonstrated, along with the capture of diverse flow regimes observed by Ref. [57] in a typical microfluidic flow-focusing setup.

III. Results and discussion
We first carry out a detailed numerical and experimental investigation of the effect of confluence angle β on the 3-D thread morphology and shape of the wetted region. Then, in Sec. III B, experimentally measured thread height ε z /h as in the case of reference flow-focusing configuration (β = 90°) is compared with the master curve to estimate the Γ e acting between the present experimental fluids.
The flow rates and rheologies of fluids in the computations and experiments are set as given earlier in Secs. II A 1 and II A 2.
A. Morphology of the thread and shape of the wetting region In this section, we compare the 3-D experimental thread shape measured with OCT for geometrically varying flowfocusing configurations described in Sec. II A, and those obtained by numerical computations for the corresponding configurations. All the numerical computations are performed with the EIT, Γ e = γ = 0.615 mN m −1 . This choice will be motivated in the upcoming Sec. III B.
Firstly, a few generic features applicable to all the geometrical configurations are noted before moving into a  measured height ε z /h (red curves) and width ε y /h of the thread (blue curves) with the numerically computed thread height and width (denoted by various color curves as per the respective configuration) depicted in Figs. 6(a) -6(e).
In view of the above observations, it is apparent that the confluence angle β has a significant influence both on the wetted region morphology and the thread development. Indeed, characterising the influence of β is an important aspect from the hydrodynamic assembly of nanofibrils point of view [34], since the shape of thread regulates the cross-section of an assembled material. This, in turn, could enhance the scope for synthesizing materials of complex shapes such as rods, ellipsoids, discs and so on. Hence, in what follows, is a systematic analysis of the effect of confluence angle β.
The flow rates, fluids rheology and channel cross-sectional width h across all the flow-focusing configurations are held fixed. So, the parameter that varies due to the effect of β, is the sheath flow momentum in the parallel and normal directions to the core flow at the confluence region. Therefore, a qualitative understanding can be gained by centering on the sheath fluid momentum and in turn, its effect on the core fluid thread topologies.
Beginning with the reference configuration (β = 90°), the sheath flows are perpendicular to the core flow. Accordingly, sheath flow momentum is acting only normal to the core flow over the confluence region 0 ≤ x/h ≤ 1, and the sheath fluid impinges the core fluid with maximum impact in the normal direction as seen in Fig. 3(a). As a result, the colloidal thread detaches with a shorter wetted length L w /h ≈ 1.8 ( Fig. 4(a)) compared to other configurations (Figs. 4(b) -4(e)). In fact, as it can be seen from Fig. 6(a), even the thread width ε y /h decreases much faster than the height ε z /h up to x/h ≈ 6, highlighting the extent of impact of sheath flow momentum along the streamwise thread development. Far downstream the width attains a constant value. However, after the thread detachment, the effect of sheath flow momentum on the height is less significant. Furthermore, since the sheath flow momentum is maximized in the direction normal to the core flow when β = 90°, all the upcoming geometrical configuration discussions will be relative to this reference configuration. In particular when the sheath flow momentum normal to the core flow is taken into consideration, and hence the choice of geometrical configurations placement in the order as shown in Fig. 1 were opted.
Continuing to the β = [60°, 120°] pair of configurations, the sheath flow momentum acts both in parallel and normal to the streamwise core flow at their respective confluence regions as seen in Figs. 3(b) and 3(c). Most importantly, both the configurations have the same length of confluence regions i.e. 0 ≤ x/h ≤ 1.15. However, in the β = 60°case, the sheath flow is impinging the core flow in the same direction as the streamwise core flow, which can be refered to as positive impingement, while in the β = 120°case, the sheath flow is impinging opposite to the direction of the streamwise core flow, i.e. a negative impingement. For the purpose of clarity, top views of both the configurations is illustrated in Figs. 1(b) and 1(c), which gives a better notion of positive and negative impingement (following the signs of arrow representing the core and sheath flows). Surprisingly, the thread morphology and shape of the wetted region for these two cases closely resembles each other as seen in Figs. 3(b), 4(b) and 3(c), 4(c) in spite of difference in the fundamental impingement direction. As it can be observed from Figs. 4(b) and 4(c), in both the cases, the wetted area increases and the wetted length extends up to L w /h ≈ 2.1, which is much longer than the reference configuration L w /h (Fig. 4(a), L w /h ≈ 1.8). This indicates that the main effect of β is through the sheath flow momentum normal to the core flow, which has been weakened by about 13.4% in relative to the reference configuration in both cases. Further, during the development of the thread width ε y /h and height ε z /h as a function of downstream positions x/h, both the configurations exhibit similar characteristics as displayed in Figs. 6(b) and 6(c). In both cases, the width and height of the thread decay slightly faster up to x/h ≈ 6 and thereafter, stays almost constant far downstream with an elliptical cross-section as visualized in Figs. 5 (panels (c),(d) and panels (e),(f)).
On the other hand, for the β = [30°, 150°] pair, the thread shape and wetted region morphologies differ substantially as viewed in Figs. 3(d), 4(d) and 3(e), 4(e). Figures 3(d) and 3(e) show the length of confluence regions 0 ≤ x/h ≤ 2 being larger than in other configurations (Figs. 3(a) or 3(b), 3(c)). Notably, in both cases, there is an expansion of the colloidal thread symmetrically in the transverse or sheath flow direction (y-direction) near the confluence region. This could be due to the decrease in the impact of sheath flow momentum normal to the core flow. Compared to the reference configuration, the magnitude of normal sheath flow momentum is reduced by around 50% due to the effect of β. In the case of postive impingement of sheath flow (β = 30°) there is a small expansion of the wetted area first and the wetted length extends up to L w /h ≈ 3.1 as depicted in Fig. 4(d). However, for the negative impingement (β = 150°), the wetted area expands considerably by curving outward in the transverse direction with a much shorter wetted length L w /h ≈ 2.1 as seen from Fig. 4(e).
The evolution of streamwise thread development in Figs. 6(d) and 6(e) shows that the width ε y /h of the thread in the confluence region 0 ≤ x/h ≤ 2 overshoots up to ε y /h ≈ 1.5 (β = 30°case) and even more for the negative impingement (β = 150°) before the decay, highlighting the expansion of the thread in the tranverse direction. Another striking behaviour is for the β = 150°configuration, where the experimental cross-section of the thread is somewhat different from the elliptical shape as illustrated in Figs. 5(i), (j). In fact, the numerical curve shows an overprediction in the upstream near the confluence junction during the expansion. This feature of experimental thread cross-section in Fig. 5(i) is not accuratley captured by the numerics.

B. Estimation of effective interfacial tension
In this section, we first present an overview of the model used to obtain a master curve that was proposed at the conclusion of our previous study in Ref [21]. Then, we fit the master curve to an exponential function, and employ it here to estimate the Γ e between the present experimental fluids. Figure 7(a) summarizes our observations performed through an extensive set of numerical computations in Ref [21]. These computations were performed with the reference flow-focusing setup for a set flow rate and the given rheology of the fluids [21]. As mentioned in the introduction, employing an immiscible fluid solver, we had demonstrated that soon after the core fluid detachment from the top and bottom channel walls of the channel at x = L w , the thread height varied from non-circular (near-ellipsoidal) cross-section to a circular shape over a range of interfacial tensions 0.024 mN m −1 < γ < 3.000 mN m −1 . The larger the value of γ, the faster the thread converges towards a circular shape for which the thread height attains a constant value of ε z /h ∼ 0.63.
Following [21], using appropriate notations, the typical time and length scales for a near-ellipsoidal thread to approach circular shape can be written as, where U , Q 1 and η 1 are the velocity, flow rate and dynamic viscosity of the core fluid. δP is the pressure gradient dependent on the thread geometry and is proportional to the interfacial tension γ. Comparing Eq. (2) and Eq. (13), η 1 Q 1 /(γh 2 ) ηQ/(Γ e h 2 ) = Ca e , the effective capillary number where the modelling interfacial tension γ ≡ Γ e (EIT) in experiments. Thus, from Eqs. (12) and (13), a scaled downstream length x * can be obtained by renormalizing the downstream length (x − L w ) with l r leading to x * = (x − L w )/(Ca e × h), respectively. As depicted in Fig 7(a), all the simulated thread heights (ε z /h) at various interfacial tensions γ collapse well on a master curve when plotted with the scaled co-ordinate x * .
The master curve depicted in Fig. 7(a) can be best fitted to an exponential function denoted by dashed-white line with the fitting parameters a, b and c given in Table I.  (Table I)  To estimate the Γ e acting between the present colloidal dispersion -solvent system, we use the spatial evolution of thread height ε z /h measured experimentally with the reference flow-focusing configuration (β = 90°) as depicted in Fig. 7(b), and compare with the exponential fit representing the master curve as shown in Fig. 7(c). By solving Eq. (14) together with the fit parameters reported in Table I, the scaled downstream length x * is obtained. As noted from above, the expression for x * is given as x * = (x − L w )/(Ca e × h). Substituting all the experimentally measured variables in x * and Ca e such as the core fluid flow rate Q 1 = 6.5 mm 3 s −1 , viscosity of the core fluid, η 1 = 4500 mPa s, the wetted length L w ∼ 1.8h along with the streamwise downstream positions x of the thread height, and the channel  representing the master curve. The estimated Γ e acting between the present experimental fluids is obtained by utilising the spatially measured experimental thread height ε z /h (β = 90°) (shown in Fig. 7(b)) and Eq. (14) along with the fit parameters, and other experimentally acquired quantities such as viscosity (η 1 ), flow rate (Q 1 ) of the core fluid and the channel width h.
0.240 ± 0.0008 2.950 ± 0.015 0.622 ± 0.0006 0.756 ± 0.0005 width h, the estimated Γ e 0.765 mN m −1 can be retreived. Furthermore, as a verification, the numerical computations were performed in close proximity to the estimated Γ e utilising the reference flow-focusing configuration (β = 90°). The difference between the numerical and experimental ( Fig. 7(b)) thread height results are plotted in Fig. 7(d) as an error, which is defined as, where i = 1 to N are the downstream positions at which thread heights ε z /h are evaluated after the thread detachment. The minimum of δ εz (Γ e ) is indicated by a filled red symbol in Fig. 7 (d) and occur at Γ e = 0.615 mN m −1 which, in turn, affirms the estimated Γ e to be good with fairly accurate order of magnitude. This value of Γ e = γ is used in the numerical simulations of all the geometrically varying flow-focusing setups as discussed in the previous Sec. III A. An interesting observation worth noting here is, that the value of Γ e (O∼10 −1 mN m −1 ) for the present study colloidal dispersion-solvent system is almost one decade higher in magnitude than our previously evaluated value (O∼10 −2 mN m −1 ) for a similar, but not identical colloidal dispersion-solvent system [21] at the same concentration. Both the colloidal dispersions exhibit non-Newtonian shear thinning behaviour, and the zero shear viscosity for the present case is η 1 = η 0 4500 mPa s while it was 1750 mPa s in Ref. [21]. The order of magnitude difference in Γ e could be attributed to the variation in length fraction of nanofibrils, interfibril interactions leading to different rheological properties of the colloidal dispersion, in turn to that on variant dispersion-solvent de-facto interface properties. Indeed, this also concurs with the observations of experimental studies performed by Refs. [8,38], where Γ e between colloidal dispersions and its own solvent were measured, and the variations in Γ e span over 5 decades for mild changes in the composition.

C. Centreline velocities
In addition to the evolution of thread shape, quantitative understanding of the flow field behaviour is vital for the modeling and prediction of hydrodynamic alignment of nanofibrils [58]. Alignment of nanofibrils is a key factor in controlling and tuning the material properties [59] of assembled materials. Figure 8 shows the velocity variation along the centreline of the colloidal thread as a function of downstream positions x/h for all the five confluence angle geometries. The agreement between numerical computations and experimental measurements is excellent. Dashed vertical lines (magenta) in all the panels (Figs. 8(a) -8(e)) mark the confluence regions. The trend of the centerline velocity (U c ) variation is more or less similar among all configurations. First, the velocity is constant in the inlet channel before the confluence region (x/h < 0). There is a slight deceleration right after the confluence point at x/h = 0 that is followed by a rapid acceleration before a high steady value at far downstream positions (x/h ≥ 10) is reached. However, a careful inspection of Fig. 8 and Fig. 9 unveils a considerable variation in the minimum and maximum centreline velocity, degree of deceleration and acceleration for different confluence angles β. These variations will now be studied in detail. Figure 9 shows the variation of minimum (U c,min ), maximum (U c,max ) centreline velocity and strain rateε ( dU /dx ) along the centreline for all the geometrical configurations. In the case of the reference configuration (β = 90°), close to the confluence point x/h 0, the deceleration is ascertained to be small as observed from Figs. 8(a) (confluence region) and 9(c) (cyan color). As a consequence, the minimum centreline velocity U c,min is higher (meaning lower deceleration) than the other configurations as seen from Fig. 9(a). Subsequent to deceleration, a substantial increase in the velocity caused by the sheath flow momentum normal to the core flow leads to acceleration of the core flow. As depicted in Fig. 9(c),ε (cyan color) continues to increase till the end of the confluence region (x/h ∼ 1), and thereafter,ε decays to 0 at around x/h ≈ 5 since the velocity is almost constant as observed in Fig. 8(a). However, far downstream, as seen in Fig. 9(b), U c,max is lowest as compared to all the other configurations.
For the β = [60°, 120°] pair, the velocity behaviour duplicate each other as observed from Figs. 8(b) and 8(c) and more clearly from the strain rateε plot in Fig. 9(c). It is worth to recall that this similarity of the flow mirrors the similarity of the thread shape evolution as seen earlier from Figs. 3(b), 3(c) to Figs. 6(b), 6(c). In both cases, the effect of β appends the streamwise flow through the contribution by sheath flow momemtum acting parallel to the core flow. As a result, the core fluid velocity increases much faster than the reference configuration velocity. For both cases, U c,min and U c,max are symmetric with respect to β = 90°as seen in Figs. 9(a) and 9(b). However, U c,min is lower for the β = 90°case (implying higher deceleration) and U c,max is highest for [60°, 120°] cases when compared to all the other configruations. As observed from Fig. 9(c), in both configurations, the variation inε (golden and green color) shows a similar trend as the β = 90°case (cyan color) but the magnitude is higher near the end of the confluence regions (x/h ∼ 1.15) and from there onwards, the decay is much slower upto x/h ≈ 5. On the contrary, Figs. 8(d) and 8(e) display a quite different behaviour for the β = [30°, 150°] pair. Even though both the configurations have the same confluence region 0 ≤ x/h ≤ 2, there is a considerable variation in the deceleration and acceleration regions as seen from Fig. 9(c). For the β = 30°case, the deceleration (x/h 0) is higher leading to the lowest U c,min among all the configurations and also lower U c,max far downstream as seen in Figs. 9(a) and 9(b), whereas in the β = 150°case, deceleration is slightly lower as compared to the β = 30°case and U c,max is also higher. This results in an asymmetry of U c,min and U c,max with respect to the β = 90°case. As seen from Fig. 9(c) for β = 30°, the strain rateε (purple color) is more steep than for the β = 150°case leading to a lower peak value close to the end of confluence region (x/h ∼ 2). For x/h > 2 onwards, the decay ofε in both cases follows a similar trend.
Thus, having affirmed the agreement between experimental measurements and numerical computations, we will now utilise purely the numerical computations for further analysis of the results in the sections that follow. In the next section, we look at the cross-sectional velocity distribution for all the geometrical configurations with varying β.  in the cross-sectional planes, albeit they are occupied by a region of high-velocities. However, the magnitudes of these high-sheath velocities seem to differ significantly for different confluence angles β. Accordingly, in the case of negative impingement configurations i.e. β = 120°and 150°, the velocity magnitudes are relatively higher and so is the length of 2-D vectors as seen in Figs. 10(c) and 10(d). They are maximum for the β = 150°case. The velocity profiles are nearly flat along the z-direction at y/h = 0 in all configurations inside the core fluid thread. Also, along the y-direction at z/h = 0, a flat profile is observed in the core with overshoots of the sheath fluid velocity close to channel walls (with an exception for the β = 150°case). Far downstream, as observed from Figs. 11(a) -11(e) for x/h = 20, the velocity profiles within the colloidal thread stays uniform in both the y and z-directions signalling the plug flow behaviour. Close examination of Figs. 11(b) and 11(c) display higher plug velocities for the β = [60°, 120°] pair amongst all the configurations, corroborating with earlier observations of maximum centreline velocity U c,max for the same cases in Fig. 9(c). Furthermore, the 2-D vectors indicate the in-plane velocities for all the configurations to be very weak at this downstream position. Another observation in all the configurations is that the velocity profiles near to the channel walls are parabolic in both the y and z-directions highlighting the sheath flows are wrapping the colloidal thread from all the four sides.

V. Replication of confluence angle effects
In this section, a strategical approach was attempted numerically to comprehend the confluence angle effects. Two computations were performed utilising the reference flow-focusing configuration (β = 90°) with varying sheath flow inlet channel widths as illustrated in Fig. 12(a). Accordingly, the side channel width (SCW) of the 'SCW-1.15h' configuration was set to 1.15h matching the confluence region (0 ≤ x/h ≤ 1.15) of the β = [60°, 120°] configurations. Similarly, the 'SCW-2h' configuration width was fixed to 2h tallying with the confluence region (0 ≤ x/h ≤ 2) of the β = [30°, 150°] configurations. Other parameters such as flow rates of the core (Q 1 ) and sheath (Q 2 ) flows, rheologies of fluids are held fixed as described in the Secs. II A 1 and II A 2. The interfacial tension γ = Γ e = 0.615 mN m −1 , and the cross-sectional width h of the central inlet and outlet channel were unchanged. This means that the sheath flow channel inlets are now rectangular in these cases.
Astonishingly, the flow in the SCW-1.15h configuration precisely replicates the morphological features of the wetted region and thread topologies of the β = [60°, 120°] configurations as observed from Figs. 12(b) and 12(d). Once again this highlights the symmetry between these two cases. Whereas, the SCW-2h configuration captures and replicates the features of only the β=150°configuration as depicted in Figs. 12(c) and 12(e) (remember the asymmetry between the β = [30°, 150°] configurations).
In the three cases, 'SCW-1.15h' and β = [60°, 120°], the sheath flow momentum normal to the core flow is attenuated by 13.4% in relative to the reference configuration (β = 90°). For the 'SCW-2h' and β = [30°, 150°] pair, the attenuation is around 50%. Therefore, inclining the sheath flow channels at various confluence angles β or setting the width of sheath flow channels of reference configuration to confluence region lengths, would essentially mean inducing the same magnitude of sheath flow momentum normal to the core flow. Thus, from these demonstrations, it is aptly clear that the magnitude of the sheath flow momentum at the confluence junction is the primary factor in controlling the morphological features of the wetted region and thread topologies except for near co-flow (β = 30°) case.

VI. A comment on effective interfacial tension and thread detachment
The comprehensive systematic demonstrations utilising geometrically varying flow-focusing setups in Sec. III A together with the alternative approach in previous Sec. V underscores the significant role of Γ e . As noticed, for a set flow rate and specified rheologies of fluids, the same value of EIT (γ = Γ e = 0.615 mN m −1 ) deduced using the reference flow-focusing configuration (β = 90°), very well reproduces the 3-D flow characteristics for the respective geometrical configurations in the numerics. All the features such as the thread evolutions, morphologies of the wetted regions, and velocity fields showed good agreement with the experimental findings in Secs. III A and III C. In addition, the replication of 3-D flow features of various confluence angle geometries captured via the modified geometries (β = 90°), numerically with the same interfacial tension, further substantiate the robustness of our numerical procedure.
An important and interesting aspect from the above illustrations is on the detachment of thread from the top and bottom channel walls. We find that, the above value of interfacial tension γ in the numerical computations not only captured the wetted region shapes and wetted lengths L w /h accurately, but also steered the thread detachment from the channel walls. Irrespective of the geometrical configurations, we also observed through a separate set of simulations that are not displayed here, that the wetted length L w /h → ∞ if γ → 0, implying no detachment of thread from the top and bottom channel walls. Therefore, presence of ultra-low interfacial tension γ (Γ e ) is a decisive factor in order for the thread detachment to occur. Nevertheless, given the presence of high viscosity contrast (χ 4500) between the colloidal dispersion (core fluid) and its solvent (sheath fluid), Cubaud and Mason [23] predict that the process of thread detachment in miscible systems having large viscosity contrast [25,26], could occur naturally through the phenomena of self-lubrication [33] depending on the fluid injection geometries and injection flow rates. According to this principle, the low-viscosity fluid enwraps the high-viscosity fluid to minimize the viscous dissipation of energy, leading to self-lubricated thread structures. However, as evident from the above observations, in high Péclet number microfluidic flows involving miscible systems with high-viscosity contrast, thread detachment is found to be clearly controlled by Γ e as well. Indeed, in order to have a finer understanding of the effects of EIT versus the ones associated with self-lubrication on the thread detachment, further detailed investigation is needed. For instance, variation of EIT, viscosity contrast and flow rate ratios, and how these parameters affect the wetted length L w /h and thereby the thread detachment is worthwhile, and such investigations are beyond the scope of this present work.

VII. Conclusions
The effective interfacial tension Γ e acting between the colloidal dispersion and its own solvent has been estimated utilising the experimentally measured 3-D spatial evolution of thread shape in the flow-focusing channel and an effective capillary number dependent master curve obtained from a simple scaling model [21]. We have verified from the 3-D numerical computations, that the Γ e which minimizes the error between computed and experimental thread heights is close to the estimated Γ e utilising the above method. The numerical computations were able to capture our experimental observations with good quantitative and qualitative agreement (both in terms of flow-topologies and flow-fields) for a range of confluence angles. By mapping the numerical observations with the experimental measurements, we gained insights into the influence of various physical mechanisms driving the process of thread detachment from the top and bottom walls in the channels. From these meticulous analyses, we find that in the case of high viscosity contrast miscible systems [25,26], the thread detachment occuring in the physical experiments of microchannels may not be entirely based on the phenomena of self-lubrication [23] related with effect of viscous dissipation of energy, but also the effect of EIT induced by composition gradients plays a major role.
Moreover, these ultra-low interfacial tensions in non-equilibrium miscible fluids are extremely difficult to evaluate in experimental measurements. The standard experimental techniques yielding reasonable results in the case of molecular miscible fluids are light scattering and spinning drop tensiometry [20]. However, these two techniques are observed to have drawbacks, in particular at the transition region due to mixing of the two fluids and difficulty to distinguish between the bulk fluids and interfacial region [20]. The recently explored new measurement strategies like examining the Saffman-Taylor instability in a Hele-Shaw cell [7,8] with miscible complex fluids, and studying the dynamics of drop shape in spinning drop tensiometry [16] tested with miscible molecular systems are noteworthy developments. In addition to the above strategies, our method is also promising and could be a suitable choice.
From a technical perspective, varying the effect of confluence angle β assists in the utilisation of sheath flows to generate effective extensional flows. This, in turn, could facilitate to achieve optimal alignment of nanofibrils in colloidal dispersions [36,37]. A detailed investigation on the orientation and alignment of nanofibrils can be pursued, utilising the velocity gradients of various flow configurations, including other relevant factors such as rotational diffusion, mobility, and rigidity of fibrils. Concluding, a thorough understanding of effects like EIT in miscible systems opens up possibilities in tuning, and controlling the material properties more efficiently via appropriate selection of microfluidic geometrical configuration based on the practical interest and application.