Contraction of a shear-thinning axisymmetric cavity

The evolution of small holes nucleated in a thin liquid sheet is critical to the integrity of bubble films and therefore important to environmental processes, such as mist formation from bursting sea bubbles, and technological processes in the food and chemical industries associated with the stability of foams and cleaning microbubbles (Laporte et al., 2016; Drenckhan and Saint-Jalmes, 2015), the stability of the falling film in curtain coating processes (Miyamoto and Katagiri, 1997; Mohammad Karim et al., 2018), and in liquid-sheet-based spray formation techniques (Villermaux, 2007; Altieri et al., 2014). Recently, the dynamics of small holes has also become relevant to nanotechnologies such as the fabrication of nanopore-based biosensors for rapid protein and genetic studies (Storm et al., 2003; 2005). The evolution of a hole nucleated in a liquid sheet is highly dependent on the initial size of the hole. Small holes (in relation to the sheet thickness) contract, driven by capillary forces, eventually restoring the integrity of the liquid phase (Taylor and Michael, 1973). Conversely, large holes expand to minimize the free energy, ultimately fragmenting the liquid sheet (Bremond and Villermaux, 2006; Lhuissier and Villermaux, 2013).

The dynamics of both large and small holes are currently well understood for holes nucleated in Newtonian liquid sheets. The dynamics of large Newtonian holes have been studied in the pioneering works of Taylor (1959) and Culick (1960). These studies demonstrated that a large hole expanding in a low-viscosity liquid sheet eventually attains a constant terminal velocity, which is currently known as the Taylor-Culick velocity. More recently, Savva and Bush (2009) extended these studies to holes expanding in a highly viscous liquid sheet. They found that a viscous hole expanding under the influence of capillarity also attains a terminal Taylor-Culick velocity after an initial transient behavior. Similarly, the dynamics of small Newtonian holes have also been studied considering both the inertial and viscous regimes. Small contracting holes in the inertial regime have been found to collapse at a velocity that increases as the hole contracts, following an inverse power law relationship with hole size (Lu et al., 2015). Conversely, holes contracting in the viscous regime have been found to collapse with a constant terminal velocity (Storm et al., 2003; Wu et al., 2006; and Asghar et al., 2011), which results from the equilibrium of viscous and surface tension forces (Lu et al., 2015). However, the dynamics of non-Newtonian liquid holes is much less understood due to the added rheological complexity of a varying viscosity with the flow near the cavity. Motivated by the ubiquity of complex fluids in natural and industrial processes, the dynamics of holes in non-Newtonian liquid sheets have been studied for the case of large expanding holes (Dalnoki-Veress et al., 1999; Roth et al., 2005; and Deyrail et al., 2007). However, a clear understanding of the dynamics of small contracting holes in non-Newtonian liquid sheets is still lacking.

Here, as a first step toward the analysis of non-Newtonian contracting holes, we study the evolution of a small hole in a highly viscous, shear-thinning liquid sheet. We use scaling arguments supported by high-fidelity simulations to study how the shear-dependent viscosity modifies the contraction dynamics, assuming that the shear-thinning properties of the sheet can be described by a Carreau-Yasuda constitutive equation. The results from the simulations demonstrate that shear-thinning effects accelerate the hole collapse due to a progressive decrease in the viscosity near the hole’s leading front where the hoop stress keeps increasing as the hole contracts. In addition, scaling arguments enable us to develop a power law to describe the evolution of the hole radius in terms of the shear-thinning properties.

In this work, the dynamics of hole contraction is studied by following the evolution of a small toroidal hole with an initial radius $r^0$ in a shear-thinning liquid sheet of density ρ, surface tension σ, and large zero-shear viscosity $μ^0$ [Fig. 1(a)]. The model is made dimensionless using the unperturbed sheet thickness ĥ as a length scale, the capillary time $μ^0ĥ/σ$ as a time scale, and the capillary stress σ/ĥ as a stress scale.

