The objective of this research paper is to relate the influence of dynamic wetting in a liquid/liquid/solid system to the breakup of emulsion droplets in capillaries. Therefore, modeling and simulation of liquid/liquid flow through a capillary constriction have been performed with varying dynamic contact angles from highly hydrophilic to highly hydrophobic. Advanced advection schemes with geometric interface reconstruction (isoAdvector) are incorporated for high interface advection accuracy. A sharp surface tension force model is used to reduce spurious currents originating from the numerical treatment and geometric reconstruction of the surface curvature at the interface. Stress singularities from the boundary condition at the three-phase contact line are removed by applying a Navier-slip boundary condition. The simulation results illustrate the strong dependency of the wettability and the contact line and interface deformation.

Liquid/liquid formulations like emulsions are widely adopted in the industry. They play an increasing role in pharmaceutical, food, and cosmetics products. For the formulation of these emulsions, the premix membrane emulsification (PME) process may be used. As a low shear stress dispersion process, it is relevant for handling shear-sensitive media, for instance, in biological systems like protein laden liquids.1 The PME process, introduced by Suzuki et al.,2 has the advantage of adjusting the final disperse phase droplet size within a narrow size distribution. This specific particle size distribution is a key parameter for tailored high-quality emulsion formulation. The influence of process parameters on the PME process, such as trans-membrane pressure,3,4 dispersed and continuous phase velocity,5 membrane porosity, and pore size,6,7 has been investigated. The influence of fluid properties, such as dispersed phase fraction, dispersed and continuous phase viscosities, and surfactant concentrations, has been quantified in Refs. 5 and 7–11. Luhede et al.12 analyzed the effect of shear stress on the stability of multiple emulsions and achieved a continuous encapsulation process by membrane functionalization. Wollborn et al. related the shear stress inside the membrane pores to the droplet breakup mechanism, utilizing idealized membrane pore geometries13 and fully resolved membranes.14 They found that the wall interaction of the droplet with the wall leads to high shear stresses for a short time.13 The importance of the wettability in the context of liquid/solid interaction with the use of surfactants was pointed out by Schroën et al.15 They concluded that most surface interactions alter the contact angle toward 90 °, either through surfactants or the disperse phase liquid (e.g., oil) that is used and, therefore, influence the breakup mechanism inside membranes.

The contact angle between the liquid and the solid phase can be described in different reference systems with different length scales. The region resp. line where the three phases (liquid disperse phase, liquid continuous phase, solid phase) meet is generally described as the contact line (CL). According to Cox,16 the microscopic contact angle (θm) is defined by molecular length scales. Outside the molecular length scale, Voinov17 stated that the static state of the contact angle is expressed by the surface energies and that the change of the contact angle with a moving liquid over a solid boundary (dynamic contact angle) depends on viscous forces. Later, Chen et al.18 stated that near the contact line, viscous stresses grow significantly to the order of the surface tension, resulting in deformation by viscous forces, setting the value of the apparent dynamic contact angle (θa). This is in line with the findings of Vecsei et al.,19 who stated that in some distance from the contact line (outside the molecular length scale), the shape is deformed by viscous stresses that leads to the apparent dynamic contact angle (θa). The concept of the microscopic (θm) and apparent dynamic contact angle (θa) is depicted in Fig. 1(a).

FIG. 1.

Sketch of the displacement regimes from (I) stable meniscus displacement (a), (II) stable cap displacement (b), and (III) unstable fingering (c), the interface is highlighted in blue. (a) Sketch of the apparent contact angle θa in relation to the microscopic contact angle θm; defining contact line and front meniscus area. (b) Sketch of the cap region: displacing disperse phase enclosed by the leading interface and the contact line. (c) Sketch of the meniscus instability (fingering at leading interface) at an early stage after its occurrence.

FIG. 1.

Sketch of the displacement regimes from (I) stable meniscus displacement (a), (II) stable cap displacement (b), and (III) unstable fingering (c), the interface is highlighted in blue. (a) Sketch of the apparent contact angle θa in relation to the microscopic contact angle θm; defining contact line and front meniscus area. (b) Sketch of the cap region: displacing disperse phase enclosed by the leading interface and the contact line. (c) Sketch of the meniscus instability (fingering at leading interface) at an early stage after its occurrence.

Close modal

With respect to the generic configuration utilized in this study, the liquid–liquid displacement process in capillaries is used to understand interface deformation and droplet formation in complex membrane geometries. For this analysis, the contact line motion describing the interface in contact with the bounding wall is of particular importance.20 The forced displacement within capillaries have been experimentally studied for instance by Hoffman21 who later derived a correlation equation to predict the dynamic contact angle.22 Fermigier and Jenffer23 performed capillary intrusion experiments which showed good agreement with Hoffmans predicted values for the dynamic contact angle. A meniscus instability residing a macroscopic film of the displaced fluid on the wall could be observed by Hoffman, where the meniscus moves at a higher velocity than the contact line. More recently, this instability is discussed in the context of the wetting transition of partially wetting fluids which has been experimentally observed.24–26 The authors reported the existence of a critical capillary number above which a finger of the displacing liquid in streamwise direction is formed, which refers to the balance between viscous forces and surface tension. Below a critical capillary number, the meniscus advances with a stable shape. Zhao et al.24 found the wetting transition to be dependent on the contact angle, which is related to the contact line motion. For the analysis of the contact line propagation and deformation this work will focus on the analysis of the intrusion of oil in a capillary initially filled with water, within a computational fluid dynamics framework. Therefore, computational fluid dynamics with a specific interface formulation is utilized in the OpenFoam software package. The term forced displacement is assigned to displacement regions from stable meniscus displacement to unstable fingering with the intermediate region of the stable cap displacement (Fig. 1) and related to the influence of the wettability. The stable cap displacement refers to an interface, which deviates from the stable meniscus with a constant curvature radius as the dynamic pressure bulges and deforms the interface from its meniscus shape, set by the apparent contact angle. The front region of the disperse phase respective to the cap propagates without the occurrence of a fingering instability.

