# Triggering the Splitting Dynamics of Low-Viscous Fingers through Surface Wettability Inside Bifurcating Channel

Authors: Akhileshwar Singh; Krishna Murari Pandey; Yogesh Singh
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This work presents the splitting dynamics of low-viscous fingers inside the single bifurcating channel through the surface wettability of daughter branches. The propagation of low-viscous fingers inside branching microchannels have importance in many applications, such as microfluidics, biofluid mechanics (pulmonary airway reopening), and biochemical testing. Several numerical simulations are performed where a water finger propagates inside the silicon oil-filled bifurcating channel, and at the bifurcating tip, it splits into two fingers and these fingers propagate into the separate daughter branches. It is noticed that the behaviour of finger splitting at the bifurcating tip depends upon numerous parameters such as surface wettability, capillary number, viscosity ratio, and surface tension. This study aims to trigger the behaviour of finger splitting through the surface wettability of daughter branchesθ1,θ2. Therefore, a series of numerical simulations are performed by considering four different surface wettability configurations of daughter branches, i.e., θ1,θ2∈78°,78°;78°,118°;78°,150°;150°,150°. According to the results obtained from numerical simulations, finger splitting may be categorized into three types based on splitting ratio λ, i.e., symmetrical splitting, nonsymmetrical splitting, and reversal (no) splitting. It is noticed that the surface wettability of both daughter branches is either hydrophilic 78°,78° or superhydrophobic 150°,150°, providing symmetrical splitting. The surface wettability of one of the daughter branches is hydrophilic and another is hydrophobic 78°,118°, providing nonsymmetrical splitting. The surface wettability of one of the daughter branches is hydrophilic and another is superhydrophobic 78°,150°, providing reversal splitting. The findings of this investigation may be incorporated in the fields of biochemical testing and occulted pulmonary airways reopening as well as respiratory diseases such as COVID-19.


---

## Body

## 1. Introduction