The contraction of the viscous toroidal hole is analyzed by solving the dimensionless continuity,

and Stokes equations,

for liquid velocity v and pressure p, where T is the dimensionless Cauchy stress tensor

where I is the identity tensor, D is the rate of strain tensor, and μ is a dimensionless viscosity to be defined shortly. Inertia forces are considered negligible in relation to viscous and surface tension forces (i.e., the Ohnesorge number $Oh≡μ^0/2ρσĥ≫1$⁠), and the surrounding gas is considered dynamically passive. Note that an alternative characteristic viscosity based on the terminal viscous-capillary velocity of the pore could also be adopted, as proposed by Thompson and Soares (2016). The shear-dependent properties of the liquid sheet are assumed to be well described by a Carreau-Yasuda constitutive equation (Barnes et al., 1989),

where Λ is the dimensionless Carreau time constant and the viscosity μ (in units of $μ^0$⁠) varies between a dimensionless zero-shear plateau viscosity μ = 1 at low strain rates $γ̇≪Λ−1$ and a dimensionless infinite-shear plateau viscosity μ = M at high strain rates $γ̇≫Λ−1M−1/(1−n)$⁠. Between these limits, the viscosity decreases as the strain rate $γ̇$ increases following a power-law with flow index n. As our work focuses on the influence of the flow index in the dynamics, the system is considered to have small infinite-shear viscosity M = 10⁻³ and moderate dimensionless time constant Λ = 10.

The gas-liquid free interface is a material surface; therefore, mass conservation is ensured by imposing the kinematic boundary condition

where $ẋ$ is the velocity of the interface and n the unit normal vector to the interface. Neglecting the dynamical effect of the outer gas, the balance of stresses at the interface is imposed through the traction boundary condition

where κ = −∇_s · n is the curvature of the interface and ∇_s ≡ ∇ − n(n · ∇) is the surface gradient operator. Computations start with a quiescent fluid and a cavity of small initial radius r₀. No penetration and vanishing tangential stress conditions are imposed along the plane of symmetry at z = 0.

The full set of axisymmetric Navier-Stokes system of equations (1) and (2) and associated boundary conditions (5) and (6) are solved simultaneously in the liquid phase for the velocity, pressure, and location of the free surface using the finite element method with an adaptive mesh and implicit time integration. To account for the movement of the free surface, the spatial derivatives were discretized using the arbitrary Lagrangian-Eulerian method of spines developed by Kistler and Scriven (1983), as extended by Xue et al. (2008) to include inelastic non-Newtonian effects. The time derivatives were discretized using a second-order trapezoidal difference method with an Adam-Bashforth predictor (Gresho et al., 1980) smoothing the first few transient solutions using a backward difference method as proposed by Kistler and Scriven (1983). The second-order Adam-Bashforth predictor was used with a trapezoidal rule to control time truncation errors, and the time steps for the time integration were adaptively chosen using a first-order continuation method (Corvalan and Saita, 1991). The domain was tessellated in a structured mesh of quadrilateral Taylor-Hood elements (Taylor and Hood, 1973), with velocity components being approximated by biquadratic shape functions and pressure by bilinear shape functions. Studies with meshes at various resolutions were carried out and meshes of $O(104)$ degrees of freedom were selected for the simulations depending on flow index. A mesh independence study was carried out until the scaling laws reported in the results remain essentially unchanged. The elements were nonuniformly spaced with higher concentration near the tip of the hole in the radial direction, and in the vicinity of the free surface in the axial direction, as exemplified in Fig. 1(b).

Here, we summarize our main results on the contraction of a small axisymmetric cavity or hole inside a generalized Newtonian liquid. We focus on the influence of shear-thinning on the hole dynamics, as characterized by the flow index n. As our focus is on the effect of the shear-dependent viscosity, it is instructive to start the discussion with the dynamics of purely Newtonian holes (n = 1) as the baseline.

