Formation of colloidal threads in geometrically varying flow-focusing channels
V. Krishne Gowda, Cecilia Rydefalk, L. Daniel Söderberg, Fredrik Lundell
FFormation of colloidal threads in geometrically varying flow-focusing channels V . Krishne Gowda , ∗ Cecilia Rydefalk , L . Daniel S¨oderberg , , and Fredrik Lundell , FLOW, Dept . of Engineering Mechanics, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden Dept . of Fibre and Polymer Technology, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden Wallenberg Wood Science Center, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden
When two miscible fluids are brought into contact with each other, the concentration gradientsinduce stresses. These are referred to as Korteweg stresses and are analogous to interfacial tensionbetween two immiscible fluids, thereby acting as an effective interfacial tension (EIT) in inhomoge-neous miscible fluid pairs. Effective interfacial tension governs the formation of a viscous thread inflow-focusing of two miscible fluids. To further investigate its significance, we have studied threadformation of a percolated colloidal dispersion focused by its own solvent. Experiments are combinedwith three-dimensional (3-D) numerical methods to systematically expand previous knowledge util-ising different flow-focusing channel setups. In the reference flow-focusing configuration, the sheathflows impinge the core flow orthogonally while in four other channel configurations, the sheath flowsimpinge the core flow at an oblique angle that is both positive and negative with respect to thereference sheath direction. As an initial estimate of the EIT between the colloidal dispersion-solventsystem, we fit the experimentally determined thread shape in the reference configuration to a mas-ter curve that depends on EIT through an effective capillary number. By numerically reproducingthese 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 mea-surements. Regardless of channel configurations, further numerical calculations performed usingthe optimal EIT evaluated from the reference configuration show good agreement with the experi-mental findings in terms of 3-D thread shapes, wetted region morphologies, and velocity flow fields.Through an additional numerical study, we demonstrate that the 3-D flow characteristics observedexperimentally in oblique channel configurations can be precisely replicated with orthogonal sheathflows if the sheath channel width is adjusted. The one-to-one comparison and analysis of numericaland experimental findings unveil the crucial role of EIT on the thread detachment from the topand bottom walls of the channel, bringing useful insights to understand the physical phenomenonsinvolved in miscible systems with a high-viscosity contrast. These key results provide the foun-dation to tune the flow-focusing for specific applications, for example in tailoring the assembly ofnanostructured materials.
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 theequilibrium 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 fluidscome into contact, a boundary between them can persist till the time the two fluids eventually form an equilibratedhomogeneous mixture due to diffusive mixing. This boundary zone can act as a de-facto interface with some of itsproperties resembling that of a distinct interface observed between immiscible fluids [1–5]. In 1901, Korteweg [6] firstproposed that when two miscible fluids are brought into contact, the composition inhomogeneties/gradients of thefluid property at the zone of contact gives rise to additional stresses (so-called Korteweg stresses). These stresseseffectively 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) formiscible fluids in non-equilibrium state, commonly denoted as Γ e , can be written as:Γ e = K ∆Φ δ , (1)where K is the Korteweg factor accounting for the relevant interaction effects (e.g. particle-solvent interactions inmiscible complex fluids) between the two fluids [7, 8], δ is the interface thickness and ∆Φ is the variation in compositionor 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. ∗ [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b The role of these transient stresses or EIT in non-equilibrium miscible fluids is evident in numerous multi-fluiddynamical processes occuring at short times. Examples are microgravity experiments [9], evolution of miscibledroplets [10, 11], modeling of hydrodynamic instabilities like viscous fingering or in Hele-Shaw flows [12, 13], sta-bilization 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 arethrough the evolution of drop shape [15, 16], examination of hydrodynamic instabilities [7, 17], and probing of cap-illary waves by light scattering [18, 19]. In spite of these attempts, measuring the EIT between miscible fluids isintrinsically difficult due to ultra-low values (Γ e ∼ − − mN m − ) 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, wehave proposed a methodology with a possibility to determine the EIT between two miscible fluids. In short, Γ e canbe estimated by measuring the spatial evolution of the thread shape formed by these two fluids, and equating it witha 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 schematicallyin Fig. 1(a). For miscible systems, the P´eclet number P e = U h/D , which measures the relative importance ofconvective 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 ofmiscible fluids is passed through the flow-focusing channel, the P´eclet number dictates the flow-patterns. At highP´eclet number under laminar flow conditions, where diffusion is almost negligible, viscous threads are formed, whileat low P´eclet number, the threads undergo diffusive instabilities leading to a wide diversity of flow-patterns [23–25]. Itis even observed that the miscible threads at high P´eclet 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´eclet number system, with a pair of miscible fluids comprising of a colloidal dispersion and its ownsolvent, we have detected and reported that the characteristics of spatial thread evolution could be accurately capturedand modelled as an immiscible fluid problem with a very weak interfacial tension γ for a set of flow rates andspecified rheology of the two fluids [21]. Ideally, at high P´eclet number, the time scale for interdiffusion between thecolloidal dispersion and its own solvent is almost negligible compared to the convective time scale of the fluids in thechannel. In such a scenario, the presence of a sharp de-facto interface due to composition gradients is expected in theexperiments [5, 7, 8]. However, it is extremely difficult to access such a sharp de-facto interface between the two fluidsexperimentally. The miscible viscous thread structures formed at high P´eclet number are in an out-of-equilibriumstate 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 numericalmodelling. In reality, the weak interfacial tension γ accounts for the Korteweg stresses induced by the compositiongradients in the experiments. The experimental observations were reproduced numerically with a very weak interfacialtension γ ( O∼ − mN m − ). Moreover, the magnitude of γ corroborated with the experimental measurements ofprevious studies conducted for miscible complex fluids involving colloidal dispersions or polymer solutions and its ownsolvent [8, 20]. Thus, it was established that the very weak interfacial tension γ acts as an effective interfacial tensionin the experiments with γ ≡ Γ e .Furthermore, in microfluidic systems, the interfacial tension and viscous dissipation of energy dominate over inertialeffects [28–30] due to the small length scales. Pertaining to miscible systems, and in the context of weakly diffusivethreads, the phenomenon of viscous dissipation of energy is pertinent as the viscosity contrast between the fluid pairsis fairly large [23, 26]. In addition to these effects, we noted that the effect of EIT due to concentration gradients alsoplay a significant role in such systems at small scales [21]. For a finite length of the channel, the streamwise spatialevolution of the high-viscosity thread cross-section from non-circular (near-ellipsoidal) to a circular shape was foundto be dependent on EIT. Indeed, it is plausible for a high-viscosity non-circular thread to evolve naturally and adapta circular shape through minimization of energy via viscous dissipation. But, such a natural process was found to beinconsequential compared to the effect of EIT. This was evident from the estimation of time and length scales derivedfrom a simple scaling model [21], wherein, the spatial evolution of the high-viscosity thread from near-ellipsoidal cross-section to circular shape could be quantified by an effective capillary number. Analogous to the capillary numberin immiscible systems, the effective capillary number based on Γ e for the non-equilibirum miscible system could bedefined as, Ca e = ηQh Γ e (2)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 threadheights at different γ (Γ e ) collapse on to a master curve.In the present work, we widen the studies of viscous thread formation at high P´eclet number in miscible environmentsinto a systematic investigation employing geometrically varying flow-focusing setups with a confluence angle β as xz y Figure 1: (a) - (e) Schematic illustration of the top view of geometrically varying flow-focusing channel configurationswith various confluence angles β . The red color represents the core fluid entering the central inlet channel with avolumetric flow rate Q . The light cyan color denotes the sheath fluid entering from the side inlet channels at anangle β with a flow rate of Q / h .illustrated in Fig. 1. The Confluence angle β refers to the angle made by the side (sheath) flow channel inlets withthe central (core) flow channel inlet. We characterise the flow in three dimensions (3-D) both experimentally andnumerically employing optical coherence tomography (OCT) and volume of fluid (VoF) method [31] implemented inthe open-source code OpenFOAM [32].There are three profound objectives for this study. The first is to investigate the effect of confluence angles onthe 3-D shape of the thread structures. Here, we aim to understand the influence of sheath fluid impingement andgeometrical channel effects on the thread formation and wetted region morphologies. Through this, we address animportant aspect: whether the process of thread detachment from the top and bottom walls of the channel in misciblesystems [25, 26] occurs (i) naturally through the self-lubrication principle [23], a phenomenon associated with theeffect of viscous dissipation of energy and originally observed in core-annular flow involving high-viscosity contrastimmiscible fluid pairs [3, 33], or (ii) through some other mechanisms with respect to high-viscosity contrast misciblefluid pairs. The effects of confluence angle β could also be potentially useful to understand how the sheath flowmomentum 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 capillarynumber dependent master curve to estimate an experimentally intractable variable, namely the Γ e between twomiscible fluids. The final goal is, to bring out the implications of EIT in microfluidic flow systems. In additionto these, a valuable feature of microfluidic flows is demonstrated, where the 3-D flow characteristics of differentconfluence 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 thecomplex microfluidic systems can be gained from a diligent analysis of numerical calculations. These insights aredifficult 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 ownsolvent as the sheath fluid. The ingredients of the colloidal dispersion consists of cellulose nanofibrils (CNF) dispersedin water. Such nanofibrils have been assembled into high-performance structural cellulose filaments via hydrodynamicfocusing [34, 35]. Understanding the underlying flow behaviour among various geometrical configurations is criticalfor controlling the hydrodynamic assembly [36, 37]. Moreover, these colloidal systems exhibit variant composition-dependent fluid properties based on particle structure and interparticle interactions, and the Γ e for such colloidaldispersion 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 briefoutline of the numerical method. In Sec. III A, we compare and discuss the numerical and experimental 3-D threadtopologies and wetted region morphologies on the top and bottom walls for various flow-focusing configurationsemphasizing 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 thevalue of Γ e obtained through an optimization procedure. In Sec. III C, numerical and experimental results of thecentreline velocity is compared. Further, a numerical comparison of strain rates along the centreline in differentgeometrical configurations is undertaken to dissect the effects of confluence angle β on the flow-field. In Sec. IV, wepresent 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 ° is presented. In Sec. VI, a remark highlighting the significance of EIT in microfluidic channels is elucidated, andfinally a brief summary of the conclusions is provided in Sec. VII. II. Experimental and numerical setupsA. Experimental setup
1. Flow-focusing geometries
In this work, we employ five distinct types of flow-focusing channel setups. The setups are planar with squarecross-sections of sidelength h = 1 mm and have the geometrical configurations as shown in Fig. 1. All the geometrieshave 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 theimpact 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 worksinvolving flow-focusing configuration with confluence angle β = 90 ° [21, 34, 36, 37]. The steel channel plate is enclosedon both sides with layers of aluminum plates and a cyclic olefin copolymer (COC) film forming an assembly of ‘alu-minum plate − COC − steel channel − COC − aluminum plate’ sandwich. The fluids are injected at constant volumetricflow rates into the core and sheath flow channel inlets by two syringe pumps (WPI, AI-4000). Furthermore, all thegeometrical 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 = 6.5 mm s − and Q = 7 . s − , 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 η = 1 mPa s. Figure 2:
Shear viscosity measurements of the colloidal dispersion represented by red dots. The solid black linerepresents the non-Newtonian Carreau model fit (see Eq. 3), and the vertical dashed-dotted line denotes the criticalshear rate ˙ γ crit for the dispersion.
2. Dispersion material and its rheological properties
The colloidal dispersion is composed of cellulose nanofibrils (CNF) suspended in water. Cellulose nanofibrils wereobtained by liberating fibrils from never dried sulfite softwood pulp (Domsj¨o Fabriker AB, Sweden) supplied by RISE(Research Institute of Sweden). Before defibrillation, never dried sulfite softwood pulp was subjected to TEMPO-mediated oxidation following the protocol described elsewhere [39, 40]. Thereafter, defibrillation of the oxidizedpulp was carried out by passing through a high pressure (1600 bars) Microfluidizer (M-110EH, microfluidics) with400 / µ m (1 pass) and 200 / µ m (4 passes) wide chambers connected in series. The resulting output is a cellulosenanofibril dispersion of 1 wt%. Finally, a transparent colloidal dispersion of concentration 3 g dm − was obtained byfurther dilution from 1 wt% to 0.3 wt% and homogenization by an Ultra-Turrax dispersing tool (IKA, Sweden) for10 min at 12000 revolutions per minute. The typical fibril lengths L vary from 100 − d is 2 . ± . η eff = η inf + ( η − η inf )[1 + ( τ ˙ γ ) ] n − , (3)where η eff is the shear viscosity, η inf the infinite shear viscosity, η 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 inFig. 2 are as follows: η inf = 12 mPa s, η = 4500 mPa s, τ = 1 .
306 s, and n = 0 . γ crit indicating the onset of shear thinning. The critical shear rate can be obtained bytaking the inverse of the orientational relaxation time of fibrils, ˙ γ crit = τ − r where τ r =1 / (6 D r ) [45]. The orientationalrelaxation time of a Brownian system can be estimated by calculating the rotational diffusion coefficient D r , which isa 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 beappoximated as [36, 37, 46, 47] D r ( L ) ≈ ˜ βk B T L ∗ ηL , (4)where ˜ β is a numerical factor (cid:39) k B = 1 . × − J K − , the temperature T = 300 K, the solvent viscosity η s = 1 mPa s. The entanglement length L ∗ is defined as [46, 47] L ∗ = (cid:18)(cid:90) + ∞ ˜ c ( L ) LdL (cid:19) − / , (5)where ˜ c is the concentration distribution dependent on fibril length L . The entanglement length L ∗ indicates ifthe fibrils are in a dilute (not entangled, cL (cid:28) c being concentration of fibrils , L < L ∗ ) or semi-dilute regime(1 (cid:28) cL (cid:28) L/d , entangled,
L > L ∗ ). More details related to L ∗ can be found in Refs. [36, 37]. For the presentcolloidal dispersion, L ∗ ≈
45 nm, and all the fibril lengths
L > L ∗ , fall in the semi-dilute regime. Substituting allthe values in Eq. (4), for an average fibril length L ∼
750 nm, the rotational diffusional coefficient close to isotropybecomes D r ≈ s − . Thus, ˙ γ crit is approximated to be around ˙ γ crit = 6 D r (cid:39) − as depicted in Fig. 2.Thus, the low shear viscosity i.e. zero shear viscosity η of the colloidal dispersion is considered near to the criticalshear rate ˙ γ crit estimated as per the Eq. (4) (see Fig. 2). The viscosity ratio between the core and sheath flow definedas χ = η /η = 4500 is arrived based on the zero shear viscosity η of the core fluid i.e. η = η = 4500 mPa s, and η = 1 mPa s, being the viscosity of sheath fluid.
3. P´eclet number estimation
In order to verify whether a de-facto interface persists in the experimental system, time scale analysis is carriedout by estimating the P´eclet number based on the characteristic length scale h of the channel system [5, 22]. Thetranslational diffusion coefficient D for Brownian fibrils of average length L ∼
750 nm and diameter d ∼ . D = k B T ln( L/d )2 πηl , (6)where the temperature T = 300 K, the Boltzmann constant k B = 1 . × − J K − , and the solvent viscosity η = 1 mPa s. Substituting all the values, D becomes approximately 5 × − m s − . The P´eclet number isestimated as P e = U h/D ≈ × for an average flow velocity of U ≈
10 mm s − . Thus, as the P´eclet numberis very large, the time scale for the interdiffusion between the colloidal dispersion and its solvent is almost negligiblecompared 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 experimentalfluids, 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 nanofibrilsform a percolating volume spanning arrested state at very low concentrations [51]. This, in turn, makes the diffusionof nanofibrils into the surrounding sheath flows very slow compared to the time scales of the dynamics in the flow-focusing channel, making the Korteweg stresses long-lived.
4. Data acquisition method
Three-dimensional core fluid thread topologies and the velocity field measurements are carried out by employing anlight-based spectral domain optical coherence tomography. Optical coherence tomography (OCT) is a non-invasivevolumetric imaging technique that uses a broadband light source, and operates based on the principle of low-coherenceinterferometry [52]. Utilising the Doppler principle, OCT can simultaneously capture the structural properties as wellas the motion of opaque and turbid media sample with micron-level resolution [53, 54]. The wavelength of the lightsource of the spectral domain OCT used in the present work is 1310 nm with a bandwidth of 270 nm, and a resolutionof ∼ µ m in both the axial and transverse directions. More details related to the working principle of OCT, subsequentdata 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 configurationsillustrated in Fig. 1 with the above mentioned fluid properties. These experimental measurements will be used inSecs. III A and III C for cross-comparison with the numerical observations, and in Sec. III B to measure the experi-mentally 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 geometricvolume of fluid (VoF) method, an interface advection algorithm called isoAdvector [31, 55]. The implementation ofthe algorithm is incorporated in interIsoFoam , a two-phase incompressible immiscible open-source flow solver whichis 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 fluidsolver was elucidated in the introduction, and will be again clarified in the upcoming Sec. III B. The set of equationsbeing solved for an immiscible system of two fluids are,the continuity equation ∇ · U = 0 , (7)the Navier-Stokes equation together with the continuum representation of an interfacial tension force F s [56] ∂ρ b U ∂t + ∇ · ( ρ b U U ) = − ∇ p + ∇ · T + F s , (8) F s = γκ ( ∇ α ) , (9)where γ is the interfacial tension and κ is the curvature of the interface κ = − ∇ · (cid:18) ∇ α | ∇ α | (cid:19) , (10)and the equation for the advection of phase/volume fraction α∂α∂t + ∇ · ( α U ) = 0 . (11)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 ρ b = ρ α + ρ (1 − α ), µ b = µ α + µ (1 − α ) based on the weighted averagedistribution of the volume fraction α of fluid where ρ , ρ , µ , µ 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 inFig. 1. At the channel walls, a no-slip velocity boundary condition and a contact angle of θ = 0 ◦ is imposed forthe phase field. A uniform velocity flow profile is prescribed at the core and sheath flow channel inlets based on theflowrates of Q = 6.5 mm s − and Q = 7 . s − . At the channel outlet, pressure is set to atmospheric, andzero gradient for the volume fraction. The non-Newtonian Carreau model depicted in Fig. 2 is implemented for therheology of the core fluid while the viscosity of water η = 1 mPa s is set for the sheath fluid.The numerical computations are performed on Cartesian meshes. The size of the computational domain comprisesof three inlet channels of length 5 h and an outlet channel of length 30 h with a square cross-section of width h = 1mm. The inlet channel domains are discretized with a unidirectional grid along the channel length and an equidistantgrid spacing across the cross-section. The outlet channel domain has an equidistant grid spacing of ∆ = 2 × − m(∆ x = ∆ y = ∆ z ) both along the channel length and across the square cross-section. An adaptive time step methodis used to achieve the stability and convergence of the computations. The computations were run on 256 processors, Figure 3: (a) - (e) 3-D view of the experimental and numerical threads for flow-focusing configurations with varyingconfluence angle β . The green curves at the top planes z/h = 0 . III. Results and discussion
We first carry out a detailed numerical and experimental investigation of the effect of confluence angle β on the 3-Dthread 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 1and 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 flow-focusing configurations described in Sec. II A, and those obtained by numerical computations for the correspondingconfigurations. All the numerical computations are performed with the EIT, Γ e = γ = 0 .
615 mN m − . This choicewill 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 Figure 4: (a) - (e) Top view of the experimental and numerical wetted region morphologies of the core fluid in theplane z/h = 0 . β . The distance from x/h = 0 to thepoint of detachment of core fluid from the top wall is referred to as ‘wetted length’ L w /h , indicated in panel (a). Thewetted region shape and length L w /h vary as per the respective geometrical configuration. These wetted boundariesare also shown in the previous Figs. 3(a) - 3(e) at the top plane z/h = 0 . β . The qualitative agreement between experiments and numericsis ascertained to be quite good. A closer view shows the core and sheath flow channel walls intersect at the point ofconfluence corresponding to x/h = 0. The region beginning from x/h = 0 to the position where the confluence of coreand sheath flow channels end (indicated by dashed magenta lines) is referred to as the confluence region . However, thelength of confluence region varies depending on the confluence angle β . The confluence region is shortest (0 ≤ x/h ≤ β = 90 ° ) and longest (0 ≤ x/h ≤
2) for β = [30 ° , ° ] pair. Meanwhile, at the topand bottom walls of the channel (i.e. at z/h = ± . x/h >
0. The dispersion continues to stay attached to the upper and lower walls upto a pinch-off point.The region originating from x/h = 0 to the colloidal thread detachment point is denoted the wetted region and thesubsequent length is called wetted length L w /h (see Fig. 4(a)). Indeed, both the wetted region and the wetted lengthvary with β as observed in Figs. 3(a) - 3(e) (marked by green curves at the top planes z/h = 0 . β in the upstream while at far downstream the thread shapes seem to be nearly elliptical. The cross-sectional viewsof the elliptical thread at downstream position x/h = 10 are displayed in Figs. 5(a) - 5(e) with its major and minoraxes denoted by ε z /h , ε y /h . Besides, quantitatively, the excellent agreement is exemplified between experimentally0 Figure 5:
Experimental (panels (a),(c),(e),(g),(i)) and numerical (panels (b),(d),(f),(h),(j)) thread cross-sections atdownstream position x/h = 10 for flow-focusing configurations with varying confluence angle β . Black color denotesthe core fluid while white color represents the sheath fluid. Thread width ε y /h (blue) and height ε z /h (red) denotingthe ellipsoid axes is pictured in panel (a).measured height ε z /h (red curves) and width ε y /h of the thread (blue curves) with the numerically computed threadheight 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 onthe wetted region morphology and the thread development. Indeed, characterising the influence of β is an importantaspect from the hydrodynamic assembly of nanofibrils point of view [34], since the shape of thread regulates thecross-section of an assembled material. This, in turn, could enhance the scope for synthesizing materials of complexshapes such as rods, ellipsoids, discs and so on. Hence, in what follows, is a systematic analysis of the effect ofconfluence angle β .The flow rates, fluids rheology and channel cross-sectional width h across all the flow-focusing configurations areheld fixed. So, the parameter that varies due to the effect of β , is the sheath flow momentum in the parallel andnormal directions to the core flow at the confluence region. Therefore, a qualitative understanding can be gained bycentering 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. Accord-ingly, sheath flow momentum is acting only normal to the core flow over the confluence region 0 ≤ x/h ≤
1, and thesheath 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 ≈ . ε y /h decreases much faster than theheight ε z /h up to x/h ≈
6, highlighting the extent of impact of sheath flow momentum along the streamwise threaddevelopment. Far downstream the width attains a constant value. However, after the thread detachment, the effect ofsheath flow momentum on the height is less significant. Furthermore, since the sheath flow momentum is maximizedin the direction normal to the core flow when β = 90 ° , all the upcoming geometrical configuration discussions willbe relative to this reference configuration. In particular when the sheath flow momentum normal to the core flow is1 Figure 6: (a) - (e) Thread width ε y /h and height ε z /h as a function of downstream positions x/h for flow-focusingconfigurations with varying confluence angle β . The experimental data are represented in blue and red color, whilethe numerical data are denoted by various solid color lines shown in the rightmost legend in panel (a) as per therespective geometrical configuration.taken into consideration, and hence the choice of geometrical configurations placement in the order as shown in Fig. 1were opted.Continuing to the β = [60 ° , ° ] pair of configurations, the sheath flow momentum acts both in parallel and normalto 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 ≤ β = 60 ° case,the sheath flow is impinging the core flow in the same direction as the streamwise core flow, which can be referedto as positive impingement, while in the β = 120 ° case, the sheath flow is impinging opposite to the direction of thestreamwise core flow, i.e. a negative impingement. For the purpose of clarity, top views of both the configurationsis illustrated in Figs. 1(b) and 1(c), which gives a better notion of positive and negative impingement (following thesigns of arrow representing the core and sheath flows). Surprisingly, the thread morphology and shape of the wettedregion 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 inthe fundamental impingement direction. As it can be observed from Figs. 4(b) and 4(c), in both the cases, the wettedarea increases and the wetted length extends up to L w /h ≈ .
1, which is much longer than the reference configuration2 L w /h (Fig. 4(a), L w /h ≈ . β is through the sheath flow momentum normalto 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 widthand height of the thread decay slightly faster up to x/h ≈ β = [30 ° , ° ] pair, the thread shape and wetted region morphologies differ substantiallyas 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 ≤ 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 thereference configuration, the magnitude of normal sheath flow momentum is reduced by around 50% due to the effectof β . In the case of postive impingement of sheath flow ( β = 30 ° ) there is a small expansion of the wetted area firstand the wetted length extends up to L w /h ≈ . β = 150 ° ), the wetted area expands considerably by curving outward in the transverse direction with a much shorterwetted length L w /h ≈ . ε y /h of the threadin the confluence region 0 ≤ x/h ≤ ε y /h ≈ . β = 30 ° case) and even more for the negativeimpingement ( β = 150 ° ) before the decay, highlighting the expansion of the thread in the tranverse direction. Anotherstriking behaviour is for the β = 150 ° configuration, where the experimental cross-section of the thread is somewhatdifferent from the elliptical shape as illustrated in Figs. 5(i), (j). In fact, the numerical curve shows an overpredictionin the upstream near the confluence junction during the expansion. This feature of experimental thread cross-sectionin 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 theconclusion of our previous study in Ref [21]. Then, we fit the master curve to an exponential function, and employ ithere 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 rheologyof the fluids [21]. As mentioned in the introduction, employing an immiscible fluid solver, we had demonstrated thatsoon after the core fluid detachment from the top and bottom channel walls of the channel at x = L w , the threadheight varied from non-circular (near-ellipsoidal) cross-section to a circular shape over a range of interfacial tensions0 .
024 mN m − < γ < .
000 mN m − . The larger the value of γ , the faster the thread converges towards a circularshape for which the thread height attains a constant value of ε z /h ∼ . τ = η δP ∝ η γ , (12) l r = U τ ∝ η Q γh . (13)where U , Q and η are the velocity, flow rate and dynamic viscosity of the core fluid. δP is the pressure gradientdependent on the thread geometry and is proportional to the interfacial tension γ . Comparing Eq. (2) and Eq. (13), η Q / ( γh ) (cid:39) ηQ/ (Γ e h ) = 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 downstreamlength ( x − L w ) with l r leading to x ∗ = ( x − L w ) / ( Ca e × h ), respectively. As depicted in Fig 7(a), all the simulatedthread heights ( ε z /h ) at various interfacial tensions γ collapse well on a master curve when plotted with the scaledco-ordinate x ∗ .The master curve depicted in Fig. 7(a) can be best fitted to an exponential function denoted by dashed-white linewith the fitting parameters a , b and c given in Table I. ε z /h = a × e − bx ∗ + c. (14)3 Figure 7: (a): The numerical thread heights ε z /h at various interfacial tensions γ plotted as a function of scaleddownstream length x ∗ (detailed in text). The data collapse on to a master curve. These simulated data are taken fromRef [21]. The master curve can be best represented by an exponential fit [see Eq. (14)] denoted by the dashed-whiteline. The fit parameters a , b and c are reported in Table I. Panel (b) shows the experimental thread height ε z /h asa function of downstream positions x/h measured in the reference configuration ( β = 90 ° ) of the present study. Thedash-dotted vertical blue line in panel (b) indicates the measured wetted length of the thread before the detachment( L w /h ∼ ε z /h is compared with the exponential fit (dashed-whiteline) representing the master curve to estimate the Γ e acting between the present experimental fluids. The light graybackdrop in panels (a) and (c) is drawn for the purpose of better visualization of the plotted data. Panel (d) showsthe percentage variation of the error (defined in text) between the numerical and experimental thread heights for arange of Γ e values near to the estimated Γ e (Table I) obtained from the master curve fit in panel (c). The red markerin panel (d) indicate the Γ e that gives the best match with the experimental measurement.The light gray backdrop in Figs. 7(a) and 7(c) is in place solely for the sake of clarity in the graphics of the plots.The fit parameters are evaluated through non-linear least-square regression performed in Matlab.To estimate the Γ e acting between the present colloidal dispersion - solvent system, we use the spatial evolutionof thread height ε z /h measured experimentally with the reference flow-focusing configuration ( β = 90 ° ) as depictedin Fig. 7(b), and compare with the exponential fit representing the master curve as shown in Fig. 7(c). By solvingEq. (14) together with the fit parameters reported in Table I, the scaled downstream length x ∗ is obtained. As notedfrom above, the expression for x ∗ is given as x ∗ = ( x − L w ) / ( Ca e × h ). Substituting all the experimentally measuredvariables in x ∗ and Ca e such as the core fluid flow rate Q = 6.5 mm s − , viscosity of the core fluid, η = 4500 mPa s,the wetted length L w ∼ . h along with the streamwise downstream positions x of the thread height, and the channel4 Table I:
Fitting parameters a , b and c of the exponential fit [Eq. (14), dashed-white line in Figs. 7(a) and 7(c)]representing the master curve. The estimated Γ e acting between the present experimental fluids is obtained byutilising the spatially measured experimental thread height ε z /h ( β = 90 ° ) (shown in Fig. 7(b)) and Eq. (14) alongwith the fit parameters, and other experimentally acquired quantities such as viscosity ( η ), flow rate ( Q ) of the corefluid and the channel width h . a b c Estimated Γ e (mN m − )0 . ± . . ± .
015 0 . ± . . ± . width h , the estimated Γ e (cid:39) .
765 mN m − 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 experimen-tal (Fig. 7(b)) thread height results are plotted in Fig. 7(d) as an error, which is defined as, δ ε z (Γ e ) = N (cid:88) i =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε numz,i − ε expz,i ε expz,i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (15)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 − which, inturn, affirms the estimated Γ e to be good with fairly accurate order of magnitude. This value of Γ e = γ is used in thenumerical 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∼ − mN m − ) for the present studycolloidal dispersion-solvent system is almost one decade higher in magnitude than our previously evaluated value( O∼ − mN m − ) 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 thepresent case is η = η (cid:39) e could be attributed to the variation in length fraction of nanofibrils, interfibril interactions leading to differentrheological properties of the colloidal dispersion, in turn to that on variant dispersion-solvent de-facto interfaceproperties. 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 formild 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 themodeling and prediction of hydrodynamic alignment of nanofibrils [58]. Alignment of nanofibrils is a key factor incontrolling 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 downstreampositions x/h for all the five confluence angle geometries. The agreement between numerical computations andexperimental measurements is excellent. Dashed vertical lines (magenta) in all the panels (Figs. 8(a) - 8(e)) mark theconfluence 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 < x/h = 0 that is followed by a rapid acceleration before a high steady value at fardownstream positions ( x/h ≥
10) is reached. However, a careful inspection of Fig. 8 and Fig. 9 unveils a considerablevariation in the minimum and maximum centreline velocity, degree of deceleration and acceleration for differentconfluence 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 ˙ ε ( d U /d x )along the centreline for all the geometrical configurations. In the case of the reference configuration ( β = 90 ° ), closeto the confluence point x/h (cid:39)
0, the deceleration is ascertained to be small as observed from Figs. 8(a) (confluenceregion) and 9(c) (cyan color). As a consequence, the minimum centreline velocity U c,min is higher (meaning lowerdeceleration) than the other configurations as seen from Fig. 9(a). Subsequent to deceleration, a substantial increasein 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 ∼ Figure 8:
Panels (a) - (e) show the experimental and numerical centreline velocity U c as a function of streamwisedownstream positions x/h for flow-focusing configurations with varying confluence angle β . The vertical dashed lines(magenta) on all the panels denote the confluence region of the respective configurations. Error bars represent thestandard deviation of experimental measurements at each position.thereafter, ˙ ε decays to 0 at around x/h ≈ U c,max is lowest as compared to all the other configurations.For the β = [60 ° , ° ] pair, the velocity behaviour duplicate each other as observed from Figs. 8(b) and 8(c) andmore clearly from the strain rate ˙ ε plot in Fig. 9(c). It is worth to recall that this similarity of the flow mirrors thesimilarity of the thread shape evolution as seen earlier from Figs. 3(b), 3(c) to Figs. 6(b), 6(c). In both cases, theeffect of β appends the streamwise flow through the contribution by sheath flow momemtum acting parallel to thecore flow. As a result, the core fluid velocity increases much faster than the reference configuration velocity. For bothcases, 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 islower for the β = 90 ° case (implying higher deceleration) and U c,max is highest for [60 ° , ° ] cases when comparedto all the other configruations. As observed from Fig. 9(c), in both configurations, the variation in ˙ ε (golden andgreen color) shows a similar trend as the β = 90 ° case (cyan color) but the magnitude is higher near the end of theconfluence regions ( x/h ∼ x/h ≈ Figure 9:
Numerical plots of variation of minimum (panel (a)) and maximum velocity (panel (b)), and strainrate ˙ ε (panel (c)) along the centreline as a function of confluence angle β . The vertical dashed lines (magenta) inpanel (c) denote the confluence regions for all the configurations. The confluence region starts at x/h = 0 for all theconfigurations.On the contrary, Figs. 8(d) and 8(e) display a quite different behaviour for the β = [30 ° , ° ] pair. Even though boththe configurations have the same confluence region 0 ≤ x/h ≤
2, there is a considerable variation in the decelerationand acceleration regions as seen from Fig. 9(c). For the β = 30 ° case, the deceleration ( x/h (cid:39)
0) is higher leadingto the lowest U c,min among all the configurations and also lower U c,max far downstream as seen in Figs. 9(a) and9(b), whereas in the β = 150 ° case, deceleration is slightly lower as compared to the β = 30 ° case and U c,max is alsohigher. 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 valueclose to the end of confluence region ( x/h ∼ x/h > ε in both cases follows a similartrend.Thus, having affirmed the agreement between experimental measurements and numerical computations, we willnow utilise purely the numerical computations for further analysis of the results in the sections that follow. In thenext section, we look at the cross-sectional velocity distribution for all the geometrical configurations with varying β . IV. Cross-sectional velocity distribution
Figures 10 and 11 show the numerical velocity maps superimposed with in-plane velocities ( V and W , 2-D vectors),together with the profiles of streamwise velocity component at y/h , z/h = 0 in different cross-sectional x/h planesfor all the geometrical configurations. The cross-sectional x/h planes in Fig. 10 corresponds to the respective threaddetachment position, i.e . the wetted length L w /h of each configuration ( x/h ∼ . β = 90 ° ; x/h ∼ .
15 for β = [60 ° , ° ]; x/h ∼ . , .
25 for β = [30 ° , ° ]) while in Fig. 11, x/h planes at far downstream ( x/h = 20) areshown.As seen in Fig. 10, velocity maps depict the core fluid thread being enclosed by high-velocity sheath flows. Further,the 2-D vector plot of in-plane velocities exhibit no evidence of appreciable recirculating/secondary flows at the corners7 -0.5 -0.25 0 0.25 0.5-0.5-0.2500.250.5 Figure 10: (a) - (e) Cross-sectional velocity maps showing streamwise velocity ( U ) superimposed with the in-planevelocity fields ( V and W , represented by black arrows) at different positions x/h of flow-focusing configurations withvarying β . The positions x/h correspond to the detachment locations of the core fluid thread from the top wallnamely the wetted lengths L w /h of respective configurations as per Figs. 4(a) - 4(e). Velocity profiles denoted bysolid and dashed lines represent the streamwise velocity component at z/h, y/h = 0 (along the dashed white lines onthe velocity maps). Note that the colorbar is different in each case (particularly for β = 150 ° ).8 ✲(cid:0)✁✂ ✲(cid:0)✁✄✂ (cid:0) (cid:0)✁✄✂ (cid:0)✁✂✲(cid:0)✁✂✲(cid:0)✁✄✂(cid:0)(cid:0)✁✄✂(cid:0)✁✂ ✵✺✶✵✶✺✷✵✷✺(cid:0)✂☎(cid:0)☎✂✄(cid:0)✄✂ (cid:0) ✂ ☎(cid:0) ☎✂ ✄(cid:0) ✄✂✲(cid:0)✁✂✲(cid:0)✁✄✂(cid:0)(cid:0)✁✄✂(cid:0)✁✂ ✵✺✶✵✶✺✷✵✷✺(cid:0)✂☎(cid:0)☎✂✄(cid:0)✄✂ ✲(cid:0)✁✂✲(cid:0)✁✄✂(cid:0)(cid:0)✁✄✂(cid:0)✁✂ ✵✺✶✵✶✺✷✵✷✺(cid:0)✂☎(cid:0)☎✂✄(cid:0)✄✂✲(cid:0)✁✂ ✲(cid:0)✁✄✂ (cid:0) (cid:0)✁✄✂ (cid:0)✁✂✲(cid:0)✁✂✲(cid:0)✁✄✂(cid:0)(cid:0)✁✄✂(cid:0)✁✂ ✵✺✶✵✶✺✷✵✷✺(cid:0)✂☎(cid:0)☎✂✄(cid:0)✄✂ (cid:0) ✂ ☎(cid:0) ☎✂ ✄(cid:0) ✄✂ ✲(cid:0)✁✂ ✲(cid:0)✁✄✂ (cid:0) (cid:0)✁✄✂ (cid:0)✁✂✲(cid:0)✁✂✲(cid:0)✁✄✂(cid:0)(cid:0)✁✄✂(cid:0)✁✂ ✵✺✶✵✶✺✷✵✷✺(cid:0)✂☎(cid:0)☎✂✄(cid:0)✄✂ (cid:0) ✂ ☎(cid:0) ☎✂ ✄(cid:0) ✄✂ Figure 11: (a) - (e) Cross-sectional velocity maps showing streamwise velocity ( U ) superimposed with in-planevelocity fields ( V and W , represented by black arrows) at far downstream position x/h = 20 for flow-focusing configu-rations with varying β . The in-plane velocities are very weak for all the configurations. Velocity profiles of streamwisecomponent are denoted by the solid and dashed lines at z/h, y/h = 0 (demarcation shown by dashed white line onthe velocity maps) show predominantly plug flow behaviour.9in the cross-sectional planes, albeit they are occupied by a region of high-velocities. However, the magnitudes of thesehigh-sheath velocities seem to differ significantly for different confluence angles β . Accordingly, in the case of negativeimpingement configurations i.e. β = 120 ° and 150 ° , the velocity magnitudes are relatively higher and so is the lengthof 2-D vectors as seen in Figs. 10(c) and 10(d). They are maximum for the β = 150 ° case. The velocity profiles arenearly flat along the z -direction at y/h = 0 in all configurations inside the core fluid thread. Also, along the y -directionat 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 threadstays 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 ° , ° ] pair amongst all the configurations, corroboratingwith earlier observations of maximum centreline velocity U c,max for the same cases in Fig. 9(c). Furthermore, the2-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 inboth 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. Twocomputations were performed utilising the reference flow-focusing configuration ( β = 90 ° ) with varying sheath flowinlet channel widths as illustrated in Fig. 12(a). Accordingly, the side channel width (SCW) of the ‘SCW-1 . h ’configuration was set to 1 . h matching the confluence region (0 ≤ x/h ≤ β = [60 ° , ° ] configurations.Similarly, the ‘SCW-2 h ’ configuration width was fixed to 2 h tallying with the confluence region (0 ≤ x/h ≤
2) of the β = [30 ° , ° ] configurations. Other parameters such as flow rates of the core ( Q ) and sheath ( Q ) flows, rheologiesof fluids are held fixed as described in the Secs. II A 1 and II A 2. The interfacial tension γ = Γ e = 0 .
615 mN m − ,and the cross-sectional width h of the central inlet and outlet channel were unchanged. This means that the sheathflow channel inlets are now rectangular in these cases.Astonishingly, the flow in the SCW-1 . h configuration precisely replicates the morphological features of the wettedregion and thread topologies of the β = [60 ° , ° ] configurations as observed from Figs. 12(b) and 12(d). Once againthis highlights the symmetry between these two cases. Whereas, the SCW-2 h configuration captures and replicatesthe features of only the β =150 ° configuration as depicted in Figs. 12(c) and 12(e) (remember the asymmetry betweenthe β = [30 ° , ° ] configurations).In the three cases, ‘SCW-1 . h ’ and β = [60 ° , ° ], the sheath flow momentum normal to the core flow is atten-uated by 13.4% in relative to the reference configuration ( β = 90 ° ). For the ‘SCW-2 h ’ and β = [30 ° , ° ] pair, theattenuation is around 50%. Therefore, inclining the sheath flow channels at various confluence angles β or settingthe width of sheath flow channels of reference configuration to confluence region lengths, would essentially meaninducing 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 incontrolling 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 Atogether with the alternative approach in previous Sec. V underscores the significant role of Γ e . As noticed, for a setflow rate and specified rheologies of fluids, the same value of EIT ( γ = Γ e = 0 .
615 mN m − ) deduced using the referenceflow-focusing configuration ( β = 90 ° ), very well reproduces the 3-D flow characteristics for the respective geometricalconfigurations 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, thereplication 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 andbottom channel walls. We find that, the above value of interfacial tension γ in the numerical computations not onlycaptured the wetted region shapes and wetted lengths L w /h accurately, but also steered the thread detachment fromthe channel walls. Irrespective of the geometrical configurations, we also observed through a separate set of simulationsthat are not displayed here, that the wetted length L w /h → ∞ if γ →
0, implying no detachment of thread from thetop and bottom channel walls. Therefore, presence of ultra-low interfacial tension γ (Γ e ) is a decisive factor in orderfor the thread detachment to occur. Nevertheless, given the presence of high viscosity contrast ( χ (cid:39) Figure 12: (a) 3-D view of the thread shapes for flow-focusing configurations with varying sheath flow channelwidths (denoted by dashed magenta lines). Width of the side channel in the ‘SCW-1 . h ’ configuration correspondsto the confluence region (0 ≤ x/h ≤ β = [60 ° , ° ] pair (see Figs. 3(b), 3(c)). Likewise, width of theside channel in ‘SCW-2 h ’ configuration corresponds to the confluence region (0 ≤ x/h ≤
2) of β = [30 ° , ° ] pair(see Figs. 3(d), 3(e)). Panels (b) and (d) shows the ‘SCW-1 . h ’ configuration wetted region, thread width ( ε y /h )and height ( ε z /h ) plotted as a function of downstream positions x/h overlapped with the wetted region and threaddimensions of β = [60 ° , ° ] pair. Similar plots of panels (b) and (d) in panels (c) and (e) for the case of ‘SCW-2 h ’and β = 150 ° configurations.phenomena of self-lubrication [33] depending on the fluid injection geometries and injection flow rates. According tothis 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´eclet numbermicrofluidic flows involving miscible systems with high-viscosity contrast, thread detachment is found to be clearlycontrolled by Γ e as well. Indeed, in order to have a finer understanding of the effects of EIT versus the ones associatedwith self-lubrication on the thread detachment, further detailed investigation is needed. For instance, variation ofEIT, viscosity contrast and flow rate ratios, and how these parameters affect the wetted length L w /h and thereby thethread detachment is worthwhile, and such investigations are beyond the scope of this present work.1 VII. Conclusions
The effective interfacial tension Γ e acting between the colloidal dispersion and its own solvent has been estimatedutilising the experimentally measured 3-D spatial evolution of thread shape in the flow-focusing channel and aneffective capillary number dependent master curve obtained from a simple scaling model [21]. We have verified fromthe 3-D numerical computations, that the Γ e which minimizes the error between computed and experimental threadheights is close to the estimated Γ e utilising the above method. The numerical computations were able to captureour experimental observations with good quantitative and qualitative agreement (both in terms of flow-topologiesand flow-fields) for a range of confluence angles. By mapping the numerical observations with the experimentalmeasurements, we gained insights into the influence of various physical mechanisms driving the process of threaddetachment from the top and bottom walls in the channels. From these meticulous analyses, we find that in thecase of high viscosity contrast miscible systems [25, 26], the thread detachment occuring in the physical experimentsof microchannels may not be entirely based on the phenomena of self-lubrication [23] related with effect of viscousdissipation 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 inexperimental measurements. The standard experimental techniques yielding reasonable results in the case of molecularmiscible fluids are light scattering and spinning drop tensiometry [20]. However, these two techniques are observedto have drawbacks, in particular at the transition region due to mixing of the two fluids and difficulty to distinguishbetween the bulk fluids and interfacial region [20]. The recently explored new measurement strategies like examiningthe Saffman–Taylor instability in a Hele-Shaw cell [7, 8] with miscible complex fluids, and studying the dynamics ofdrop 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 togenerate effective extensional flows. This, in turn, could facilitate to achieve optimal alignment of nanofibrils incolloidal 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 dif-fusion, mobility, and rigidity of fibrils. Concluding, a thorough understanding of effects like EIT in miscible systemsopens up possibilities in tuning, and controlling the material properties more efficiently via appropriate selection ofmicrofluidic geometrical configuration based on the practical interest and application. Acknowledgments
Financial support by the Swedish Research Council for Environment, Agricultural Sciences and Spatial Plan-ning (FORMAS) and Swedish Research Council (VR) are gratefully acknowledged. The simulations were performedon computer resources provided by the Swedish National Infrastructure for Computing (SNIC) at National Super-computer Centre, Link¨oping and at the Centre for High Performance Computing (PDC), KTH Royal Institute ofTechnology, Stockholm. Dr. Korneliya Gordeyeva is acknowledged for experimental assistance. The authors thankDr. Christophe Brouzet for helpful discussions. [1] D. D. Joseph, Fluid dynamics of two miscible liquids with diffusion and gradient stresses, Eur. J. Mech. B Fluids , 565(1990).[2] P. Garik, J. Hetrick, B. Orr, D. Barkey, and E. Ben-Jacob, Interfacial cellular mixing and a conjecture on global depositmorphology, Phys. Rev. Lett. , 1606 (1991).[3] D. D. Joseph and Y. Y. Renardy, Fundamentals of Two-Fluid Dynamics. Part II: Lubricated Transport, Drops and MiscibleLiquids (Springer, New York, 1993).[4] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annual Review ofFluid Mechanics , 139 (1998).[5] J. Atencia and D. J. Beebe, Controlled microfluidic interfaces, Nature , 648 (2005).[6] D. J. Korteweg, Sur la forme que prennent les ´equations du mouvement des fluides si l’on tient compte des forces capillairescaus´ees par des variations de densit´e, Arch. N´eerlandaises Sci. Exactes Naturelles (1901).[7] D. Truzzolillo, S. Mora, C. Dupas, and L. Cipelletti, Off-equilibrium surface tension in colloidal suspensions, Phys. Rev.Lett. , 128303 (2014).[8] D. Truzzolillo, S. Mora, C. Dupas, and L. Cipelletti, Nonequilibrium interfacial tension in simple and complex fluids, Phys.Rev. X , 041057 (2016). [9] N. Bessonov, V. Volpert, J. A. Pojman, and B. Zoltowski, Numerical simulations of convection induced by Kortewegstresses in miscible polymer-monomer systems, Microgravity-Science and Technology , 8 (2005).[10] C. Chen, L. Wang, and E. Meiburg, Miscible droplets in a porous medium and the effects of Korteweg stresses, Phys.Fluids , 2447 (2001).[11] C. Chen and E. Meiburg, Miscible displacements in capillary tubes: Influence of Korteweg stresses and divergence effects,Phys. Fluids , 2052 (2002).[12] C. Chen, C. Huang, H. Gadˆelha, and J. Miranda, Radial viscous fingering in miscible Hele-Shaw flows: A numerical study,Phys. Rev. E , 016306 (2008).[13] E. O. Dias and J. A. Miranda, Wavelength selection in Hele-Shaw flows: A maximum-amplitude criterion, Phys. Rev. E , 013016 (2013).[14] E. Mossige, C. S.V., M. Islamov, S. Wheeler, and G. G. Fuller, Evaporation-induced Rayleigh–Taylor instabilities inpolymer solutions, Philosophical Transactions of the Royal Society A , 20190533 (2020).[15] J. A. Pojmanm, C. Whitmore, L. Turco, L. Maria, R. Lombardo, J. Marszalek, R. Parker, and B. Zoltowski, Evidence forthe existence of an effective interfacial tension between miscible fluids: Isobutyric acid- water and 1-butanol- water in aspinning-drop tensiometer, Langmuir , 2569 (2006).[16] A. Carbonaro, L. Cipelletti, and D. Truzzolillo, Ultralow effective interfacial tension between miscible molecular fluids,Phys. Rev. Fluids , 074001 (2020).[17] V. Shevtsova, Y. Gaponenko, V. Yasnou, A. Mialdun, and A. Nepomnyashchy, Two-scale wave patterns on a periodicallyexcited miscible liquid–liquid interface, J. Fluid. Mech. , 409 (2016).[18] S. May and J. Maher, Capillary-wave relaxation for a meniscus between miscible liquids, Phys. Rev. Lett , 2013 (1991).[19] P. Cicuta, A. Vailati, and M. Giglio, Capillary-to-bulk crossover of nonequilibrium fluctuations in the free diffusion of anear-critical binary liquid mixture, Applied Optics , 4140 (2001).[20] D. Truzzolillo and L. Cipelletti, Off-equilibrium surface tension in miscible fluids, Soft Matter , 13 (2017).[21] K. Gowda, C. Brouzet, T. Lefranc, L. S¨oderberg, and F. Lundell, Effective interfacial tension in flow-focusing of colloidaldispersions: 3-D numerical simulations and experiments, J. Fluid. Mech. , 1052 (2019).[22] T. Cubaud and T. Mason, Folding of viscous threads in diverging microchannels, Phys. Rev. Lett. , 114501 (2006).[23] T. Cubaud and T. Mason, Interacting viscous instabilities in microfluidic systems, Soft Matter , 10573 (2012).[24] O. Bonhomme, J. Leng, and A. Colin, Microfluidic wet-spinning of alginate microfibers: A theoretical analysis of fiberformation, Soft Matter , 10641 (2012).[25] T. Cubaud and S. Notaro, Regimes of miscible fluid thread formation in microfluidic focusing sections, Phys. Fluids ,122005 (2014).[26] T. Cubaud and T. G. Mason, High-viscosity fluid threads in weakly diffusive microfluidic systems, New J. Phys. , 075029(2009).[27] T. Cubaud, Swelling of diffusive fluid threads in microchannels, Phys. Rev. Lett. , 174502 (2020).[28] S. L. Anna, N. Bontoux, and H. A. Stone, Formation of dispersions using flow focusing in microchannels, App. Phys. Lett. , 364 (2003).[29] P. Garstecki, I. Gitlin, W. DiLuzio, G. M. Whitesides, E. Kumacheva, and H. A. Stone, Formation of monodisperse bubblesin a microfluidic flow-focusing device, App. Phys. Lett. , 2649 (2004).[30] G. M. Whitesides, The origins and the future of microfluidics, Nature , 368 (2006).[31] C. W. Hirt and B. D. Nichols, Volume of fluid (VoF) method for the dynamics of free boundaries, J. Comput. Phys. ,201 (1981).[32] E. OpenCFD, OpenFOAM v1806 Userguide (ESI-OpenCFD, 2018).[33] D. Joseph, K. Nguyen, and G. Beavers, Non-uniqueness and stability of the configuration of flow of immiscible fluids withdifferent viscosities, J. Fluid. Mech. , 319 (1984).[34] K. M. O. H˚akansson, A. B. Fall, F. Lundell, S. Yu, C. Krywka, S. V. Roth, G. Santoro, M. Kvick, L. Prahl-Wittberg,L. W˚agberg, and L. D. S¨oderberg, Hydrodynamic alignment and assembly of nanofibrils resulting in strong cellulosefilaments, Nat. Commun. , 4018 (2014).[35] N. Mittal, F. Ansari, G. V. K., C. Brouzet, P. Chen, P. T. Larsson, S. V. Roth, F. Lundell, L. Wagberg, N. A. Kotov,and L. D. S¨oderberg, Multiscale control of nanocellulose assembly: Transferring remarkable nanoscale fibril mechanics tomacroscale fibers, ACS Nano , 6378 (2018).[36] C. Brouzet, N. Mittal, L. D. S¨oderberg, and F. Lundell, Size-dependent orientational dynamics of Brownian nanorods,ACS Macro Letters , 1022 (2018).[37] C. Brouzet, N. Mittal, , F. Lundell, and L. D. S¨oderberg, Characterizing the orientational and network dynamics ofpolydisperse nanofibers on the nanoscale, Macromolecules , 2286 (2019).[38] D. Truzzolillo, V. Roger, C. Dupas, S. Mora, and L. Cipelletti, Bulk and interfacial stresses in suspensions of soft and hardcolloids, Journal of Physics: Condensed Matter , 194103 (2015).[39] T. Saito, S. Kimura, Y. Nishiyama, and A. Isogai, Cellulose nanofibers prepared by Tempo-mediated oxidation of nativecellulose, Biomacromolecules , 2485 (2007).[40] A. Isogai, T. Saito, and H. Fukuzumi, Tempo-oxidized cellulose nanofibers, Nanoscale , 71 (2011).[41] F. Marto¨ıa, P. J. J. Dumont, L. Org´eas, M. N. Belgacem, and J. L. Putaux, Micro-mechanics of electrostatically stabilizedsuspensions of cellulose nanofibrils under steady state shear flow, Soft matter , 1721 (2016).[42] O. Nechyporchuk, M. N. Belgacem, and F. Pignon, Current progress in rheology of cellulose nanofibril suspensions,Biomacromolecules , 2311 (2016).[43] L. Geng, N. Mittal, C. Zhan, F. Ansari, P. R. Sharma, X. Peng, B. S. Hsiao, and L. S¨oderberg, Understanding the mechanistic behavior of highly charged cellulose nanofibers in aqueous systems, Macromolecules , 1498 (2018).[44] L. H. Switzer and D. J. Klingenberg, Simulations of fiber floc dispersion in linear flow fields, Nordic Pulp & Paper ResearchJournal , 141 (2003).[45] M. Doi and S. F. Edwards, The theory of polymer dynamics , Vol. 73 (Oxford University Press, 1988).[46] G. Marrucci and N. Grizzuti, The effect of polydispersity on rotational diffusivity and shear viscosity of rodlike polymersin concentrated solutions, Journal of Polymer Science: Polymer Letters Edition , 83 (1983).[47] G. Marrucci and N. Grizzuti, Predicted effect of polydispersity on rodlike polymer behaviour in concentrated solutions,Journal of Non-Newtonian Fluid Mechanics , 103 (1984).[48] A. W. Chow, G. G. Fuller, D. G. Wallace, and J. A. Madri, Rheo-optical response of rodlike, shortened collagen proteinto transient shear flow, Macromolecules , 805 (1985).[49] R. Pecora, Dynamics of rodlike macromolecules in semidilute solutions, Journal of Polymer Science: Polymer Symposia , 83 (1985).[50] S. S. Rogers, P. Venema, L. M. C. Sagis, E. van der Linden, and A. M. Donald, Measuring the length distribution of afibril system: a flow birefringence technique applied to amyloid fibrils, Macromolecules , 2948 (2005).[51] T. Ros´en, B. S. Hsiao, and L. D. S¨oderberg, Elucidating the opportunities and challenges for nanocellulose spinning,Advanced Materials , 2001238 (2020).[52] D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, and C. A. Puliafito,Optical coherence tomography, Science , 1178 (1991).[53] W. Drexler and J. G. Fujimoto, Optical coherence tomography: technology and applications (Springer Science & BusinessMedia, 2008).[54] R. A. Leitgeb, R. M. Werkmeister, C. Blatter, and L. Schmetterer, Doppler optical coherence tomography, Progress inRetinal and Eye Research , 26 (2014).[55] J. Roenby, H. Bredmose, and H. Jasak, A computational method for sharp interface advection, R. Soc. Open Sci. , 160405(2016).[56] J. U. Brackbill, D. B. Kothe, and C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. ,335 (1992).[57] T. Cubaud and T. G. Mason, Capillary threads and viscous droplets in square microchannels, Phys. Fluids , 053302(2008).[58] G. B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proceedings of the Royal Society of London.Series A, Containing papers of a mathematical and physical character , 161 (1922).[59] K. M. O. H˚akansson, F. Lundell, L. Prahl-Wittberg, and L. D. S¨oderberg, Nanofibril alignment in flow focusing: Measure-ments and calculations, J. Phys. Chem. B120