Viscous fingering (VF) is a hydrodynamic interfacial instability that is inherently evolved at the interface of two fluids when a low-viscous fluid displaces another more viscous fluid in the porous medium or Hele-Shaw cell [1–3]. These two fluids may be miscible or immiscible. The low-viscous fluid may be invaded radially or linearly in the results; radial or rectilinear linear displacement of more viscous fluid is observed. According to these displacements, VF is categorized into two groups, radial VF [4–7] and rectilinear VF[8–11]. The displacement of more viscous fluid by less viscous fluid leads to be unstable due to viscosity differences (viscous force) and density differences (gravity force) so that finger-shaped instability appears at the interface of two fluids. It is so-called viscous fingering. Viscous fingering is seen in many natural and industrial processes such as oil recovery [12], mixing of fluids in microfluidic devices [13], drugs delivery [14], pulmonary biomechanics [15], and complicated branching channels of biochemical testing [16]. In two-phase displacement flows, where fluid-1 and fluid-2 are represented displaced (defending) fluid and displacing (invading) fluid, respectively. The physical properties of fluid-1 and fluid-2, i.e., densities and viscosities, are represented by ρ1,μ1 and ρ2,μ2, respectively. It is assumed that the viscosity of fluid-1 is more than fluid-2 μ1>μ2. Saffman and Taylor [17] conducted first experiments in 1958 on viscous fingering instability. They performed experiments on Hele-Shaw cell (HSC) by considering the less viscous fluid-2 displaces more viscous fluid-1. The fluid-2 was invaded linearly with constant velocity in HSC, after that fluid-2 was penetrated in the fluid-1 in the form of finger shape. Again, Taylor [18] performed experiments on the fluid-fluid displacement in the microtube. It was noticed that more viscous fluid was displaced by less viscosity fluid from the middle portion of the microtube but near the wall was not displaced by it. Therefore, the displacement front looked finger-shaped. Reinelt et al. [19] numerically investigated linear penetration of low-viscous fingers into the more viscous fluid inside the channel and tube. It has been found that viscous fingering instabilities are driven by three forces, i.e., viscous force (driven by viscosity differences), gravity force (driven by density differences), and surface tension force. But gravity force may be neglected in the horizontal channel and surface tension force may also be neglected in miscible fluids. It is observed that viscous fingering instability appears on the fluids interface inherently when the viscosity ratio μr=μ2/μ1 is less than unity and it is called an unfavourable or an adverse viscosity ratioμr<1.VF is again classified into two groups according to the nature of fluids, i.e., miscible VF [20–22] and immiscible VF [23–25]. In the miscible VF, fluid pair is miscible to each other and the effect of surface tension is neglected at the interface of two fluids because surface tension is considered to be zero in this VF. In immiscible VF, fluid pair is immiscible to each other, and pressure jump at the interface of two-fluid interface occurs due to surface tension. Sahu et al. [20] examined the linear stability analysis of miscible two-phase flow inside the horizontal rectangular channel. They used both the direct numerical simulation method and linear stability theory for investigation. The convective-diffusion equation was coupled with transport equations because the viscosity was a function of the concentration gradient. They performed numerical simulations at high Reynolds numbers and they assumed that the density of fluids is constant (neutrally buoyance fluids). The operating parameters of investigations were Reynolds number, viscosity ratio, and Schmidt number. They explored the roll-up phenomena as well as intense convective mixing. They found these phenomena are very noticeable at high flow rates and a large viscosity ratio. Sahu et al. [21] investigated the effect of buoyancy on the miscible displacement flow inside an inclined channel. They adopted the direct numerical simulations method for the study. They assumed that both density and viscosity are dependent upon the concentration gradients. The miscible displacement flows were driven by viscous force (due to viscosity gradient) and gravity force (due to density gradient). Numerical experiments were conducted at moderate Reynolds numbers and viscosity ratio. This investigation was based on these parameters-density ratio, channel inclination, and Froude number. They found that the rate of mixing and displacement may be enhanced by increasing the density ratio and Froude number as well as inclination angle. Mishra et al. [22] investigated the impact of double-diffusion on miscible displacement flow in a horizontal channel. If a liquid solution has two solutes, then the viscosity of the liquid solution depends upon the diffusing rate of solutes. They considered that viscosity is an exponential function of the concentration gradient of solutes. They also assumed that density of both fluids is equal. They reported different types of instability patterns formation on the fluid interface when two solutes diffuse at various rates. They found instability increases with increasing the diffusivity ratio between the faster-diffusing and the slower-diffusing solutes. They also studied the effect of Reynolds number and Schmidt number on instability patterns formation. Several experimental, numerical, and theoretical studies have been performed on controlling viscous fingering. Because viscous fingering is beneficial in some engineering applications (mixing of two fluids in microchannel) and detrimental in some engineering applications (enhance oil recovery). Controlling (enhancing or damping) the viscous fingering is studied into two parts: active control and passive control. In the active control [23–26], external surface forces such as electric field force and magnetic field force are imposed against hydrodynamic interfacial instabilities to control the viscous fingering. Flament et al. [23] investigated the effect of the magnetic force on radial viscous fingering. They observed that hydrodynamic instability between magnetic fluid and air may be suppressed by external magnetic force. Sinzato and Cunha [24] imposed an external magnetic field against the rectilinear viscous fingering. They found that finger-shaped instability may be damped by magnetic force. Mirzadeh and Bazant [25] investigated theoretically the effect of external electric fields on linear viscous fingering. They used two immiscible electrolytes and performed a theoretical study. They found that viscous fingering may be hindered by the electric field. They also found that after the critical injection ratio, interfacial stability is a function of the relative direction of flow and current, regardless of the viscosity ratio. Gao et al. [26] performed experiments and found that radial viscous fingering may be suppressed or enhanced by the external electric force. They also found that radial viscous fingering may be generated by an electric field when the viscosity ratio is favourable μr>1.In the passive control, the external surface forces are absent and it is achieved through varying the operating parameters such as injection rate [27, 28], surface tension [29, 30], viscosity ratio [31, 32], surface wettability [33, 34], and geometry [35–40]. Kozlov et al. [27] performed experiments on controlling immiscible radial viscous fingering through flow rate. They were injected low-viscous fluid radially with harmonic flow rate. They found that stable radial displacement may be achieved only at small amplitudes of oscillations. Li et al. [28] reported damping of radial viscous fingering through the injection rate of the less viscous fluid. They performed a numerical and experimental investigation and found that the morphology of viscous fingering may be precisely controlled by injection rate. McLean and Saffman [29] reported the effect of surface tension on the shape of immiscible rectilinear viscous fingers. They included surface tension in numerical modeling. After the linear stability analysis, they found that surface tension does not stabilize two-dimensional fingers. Saffman [30] reported a theoretical study on the shape and stability of immiscible fingers in HSC cells. They explained various possible shapes of immiscible rectilinear fingers. They included surface tension in this theoretical study. Bischofberger et al. [31] reported morphology of radial miscible viscous fingering at various viscosity ratios. They performed experiments and found that miscible VF may be damped by increasing the viscosity ratio. They developed empirical relation between viscosity ratio and finger lengths. The surface tension was considered to be zero. Varges et al. [32] performed experiments on the rectilinear immiscible displacement. They found that immiscible displacement strongly depends upon the viscosity ratio. They observed that unstable displacement is sifted into stable displacement when the viscosity ratio is approaching unit value. It means that stable displacement (no fingering) occurs inherently when the viscosity ratio is more than one. They also found that a more stable interface appears when the Newtonian liquids displace the shear-thinning liquids. Jackson et al. [33] analyzed the effect of wettability on immiscible radial VF. They noticed that wettability may be affected by interfacial pressure drop dynamically. They found that wettability may delay finger splitting at high capillary numbers and finger shielding may be also reduced by it. Hu et al. [41] reported the effect of wettability on immiscible fluid-fluid displacement in porous media. They used a theoretical model for study. They described briefly the impact of contact angle and capillary number on stable and unstable displacements. Primkulov et al. [34] examined the effect of wettability and capillary number on the dynamics of immiscible displacement in porous media. They used the moving-capacitor dynamic network model which captured the motion of each fluid interface by pore geometry. They analyzed viscous-dominated and capillary-dominated flow domains for a wide range of capillary numbers. Controlling the viscous fingering through HSC geometry has been extensively studied. Al-Housseiny et al. [35] performed experiments on hindering the immiscible rectilinear VF by tapering the geometry. They provided decreasing constant depth gradient along the length of HSC (converging HSC). They observed that the interface of oil and air becomes stable (no finger) when HSC was converged along the length. They concluded that stable immiscible displacement may be achieved by tapering the HSC. Again, Al-Housseiny et al. [36] reported a theoretical study on controlling VF in radial and rectilinear displacements by tapering the HSC in the direction of displacements. They collectively investigated the effect of depth gradient and capillary number. Anjos et al. [37] examined the role of capillary numbers on fingering instability in tapered radial HSC. They observed that stable displacement may be achieved at low capillary numbers in converging HSC, and at high capillary numbers, it may be achieved in diverging HSC. The results are altered at high capillary numbers.Bongrand and Tsai [38] tuned the radial VF in a nonuniform channel by injection rate. They experimentally investigated the collective effect of gap gradient and injection rate on VF. They found that radially converging cells may damp VF viscous fingering. Morrow et al. [39] studied suppressing the growth of the fingering instability in nonstandard HSC. They injected inviscid fluid with time-dependent injection rates in various configurations of HSC. They recorded various morphologies of fingering instability. Lu et al. [40] reported theoretical and numerical studies on immiscible fingering instability in tapered HSC. They adopted two types of tapered HSC, i.e., diverging cells with a positive depth gradient and converging cells with a negative depth gradient. They discussed various regions of displacement accordingly critical capillary number.Pihler-Puzovic et al. [42] performed experiments on suppressing the radial fingering instability in the narrow gap between two parallel plates. They were replaced one of the rigid plates with an elastic membrane and performed the experiments. They observed that fingering instability may be suppressed by replacing one of the rigid plates with elastic membrane. Again, Pihler-Puzovic et al. [43] reported radial VF in HSC where the top boundary was replaced with a thin elastic sheet. They observed that elastic sheet provides stability. Redapangu et al. [44] reported the displacement flow of two immiscible fluids inside an inclined channel. They used the lattice Boltzmann method for investigation. Here, displacement flow was driven by both viscous force (due to viscosity differences) and gravity force (due to density differences). They analyzed the effect of viscosity ratio, Froude number, capillary number, Atwood number, and inclination angle on the interfacial instability patterns formation. They found that by increasing the viscosity ratio or decreasing the inclination angle, the displacement rate decreases. They also observed lock-exchange flow. Kang et al. [45] numerically investigated the effects of viscosity ratio, capillary number, and wettability on immiscible fluid displacement in a 2D rectangular channel. They used the lattice Boltzmann model for interface tracking. They found that slip distance, as well as the width of a finger, decreases, but finger length increases with increasing the capillary number or viscosity ratio. They also found that slip distance decreases, but finger length increases with decreasing wettability of the displacing fluid.It is found that from the above literature, immiscible fingers propagation and splitting dynamics depend upon the numerous parameters such as surface tension (σ), surface wettability θ, viscosity ratio μr, and capillary number Ca. The splitting of the fingers is more fascinating in pulmonary biomechanics because occulted pulmonary airways are reopened by air fingers propagation into pulmonary airways. Numerous respiratory diseases such as respiratory distress syndrome as well as COVID-19 occur due to blockage of the pulmonary airways. The blockage pulmonary airways may be reopened through low-viscous finger propagation. Gaver et al. [46] experimented first time on the reopening of the occulted airway through the air finger. The limitation of this work was that they investigated the reopening of the occulted pulmonary airway by considering the single-channel pulmonary airway. But pulmonary airways exist inherently in the form of branching networks. Later on, several authors investigated the reopening of occulted bifurcation or branching type pulmonary airways through the low-viscous finger propagation [47]. The low-viscous finger may be also relevant for starving tumour cells in cancer treatment methods such as gas embolotherapy [48, 49].Numerous experimental and numerical research works [50–53] have been done on the splitting dynamics of the droplet into the bifurcating channels. Recently, Dekha et al. have performed computation work on the splitting dynamics of the droplet in a Y-shaped channel. They tuned the splitting dynamic of the droplet through surface wettability. They used the volume of fluid model (VOF) for the tacking of the interface. VOF model follows the conservation of mass and it is used to capture the interface of both continuous-continuous phases and continuous-discrete phases [54, 55]. VOF model provides the sharp and accurate interface of two fluids. But, the limitation of this model is that both fluids must be immiscible. Very few research works have been done on the low-viscous fingers propagation in the bifurcating channels or displacement flows in the branching channels. To our best of knowledge, we first time tune the spitting dynamics of the low-viscous finger inside a single bifurcating channel through the surface wettability of daughter branches. The motivation of this work has come from Baroud et al.’s [56] work.In this paper, splitting behaviour of a water finger is investigated inside the two-dimensional single bifurcating (Y-shaped) channel. Here, fluid-fluid displacement is driven by shear or pressure-driven flow in the absence of gravity force and inertial force. The entire Y-shaped channel is saturated with silicon oil after that water is invaded with a constant velocity V into the Y-shaped saturated channel. The parent channel walls are considered perfectly nonwetting in each simulation. But daughter branches are considered a combination of different wall surface wettability configurations such as hydrophilic, hydrophobic, and superhydrophobic. A low-viscous fluid displaces another much more viscous fluid in the microchannel, then the interface of fluids moves in the displacement direction in the form of finger shape due to viscosity differences. We report the propagation of water fingers in the silicone oil filled Y-shaped channel and splitting of a water finger at bifurcating tip is triggered by the surface wettability of daughter branches.


## 2. Numerical Method

### 2.1. Problem Formulation