We first consider the results from the work of Savva and Bush (2009) for an expanding Newtonian hole. The dynamics of expanding Newtonian holes was characterized independently by Taylor (1959) and Culick (1960) for inertial fluids, and more recently, Savva and Bush (2009) developed a theoretical and numerical lubrication model including the influence of the fluid viscosity. Figure 2 shows the results from the work of Savva and Bush (2009) for an expanding hole with large initial size r₀ = 50 (in units of the unperturbed sheet thickness) in a moderately viscous fluid sheet with Oh = 10. This figure confirms that as the hole expands driven by capillary forces, the hole velocity grows rapidly but then slows down and approaches the terminal Taylor-Culick velocity û ≡ (2σ/(ρĥ))^1/2 expected on the basis of momentum conservation. In addition, this figure shows excellent agreement between the results from the lubrication approximation (symbols) and the solution of the full Navier-Stokes system (solid line) for both the hole velocity [Fig. 2(a)] and a typical interfacial shape [Fig. 2(b)].

Similarly, the speed of contracting Newtonian holes also become constant at later times provided that the viscosity of the liquid sheet is sufficiently large. A constant velocity of contraction has been observed in highly viscous holes during experiments motivated by the fabrication of nanopore based sensors for DNA sequencing (Storm et al., 2003; 2005; and Wu et al., 2006) and during the collapse of the cavity formed during bubble pinch-off (Thoroddsen et al., 2007; Burton et al., 2005). Figure 3 shows simulations results for a Newtonian contracting hole with small initial size r₀ = 0.2 in a purely viscous fluid sheet (blue line). This figure confirms that after an initial transient, the contracting hole attains a terminal viscous-capillary velocity $ŵ=σ/(2μ^0)$ (or $w$ = 1/2 in dimensionless terms), which emerges from the eventual equilibrium of viscous and surface tension forces (Lu et al., 2015; Burton et al., 2005).

The terminal viscous-capillary velocity $ŵ=σ/(2μ^0)$ for viscous contracting holes is a function of the fluid viscosity and therefore likely to be influenced by potential shear-thinning effects. In addition, shear-thinning effects are expected to become important in collapsing cavities because viscous stresses grow rapidly near the cavity’s leading front as the minimum hole radius r_m → 0 (Lu et al., 2015; Bolanos-Jiménez et al., 2009).

To gain preliminary insight into the influence of shear-thinning on the hole dynamics, Fig. 3 compares the dynamics of the viscous Newtonian hole to that of the corresponding shear-thinning hole with a small Carreau flow index n = 0.7. The results show that the behavior of the shear-thinning hole (solid black line) is markedly different from that of the Newtonian one (blue line). Critically, the advancing front of the shear-thinning cavity no longer shows a constant terminal velocity. Instead, the velocity of contraction $v$_m grows following a power-law relationship with the hole size

with α ≈ −3/7 (dashed line), indicating that shear thinning is a critical feature modulating the contraction rate of non-Newtonian holes (Jiang et al., 2017).

Despite its important influence on the dynamics, we observe that shear-thinning is localized in a neighborhood of the cavity’s leading front. We calculate the local viscosity in the liquid sheet to identify how the shear-thinning effect contributes to the flow dynamics. Figure 4 illustrates in detail the evolution of the cross-sectional viscosity field for the shear-thinning hole of Fig. 3. The results show that the overall shear-thinning influence extends a radial distance of the order of the film thickness and that the dominant shear-thinning effect is highly localized in a small neighborhood of the hole tip. Thus, shear-thinning accelerates the hole collapse due to the progressive decrease in the shear-dependent viscosity in the vicinity of the tip of the pore.

A simple scaling argument leads to a relationship that supports the value of the power-law exponent estimated in Fig. 3 and predicts the exponent α for other values of the Carreau flow index in the high Oh limit.