Membranes used in premix emulsification have, in general, a complex geometric pore structure with a branched capillary network in which breakup and coalescence events simultaneously occur. This fact makes it challenging to observe single droplet breakup mechanisms directly. To address the complexity of two-phase flow through a complex shaped pore system, a straight microcapillary is used as geometric constraint for the present simulations. As the studied flow case is assumed to be axisymmetric, a two dimensional axisymmetric domain is used as shown in Fig. 2(a). The capillary is initially filled with water (liquid 1 cont. phase). At the start of the simulation, the oil (liquid 2 disp. phase) is injected into the domain from the inlet boundary as shown in Fig. 2(b). Consequently, the three phase (liquid 1/solid wall/liquid 2) contact line develops from the inlet boundary and propagates toward the outlet. The axisymmetric assumption leads to a wedge shaped domain geometry that is used to describe rotational symmetry in two dimensions. This allows an isolated insight into the interface front cap- and finger-formation mechanisms and associated wall detachment at the droplet front within the pore. The full momentum equation is solved in a one-field approach with spatially varying physical properties with respect to liquid 1 and liquid 2. Mixture values are taken in cells at the interface as described in the volume-of-fluid method. A constant velocity is initialized at the inlet boundary, and the constant pressure boundary condition is imposed at the outlet boundary. On the side walls, a Navier-slip boundary condition is imposed for the velocity field to remove stress singularities originating at the three-phase contact line as proposed by Huh and Scriven27 as
u | r = r p = u w = l λ u r | r = r p ,
(1)
with the velocity u, the radius of the capillary rp, the wall velocity u w, and the slip length l λ, depicted in Fig. 2(c). A constant slip length of l λ = 1 μ m has been used in the final simulation runs. The relevance of a proper boundary condition at the contact line within wetting processes has been highlighted by Gründing et al.28 The implementation of the Navier-slip boundary condition in OpenFoam has been done by Gründing.29 As we have a moving interface over the solid boundary, surface tension is present at the wall boundary. To prevent any mass flux through the wall boundaries, the fluxes are corrected within the pressure correction step of the solver during the velocity pressure coupling. A sharp surface tension force model by Raeini et al.30 is used to reduce spurious currents originating from the interfacial tension implementation at high surface curvatures at the discontinuity of the interface. It allows for the calculation of the capillary pressure pc, which arises from the interface curvature due to the Laplace pressure. The open source software OpenFoam is used with the interFlow31 solver. This solver is based on the Volume Of Fluid solver interFoam with the advantage of a geometric interface reconstruction and sharp interface advection between two incompressible fluids known as IsoAdvector.32 In the context of capillary flows with the introduction of wettability and consequent high interface curvatures in the vicinity of the wall, it is of great importance to keep curvature calculations precise with a low advection error to guarantee bounding of the scalar fields. The multiphase field is represented by an indicator function α, where α = 0 is equivalent to the continuous phase and α = 1 describes the disperse phase. In this interface capturing approach, the interface is represented by α = 0.5.
FIG. 2.

Numerical setup. (a) Numerical domain as an axisymmetric wedge. (b) Initial setup of liquid/liquid displacement in a capillary with geometrical details. (c) Illustration of the slip length lλ based on a flow profile in a channel.

FIG. 2.

Numerical setup. (a) Numerical domain as an axisymmetric wedge. (b) Initial setup of liquid/liquid displacement in a capillary with geometrical details. (c) Illustration of the slip length lλ based on a flow profile in a channel.

Close modal
The apparent dynamic contact angle θa is evaluated as a function of the contact line velocity and the equilibrium contact angle θe. The equilibrium contact angle θe describes the contact angle at rest with a zero velocity at the contact line. The wettability is implemented as a boundary condition with a dynamic contact angle model according to the Kistler correlation33 in OpenFOAM by Nabil and Rattner34 and validated by Sykes et al.35 as
θ a = f H ( C a c l + f H 1 ( θ e ) ) ,
(2)
utilizing the Hoffmann function fH,
f H = arccos [ 1 2 tanh [ 5.16 ( C a c l 1 + 1.31 C a c l 0.99 ) 0.706 ] ] .
(3)
The capillary number at the contact line Cacl is calculated as
C a c l = μ int | U c l | σ ,
(4)
where Ucl is resolved normal to the interface, parallel to the wall in the first cell above the wall. The dynamic viscosity at the interface is described by μint and is calculated as a mixture of the alpha field. The Kistler correlation has been tested against experiments36,37 and has proven to accurately predict spreading and recoiling dynamics. The apparent contact angle θa is defined from the continuous phase; hence, greater contact angles refer to disperse phase (oil) wetting as depicted in Fig. 1(a). The normal vector at the interface is calculated with the alpha field gradient as
n int = ( a ) | ( a ) | .
(5)

A hexahedral mesh with a grid size of 2.9 μ m ( d cap / Δ x 70 ) resulting in 35 cells across the radius of the capillary (rp) is used. Grid independence of the continuous flow profile with its converging maximum velocity in the capillary is achieved at 30 cells per radius as shown in Fig. 3. To evaluate the independence of the mesh from the applied contact angle, equilibrium interfacial positions are evaluated for different grid sizes, shown in Fig. 4(a). Since the meniscus shape and the apparent contact angle does not change significantly above 35 cells per radius, this grid size is chosen. Furthermore, the cap deformation of Taylor droplets from experimental literature38 is compared to the numerical results at a resolution of 35 cells per radius as shown in Fig. 4(b). Helmers et al.38 performed microfluidic experiments for Taylor slug flow in channels to determine the cap deformation of the Taylor droplets with varying capillary numbers and derived a correlation describing the cap deformation factor k c , f with respect to the capillary number Cac (based on continuous phase properties) from it. For these simulations, the numerical domain in Fig. 1 is used. Agreement of the Taylor droplet simulations with the experimental results and the correlation proves the correct behavior of the numerical setup in terms of interfacial deformations and curvature calculations in a multiphase environment.

FIG. 3.

Validation of the momentum transport in the capillary at dimensionless axial position x d cap = 15 , C a d = 0.24 , W e d = 0.7, and Red = 3. (a) Dimensionless velocity over dimensionless radial position, in a single phase flow through the capillary. (b) The deviation of maximum velocity from the analytical solution with respect to on the number of grid cells.

FIG. 3.