In this work, we investigate the effect of surface wettability on the splitting dynamics of a low-viscous finger in the single bifurcating channel. A schematic diagram of a two-dimensional (2D) single bifurcating channel with inlet and outlet boundaries conditions is shown in Figure1(a). The parent channel is bifurcated into two identical daughter branches which are shown in Figure 1(a) as branch-1 and branch-2. The size of each daughter branch Ld1×Wd1=Ld2×Wd2 is equal and it is taken as 5mm×0.75mm. The length ratio Lr=Ld/Lp and width ratio Wr=Wd/Wp of the bifurcating channel are taken as 1 and 0.75, respectively. The bifurcation angle (β) is taken as 39°. The snapshot of meshed computational domain with proper dimensions is shown in Figure 1(b). The meshing of the computational domain has been done by the meshing module of ANSYS FLUENT software. The fine structural meshing has been adopted for precise results. Initially, the entire Y-shaped channel is filled with high viscous fluid-1 (red) after that fluid-2 is invaded through the inlet with the velocity, V. The orientation of the Y-shaped channel is considered as horizontal. Therefore, the effect of gravity force due to density difference may be neglected in this investigation. The nature of both fluids is considered immiscible and incompressible. The fluid-2 is injected into the bifurcating channel with a constant velocity of 0.015 m/s. The pressures at the outlets of the daughter branches are considered as atmospheric pressure (P = 0 Pa). The wall surface wettability is governed by the wetting contact angle and it is defined by Young’s equation, i.e., cosθ=σw1−σw2/σ, where θ denotes the wall wetting contact angle and σw1 denotes interfacial tension between wall and fluid-1, σw2 denotes interfacial tension between wall and fluid-2, and σ denotes interfacial tension between fluid-1 and fluid-2. After several investigations, four different wall surface wettability configurations of daughter branches are taken, keeping the parent channel as strongly hydrophobic or nonwetting θp=180°. These surface wettability configurations are created by incorporating hydrophilic, hydrophobic, and superhydrophobic wettability conditions. All surface wettability configurations are shown in Figure 2. In these figures, θp,θ1, and θ2 are representing wall wetting contact angle of parent channel, daughter branch-1, and daughter branch-2, respectively. The surface wettability configurations of daughter branches are as follows:(i)
The surfaces of both daughter branches are taken as hydrophilicθ1,θ2=78°,78°(ii)
The surfaces of branch-1 and branch-2 are taken as hydrophilic and hydrophobicθ1,θ2=78°,118°, respectively(iii)
The surfaces of branch-1 and branch-2 are taken as hydrophilic and superhydrophobicθ1,θ2=78°,150°, respectively(iv)
The surfaces of both daughter branches are taken as superhydrophobicθ1,θ2=150°,150°Figure 1
(a) Schematic diagram of the two-dimensional single bifurcating channel. (b) Snapshot of the meshed computational domain.
(a)(b)Figure 2
Different wall surface wettability configurations of daughter branches keeping parent channel as the nonwettingθp=180°.
(a)(b)(c)(d)


### 2.2. Mathematical Model

ANSYS FLUENT commercial CFD (computational fluid dynamic) software has been used for numerical simulations. The numerical simulations have been performed by considering a pressure-based solver with unsteady fluid displacement in the two-dimensional Y-shaped channel. The shape of the moving interface is tracked by using the volume of fluid (VOF) method. The VOF method is one of the numerical tools of interface tracking which is used to track the interface of two or more immiscible fluids in the multiphase flows [57]. It is more renowned among interface tracking methods because it follows the conservation of mass. In the VOF method, the transport equations are solved simultaneously for each phase with a single set of continuity and momentum equations in the entire flow field. In the present two-phase displacement flow, the volume fraction of primary and secondary phases is represented by α1 and α2, respectively. The sum of the volume fractions of two phases must be unity α1+α2=1 in each cell as well as the entire computational domain. The volume fractions are a function of space and time. The volume fraction equation is solved only for the second phase in the VOF model and the result; we get the value of α2. The value of α1 is calculated from the conservation equation of volume fraction, i.e., α1=1−α2. The moving interface is constructed in the computational domain on the value of volume fractions. The fluids inside the computational domain are defined by values of α2 in the following manner: α2=0⇒fluid1;1⇒fluid2;0<α2<1⇒interface.We assume that unsteady laminar flow occurs in the entire investigation and primary and secondary phases, i.e., fluid-1 (silicon oil) and fluid-2 (water) are immiscible and incompressible.The governing equations of immiscible and incompressible fluid flows can be expressed by the following mathematical expressions:Continuity equation:(1)∇⋅V=0.Momentum conservation equation (equation (2)):A single set of momentum equations is solved in entire the computational domain, and the obtaining velocity is shared between the immiscible fluids. The mathematical form of the momentum equation can be expressed in the following manner:(2)ρ∂V∂t+V⋅∇V=−∇p+∇⋅τ+ρg+Fs,where V denotes the velocity vector field and t denotes the time. p denotes the pressure field and g denotes the gravitational acceleration vector and τ=µ∇V+∇VT denotes the stress tensor. Fs denotes the volumetric surface tension force, and it is defined by the continuum surface force (CSF) model [58]. The surface tension force is considered as volumetric body force in the CSF model and pressure jump on the fluid interfaces due to surface tension is included in this model. It is defined by the following mathematical expression in the CFS model:(3)Fs=σk∇α2,where σ denotes the surface tension and k=∇⋅∇α2/∇α2 denotes the curvature of the interface.Volume fraction equation:In the VOF model, the volume fraction of the secondary phaseα2 is a function of space and time. It can be expressed by the following mathematical expressions:(4)∂α2∂t+∇⋅α2V=0.Mixture property equations:The average mixture properties (i.e., density,ρ, and viscosity, μ) of two phases can be calculated by the following mathematical expressions:(5)µ=α2µ2+1−α2µ1,(6)ρ=α2ρ2+1−α2ρ1.The average mixture properties, i.e., apparent viscosity and densityμ, ρ, are calculated by the above equations (5) and (6), respectively. These properties are used to solving the transport equations.


### 2.3. Solution Methodology

The VOF method is adopted to capture the moving interface of the two immiscible fluids by solving a single set of momentum equations in the entire flowing field. VOF method calculates the value of volume fractions of each fluid in each cell in the computational domain and fluid interface in each cell is defined through the value of volume fractions of fluids. The transport equations (1) and (2) are coupled with the volume fraction equation (equation (4)). The pressure implicit with the splitting of operators (PISO) scheme is adapted for coupling the pressure and velocity. The momentum equation is discretized by a second-order upwind scheme. The volume fraction equation is discretized by a geometrical-reconstruction scheme. It works on linear interpolation and the precise interface of immiscible fluids are traced by this scheme. The gradients present in the governing equations are spatially discretized by the least-square scheme and the pressure terms are discretized by the PRESTO scheme. The dynamic behaviour of the two-phase fluids flow is furnished within a time step 10−4 s. It is taken fixed for each simulation. Each numerical simulation is initialized from the inlet. The transient formulations are handled by a first-order implicit scheme. The convergence of numerical simulation is based on 10−6 relative error criteria. Numerical solutions are stabilized by using Courant-Friedrichs-Lewy (CFL) criterion, i.e., courant number must be less than one (0.5) in each iteration.


#### 2.3.1. Boundary Conditions

The results of numerical simulations are based on initial boundary conditions. Here, velocity inlet and pressure outlet are selected as inlet and outer boundary conditions. The stationary walls with no-slip criteria are adopted as walls boundary conditions. Initially, the entire computational domain is filled with more viscous fluid-1 (i.e., silicon oil). After that, fluid-2 (i.e., water) is invaded with the velocity of 0.015 m/s through the inlet into the silicon oil-saturated Y-shaped channel. The pressure at the outlets is considered as constant ambient pressure.


#### 2.3.2. Physical Properties of Fluids

The high viscous fluid-1, i.e., silicon oil, is considered as displaced (defending) fluid because it is used in vitreoretinal surgery and many more clinical applications [59]. The low-viscous fluid-2, i.e., water, is considered as displacing fluid. Silicon oil is more viscous than water and both are immiscible. The physical properties of working fluids are given in Table 1.Table 1
Physical properties of fluids.
FluidsDensity (kg/m3)Viscosity (Pa-s)Surface tension (N/m)Viscosity ratio,μr=μ2/μ1Fluid-1 (silicone oil)8500.0500.0300.02Fluid-2 (water)10000.001


### 2.4. Assumptions

The following assumptions have been taken to perform numerical simulations:(i)
The fluid pair should be Newtonian(ii)
The fluid pair should be immiscible and incompressible(iii)
Flow should be unsteady and laminar(iv)
Reynolds’s number should be very low(v)
The viscosity ratio should be less than one(vi)
Inertia and gravity forces should be absent


### 2.5. Grid-Independent Test

Before the primary investigation, a grid-independent test has been performed to ensure that the numerical results do not depend on the grid size or the number of grids. The optimum number of grids at which the numerical results do not alter is found to be 98600. All numerical simulations have been performed at the optimum number of grids.


## 2.1. Problem Formulation