As shown in Fig. 4, the scaling analysis is based on the interfacial force balance in the vicinity of the hole tip r ≈ r_m. As the cavity approaches collapse, the tip curvature becomes dominated by the radial curvature, and thus, the interfacial force balance in the vicinity of the tip is well approximated by [Eqs. (3) and (6)]

since the tip curvature κ → 1/r_m as r_m → 0.

In addition, as the hole approaches collapse the flow field near the tip becomes essentially one dimensional (radial), and therefore, the strain rate is largely determined by the strain rate in the radial direction $γ̇≈∂vr/∂r∼vm/rm$⁠. Consequently, the apparent viscosity near the tip of the hole scales as [Eq. (4)]

provided that the infinite-shear viscosity is sufficiently low. In the high Oh limit, the momentum balance in the one-dimensional flow near the tip balances fluid pressure and viscous stress

where we use Eq. (9) to estimate the shear-dependent viscosity.

The scalings described above are exemplified quantitatively in Fig. 5. This figure shows both the fluid viscosity μ_m [Fig. 5(a)] and the pressure p_m [Fig. 5(b)] evaluated at the tip of the hole (r = r_m) as a function of $v$_m/r_m for the shear-thinning hole of Fig. 3. The results show that the evolution is in general nonlinear and, as the hole collapses and $v$_m/r_m increases, the tip viscosity and pressure eventually attain the scaling of Eqs. (9) and (10), respectively. The substitution of these scalings into the stress balance at the tip of the hole [Eq. (8)] yields

where the scaling exponent α = (n − 1)/n is in excellent agreement with the value α ≈ −3/7 found in Fig. 3 for the liquid sheet with n = 0.7.

Broadening the analysis, we now use the full numerical simulations to examine the degree to which the scaling predicted by Eq. (11) is able to describe the closing dynamics for different n values. To this end, we show in Fig. 6(a) the velocity of collapse $v$_m for five Carreau indices ranging from n = 0.6 to n = 1. Near the singularity, the solutions of the full governing equations (continuous lines) are in good agreement with the exponent α = (n − 1)/n estimated from scaling analysis (dashed lines) for the different Carreau indices. Moreover, the simulations in Fig. 6(a) can be used to calculate the instantaneous rate of the change of the velocity, defined as α^*(r_m) ≡ d(log $v$_m)/d(log r_m). The results summarized in Fig. 6(b) show that the instantaneous logarithmic rate α^*(r_m) tends to α = (n − 1)/n as r_m → 0, further confirming the scaling of Eq. (11).

In conclusion, we show that the micropores nucleated in shear-thinning liquid sheets contract at a faster speed than the corresponding Newtonian holes due to the localized decrease in the apparent viscosity in the vicinity of the hole front. Scaling arguments show that the shear thinning effect is reflected in a power-law relationship $vm∼rm(n−1)/n$ between the rate of contraction $v$_m and the hole radius r_m for purely viscous liquids. The results from the scaling arguments are supported and fully characterized using high-fidelity simulations.

This work provides a promising starting point for investigations into the collapse of micropores that incorporate the full non-Newtonian physics of the system into the studies of hole collapse. This analysis can be extended in many ways by addressing the role of the fluid’s elasticity, as has been done for opening holes (Dusterhoft and Penno, 2011; Deka et al., 2019), for example. Moreover, our analysis has been limited to purely viscous fluids and should be extended to address the influence of the Ohnesorge number. In this regard, it is instructive to estimate how the local Reynolds number $Rem≡ρv^mr^m/μ^m$ scales with the radius of the orifice. According to Eqs. (9) and (11), the local Reynolds varies with the hole radius as

which shows that Re_m decreases during collapse for n > 2/3 and increases for n < 2/3. The latter case suggests that for real fluids with finite values of Oh, the effect of inertia would become relevant at some point during collapse, regardless of Oh.

This project was partially supported by the USDA National Institute of Food and Agriculture, AFRI Project No. 2017-05024 and Hatch Project No. 1017342. S.U. acknowledges support from CONICET Grant No. PIP-2017-843.