Validation of the momentum transport in the capillary at dimensionless axial position x d cap = 15 , C a d = 0.24 , W e d = 0.7, and Red = 3. (a) Dimensionless velocity over dimensionless radial position, in a single phase flow through the capillary. (b) The deviation of maximum velocity from the analytical solution with respect to on the number of grid cells.

Close modal
FIG. 4.

Validation cases for numerical setup. (a) Equilibrium interfacial position, θa = 150°, Cad = 0.24, Wed = 0.7, Red = 3. (b) Front cap deformation ratio with respect to Cac (cont. phase); experimental data from Helmers et al.,38 correlation for Helmers and Taylor droplet simulations performed for this paper: a fixed volume Taylor droplet with disperse phase non wetting condition on the wall, dcap = 200 μm.

FIG. 4.

Validation cases for numerical setup. (a) Equilibrium interfacial position, θa = 150°, Cad = 0.24, Wed = 0.7, Red = 3. (b) Front cap deformation ratio with respect to Cac (cont. phase); experimental data from Helmers et al.,38 correlation for Helmers and Taylor droplet simulations performed for this paper: a fixed volume Taylor droplet with disperse phase non wetting condition on the wall, dcap = 200 μm.

Close modal
The evaluation of shear stress τint at the liquid/liquid interface is based on the work of Wollborn et al.13 and calculated as
τ int = μ int n int D int .
(6)
With the interfacial normal vector n int, the dynamic viscosity at the interface μint. Dint is the velocity gradient tensor, generally described by D = grad ( u ) + grad ( u ) T. For all simulations, the time step Δ t is chosen to satisfy the Courant–Friedrichs–Lewy (CFL) criterion, for a maximum CFL number of 0.5 and time step restrictions to resolve capillary waves39 as follows:
Δ t < Δ t CLF ,
(7)
Δ t < Δ t σ = 0.5 ( ρ c + ρ d ) * Δ x 3 2 π σ ,
(8)
with the continuous and disperse liquid density ρ c , ρ d, respectively, the grid spacing Δ x, and the surface tension σ.
The input parameters for the simulations have been chosen according to a range of dimensionless numbers presented in Table I, to cover cases from stable displacement to unstable displacement, governed by interfacial instabilities. Two-phase flow in capillaries can be described by the interplay of viscous-interfacial- and inertial- forces. Therefore, the following dimensionless numbers are chosen. The capillary number with respect to the disperse phase properties is defined as
C a d = μ d * u ¯ σ .
(9)
TABLE I.

Simulation parameters.

Parameter Range of values
Capillary number (-)  0.04–1 
Weber number (-)  0.02–13.3 
Reynolds number (-)  0.5–13 
Apparent contact angle, θa (°)  0–170 
Capillary diameter, dcap (μm)  200 
Capillary length, l (mm)  6–9 
Viscosity ratio μ c μ d (-)  0.05 
Parameter Range of values
Capillary number (-)  0.04–1 
Weber number (-)  0.02–13.3 
Reynolds number (-)  0.5–13 
Apparent contact angle, θa (°)  0–170 
Capillary diameter, dcap (μm)  200 
Capillary length, l (mm)  6–9 
Viscosity ratio μ c μ d (-)  0.05 
The capillary number accounts for the ratio between viscous forces and interfacial forces. Here, μd describes the disperse phase dynamic viscosity, u ¯ the mean velocity inside the capillary, and σ the interfacial tension between the continuous and disperse phase. The Reynolds number Re describes the ratio between inertial forces to viscous forces and is written as
R e d = ρ d * u ¯ * d cap μ d ,
(10)
with ρd being the disperse phase density and dcap the capillary diameter. The Weber number We defined as
W e d = ρ d * u ¯ 2 * d cap σ
(11)
accounts for the ratio of inertial forces to interfacial forces.
The dimensionless time τ is calculated as
τ = t t ¯ ,
(12)
with the time t and a characteristic reference time t ¯ as
t ¯ = d p u ¯ ,
(13)
with the pore diameter dp and the mean velocity u ¯.
For laminar pipe flows, the hydraulic entrance length xh can be calculated as
x h = 0.05 * R e c * d cap ,
(14)
which describes the length necessary for the transition from a uniform plug flow profile to a Poiseuille-profile.40 Inside the continuous phase, it is calculated with Rec being the Reynolds number of the continuous phase with the continuous phase properties.

As the disperse phase enters the capillary, which is initially filled with the continuous phase, depending on the capillary number, three different scenarios of interfacial propagation are found as illustrated in Fig. 5: (I) If the mean velocity is the same as the contact line (CL) velocity and the dynamic pressure is low, a spherical meniscus is formed, which stably displaces the water ( C a d < 0.1). (II) If the mean velocity is the same as the contact line velocity and the dynamic pressure significantly contributes to the total pressure acting on the interface ( p t = p c + p d), an interface deformation in the form of a cap propagating at a constant velocity is present ( 0.1 < C a d < 0.25). (III) If the mean velocity exceeds the contact line (CL) velocity, the cap undergoes further dynamic deformation and a finger in the center of the capillary is formed, which penetrates the continuous phase ( C a d > 0.25). This finger dynamically grows in the flow direction and may break into dispersed droplets further downstream, the interface contour in Fig. 5 is, therefore, a snapshot. The different meniscus shapes are depicted in Fig. 5. The transition of the scenarios can be related to the balance between viscous and surface tension forces. When the capillary number is low, the interfacial forces exceed the viscous forces; consequently, a meniscus or cap propagating with a constant velocity is present. As the viscous force exceeds the interfacial force, the interface deforms under the influence of viscous stress. Different regimes of the stable liquid/liquid displacement to unstable displacement with the onset of the meniscus instability were experimentally reported in the literature23,24 and related to a critical capillary number. Furthermore, it can be shown from slip-length sensitivity studies that the difference between the CL velocity and the mean pore velocity is the decisive factor. Therefore, in a series of simulation runs, the slip length l λ = 1 μ m has been varied between 0.01 μ m and 10 mm for capillary numbers of C a d = 0.24 0.315 The slip length refers to the extrapolated distance from the wall boundary where the tangential velocity component becomes zero. Increasing the slip length results in a higher CL velocity, finger formation can be artificially prevented for any contact angle within the range of capillary numbers. It has been shown that the balance between the mean pore velocity and the contact line velocity is impacted by the slip length and that this balance is important for the onset of the CL-instabilities. This was shown experimentally in the study by Zhao et al.24 who related the forced wetting transition to this balance. Therefore, the shear stress at the three-phase CL is crucial for the CL propagation and, consequently, for the interface deformation. However, a constant slip length of 1 μ m according to Li et al.41 is chosen in this study.