In this work, we investigate the effect of surface wettability on the splitting dynamics of a low-viscous finger in the single bifurcating channel. A schematic diagram of a two-dimensional (2D) single bifurcating channel with inlet and outlet boundaries conditions is shown in Figure1(a). The parent channel is bifurcated into two identical daughter branches which are shown in Figure 1(a) as branch-1 and branch-2. The size of each daughter branch Ld1×Wd1=Ld2×Wd2 is equal and it is taken as 5mm×0.75mm. The length ratio Lr=Ld/Lp and width ratio Wr=Wd/Wp of the bifurcating channel are taken as 1 and 0.75, respectively. The bifurcation angle (β) is taken as 39°. The snapshot of meshed computational domain with proper dimensions is shown in Figure 1(b). The meshing of the computational domain has been done by the meshing module of ANSYS FLUENT software. The fine structural meshing has been adopted for precise results. Initially, the entire Y-shaped channel is filled with high viscous fluid-1 (red) after that fluid-2 is invaded through the inlet with the velocity, V. The orientation of the Y-shaped channel is considered as horizontal. Therefore, the effect of gravity force due to density difference may be neglected in this investigation. The nature of both fluids is considered immiscible and incompressible. The fluid-2 is injected into the bifurcating channel with a constant velocity of 0.015 m/s. The pressures at the outlets of the daughter branches are considered as atmospheric pressure (P = 0 Pa). The wall surface wettability is governed by the wetting contact angle and it is defined by Young’s equation, i.e., cosθ=σw1−σw2/σ, where θ denotes the wall wetting contact angle and σw1 denotes interfacial tension between wall and fluid-1, σw2 denotes interfacial tension between wall and fluid-2, and σ denotes interfacial tension between fluid-1 and fluid-2. After several investigations, four different wall surface wettability configurations of daughter branches are taken, keeping the parent channel as strongly hydrophobic or nonwetting θp=180°. These surface wettability configurations are created by incorporating hydrophilic, hydrophobic, and superhydrophobic wettability conditions. All surface wettability configurations are shown in Figure 2. In these figures, θp,θ1, and θ2 are representing wall wetting contact angle of parent channel, daughter branch-1, and daughter branch-2, respectively. The surface wettability configurations of daughter branches are as follows:(i)
The surfaces of both daughter branches are taken as hydrophilicθ1,θ2=78°,78°(ii)
The surfaces of branch-1 and branch-2 are taken as hydrophilic and hydrophobicθ1,θ2=78°,118°, respectively(iii)
The surfaces of branch-1 and branch-2 are taken as hydrophilic and superhydrophobicθ1,θ2=78°,150°, respectively(iv)
The surfaces of both daughter branches are taken as superhydrophobicθ1,θ2=150°,150°Figure 1
(a) Schematic diagram of the two-dimensional single bifurcating channel. (b) Snapshot of the meshed computational domain.
(a)(b)Figure 2
Different wall surface wettability configurations of daughter branches keeping parent channel as the nonwettingθp=180°.
(a)(b)(c)(d)


## 2.2. Mathematical Model

ANSYS FLUENT commercial CFD (computational fluid dynamic) software has been used for numerical simulations. The numerical simulations have been performed by considering a pressure-based solver with unsteady fluid displacement in the two-dimensional Y-shaped channel. The shape of the moving interface is tracked by using the volume of fluid (VOF) method. The VOF method is one of the numerical tools of interface tracking which is used to track the interface of two or more immiscible fluids in the multiphase flows [57]. It is more renowned among interface tracking methods because it follows the conservation of mass. In the VOF method, the transport equations are solved simultaneously for each phase with a single set of continuity and momentum equations in the entire flow field. In the present two-phase displacement flow, the volume fraction of primary and secondary phases is represented by α1 and α2, respectively. The sum of the volume fractions of two phases must be unity α1+α2=1 in each cell as well as the entire computational domain. The volume fractions are a function of space and time. The volume fraction equation is solved only for the second phase in the VOF model and the result; we get the value of α2. The value of α1 is calculated from the conservation equation of volume fraction, i.e., α1=1−α2. The moving interface is constructed in the computational domain on the value of volume fractions. The fluids inside the computational domain are defined by values of α2 in the following manner: α2=0⇒fluid1;1⇒fluid2;0<α2<1⇒interface.We assume that unsteady laminar flow occurs in the entire investigation and primary and secondary phases, i.e., fluid-1 (silicon oil) and fluid-2 (water) are immiscible and incompressible.The governing equations of immiscible and incompressible fluid flows can be expressed by the following mathematical expressions:Continuity equation:(1)∇⋅V=0.Momentum conservation equation (equation (2)):A single set of momentum equations is solved in entire the computational domain, and the obtaining velocity is shared between the immiscible fluids. The mathematical form of the momentum equation can be expressed in the following manner:(2)ρ∂V∂t+V⋅∇V=−∇p+∇⋅τ+ρg+Fs,where V denotes the velocity vector field and t denotes the time. p denotes the pressure field and g denotes the gravitational acceleration vector and τ=µ∇V+∇VT denotes the stress tensor. Fs denotes the volumetric surface tension force, and it is defined by the continuum surface force (CSF) model [58]. The surface tension force is considered as volumetric body force in the CSF model and pressure jump on the fluid interfaces due to surface tension is included in this model. It is defined by the following mathematical expression in the CFS model:(3)Fs=σk∇α2,where σ denotes the surface tension and k=∇⋅∇α2/∇α2 denotes the curvature of the interface.Volume fraction equation:In the VOF model, the volume fraction of the secondary phaseα2 is a function of space and time. It can be expressed by the following mathematical expressions:(4)∂α2∂t+∇⋅α2V=0.Mixture property equations:The average mixture properties (i.e., density,ρ, and viscosity, μ) of two phases can be calculated by the following mathematical expressions:(5)µ=α2µ2+1−α2µ1,(6)ρ=α2ρ2+1−α2ρ1.The average mixture properties, i.e., apparent viscosity and densityμ, ρ, are calculated by the above equations (5) and (6), respectively. These properties are used to solving the transport equations.


## 2.3. Solution Methodology

The VOF method is adopted to capture the moving interface of the two immiscible fluids by solving a single set of momentum equations in the entire flowing field. VOF method calculates the value of volume fractions of each fluid in each cell in the computational domain and fluid interface in each cell is defined through the value of volume fractions of fluids. The transport equations (1) and (2) are coupled with the volume fraction equation (equation (4)). The pressure implicit with the splitting of operators (PISO) scheme is adapted for coupling the pressure and velocity. The momentum equation is discretized by a second-order upwind scheme. The volume fraction equation is discretized by a geometrical-reconstruction scheme. It works on linear interpolation and the precise interface of immiscible fluids are traced by this scheme. The gradients present in the governing equations are spatially discretized by the least-square scheme and the pressure terms are discretized by the PRESTO scheme. The dynamic behaviour of the two-phase fluids flow is furnished within a time step 10−4 s. It is taken fixed for each simulation. Each numerical simulation is initialized from the inlet. The transient formulations are handled by a first-order implicit scheme. The convergence of numerical simulation is based on 10−6 relative error criteria. Numerical solutions are stabilized by using Courant-Friedrichs-Lewy (CFL) criterion, i.e., courant number must be less than one (0.5) in each iteration.


### 2.3.1. Boundary Conditions

The results of numerical simulations are based on initial boundary conditions. Here, velocity inlet and pressure outlet are selected as inlet and outer boundary conditions. The stationary walls with no-slip criteria are adopted as walls boundary conditions. Initially, the entire computational domain is filled with more viscous fluid-1 (i.e., silicon oil). After that, fluid-2 (i.e., water) is invaded with the velocity of 0.015 m/s through the inlet into the silicon oil-saturated Y-shaped channel. The pressure at the outlets is considered as constant ambient pressure.


### 2.3.2. Physical Properties of Fluids

The high viscous fluid-1, i.e., silicon oil, is considered as displaced (defending) fluid because it is used in vitreoretinal surgery and many more clinical applications [59]. The low-viscous fluid-2, i.e., water, is considered as displacing fluid. Silicon oil is more viscous than water and both are immiscible. The physical properties of working fluids are given in Table 1.Table 1
Physical properties of fluids.
FluidsDensity (kg/m3)Viscosity (Pa-s)Surface tension (N/m)Viscosity ratio,μr=μ2/μ1Fluid-1 (silicone oil)8500.0500.0300.02Fluid-2 (water)10000.001


## 2.3.1. Boundary Conditions

The results of numerical simulations are based on initial boundary conditions. Here, velocity inlet and pressure outlet are selected as inlet and outer boundary conditions. The stationary walls with no-slip criteria are adopted as walls boundary conditions. Initially, the entire computational domain is filled with more viscous fluid-1 (i.e., silicon oil). After that, fluid-2 (i.e., water) is invaded with the velocity of 0.015 m/s through the inlet into the silicon oil-saturated Y-shaped channel. The pressure at the outlets is considered as constant ambient pressure.


## 2.3.2. Physical Properties of Fluids

The high viscous fluid-1, i.e., silicon oil, is considered as displaced (defending) fluid because it is used in vitreoretinal surgery and many more clinical applications [59]. The low-viscous fluid-2, i.e., water, is considered as displacing fluid. Silicon oil is more viscous than water and both are immiscible. The physical properties of working fluids are given in Table 1.Table 1
Physical properties of fluids.
FluidsDensity (kg/m3)Viscosity (Pa-s)Surface tension (N/m)Viscosity ratio,μr=μ2/μ1Fluid-1 (silicone oil)8500.0500.0300.02Fluid-2 (water)10000.001


## 2.4. Assumptions

The following assumptions have been taken to perform numerical simulations:(i)
The fluid pair should be Newtonian(ii)
The fluid pair should be immiscible and incompressible(iii)
Flow should be unsteady and laminar(iv)
Reynolds’s number should be very low(v)
The viscosity ratio should be less than one(vi)
Inertia and gravity forces should be absent


## 2.5. Grid-Independent Test

Before the primary investigation, a grid-independent test has been performed to ensure that the numerical results do not depend on the grid size or the number of grids. The optimum number of grids at which the numerical results do not alter is found to be 98600. All numerical simulations have been performed at the optimum number of grids.


## 3. Results and Discussion

### 3.1. Validation of Numerical Method

Before presenting numerical results, validation of the numerical method is needed to ensure that the numerical solver is valid for the current investigation. Here, the numerical method is validated with the experimental study of Baroud et al. [56]. We have numerically replicated the experimental setup of Baroud et al. by assuming perfluorodecalin (PFD) as a displacing fluid and silicone oil as a displaced fluid inside a single bifurcating channel. They considered that fluids were immiscible and incompressible. They took the viscosity of PFD and silicon oil as 2.9 cP and 100 cP, respectively. They took surface tension at the fluid interface, 50 mN/m. They took the size of the single bifurcating channel as 200 μm × 100 μm. The entire Y-shaped channel was filled with silicon oil, after that displacing fluid was injected through a syringe pump with the range of flow rates, 1–60 nls−1. The image of the experimental result, i.e., propagation of low-viscous finger inside the single bifurcation channel, is shown in Figure 3(a). Numerical simulation has been performed by considering the same setups as well as the same working fluids which were used in the experimental investigation. We have used the VOF model for tracking the interface. The numerical result of the propagation of the low-viscous finger inside the single bifurcation channel is shown in Figure 3(b). The numerical result is compared with the experimental result and a good agreement is to be found. Figure 3 shows a good agreement between the numerical result and experiment result so that the numerical method is valid for the current investigation and the accuracy of this method is also appropriate.Figure 3
Validation: comparison of propagation of low-viscous finger inside the single bifurcating channel: (a) experimental snapshot [56] and (b) computational snapshot of volume fraction.


### 3.2. Mechanism of Low-Viscous Fingers Formation in the Horizontal Channel

We have assumed that inertia and gravity forces are absent (tend to be zero) because quasistatic flow (creep flow) occurs in the horizontal channel. Therefore, the momentum equation for one-dimensional flow in the parent channel can be written in the following manner [60]:(7a)V=−h212µdpdx=−kµdpdx,where h2/12 is equal to a constant k and h denotes channel height.We can rewrite the momentum equation in the form of a pressure gradient in the following manner:(7b)dpdx=−µVkThe net pressure force,p2−p1, is exerted by displacing fluid-2 on displaced fluid-1. Then, stationary fluid-1 is virtually displaced, dx, from the initial position in the parent channel.The net pressure force exerted by fluid-2 on fluid-1 can be expressed by the following mathematical expression:(7c)p2−p1=µ1−µ2Vkdx.Here, if the net pressure force is positive, then any small displacement of fluid-1 will be created viscous fingering instability on the interface of fluid-1 and fluid-2. We have already assumed that fluid-1 is more viscous than fluid-2 so thatµ1−µ2 will be positive. In equation (7c), each term of right hand side is positive. Therefore, the net pressure force must be positive. We can say that viscous fingering instability is always evolved on the fluid interface at any velocity when less viscous fluid displaces another more viscous fluid in the horizontal channel (gravity force absent). In other words, if the more viscous fluid is being driven out by a less viscous fluid in horizontal channel, then fingering instability always appears on the fluid interface at any velocity. If the viscosity of both fluids is equal or fluid-2 is more viscous than fluid-1, then the interface will be always flat or stable (no fingering instability). Note down, the above information is only valid for displacement flow in horizontal channel (not for inclined channel), because in this flow, density changes do not play any role in fingering instability formation. The physical explanation of viscous fingering instability is that the effective viscosity is always lower in a thin band, hence it moves faster; therefore, low-viscous finger-shaped instability appears on fluid interface in the displacement flow.In the present study, fluid-fluid displacement flows are quasistatic because the Reynolds numberRe=ρ1VWp/μ1 in the entire flows domain is considered to be less than one Re=0.25. In quasistatic flows, the effect of the inertia force may be neglected. The impact of gravity force may be ignored when displacement flow occurs in the horizontal microchannel. Here, the effects of inertia force and the gravity force are neglected because quasistatic flow occurs in the horizontal Y-shaped channel. Immiscible displacement flow is driven by two forces, i.e., viscous force and surface tension force. In another way, viscous force is responsible for the overall pressure gradient, but surface tension force is responsible for pressure jump on the interfaces. The ratio of viscous force and surface tension force is called the capillary number Ca=µ1V/σ. It measures the relative importance of viscous force as compared to surface tension force.Here, we investigate the dynamics of water finger propagation and the behaviour of finger splitting inside the single bifurcating channel using the VOF model. Numerical simulations have been performed by considering four different surface wettability configurations of daughter branches, i.e., (i)θ1,θ2=78°,78°, (ii) θ1,θ2=78°,118°, (iii) θ1,θ2=78°,150°, and (iv) θ1,θ2=150°,150°. The capillary number is taken, 0.025, in each simulation. The primary phase (silicon oil) is initially filled into the entire Y-shaped channel by patching the volume fraction of the secondary phase (water) which is zero α2=0. It is shown in Figure 4(a) and the red color represents the silicon oil in this figure. After that, water is injected through the inlet with a velocity of 0.015 m/s. The viscosity of displacing fluid, i.e., water, is less than the viscosity of the displaced fluid, i.e., silicon oil, so that unstable displacement flow occurs, and as a result, a finger of water appears on the fluid interface. The water finger propagates towards the displacement direction inside the parent channel. Displacement of silicon oil by water from the parent channel or propagation of water finger in the silicon oil-saturated parent channel is shown in Figure 4. Propagation of water finger is captured at equal time steps. The unit of time steps (t) is taken in seconds (s). The shape of blue color in each figure represents water finger.Figure 4
Water finger (blue) propagation inside the parent channel at equal time steps (seconds): (a)t = 0.00, (b) t = 0.05, (c) t = 0.10, (d) t = 0.15, (e) t = 0.20, and (f) t = 0.25.
(a)(b)(c)(d)(e)(f)


### 3.3. Splitting Dynamics

From the numerical results, finger splitting is categorized into three types based on splitting ratio: symmetrical splitting, nonsymmetrical splitting, and reversal splitting (no splitting). The splitting ratio (λ) is defined as the ratio of finger lengths in branch-2 and branch-1.(8)λ=Lf2/Lf1,where, λ denotes the splitting ratio. Lf1 and Lf2 denote the finger lengths in branch-1 and branch-2. According to the splitting ratio, fingers splitting dynamics are classified in the following manner:(9)λ=1,⇒Symmetrical splitting,0<λ<1,⇒Nonsymmetrical splitting,0,⇒Reversalnosplitting.The physical significance of the splitting ratio is that whenever finger lengths in branch-1 and branch-2 are the same, then it looks symmetrical so that it is called symmetrical splitting. It has one splitting ratio valueλ=1. Finger lengths are not the same in both the branches, i.e. finger length in branch-2 is less compared to the finger length in branch-1. It is called nonsymmetrical splitting. It has less than one splitting ratio value λ<1. The finger length in branch-2 is zero or the finger length in branch-2 is zero compared to finger length in branch-1. It is called reversal (no) splitting. It has zero splitting ratio value λ=0.


#### 3.3.1. Symmetrical Splitting