FIG. 5.

Meniscus shapes for three different interfacial propagation scenarios at τ = 6, θ a = 30 ° from (I) stable meniscus displacement ( C a d = 0.04), (II) stable cap propagation ( C a d = 0.24), and (III) unstable finger formation ( C a d = 0.26).

FIG. 5.

Meniscus shapes for three different interfacial propagation scenarios at τ = 6, θ a = 30 ° from (I) stable meniscus displacement ( C a d = 0.04), (II) stable cap propagation ( C a d = 0.24), and (III) unstable finger formation ( C a d = 0.26).

Close modal

In Fig. 6, the dependence of the capillary number and the apparent contact angle on the maximum shear rate at the three-phase contact line is shown. Overall, the shear rate increases as expected with the velocity and, therefore, the capillary number. Starting from low capillary numbers and low contact angles (continuous phase wetting), the shear rate decreases toward θ a = 90 ° and finds its global minimum at θ a = 90 °. As the meniscus recedes behind the contact line for θ a > 90 °, the shear rate increases. The meniscus is dragged behind the contact line, resulting in additional shear at the contact line. Up to a capillary number C a d < 0.25, the shear rate increases linearly with the velocity for any given contact angle. For contact angles θ a 60 °, the shear rate decreases for capillary numbers from C a d = 0.26 to C a d = 0.32 as a meniscus instability in the shape of a finger is formed. The deformation of the interface and the unstable displacement of the continuous phase refers to the regime (III) unstable finger formation, as described in Fig. 5. As the dynamic pressure bulges the leading interface, the contact line velocity decreases under the mean velocity, therefore resulting in lower shear rates as depicted in Fig. 7. This can be attributed to a decoupling of the contact line from the leading interface. With increasing capillary numbers, a central finger is formed for all contact angles, consequently reducing the slope of the shear rate with the capillary number. Furthermore, the shear rate increases for C a d 0.25 with the increasing contact angle (disperse phase wetting), which can be attributed to higher contact line velocities. Zhao et al.24 reported a decoupling of the dynamics close to the contact line and the fingertip, beyond the wetting transition, which led to slightly lower contact line velocities. The shear stress affects the interface deformation in regions away from the wall and the CL, where microscopic effects are negligible. It could be shown that the CL shear is related to the meniscus position in relation to the CL. With the onset of the meniscus instability, the CL velocity and therefore, the shear rate is reduced.

FIG. 6.

Shear rate at the contact line over capillary number and apparent contact angle, the surface is generated by interpolation of data points.

FIG. 6.

Shear rate at the contact line over capillary number and apparent contact angle, the surface is generated by interpolation of data points.

Close modal
FIG. 7.

Dimensionless contact line velocity for θ a = 0 ° and 30 ° (instability) and θ a = 60 ° 170 ° (no instability), C a d = 0.26 , W e d = 0.86 , R e d = 3.3.

FIG. 7.

Dimensionless contact line velocity for θ a = 0 ° and 30 ° (instability) and θ a = 60 ° 170 ° (no instability), C a d = 0.26 , W e d = 0.86 , R e d = 3.3.

Close modal

The secondary flow field is important for describing the interface deformation in two-phase microcapillary flows in addition to the primary flow field. Figure 8 shows the streamlines of the secondary flow field for C a d = 0.24 and varying θa in a Surface Line Integral Convolution (LIC) representation that has been achieved by subtracting the mean velocity from the velocity field. There is a secondary vortex forming for θ a 90 ° in the disperse phase (oil), in the cap region. In the continuous phase a recirculation flow in front of the interface is present, as a result of the transition from the plug-like velocity profile at the interface to a parabolic Poiseuille flow profile further downstream. The same recirculating flow can be observed in the disperse phase, originating from the parabolic Poiseuille flow profile in the disperse phase to the plug-like profile at the interface. Utilizing Eq. (14), the transition length in this study at Rec = 60, Red = 3, is three times the capillary diameter in the continuous phase and 0.15 times the diameter in the disperse phase. Therefore, most of the transition zones are visible in Fig. 8. The secondary vortex is located inside the drop limited by the concave side of the interface. For contact angles θ a 90 °, the secondary vortex vanishes inside the droplet. For the vortex to form inside the disperse phase, the interface needs to be mobile. As shown by Li et al.,42 who resolved the flow around bubbles with different surface mobility, an immobile bubble surface results in vortices in the wake of the bubble outside the bubble. These structures could not be observed with the mobile bubble surface. With reference to Mießner et al.,43 it can be discussed that these vortices located in the wake of the bubble for immobile surfaces move into the bubble for mobile surfaces. Therefore, the reduced mobility of the interface in the case of disperse phase wetting in this study affects the existence of the secondary vortices within the displacing phase. For θ a 90 °, secondary vortices in front of the interface originate from the contact line become prominent. This vortex at high local curvatures in the vicinity of the wall has been also reported by Sui and Spelt.44 Furthermore, there is a stagnation point for θ a 90 ° in the vicinity of the interface between the recirculating flow and the secondary vortex (highlighted). The stagnation point moves further into the center of the pore as the contact angle increases. In the context of liquid film deposition Mayer and Krechetnikov45 reported the film thickening phenomenon of surfactant–laden interfaces with an interior stagnation point. Also, Mießner et al.46 reported an annular stagnation region close to the walls and define it as the region where droplet deformation occurs. Consequently increased liquid film deposition with increased disperse phase wetting (greater contact angles) can be attributed to this stagnation point, as the interface bulges at this stagnation point at sufficiently high capillary numbers, creating a liquid film of continuous phase at the wall, as depicted in Fig. 9. The liquid film thickness of continuous phase at the wall is approximately 1/5 of the channel width if the disperse phase wets the wall ( θ a = 150 °), while it is approximately 1/8 of the channel width in the inverse case. Vécsei et al.19 found that interfacial instabilities of liquid films coating the walls are highly affected by the wall film thickness when sheared with a gas flow. For small film thicknesses ( 1 / 8 of the channel width), the interface is mainly deformed by viscous stress, while for thicker films, the instabilities are driven by pressure gradients. The instabilities due to viscous stresses are driven by the counter-rotating vortices, whereas the pressure driven instability occurs due to the widening of the film layer and an increase in the dynamic pressure analogous to the Venturi effect as shown by Vécsei et al.19 