Figures5 and 6 depict symmetrical splitting. An immiscible single finger of water propagates into the parent channel. At the bifurcating tip, it splits into two fingers and each finger propagates into each daughter branch. It is observed that finger splitting behaviour depends upon the surface wettability of the daughter branches. Figures 5 and 6 show the symmetrical finger splitting behaviour inside the single bifurcating channel when both the daughter branches are hydrophilic, i.e., (i) θ1,θ2=78°,78° and superhydrophobic, i.e., (iv) θ1,θ2=150°,150°, respectively. It is noticed that in both cases, symmetrical splitting occurs. The splitting ratio in both cases is unity (λ = 1). In the symmetrical splitting, fingers propagate into the daughter branches with the same velocity because the surface wettability of daughter branches is the same. It is concluded that symmetrical splitting may be achieved when the surface wettability of both the daughter branches is either hydrophilic or superhydrophobic. The tips of the fingers are slightly concave when the surfaces of daughter branches are hydrophilic and convex when the surfaces of daughter branches are superhydrophobic. The velocities of fingers are approximately the same in both cases.Figure 5
Symmetrical finger splitting (λ = 1) inside the single bifurcating channel when the surfaces of both the daughter branches are hydrophilic θ1,θ2=78°,78°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
Symmetrical finger splitting (λ = 1) inside the single bifurcating channel when the surfaces of both the daughter branches are superhydrophobic θ1,θ2=150°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)


#### 3.3.2. Nonsymmetrical Splitting

Figure7 depicts nonsymmetrical splitting. Nonsymmetrical splitting occurs when the finger velocities and finger lengths in both the daughter branches are different. If the surface wettability of both daughter branches is different, then resistance offered by the displaced fluid against the displacing fluid inside daughter branches is also different. Therefore, finger velocities in both the daughter branches are not the same and finger lengths in both daughter branches are also not equal. It can be seen in Figure 7. Here, the surface of branch-1 is hydrophilic and branch-2 is hydrophobic, i.e., (ii) θ1,θ2=78°,118°. The finger velocity in branch-1 is more than the finger velocity in branch-2, so the finger length in branch-1 is more than the finger length in branch-2. The reason is that the net effect of hydraulic resistance and intermolecular resistance is less in branch-1 with respect to branch-2. The intermolecular force appears due to surface wettability and it provides additional resistance against finger propagation. The splitting ratio is calculated less than one (λ < 0.20) in nonsymmetrical splitting. It is observed that the finger velocity and finger length are more in hydrophilic daughter branch-1. Nonsymmetrical splitting is found when the surface wettability of one of the daughter branches is hydrophilic and another is hydrophobic.Figure 7
Nonsymmetrical finger splitting (λ < 1) inside the single bifurcating channel when the surfaces of daughter branches are hydrophilic and hydrophobic θ1,θ2=78°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)


#### 3.3.3. Reversal Splitting

Figure8 depicts the reversal (no) splitting. Here, the surface of daughter branch-1 is hydrophilic and daughter branch-2 is superhydrophobic, i.e., (iii) θ1,θ2=78°,150°. It can be seen that the water finger propagates toward the bifurcating tip but it does not split at bifurcating tip. It passes completely through the daughter branch-1 without splitting. Hydrophilic surface provides more affinity to fluid-wall interaction as compared to superhydrophobic. The net effect of hydraulic resistance and intermolecular resistance is very large in branch-2 compared to branch-1 because superhydrophobic surface offers strong additional resistance against water finger propagation. Therefore, the water finger passes completely through hydrophilic branch-1 without splitting. It can be concluded that reversal splitting evolves when the surface wettability of one of the daughter branches is hydrophilic and another is superhydrophobic. The low-viscous finger always passes through the hydrophilic daughter branch. The splitting ratio is zero in this splitting behaviour.Figure 8
Reversal finger splitting (λ = 0) inside the single bifurcating channel when the surfaces of daughter branches are hydrophilic and superhydrophobic θ1,θ2=78°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)In Figure9, all finger splitting behaviours are incorporated. Numerical simulations are performed at Ca = 0.025, Re = 0.25, and μr=0.02. It is seen that the splitting ratio is one when the daughter branches have surface wettability configurations, i.e., θ1,θ2=78°,78° and θ1,θ2=150°,150°. This splitting behaviour is called symmetrical splitting. The splitting ratio is 0.20 when the daughter branches have surface wettability configurations, i.e., θ1,θ2=78°,118°. This splitting behaviour is called nonsymmetrical spitting. It is noticed that the splitting ratio is zero when the daughter branches have surface wettability configurations θ1,θ2=78°,150°. This splitting behaviour is called reversal splitting.Figure 9
Different types of splitting behaviour of a water finger through the surface wettability of daughter branches at Ca = 0.025, Re = 0.25, and, V.R = 0.02.


## 3.1. Validation of Numerical Method

Before presenting numerical results, validation of the numerical method is needed to ensure that the numerical solver is valid for the current investigation. Here, the numerical method is validated with the experimental study of Baroud et al. [56]. We have numerically replicated the experimental setup of Baroud et al. by assuming perfluorodecalin (PFD) as a displacing fluid and silicone oil as a displaced fluid inside a single bifurcating channel. They considered that fluids were immiscible and incompressible. They took the viscosity of PFD and silicon oil as 2.9 cP and 100 cP, respectively. They took surface tension at the fluid interface, 50 mN/m. They took the size of the single bifurcating channel as 200 μm × 100 μm. The entire Y-shaped channel was filled with silicon oil, after that displacing fluid was injected through a syringe pump with the range of flow rates, 1–60 nls−1. The image of the experimental result, i.e., propagation of low-viscous finger inside the single bifurcation channel, is shown in Figure 3(a). Numerical simulation has been performed by considering the same setups as well as the same working fluids which were used in the experimental investigation. We have used the VOF model for tracking the interface. The numerical result of the propagation of the low-viscous finger inside the single bifurcation channel is shown in Figure 3(b). The numerical result is compared with the experimental result and a good agreement is to be found. Figure 3 shows a good agreement between the numerical result and experiment result so that the numerical method is valid for the current investigation and the accuracy of this method is also appropriate.Figure 3
Validation: comparison of propagation of low-viscous finger inside the single bifurcating channel: (a) experimental snapshot [56] and (b) computational snapshot of volume fraction.


## 3.2. Mechanism of Low-Viscous Fingers Formation in the Horizontal Channel

We have assumed that inertia and gravity forces are absent (tend to be zero) because quasistatic flow (creep flow) occurs in the horizontal channel. Therefore, the momentum equation for one-dimensional flow in the parent channel can be written in the following manner [60]:(7a)V=−h212µdpdx=−kµdpdx,where h2/12 is equal to a constant k and h denotes channel height.We can rewrite the momentum equation in the form of a pressure gradient in the following manner:(7b)dpdx=−µVkThe net pressure force,p2−p1, is exerted by displacing fluid-2 on displaced fluid-1. Then, stationary fluid-1 is virtually displaced, dx, from the initial position in the parent channel.The net pressure force exerted by fluid-2 on fluid-1 can be expressed by the following mathematical expression:(7c)p2−p1=µ1−µ2Vkdx.Here, if the net pressure force is positive, then any small displacement of fluid-1 will be created viscous fingering instability on the interface of fluid-1 and fluid-2. We have already assumed that fluid-1 is more viscous than fluid-2 so thatµ1−µ2 will be positive. In equation (7c), each term of right hand side is positive. Therefore, the net pressure force must be positive. We can say that viscous fingering instability is always evolved on the fluid interface at any velocity when less viscous fluid displaces another more viscous fluid in the horizontal channel (gravity force absent). In other words, if the more viscous fluid is being driven out by a less viscous fluid in horizontal channel, then fingering instability always appears on the fluid interface at any velocity. If the viscosity of both fluids is equal or fluid-2 is more viscous than fluid-1, then the interface will be always flat or stable (no fingering instability). Note down, the above information is only valid for displacement flow in horizontal channel (not for inclined channel), because in this flow, density changes do not play any role in fingering instability formation. The physical explanation of viscous fingering instability is that the effective viscosity is always lower in a thin band, hence it moves faster; therefore, low-viscous finger-shaped instability appears on fluid interface in the displacement flow.In the present study, fluid-fluid displacement flows are quasistatic because the Reynolds numberRe=ρ1VWp/μ1 in the entire flows domain is considered to be less than one Re=0.25. In quasistatic flows, the effect of the inertia force may be neglected. The impact of gravity force may be ignored when displacement flow occurs in the horizontal microchannel. Here, the effects of inertia force and the gravity force are neglected because quasistatic flow occurs in the horizontal Y-shaped channel. Immiscible displacement flow is driven by two forces, i.e., viscous force and surface tension force. In another way, viscous force is responsible for the overall pressure gradient, but surface tension force is responsible for pressure jump on the interfaces. The ratio of viscous force and surface tension force is called the capillary number Ca=µ1V/σ. It measures the relative importance of viscous force as compared to surface tension force.Here, we investigate the dynamics of water finger propagation and the behaviour of finger splitting inside the single bifurcating channel using the VOF model. Numerical simulations have been performed by considering four different surface wettability configurations of daughter branches, i.e., (i)θ1,θ2=78°,78°, (ii) θ1,θ2=78°,118°, (iii) θ1,θ2=78°,150°, and (iv) θ1,θ2=150°,150°. The capillary number is taken, 0.025, in each simulation. The primary phase (silicon oil) is initially filled into the entire Y-shaped channel by patching the volume fraction of the secondary phase (water) which is zero α2=0. It is shown in Figure 4(a) and the red color represents the silicon oil in this figure. After that, water is injected through the inlet with a velocity of 0.015 m/s. The viscosity of displacing fluid, i.e., water, is less than the viscosity of the displaced fluid, i.e., silicon oil, so that unstable displacement flow occurs, and as a result, a finger of water appears on the fluid interface. The water finger propagates towards the displacement direction inside the parent channel. Displacement of silicon oil by water from the parent channel or propagation of water finger in the silicon oil-saturated parent channel is shown in Figure 4. Propagation of water finger is captured at equal time steps. The unit of time steps (t) is taken in seconds (s). The shape of blue color in each figure represents water finger.Figure 4
Water finger (blue) propagation inside the parent channel at equal time steps (seconds): (a)t = 0.00, (b) t = 0.05, (c) t = 0.10, (d) t = 0.15, (e) t = 0.20, and (f) t = 0.25.
(a)(b)(c)(d)(e)(f)