FIG. 8.

Radial velocity and streamlines of secondary flow in LIC representation for varying contact angles without finger formation, the location of the interface is highlighted in white, the stagnation point between the recirculating flow and the secondary vortex highlighted with the orange circle, τ = 15, C a d = 0.24 , W e d = 0.7, Red = 3. (a) θa = 30°, (b) θa = 60°, (c) θa = 90°, (d) θa = 120°, (e) θa = 150°, and (f) θa = 170°.

FIG. 8.

Radial velocity and streamlines of secondary flow in LIC representation for varying contact angles without finger formation, the location of the interface is highlighted in white, the stagnation point between the recirculating flow and the secondary vortex highlighted with the orange circle, τ = 15, C a d = 0.24 , W e d = 0.7, Red = 3. (a) θa = 30°, (b) θa = 60°, (c) θa = 90°, (d) θa = 120°, (e) θa = 150°, and (f) θa = 170°.

Close modal
FIG. 9.

Snapshot of the interface contour for four contact angles after the formation of the meniscus instability, the continuous phase liquid film deposition is highlighted with arrows, τ = 6, C a d = 0.47 , W e d = 2.8, Red = 6.

FIG. 9.

Snapshot of the interface contour for four contact angles after the formation of the meniscus instability, the continuous phase liquid film deposition is highlighted with arrows, τ = 6, C a d = 0.47 , W e d = 2.8, Red = 6.

Close modal

With respect to these results the wettability is suspect to change the nature of the interface instabilities from shear to pressure induced that eventually breakup the interface. Due to the different radial locations of the bulged interface within the capillary, the interfacial shear is, therefore, dependent on the wettability, resulting in different shear histories for wetting and non-wetting systems. However, this only applies to a certain range of capillary numbers, since the dynamic contact angle approaches full wetting at Ca = 1 according to Kistler,33 therefore eliminating any dependence of the final film height with the equilibrium contact angle. As the capillary number is further increased into the regime of unstable displacement, the velocity magnitude in the region of the cap, normal to the interface exceeds the mean velocity, and the secondary vortex vanishes with the interfacial growth. The origin of the secondary vortex is defined by the mean velocity of the flow field. When the interface bulges due to the interfacial growth, the velocity in the cap region exceeds the mean velocity and the vortex vanishes. Experimental observations for this transition are given by Vagner et al.47 For growing microdroplets during co-flow, they used micro-particle image velocimetry ( μ PIV ) and showed that vortices inside the droplet vanish, if the droplet growth velocity dominates over the velocity of the continuous phase in the gap between the wall and the droplet.

Figure 10(a) shows the distribution of the total pressure p and the capillary pressure pc that arises from the interfacial curvature along the centerline of the capillary for two different contact angles, for a time step where no meniscus instability is present. Slightly lower pressures at the entrance for the disperse phase wetting system can be noted and is in line with classic lubrication theory. From the Young–Laplace equation,48 the pressure difference required for a liquid to enter the capillary (breakthrough pressure) decreases with the decreasing contact angle of the displacing (disperse) liquid. The total pressure decreases linearly with the same gradient for both systems, until it reaches the area limited by the interface. For high interfacial curvatures ( θ a < 90 °), the total pressure rises toward the interface, while for lower interfacial curvatures ( θ a 90 °), the gradient flattens toward the interface. This is followed by the pressure jump across the interface (Laplace pressure) and subsequent linear decrease in pressure due to viscous dissipation in the continuous phase. The dashed lines indicate that the observed increasing pressure arises from the capillary pressure and, therefore, the interfacial curvature. The increase in the capillary pressure for θ a < 90 ° under interfacial deformation can be attributed to the effect of the dynamic capillary pressure caused by the viscous forces exerted by the secondary flow vortices. The secondary vortices create a stagnation point on the interface (at the centerline) as described by Mießner et al.46 for Taylor droplets, which contributes to the increase in total pressure toward the interface by viscous displacement. As the droplet propagates into the pore eventually, a meniscus instability is formed for θ a = 30 °. The corresponding pressure plot is depicted in Fig. 10(b). As the droplet further penetrates into the capillary, the back pressure at the inlet increases. When the interface bulges into a finger the back pressure decreases, resulting in lower pressure at the entrance of the pore for θ a = 30 ° compared to θ a = 150 ° where no instability is present. Inside the bulged interface (highlighted), the pressure gradient reduces to the level of the viscous dissipation inside the continuous phase in front of the droplet. Due to the finger formation and the decoupling of the contact line, the contact line velocity decreases; thus, the shear induced pressure in Fig. 10(b) at x / l = 0 decreases. When the interface bulges into a finger, the curvature of the interface on the centerline decreases; therefore, the increase in pressure toward the interface decreases, as depicted in Fig. 10(a) vs Fig. 10(b).

FIG. 10.

Center line pressure and capillary pressure for θ a = 30 ° and θ a = 150 °, (a) τ = 3 and (b) τ = 9, C a d = 0.26 , W e d = 0.86 , R e d = 3.3. (a) τ = 3. (b) τ = 9, meniscus instability highlighted.

FIG. 10.

Center line pressure and capillary pressure for θ a = 30 ° and θ a = 150 °, (a) τ = 3 and (b) τ = 9, C a d = 0.26 , W e d = 0.86 , R e d = 3.3. (a) τ = 3. (b) τ = 9, meniscus instability highlighted.

Close modal