## 3.3. Splitting Dynamics

From the numerical results, finger splitting is categorized into three types based on splitting ratio: symmetrical splitting, nonsymmetrical splitting, and reversal splitting (no splitting). The splitting ratio (λ) is defined as the ratio of finger lengths in branch-2 and branch-1.(8)λ=Lf2/Lf1,where, λ denotes the splitting ratio. Lf1 and Lf2 denote the finger lengths in branch-1 and branch-2. According to the splitting ratio, fingers splitting dynamics are classified in the following manner:(9)λ=1,⇒Symmetrical splitting,0<λ<1,⇒Nonsymmetrical splitting,0,⇒Reversalnosplitting.The physical significance of the splitting ratio is that whenever finger lengths in branch-1 and branch-2 are the same, then it looks symmetrical so that it is called symmetrical splitting. It has one splitting ratio valueλ=1. Finger lengths are not the same in both the branches, i.e. finger length in branch-2 is less compared to the finger length in branch-1. It is called nonsymmetrical splitting. It has less than one splitting ratio value λ<1. The finger length in branch-2 is zero or the finger length in branch-2 is zero compared to finger length in branch-1. It is called reversal (no) splitting. It has zero splitting ratio value λ=0.


### 3.3.1. Symmetrical Splitting

Figures5 and 6 depict symmetrical splitting. An immiscible single finger of water propagates into the parent channel. At the bifurcating tip, it splits into two fingers and each finger propagates into each daughter branch. It is observed that finger splitting behaviour depends upon the surface wettability of the daughter branches. Figures 5 and 6 show the symmetrical finger splitting behaviour inside the single bifurcating channel when both the daughter branches are hydrophilic, i.e., (i) θ1,θ2=78°,78° and superhydrophobic, i.e., (iv) θ1,θ2=150°,150°, respectively. It is noticed that in both cases, symmetrical splitting occurs. The splitting ratio in both cases is unity (λ = 1). In the symmetrical splitting, fingers propagate into the daughter branches with the same velocity because the surface wettability of daughter branches is the same. It is concluded that symmetrical splitting may be achieved when the surface wettability of both the daughter branches is either hydrophilic or superhydrophobic. The tips of the fingers are slightly concave when the surfaces of daughter branches are hydrophilic and convex when the surfaces of daughter branches are superhydrophobic. The velocities of fingers are approximately the same in both cases.Figure 5
Symmetrical finger splitting (λ = 1) inside the single bifurcating channel when the surfaces of both the daughter branches are hydrophilic θ1,θ2=78°,78°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
Symmetrical finger splitting (λ = 1) inside the single bifurcating channel when the surfaces of both the daughter branches are superhydrophobic θ1,θ2=150°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)


### 3.3.2. Nonsymmetrical Splitting

Figure7 depicts nonsymmetrical splitting. Nonsymmetrical splitting occurs when the finger velocities and finger lengths in both the daughter branches are different. If the surface wettability of both daughter branches is different, then resistance offered by the displaced fluid against the displacing fluid inside daughter branches is also different. Therefore, finger velocities in both the daughter branches are not the same and finger lengths in both daughter branches are also not equal. It can be seen in Figure 7. Here, the surface of branch-1 is hydrophilic and branch-2 is hydrophobic, i.e., (ii) θ1,θ2=78°,118°. The finger velocity in branch-1 is more than the finger velocity in branch-2, so the finger length in branch-1 is more than the finger length in branch-2. The reason is that the net effect of hydraulic resistance and intermolecular resistance is less in branch-1 with respect to branch-2. The intermolecular force appears due to surface wettability and it provides additional resistance against finger propagation. The splitting ratio is calculated less than one (λ < 0.20) in nonsymmetrical splitting. It is observed that the finger velocity and finger length are more in hydrophilic daughter branch-1. Nonsymmetrical splitting is found when the surface wettability of one of the daughter branches is hydrophilic and another is hydrophobic.Figure 7
Nonsymmetrical finger splitting (λ < 1) inside the single bifurcating channel when the surfaces of daughter branches are hydrophilic and hydrophobic θ1,θ2=78°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)


### 3.3.3. Reversal Splitting

Figure8 depicts the reversal (no) splitting. Here, the surface of daughter branch-1 is hydrophilic and daughter branch-2 is superhydrophobic, i.e., (iii) θ1,θ2=78°,150°. It can be seen that the water finger propagates toward the bifurcating tip but it does not split at bifurcating tip. It passes completely through the daughter branch-1 without splitting. Hydrophilic surface provides more affinity to fluid-wall interaction as compared to superhydrophobic. The net effect of hydraulic resistance and intermolecular resistance is very large in branch-2 compared to branch-1 because superhydrophobic surface offers strong additional resistance against water finger propagation. Therefore, the water finger passes completely through hydrophilic branch-1 without splitting. It can be concluded that reversal splitting evolves when the surface wettability of one of the daughter branches is hydrophilic and another is superhydrophobic. The low-viscous finger always passes through the hydrophilic daughter branch. The splitting ratio is zero in this splitting behaviour.Figure 8
Reversal finger splitting (λ = 0) inside the single bifurcating channel when the surfaces of daughter branches are hydrophilic and superhydrophobic θ1,θ2=78°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)In Figure9, all finger splitting behaviours are incorporated. Numerical simulations are performed at Ca = 0.025, Re = 0.25, and μr=0.02. It is seen that the splitting ratio is one when the daughter branches have surface wettability configurations, i.e., θ1,θ2=78°,78° and θ1,θ2=150°,150°. This splitting behaviour is called symmetrical splitting. The splitting ratio is 0.20 when the daughter branches have surface wettability configurations, i.e., θ1,θ2=78°,118°. This splitting behaviour is called nonsymmetrical spitting. It is noticed that the splitting ratio is zero when the daughter branches have surface wettability configurations θ1,θ2=78°,150°. This splitting behaviour is called reversal splitting.Figure 9
Different types of splitting behaviour of a water finger through the surface wettability of daughter branches at Ca = 0.025, Re = 0.25, and, V.R = 0.02.


## 3.3.1. Symmetrical Splitting

Figures5 and 6 depict symmetrical splitting. An immiscible single finger of water propagates into the parent channel. At the bifurcating tip, it splits into two fingers and each finger propagates into each daughter branch. It is observed that finger splitting behaviour depends upon the surface wettability of the daughter branches. Figures 5 and 6 show the symmetrical finger splitting behaviour inside the single bifurcating channel when both the daughter branches are hydrophilic, i.e., (i) θ1,θ2=78°,78° and superhydrophobic, i.e., (iv) θ1,θ2=150°,150°, respectively. It is noticed that in both cases, symmetrical splitting occurs. The splitting ratio in both cases is unity (λ = 1). In the symmetrical splitting, fingers propagate into the daughter branches with the same velocity because the surface wettability of daughter branches is the same. It is concluded that symmetrical splitting may be achieved when the surface wettability of both the daughter branches is either hydrophilic or superhydrophobic. The tips of the fingers are slightly concave when the surfaces of daughter branches are hydrophilic and convex when the surfaces of daughter branches are superhydrophobic. The velocities of fingers are approximately the same in both cases.Figure 5
Symmetrical finger splitting (λ = 1) inside the single bifurcating channel when the surfaces of both the daughter branches are hydrophilic θ1,θ2=78°,78°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
Symmetrical finger splitting (λ = 1) inside the single bifurcating channel when the surfaces of both the daughter branches are superhydrophobic θ1,θ2=150°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)