To further evaluate the influence of the secondary vortices on the droplet deformation, a case with zero continuous phase viscosity μ c = 0 was conducted, hence keeping C a d = const . and eliminating fluid/fluid shear at the interface. The simulation results indicate that no secondary vortex in the cap region exists [Fig. 11(c)]. Consequently, the secondary vortex limited by the concave side of the interface is driven by fluid/fluid shear, thus driven by the continuous phase. This is, in line, with the findings of Mießner et al.46 who stated that the momentum coupling across the interface transfers kinetic energy into the droplet phase. Furthermore, the interfacial deformation is reduced significantly resulting in lower capillary pressures [Fig. 11(a)]. The overall lower pressure in this case is attributed to the missing viscous pressure drop in the continuous phase. In this case, no instability is formed, emphasizing the importance of fluid/fluid shear at the interface in this context and the transition from capillary- to viscosity-dominated flow dynamics. Therefore, the meniscus instability can be attributed to the co-existence of the recirculating flow in the continuous phase and the secondary vortex in the disperse phase. If there is no secondary vortex, as for instance depicted in Fig. 11(c), there is no secondary flow and thus no stagnation point. Consequently, no meniscus instability is present (proven at capillary numbers up to Cad = 1).

FIG. 11.

Comparison of standard case and zero continuous viscosity case for stable cap displacement, τ = 3, C a d = 0.26 , W e d = 0.86 , R e d = 3.3. (a) Centerline pressure and capillary pressure for θa = 30° and θa = 150°. (b) Interface contour for θa = 30° for standard case and μc = 0. (c) Radial velocity and streamlines of secondary flow in LIC representation for θa = 30°, μc = 0 with no secondary vortex in the cap region.

FIG. 11.

Comparison of standard case and zero continuous viscosity case for stable cap displacement, τ = 3, C a d = 0.26 , W e d = 0.86 , R e d = 3.3. (a) Centerline pressure and capillary pressure for θa = 30° and θa = 150°. (b) Interface contour for θa = 30° for standard case and μc = 0. (c) Radial velocity and streamlines of secondary flow in LIC representation for θa = 30°, μc = 0 with no secondary vortex in the cap region.

Close modal

The relationship between shear and capillary pressure can be described as follows: (I) The interface is sheared by fluid/fluid interactions; (II) A secondary vortex arises, limited by the concave side of the interface; (III) The viscous displacement results in a stagnation point at the droplet front; (IV) The normal force exerted at the stagnation point on the centerline deforms the interface; and (V) The capillary pressure rises. It can be discussed that the increased pressure gradient in the cap of the droplet further accelerates the secondary vortex, hence accelerating the occurrence of the meniscus instability. Comparing the contour of the interface over time for two contact angles as shown in Fig. 12, it is evident that the interface undergoes axial movement in the disperse phase non-wetting case [Fig. 12(a)], whereas the position of the interface in the disperse phase wetting case [Fig. 12(b)] is pinned to its location. This axial movement follows a dampened oscillation, originating from an “overshooting” of the interfacial position during the entrance in the pore; therefore, high local curvatures occur. The overshooting phenomenon describes an initial increase in the degree of interface deformation due to oscillatory motion during the disperse (second) phase entering the capillary. As the contact line forms at the inlet, the displacing phase may be pinned to the entrance until the breakthrough pressure is reached, resulting in a temporal overshoot and interface oscillation. The overshooting phenomenon contributes to the earlier transition from the stable interface translation state to the capillary-fingering regime for non-wetting displacing fluids. Furthermore, the presence of a meniscus bends streamlines close to the interface by capillary force, which results in increased shear and viscous dissipation in the continuous phase film respective to the wedge between the interface and the wall, compared to Poiseuille flow.49 Since the capillary pressure is proportional to the interface curvature, which is greater for smaller advancing contact angles as shown in this study, the viscous dissipation in the case for smaller contact angles is greater in the wedge. The higher velocity associated with greater viscous dissipation in the continuous phases transfers kinetic energy via momentum coupling over the interface into the disperse phase. This accelerates the secondary vortices and hence increases the pressure at the stagnation point at the centerline. Capillary oscillations have been reported by Sui and Spelt44 at significantly higher Reynolds numbers, they argue, however, that realistically small slip length values would shift the oscillatory regime to lower Reynolds numbers, which applies to this study. However, while we could observe oscillations of the meniscus, no oscillations at the contact line were present, which is fundamentally different. These oscillations correspond to a change in the two-dimensional origin of the secondary vortex with a steady contact line velocity. These oscillations were confirmed with a parabolic Poiseuille flow profile at the inlet, as used in the work of Sui and Spelt.44 

FIG. 12.

Interface contour over time for hydrophilic and hydrophobic contact angles, C a d = 0.24 , W e d = 0.71, Red = 3: (a) θa = 30° and (b) θa = 150°.

FIG. 12.

Interface contour over time for hydrophilic and hydrophobic contact angles, C a d = 0.24 , W e d = 0.71, Red = 3: (a) θa = 30° and (b) θa = 150°.

Close modal

We have investigated the effect of wettability in a microcapillary tube in a liquid/liquid system where a fluid is displacing another fluid. Three regimes of interfacial propagation are identified: (I) stable meniscus displacement, (II) stable cap displacement, and (III) unstable fingering and breakup. The regime transitions are impacted by the chosen slip length and, therefore, dependent on the balance between the contact line velocity and the mean pore velocity. These regimes can be assigned to the balance of viscous forces to interfacial forces. If the interfacial forces dominate, stable displacement occurs. If the viscous forces exceed the interfacial forces, unstable displacement occurs and the contact line decouples from the leading interface. The wettability smears the regime transitions, resulting in earlier unstable fingering for non-wetting displacing fluids. The dependence of the wettability on the regime transitions is in line with experimental data in the literature.24 In particular, the formation of secondary vortices driven by fluid/fluid shear is greatly dependent on the wettability and affects the interfacial deformation due to viscous displacement, thereby the capillary pressure. The capillary pressure rises toward the curved meniscus and can be seen as the driving force of capillary fingering. Furthermore, the displaced liquid film deposition strongly depends on the wettability and can be assigned to a stagnation point in the secondary flow field. For non-wetting displacing fluids oscillations of the meniscus are observed, originating from an overshooting of the meniscus at the entrance of the capillary and not from contact line oscillations as reported in the literature for higher Reynolds numbers. Lower shear rates at the liquid/liquid interface for wetting displacing fluids can be tuned by the wettability, thereby affecting the radial position of the interface in the capillary. This understanding of the interface propagation in microcapillary liquid/liquid flow is fundamental for proper handling of shear-sensitive media like interfacial active proteins or surfactants for instance in emulsification processes.

This project has been funded by the German Research Foundation (DFG) within the SPP 1934 “DiSPBiotech” under Grant No. FR 912/40. Computational resources were provided by The North-German Supercomputing Alliance (HLRN) under Grant No. hbi00036.

The authors have no conflicts to disclose.

Patrick Giefer: Conceptualization (lead); Investigation (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (equal). Apostolos Kyrloglou: Software (equal); Writing – review & editing (equal). Udo Fritsching: Funding acquisition (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.7551242, Ref. 50.

1.
T.
Wollborn
,
M.
Michaelis
,
L. C.
Ciacchi
, and
U.
Fritsching
, “
Protein conformational changes at the oil/water-interface induced by premix membrane emulsification
,”
J. Colloid Interface Sci.
628
,
72
81
(
2022
).
2.
K.
Suzuki
,
K.
Hayakawa
, and
Y.
Hagura
, “
Preparation of high concentration O/W and W/O emulsions by the membrane phase inversion emulsification using PTFE membranes
,”
Food Sci. Technol. Res.
5
,
234
238
(
1999
).
3.
J.
Wu
,
W.
Jing
,
W.
Xing
, and
N.
Xu
, “
Preparation of W/O emulsions by membrane emulsification with a mullite ceramic membrane
,”
Desalination
193
,
381
386
(
2006
).
4.
G.
Vladisavljevic
, “
Influence of process parameters on droplet size distribution in SPG membrane emulsification and stability of prepared emulsion droplets
,”
J. Membr. Sci.
225
,
15
23
(
2003
).
5.
I.
Kobayashi
,
M.
Yasuno
,
S.
Iwamoto
,
A.
Shono
,
K.
Satoh
, and
M.
Nakajima
, “
Microscopic observation of emulsion droplet formation from a polycarbonate membrane
,”
Colloids Surf., A
207
,
185
196
(
2002
).
6.
S.
Omi
, “
Preparation of monodisperse microspheres using the Shirasu porous glass emulsification technique
,”
Colloids Surf., A
109
,
97
107
(
1996
).
7.
K.
Kandori
,
K.
Kishi
, and
T.
Ishikawa
, “
Formation mechanisms of monodispersed W/O emulsions by SPG filter emulsification method
,”
Colloids Surf.
61
,
269
279
(
1991
).
8.
Z.
Liu
,
M.
Chai
,
X.
Chen
,
S. H.
Hejazi
, and
Y.
Li
, “
Emulsification in a microfluidic flow-focusing device: Effect of the dispersed phase viscosity
,”
Fuel
283
,
119229
(
2021
).
9.
G.
Muschiolik
,
S.
Dräger
,
I.
Scherze
,
H. M.
Rawel
, and
M.
Stang
, “
Protein-stabilized emulsions prepared by the micro-porous glass method
,” in
Food Colloids
(
Elsevier
,
2004
), pp.
393
400
.
10.
K.
van Dijke
,
I.
Kobayashi
,
K.
Schroën
,
K.
Uemura
,
M.
Nakajima
, and
R.
Boom
, “
Effect of viscosities of dispersed and continuous phases in microchannel oil-in-water emulsification
,”
Microfluid. Nanofluid.
9
,
77
85
(
2010
).
11.
K.
Kandori
,
K.
Kishi
, and
T.
Ishikawa
, “
Preparation of monodispersed W/O emulsions by Shirasu-porous-glass filter emulsification technique
,”
Colloids Surf.
55
,
73
78
(
1991
).
12.
L.
Luhede
,
T.
Wollborn
, and
U.
Fritsching
, “
Stability of multiple emulsions under shear stress
,”
Can. J. Chem. Eng.
98
,
186
193
(
2020
).
13.
T.
Wollborn
,
L.
Luhede
, and
U.
Fritsching
, “
Evaluating interfacial shear and strain stress during droplet deformation in micro-pores
,”
Phys. Fluids
31
,
012109
(
2019
).
14.
T.
Wollborn
,
P.
Giefer
,
H.
Kieserling
,
A. M.
Wagemans
,
S.
Drusch
, and
U.
Fritsching
, “
Investigation of local and temporal interfacial shear stress distribution during membrane emulsification
,”
Can. J. Chem. Eng.
100
,
1061
1078
(
2022
).
15.
K.
Schroën
,
M.
Ferrando
,
S.
de Lamo-Castellví
,
S.
Sahin
, and
C.
Güell
, “
Linking findings in microfluidics to membrane emulsification process design: The importance of wettability and component interactions with interfaces
,”
Membranes
6
,
26
(
2016
).
16.
R. G.
Cox
, “
The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow
,”
J. Fluid Mech.
168
,
169
(
1986
).
17.
O. V.
Voinov
, “
Inclination angles of the boundary in moving liquid layers
,”
J. Appl. Mech. Tech. Phys.
18
,
216
222
(
1977
).
18.
X.
Chen
,
E.
Ramé
, and
S.
Garoff
, “
The effects of thin and ultrathin liquid films on dynamic wetting
,”
Phys. Fluids
16
,
287
297
(
2004
).
19.
M.
Vécsei
,
M.
Dietzel
, and
S.
Hardt
, “
Interfacial instability of liquid films coating the walls of a parallel-plate channel and sheared by a gas flow
,”
Microfluid. Nanofluid.
22
,
91
(
2018
).
20.
D.
Bonn
,
J.
Eggers
,
J.
Indekeu
,
J.
Meunier
, and
E.
Rolley
, “
Wetting and spreading
,”
Rev. Mod. Phys.
81
,
739
805
(
2009
).
21.
R. L.
Hoffman
, “
A study of the advancing interface. I. Interface shape in liquid–gas systems
,”
J. Colloid Interface Sci.
50
,
228
241
(
1975
).
22.
R. L.
Hoffman
, “
A study of the advancing interface
,”
J. Colloid Interface Sci.
94
,
470
486
(
1983
).
23.
M.
Fermigier
and
P.
Jenffer
, “
An experimental investigation of the dynamic contact angle in liquid-liquid systems
,”
J. Colloid Interface Sci.
146
,
226
241
(
1991
).
24.
B.
Zhao
,
A. A.
Pahlavan
,
L.
Cueto-Felgueroso
, and
R.
Juanes
, “
Forced wetting transition and bubble pinch-off in a capillary tube
,”
Phys. Rev. Lett.
120
,
084501
(
2018
).
25.
A. A.
Pahlavan
,
H. A.
Stone
,
G. H.
McKinley
, and
R.
Juanes
, “
Restoring universality to the pinch-off of a bubble
,”
Proc. Natl. Acad. Sci.
116
,
13780
13784
(
2019
).
26.
C.-E.
Wu
,
H.-W.
Du
,
J.
Qin
,
E.-Q.
Li
, and
P.
Gao
, “
Experiment on bubble formation through dynamical wetting transition in a square capillary
,”
AIP Adv.
11
,
075107
(
2021
).
27.
C.
Huh
and
L.
Scriven
, “
Hydrodynamic model of steady movement of a solid/liquid/fluid contact line
,”
J. Colloid Interface Sci.
35
,
85
101
(
1971
).
28.
D.
Gründing
,
M.
Smuda
,
T.
Antritter
,
M.
Fricke
,
D.
Rettenmaier
,
F.
Kummer
,
P.
Stephan
,
H.
Marschall
, and
D.
Bothe
, “
A comparative study of transient capillary rise using direct numerical simulations
,”
Appl. Math. Modell.
86
,
142
165
(
2020
).
29.
D.
Gründing
, “
An arbitrary Lagrangian-Eulerian method for the direct numerical simulation of wetting processes
,” Ph.D. thesis (
Technische Universität Darmstadt
,
2020
).
30.
A. Q.
Raeini
,
M. J.
Blunt
, and
B.
Bijeljic
, “
Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method
,”
J. Comput. Phys.
231
,
5653
5668
(
2012
).
31.
J.
Roenby
,
H.
Bredmose
, and
H.
Jasak
, “
A computational method for sharp interface advection
,”
R. Soc. Open Sci.
3
,
160405
(
2016
).
32.
J.
Roenby
,
H.
Bredmose
, and
H.
Jasak
, “
IsoAdvector: Geometric VOF on general meshes
,” in
OpenFOAM®
(
Springer International Publishing
,
2019
), pp.
281
296
.
33.
S.
Kistler
, “
Hydrodynamics of wetting
,” in
Wettability
(
Taylor & Francis
,
1993
).
34.
M.
Nabil
and
A. S.
Rattner
, “
interThermalPhaseChangeFoam—A framework for two-phase flow simulations with thermally driven phase change
,”
SoftwareX
5
,
216
226
(
2016
).
35.
T. C.
Sykes
,
D.
Harbottle
,
Z.
Khatir
,
H. M.
Thompson
, and
M. C. T.
Wilson
, “
Substrate wettability influences internal jet formation and mixing during droplet coalescence
,”
Langmuir
36
,
9596
9607
(
2020
).
36.
K.
Vontas
,
C.
Boscariol
,
M.
Andredaki
,
A.
Georgoulas
,
C.
Crua
,
J. H.
Walther
, and
M.
Marengo
, “
Droplet impact on suspended metallic meshes: Effects of wettability, Reynolds and Weber numbers
,”
Fluids
5
,
81
(
2020
).
37.
J.
Göhl
,
A.
Mark
,
S.
Sasic
, and
F.
Edelvik
, “
An immersed boundary based dynamic contact angle framework for handling complex surfaces of mixed wettabilities
,”
Int. J. Multiphase Flow
109
,
164
177
(
2018
).
38.
T.
Helmers
,
P.
Kemper
,
J.
Thöming
, and
U.
Mießner
, “
Determining the flow-related cap deformation of Taylor droplets at low Ca numbers using ensemble-averaged high-speed images
,”
Exp. Fluids
60
,
113
(
2019
).
39.
J.
Brackbill
,
D.
Kothe
, and
C.
Zemach
, “
A continuum method for modeling surface tension
,”
J. Comput. Phys.
100
,
335
354
(
1992
).
40.
T. L.
Bergman
,
Fundamentals of Heat and Mass Transfer
(
Wiley
,
2011
), p.
1048
.
41.
X.
Li
,
L.
He
,
S.
Lv
,
C.
Xu
,
P.
Qian
,
F.
Xie
, and
M.
Liu
, “
Effects of wall velocity slip on droplet generation in microfluidic T-junctions
,”
RSC Adv.
9
,
23229
23240
(
2019
).
42.
S.
Li
,
K.
Jue
, and
C.
Sun
, “
Effect of bubble surface properties on bubble–particle collision efficiency in froth flotation
,”
Minerals
10
,
367
(
2020
).
43.
U.
Mießner
,
T.
Helmers
,
R.
Lindken
, and
J.
Westerweel
, “
Reconstruction of the 3D pressure field and energy dissipation of a Taylor droplet from a μPIV measurement
,”
Exp. Fluids
62
,
83
(
2021
).
44.
Y.
Sui
and
P. D. M.
Spelt
, “
Sustained inertial-capillary oscillations and jet formation in displacement flow in a tube
,”
Phys. Fluids
23
,
122104
(
2011
).
45.
H. C.
Mayer
and
R.
Krechetnikov
, “
Landau-Levich flow visualization: Revealing the flow topology responsible for the film thickening phenomena
,”
Phys. Fluids
24
,
052103
(
2012
).
46.
U.
Mießner
,
T.
Helmers
,
R.
Lindken
, and
J.
Westerweel
, “
μPIV measurement of the 3D velocity distribution of Taylor droplets moving in a square horizontal channel
,”
Exp. Fluids
61
,
125
(
2020
).
47.
S. A.
Vagner
,
S. A.
Patlazhan
,
C. A.
Serra
, and
D.
Funfschilling
, “
Vortex flow evolution in a growing microdroplet during co-flow in coaxial capillaries
,”
Phys. Fluids
33
,
072010
(
2021
).
48.
T.
Young
, “
III. An essay on the cohesion of fluids
,”
Philos. Trans. R. Soc. London
95
,
65
87
(
1805
).
49.
I.
Lunati
and
D.
Or
, “
Gravity-driven slug motion in capillary tubes
,”
Phys. Fluids
21
,
052003
(
2009
).
50.
P.
Giefer
,
A.
Kyrloglou
, and
U.
Fritsching
, “Leibniz-IWT/capillaryWettingTransition: Initial Release,”
Zenodo
.