## 3.3.2. Nonsymmetrical Splitting

Figure7 depicts nonsymmetrical splitting. Nonsymmetrical splitting occurs when the finger velocities and finger lengths in both the daughter branches are different. If the surface wettability of both daughter branches is different, then resistance offered by the displaced fluid against the displacing fluid inside daughter branches is also different. Therefore, finger velocities in both the daughter branches are not the same and finger lengths in both daughter branches are also not equal. It can be seen in Figure 7. Here, the surface of branch-1 is hydrophilic and branch-2 is hydrophobic, i.e., (ii) θ1,θ2=78°,118°. The finger velocity in branch-1 is more than the finger velocity in branch-2, so the finger length in branch-1 is more than the finger length in branch-2. The reason is that the net effect of hydraulic resistance and intermolecular resistance is less in branch-1 with respect to branch-2. The intermolecular force appears due to surface wettability and it provides additional resistance against finger propagation. The splitting ratio is calculated less than one (λ < 0.20) in nonsymmetrical splitting. It is observed that the finger velocity and finger length are more in hydrophilic daughter branch-1. Nonsymmetrical splitting is found when the surface wettability of one of the daughter branches is hydrophilic and another is hydrophobic.Figure 7
Nonsymmetrical finger splitting (λ < 1) inside the single bifurcating channel when the surfaces of daughter branches are hydrophilic and hydrophobic θ1,θ2=78°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)


## 3.3.3. Reversal Splitting

Figure8 depicts the reversal (no) splitting. Here, the surface of daughter branch-1 is hydrophilic and daughter branch-2 is superhydrophobic, i.e., (iii) θ1,θ2=78°,150°. It can be seen that the water finger propagates toward the bifurcating tip but it does not split at bifurcating tip. It passes completely through the daughter branch-1 without splitting. Hydrophilic surface provides more affinity to fluid-wall interaction as compared to superhydrophobic. The net effect of hydraulic resistance and intermolecular resistance is very large in branch-2 compared to branch-1 because superhydrophobic surface offers strong additional resistance against water finger propagation. Therefore, the water finger passes completely through hydrophilic branch-1 without splitting. It can be concluded that reversal splitting evolves when the surface wettability of one of the daughter branches is hydrophilic and another is superhydrophobic. The low-viscous finger always passes through the hydrophilic daughter branch. The splitting ratio is zero in this splitting behaviour.Figure 8
Reversal finger splitting (λ = 0) inside the single bifurcating channel when the surfaces of daughter branches are hydrophilic and superhydrophobic θ1,θ2=78°,150°. (a) t=0.30, (b) t=0.35, (c) t=0.38, (d) t=0.40, (e) t=0.45, (f) t=0.50, (g) t=0.52, (h) t=0.55.
(a)(b)(c)(d)(e)(f)(g)(h)In Figure9, all finger splitting behaviours are incorporated. Numerical simulations are performed at Ca = 0.025, Re = 0.25, and μr=0.02. It is seen that the splitting ratio is one when the daughter branches have surface wettability configurations, i.e., θ1,θ2=78°,78° and θ1,θ2=150°,150°. This splitting behaviour is called symmetrical splitting. The splitting ratio is 0.20 when the daughter branches have surface wettability configurations, i.e., θ1,θ2=78°,118°. This splitting behaviour is called nonsymmetrical spitting. It is noticed that the splitting ratio is zero when the daughter branches have surface wettability configurations θ1,θ2=78°,150°. This splitting behaviour is called reversal splitting.Figure 9
Different types of splitting behaviour of a water finger through the surface wettability of daughter branches at Ca = 0.025, Re = 0.25, and, V.R = 0.02.


## 4. Conclusions

In the present work, the splitting dynamics of a water finger is tuned by the surface wettability of daughter branches inside the 2D single bifurcating channel. The splitting of a water finger is captured by using the VOF method. This method is applicable for only tracking the interfaces of immiscible fluids. In this numerical study, the effects of inertia force and the gravity force are neglected because quasistatic or creep flow occurs in the horizontal channel. Here, low-viscous fluid-2 (water) drives more viscous fluid-1 (silicon oil) with a velocity of 0.015 m/s in the parent channel. Due to an unfavourable viscosity ratio (0.020), a finger-like interfacial instability evolves on the fluids interface. It is called the viscous fingering phenomenon. It is driven by mainly viscous force (due to viscosity differences) and surface tension force. The ratio of these forces is called the capillary number. A fixed capillary number, 0.025, and a fixed Reynolds number, 0.25, are taken in each numerical simulation.Initially, bifurcating channel is filled with a stationary more viscous fluid-1 (silicon oil). After that, low-viscous fluid-2 (water) is invaded through the inlet. It penetrates into silicon oil in the form of finger shape. The water finger continuously moves toward the bifurcating tip and displaces silicon oil from the parent channel. But, silicon oil near the wall offered strong hydraulic resistance, therefore it is not displaced by water. It is shown in Figure4(f). Generally, a single finger splits into two fingers at the bifurcating tip and each finger propagates into each daughter’s branches, but we found finger splitting may be triggered by surface wettability of daughter branches. The numerical investigation has been carried out by considering four different surface wettability configurations of daughter branches θ1,θ2, i.e., (i) θ1,θ2=78°,78°, (ii) θ1,θ2=78°,118°, (iii) θ1,θ2=78°,150°, and (iv) θ1,θ2=150°,150. From the numerical results, we found three distinct behaviours of finger splitting. These finger splitting behaviours are classified into three types based on splitting ratio: symmetrical splitting, nonsymmetrical splitting, and reversal (no) splitting. It is observed that a water finger splits at bifurcating tip into two similar fingers when the surface wettability of both daughter branches is either hydrophilic θ1,θ2=78°,78° or superhydrophobic θ1,θ2=150°,150°. These fingers propagate into the daughter branch-1 and daughter branch-2. The lengths of both the fingers are the same so that the value of the spitting ratio is one. This type of splitting behaviour is called symmetrical splitting. It is noticed that a water finger splits into two dissimilar fingers at bifurcating tip when the surface wettability of branch-1 is hydrophilic and branch-2 is hydrophobic θ1,θ2=78°,118°. The lengths of both the fingers are not the same. The length of a finger in branch-1 is more than the length of a finger in branch-2 so that the value of the splitting ratio is less than one (0.20). This type of splitting behaviour is called nonsymmetrical splitting. It is observed that a water finger does not split at bifurcating tip when the surface wettability of branch-1 is hydrophilic and branch-2 is superhydrophobic θ1,θ2=78°,150°. The water finger passes through the hydrophilic branch-1 without splitting so that value of the spitting ratio is zero. This type of splitting behaviour is called reversal splitting. We have summarized the outcomes of finger splitting dynamics of four distinct wettability configurations of daughter branches in Table 2.Table 2
Finger splitting dynamics with four distinct wettability configurations of daughter branches at Ca = 0.025, Re = 0.25, andμr=0.02.
S. no.Wettability configuration of daughter branches,θ1,θ2Splitting ratio,λType of splitting(i)78°,78°1Symmetrical splitting(ii)78°,118°0.20Nonsymmetrical splitting(iii)78°,150°0Reversal (no) splitting(iv)150°,150°1Symmetrical splittingThe outcomes of this work may be incorporated in many engineering and biomedical applications where fluid displacement flows occur inside branching microchannels such as biochemical testing, occulted pulmonary airways reopening, starving the tumour cells in embolotherapy cancer treatment, and respiratory diseases such as COVID-19 [61].


## 5. Future Scope

Further investigation on splitting dynamics of the low-viscous finger may be performed by considering multibifurcating (networking) channels instead of single bifurcating channels. This work may be also extended by considering distinct viscosity ratio and surface tension of fluid pairs as well as varying the inlet velocity. The dimensionless terms such as length ratio, width ratio, and bifurcating angle may be used in the future for optimizing the dimensions as well as the effect of these terms on splitting ratio.



---

*Source: 3462844-2022-02-10.xml*