Mode-coupling approach to non-Newtonian Hele-Shaw flow

Magdalena Constantin

Department of Physics, Theoretical Condensed Matter Center, University of Maryland, College Park, MD 20742

Michael Widom

Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213

José A. Miranda

Laboratório de Física Teórica e Computacional, Departamento de Física,

Universidade Federal de Pernambuco, Recife, PE 50670-901 Brazil (September 23, 2002)

Abstract

The Saffman-Taylor viscous fingering problem is investigated for the displacement of a non-Newtonian fluid by a Newtonian one in a radial Hele-Shaw cell. We execute a mode-coupling approach to the problem and examine the morphology of the fluid-fluid interface in the weak shear limit. A differential equation describing the early nonlinear evolution of the interface modes is derived in detail. Owing to vorticity arising from our modified Darcy’s law, we introduce a vector potential for the velocity in contrast to the conventional scalar potential. Our analytical results show how mode-coupling dynamics leads to suppression or enhancement of tip-splitting and side branching in both shear thinning and shear thickening cases. The development of non-Newtonian interfacial patterns in rectangular Hele-Shaw cells is also analyzed.

Typeset using REVTEX

I. INTRODUCTION

The celebrated Saffman-Taylor instability [1] arises when a less viscous fluid pushes a more viscous one in the thin gap of a Hele-Shaw cell. The less viscous fluid can be either injected at an end of a channel-shaped cell (rectangular geometry) [1], or from the center of the cell (radial geometry) [2]. In both geometries the interface separating the fluids may deform, leading to the formation of fingerlike patterns, commonly known as viscous fingers. During the last four decades, the viscous fingering instability has been extensively studied, both theoretically [3] and experimentally [4]. Much of the research in this area has examined the case in which the fluids involved are Newtonian. For Newtonian fluids, it is observed that the fingers grow and compete dynamically, resulting in a single stable finger in the rectangular geometry, and in patterned structures markedly characterized by the spreading, and subsequent splitting of the finger tips in the radial setup.

A whole different class of interfacial patterns arise when the Saffman-Taylor instability is studied by taking the displaced fluid as non-Newtonian [4,5]. In contrast to most Newtonian fluids, non-Newtonian fluids differ widely in their hydrodynamic properties, with different fluids exhibiting a range of effects from elasticity and plasticity to shear thinning and shear thickening. Experiments using non-Newtonian fluids in radial and rectangular Hele-Shaw cells have revealed a wide variety of new patterns, showing snowflake-like shapes [4,6] and fracturelike structures [7,8]. Instead of the traditional, tip-splitting-dominated Newtonian patterns, these experiments [4,6–8] exhibit dendritic fingers and side branching. This morphological diversity and rich dynamical behavior motivated a number of theoretical studies of the problem [9–13]. One major difficulty faced by researchers is that, as opposed to the Newtonian case, the pressure field is no longer Laplacian. This implies a serious theoretical challenge, since a priori one would not be allowed to directly apply a Darcy’s law approach to attack the problem.

The instability of radial Hele-Shaw flows involving non-Newtonian fluids has been studied theoretically by Sader et al. [9]. They considered power law fluids, and performed a linear stability analysis without making use of Darcy’s law. Essentially, they showed that decreasing the power law index dramatically increases the growth rates leading to a more rapid development of the fingering patterns. A Darcy’s law-type approach has been proposed by Bonn and collaborators [10]. Bonn et al. suggested a modified Darcy’s law including a shear rate dependent viscosity, and showed that, within their approach, the pressure field remains Laplacian. In a series of interesting papers Kondic et al. [11,12], and subsequently Fast et al. [13] extended reference [10] ideas and derived a generalized Darcy’s law from first principles, where viscosity depends upon the squared pressure gradient. It turns out that the Darcy’s law formula proposed by Bonn et al. [10] follows from the more basic version rigorously derived by Kondic et al. [11,12]. Efficient numerical simulations performed in references [12,13] have shown that shear thinning can suppress tip-splitting and leads to the formation of dendritic structures, presenting a clear side branching behavior.

Theoretical studies of the fully nonlinear stages of Hele-Shaw flow with non-Newtonian fluids rely heavily on intensive numerical simulations [12,13]. On the analytical side, the structure of the fingering dynamics in such complex fluids is largely restricted to linear stability investigations [9,11]. Much less attention has been paid to the analytical investigation of the dynamics that bridges the initial (purely linear) and final (fully nonlinear) time regimes. In addition, theoretical as well as experimental analyses of flow of shear thickening fluids in Hele-Shaw cells still need to be addressed. In this paper we carry out the analytical weakly nonlinear analysis for the intermediate stages of evolution, and examine both shear-thinning and shear thickening cases. We adapt a weakly nonlinear approach originally developed to study Newtonian Hele-Shaw flows [14,15], to the non-Newtonian situation. We focus on the onset of the nonlinear effects, and try to understand how mode-coupling dynamics leads to basic morphological features and behaviors observed in non-Newtonian Hele-Shaw flows.

The article is organized as follows: section II formulates our theoretical approach. We perform a Fourier decomposition of the interface shape, and from an alternative form of Darcy’s law study the influence of weak shear effects on the development of interfacial patterns. Coupled, nonlinear, ordinary differential equations governing the time evolution of Fourier amplitudes are derived in detail. Section III discusses both linear and weakly nonlinear dynamics. It concentrates on the effect of shear thinning and shear thickening on finger tip-splitting and side branching. Section III A briefly discusses our linear stability results. Linear results are useful and instructive, but do not allow accurate predictions of important interfacial features. In section III B we show that some of such features can indeed be predicted and more quantitatively explained by our analytical mode-coupling approach. At second order we find a mechanism for finger tip-splitting in non-Newtonian Hele-Shaw flow: it is suppressed (favored) for shear thinning (thickening) fluids. Our results indicate absence of side branching in the weak shear limit and early flow stages, but suggest that it would be enhanced (inhibited) for shear thinning (thickening) fluids. Section III C discusses mode coupling in rectangular flow geometry. Our chief conclusions are summarized in Sec. IV.

II. THE MODE COUPLING DIFFERENTIAL EQUATION

The Hele-Shaw cell (Fig.1) consists of two parallel plates separated by a small distance. The cell thickness b is considered to be much smaller that any other length scale in the problem, so that the system is essentially two-dimensional. Consider the displacement of a viscous, non-Newtonian fluid by a Newtonian fluid of negligible viscosity in such confined geometry. The surface tension between the fluids is denoted by σ. The Newtonian fluid is injected at a constant areal flow rate Q at the center of the cell.

The initial circular fluid-fluid interface is slightly perturbed, R = R(t) + ζ (θ, t) (ζ/R ≪ 1), where the time dependent unperturbed radius is given by

R(t) = \sqrt{R_{0}^{2} + \frac{Qt}{π}},  (1)

R_{0} being the unperturbed radius at t = 0. The interface perturbation is written in the form of a Fourier expansion

ζ(θ, t) = ∑_{n=−∞}^{+∞} ζ_{n}(t) exp (inθ),  (2)

where ζ_{n}(t) = (1/2π) ∫_{0}^{2π} ζ(θ, t) exp (−inθ) dθ denotes the complex Fourier mode amplitudes and n=0, ±1, ±2, ... is the discrete polar wave number.

The relevant hydrodynamic equation for Newtonian Hele-Shaw flows, Darcy’s law [1,3], states that the gradient of the pressure is proportional to the fluid velocity, and oriented in the opposite direction with respect to the fluid flow

∇p = −\frac{12μ}{b^{2}} v,  (3)

where v = v(r, θ) and p = p(r, θ) are, respectively, the gap-averaged velocity and pressure in the fluid. The viscosity of the fluid is represented by μ. To model non-Newtonian Hele-Shaw flows, we use a suitable form of the Darcy’s law in the weak shear thinning or thickening limit. We follow Bonn et al. [10] and consider their shear rate dependent viscosity in the weak non-Newtonian limit (ωτ)^{2} ≪ 1

μ(ω^{2}) ≈ μ_{0} [1 − (1 − α) ω^{2}τ^{2}],  (4)

where ω ≈ v/b is the shear rate, v = |v|, τ denotes the characteristic relaxation time of the fluid, and μ_{0} is a constant, zero-shear viscosity. The parameter α measures the shear dependence: α = 1 corresponds to the Newtonian fluids, and α < 1 (α > 1) gives the shear thinning (thickening) case.

By substituing Eq. (4) into Eq. (3) we obtain an alternative form of Darcy’s law ideally suited to describe weak non-Newtonian effects

∇P = −v + δv^{2}v.  (5)

Here P = [b^{2}/(12μ_{0})] p represents a generalized pressure field, and δ = (1−α)(τ/b)^{2} is a small parameter that expresses the non-Newtonian nature of the displaced fluid: δ = 0 corresponds to the Newtonian case, while δ > 0 (δ < 0) describes the shear thinning (thickening) case. Our pertubative approach keeps terms up to the second order in ζ and up to first order in δ. Considering the fact that the viscosity of the Newtonian fluid is negligible, the generalized pressure jump condition at the interface can be written as [3]

P |_{R} = −γ κ_{∥}|_{R},  (6)

where γ = b^{2}σ/(12μ_{0}) and κ_{∥} is the curvature in the direction parallel to the plates.

The weakly nonlinear approach to radial, Newtonian Hele-Shaw flow developed in reference [15], related the fluid velocity to a scalar velocity potential v = −∇ϕ, this replacement made possible by the irrotational nature of the flow for Newtonian fluids. For non-Newtonian fluids, in contrast, flows governed by the modified Darcy’s law (5) exhibit vorticity. Hence we perform our calculations using a vector potential v = ∇ × A. The most general form of the vector potential can be written as

A = [\frac{Q}{2π} θ + ∑_{m,n≠0} A_{mn} (\frac{R}{r})^{m} exp (inθ)] ẑ,  (7)

where A_{mn} are the Fourier coefficients of the velocity vector potential and ẑ is the outward unit-normal to the upper cell plate. The radial and polar components of the fluids velocities are

v_{r} = \frac{Q}{2πr} + ∑_{m,n≠0} inA_{mn} (\frac{R^{m}}{r^{m+1}} ) exp (inθ),  (8)

and

v_{θ} = ∑_{m,n≠0} mA_{mn} (\frac{R^{m}}{r^{m+1}} ) exp (inθ).  (9)

Note that the vector potential reduces to the unperturbed steady flow with a circular interface (v_{r} = Q/2πr, v_{θ} = 0) when R → 0 and also when r → ∞.

We exploit the fact that ∇P must be curl free, and impose the so-called solvability condition ∇×∇P = 0. It simplifies the general form of the vector potential expansion given in Eq. (7). The solvability condition reveals that, without loss of generality, one can rewrite the vector potential as

A = \{\frac{Q}{2π} θ + ∑_{n≠0} A_{n} (\frac{R}{r})^{|n|} exp (inθ) + δ [∑_{n≠0} B_{n} (\frac{R}{r})^{|n|} \frac{1}{r^{2}} exp (inθ)]\} ẑ,  (10)

replacing the array of coefficients A_{mn} with the simpler set of A_{n} and B_{n}. Observe that the vector potential (10) is simply a superposition of a purely Newtonian term (∝ δ^{0}, coefficients A_{n}) and a non-Newtonian contribution (∝ δ^{1}, coefficients B_{n})

A = A_{N} + A_{NN}.  (11)

The flow described by A_{N} is irrotational, while A_{NN} has a curl.

Similarly, we express the pressure of the outer fluid as a sum of Newtonian and non-Newtonian pressures, and propose a general form for their Fourier expansion

P = P_{N} + P_{NN},  (12)

where

P_{N} = −\frac{Q}{2π} log (\frac{r}{R}) + ∑_{n≠0} p_{n} (\frac{R}{r})^{|n|} exp (inθ),  (13)

and

P_{NN} = δ [− (\frac{Q}{2π} )^{3} \frac{1}{2r^{2}} + ∑_{n≠0} q_{n} (\frac{R}{r})^{|n|} \frac{1}{r^{2}} exp (inθ)] .  (14)

The gradient of the complex pressure field (12) must satisfy the non-Newtonian Darcy’s law given by Eq. (5). By inspecting the r and θ components of (5), and by examining the Newtonian and non-Newtonian components of it, we can express the Fourier coefficients of P_{N}, P_{NN} and A_{NN} in terms of the Fourier coefficients of A_{N},

p_{n} = sgn(n) (iA_{n}) ,  (15)
q_{n} = −iA_{n}β(n)\frac{n}{|n|} (\frac{Q}{2π} )^{2} + \frac{Q}{4π (|n| + 1)} ∑_{m≠0,m≠n} mA_{m} (n − m) A_{n−m} k (n, m) , (16)
B_{n} = −A_{n}α(n) (\frac{Q}{2π} )^{2} + \frac{Q}{4π (|n| + 1)} ∑_{m≠0,m≠n} m (iA_{m}) (n − m) A_{n−m} h (n, m) , (17)

where in order to keep the results in a compact form, we introduced the coefficients

α (n) = \frac{|n| (|n| − 1)}{2 (|n| + 1)},  (18)
β (n) = \frac{|n| (|n| + 3)}{2 (|n| + 1)},  (19)
h (n, m) = (|n| + 2) sgn(n − m) − \frac{n}{2} \{3 − sgn [m (n − m)]\} ,  (20)

and

k (n, m) = −n sgn(n − m) + \frac{1}{2} (|n| + 2) \{3 − sgn [m (n − m)]\} .  (21)

Note that sgn(n) = 1 if n > 0 and sgn(n) = −1 if n < 0.

Using Eqs. (15)-(17), which are consistent with the solvability condition and Darcy’s law (5), we can derive the general expression of the vector potential Fourier coefficients in terms of the perturbation amplitudes. To fulfill this goal, consider the kinematic boundary condition [3,4]

\frac{∂R}{∂t} = [\frac{1}{r} \frac{∂R}{∂θ} (−v_{θ}) + v_{r}]_{|_{R}} ,  (22)

which states that the normal components of each fluid’s velocity at the interface equals the velocity of the interface itself. By expanding Eq. (22) up to the second order in ζ and up to first order in δ we find the coefficient of the vector potential corresponding to the n-th evolution mode, A_{n}^{(k)}, in terms of δ and the k-th order in ζ (k = 1, 2)

iA_{n}^{(1)}(t) = [\frac{R}{n}ζ̇_{n} + \frac{Ṙ}{n}ζ_{n}] [1 + δṘ^{2}α (n)] ,  (23)
iA_{n}^{(2)}(t) = \frac{Ṙ}{R} ∑_{m≠0,n} [\frac{|m|}{m} + δṘ^{2}u(n, m)] ζ_{m}ζ_{n−m} + ∑_{m≠0,n} [\frac{1}{n} + \frac{|m|}{m} + δṘ^{2}v(n, m)] ζ̇_{m}ζ_{n−m}
− \frac{δṘR}{2 (|n| + 1)} ∑_{m≠0,n} h(n, m) ζ̇_{m}ζ̇_{n−m} − \frac{δṘ^{2}}{2 (|n| + 1)} ∑_{m≠0,n} h(n, m) ζ_{m}ζ̇_{n−m},

(24)

where the overdot denotes total time derivative, Ṙ = Q/(2πR) represents the unperturbed interface velocity, and the coefficients

u (n, m) = α(n)\frac{|m|}{m} − 2\frac{α(m)}{m} − \frac{h(n, m)}{2 (|n| + 1)},  (25)
v (n, m) = α(n) (\frac{1}{n} + \frac{|m|}{m}) − 2\frac{α(m)}{m} − \frac{h(n, m)}{2 (|n| + 1)}.  (26)

To conclude our derivation we need one more step. The vector potential coefficients can be introduced into the pressure jump condition (Eq.(6)) and using Darcy’s law (Eq.(5)) one can finally find the equation of motion for perturbation amplitudes ζ_{n}. We present the evolution of the perturbation amplitudes in terms of δ and the k-th order in the perturbation amplitude ζ

ζ̇_{n} = ζ̇_{n}^{(1)} + ζ̇_{n}^{(2)},  (27)

where

ζ̇_{n}^{(1)} = λ (n) ζ_{n},  (28)

λ(n) = \frac{Ṙ}{R} (|n| − 1) − \frac{γ}{R^{3}} |n| (n^{2} − 1) + δṘ^{2} \frac{|n|}{R} [Ṙ\frac{|n| − 1}{|n| + 1} − \frac{2γ}{R^{2}} \frac{|n| (n^{2} − 1)}{|n| + 1}] ,  (29)

is the non-Newtonian linear growth rate first reported in Kondic, et al. [11], and

ζ̇_{n}^{(2)} = ∑_{m≠n,0} [F_{N}(n, m) + δ F_{NN}(n, m)] ζ_{m}ζ_{n−m} + ∑_{m≠n,0} [G_{N}(n, m) + δ G_{NN}(n, m)] ζ̇_{m}ζ_{n−m}
+ δ ∑_{m≠n,0} H_{NN} (n, m) ζ_{m}ζ̇_{n−m} + δ ∑_{m≠n,0} J_{NN}(n, m) ζ̇_{m}ζ̇_{n−m}.  (30)

In Eq. (30) the coefficients F_{N}, F_{NN}, G_{N}, G_{NN}, H_{NN}, and J_{NN} represent the second order Newtonian (N) and non-Newtonian (NN) terms. Their detailed functional form are presented in the appendix. An important feature of the second order coefficients is that they present special reflection symmetries

C (n, −m) = C (−n, m)
C (−n, −m) = C (n, m)  (31)

where C = F_{N}, F_{NN}, G_{N}, G_{NN}, H_{NN}, and J_{NN}, respectively. Equation (27) is the mode coupling equation of the non-Newtonian Saffman-Taylor problem in radial geometry. It gives us the time evolution of the perturbation amplitudes ζ_{n}, accurate to second order, in the weak shear limit. The rest of the paper uses Eq. (27) to study the development of interfacial instabilities, and to examine how the non-Newtonian parameter δ affects pattern morphology.

III. DISCUSSION

In the next three subsections we use our mode-coupling approach to investigate the interface evolution at first and second order. To simplify our discussion it is convenient to rewrite the net perturbation (2) in terms of cosine and sine modes

ζ(θ, t) = ζ_{0} + ∑_{n=1}^{∞} [a_{n}(t) cos(nθ) + b_{n}(t) sin(nθ)] ,  (32)

where a_{n} = ζ_{n} + ζ_{−n} and b_{n} = i (ζ_{n} − ζ_{−n}) are real-valued. Without loss of generality we may choose the phase of the fundamental mode so that a_{n} > 0 and b_{n} = 0.

A. First order

Although at the level of linear analysis we do not expect to detect or rigorously predict important nonlinear effects such as tip-splitting and side branching, some useful information may still be extracted. A nice example of how purely linear results can help to understand complicated morphological features appearing in non-Newtonian radial Hele-Shaw flows is found in references [11,12]. By studying their linear growth rate, Kondic et al. [11,12] found that shear thinning decreases the wave number of maximal growth, increases the maximum growth rate, and tightens the band of unstable modes. Based on this increased selectivity of wavelengths, they postulated that shear thinning can lead to suppression of tip-splitting. Their speculations have not been further investigated analytically, but instead have been supported by their own intensive numerical simulations [12,13].

Our first order results, namely our linear growth rate expression (29), agrees with Kondic et al. [11] equivalent formula, if we set their constant C equal to Ṙ^{2}. We use growth rate (29) to gain insight into shear thinning/thickening behavior. Figure 2 plots λ(n) as a function of mode number n for three different values of the non-Newtoniam parameter: a) δ = 0; b) δ = 0.05; and c) δ = −0.05. By inspecting Fig. 2, we notice that, unlike the shear thinning case (δ > 0) discussed in references [11,12], shear thickening (δ < 0) widens the band of unstable modes, and decreases growth rate for the wave number of maximal growth. So, by applying similar arguments as those used by Kondic et al. [11,12] for the weak shear thinning case, we postulate that shear thickening would lead to enhanced tip-splitting. Unfortunately, for the shear thickening radial Hele-Shaw flow, both numerical simulations and experiments are not available in the literature to confirm this claim. Of course, it is also of interest to study such a possibility analytically. This is one of the topics we examine in Sec. III B.

B. Second order

On the analytical side, one must go beyond linear analysis in order to investigate in more detail the main morphological features found in Hele-Shaw flows with both shear thinning and shear thickening fluids. To this day, besides a few first order linear growth studies [11–13], there are no other analytical results concerning non-Newtonian fluid evolution in radial Hele-Shaw geometry. If on one hand the non-Newtonian problem is vastly richer, it is also true it is much less amenable to analytic attack than its Newtonian counterpart. Our mode-coupling approach intends to provide a better and more complete analytical understanding of the complex non-Newtonian time evolution dynamics in Hele-Shaw cells.

1. Action of the parameter δ on finger tip-splitting

We use the mode coupling equation (27) to investigate the coupling of a small number of modes. At second order the most noteworthy effect refers to the action of the non-Newtonian parameter δ on finger tip-splitting. Tip-splitting is related to the influence of a fundamental mode n on the growth of its harmonic 2n [15]. For consistent second order expressions, we replace the time derivative terms ȧ_{n} and ḃ_{n} by λ(n) a_{n} and λ(n) b_{n}, respectively. Under these circumstances the equation of motion for the cosine mode 2n is

ȧ_{2n} = λ(2n) a_{2n} + \frac{1}{2} T (2n, n) a_{n}^{2},  (33)

where the function that multiplies a_{n}^{2}, T (2n, n), is called the tip-splitting function and its general expression is

T (n, m) = F_{N}(n, m) + λ(m) G_{N}(n, m) + δ [F_{NN}(n, m) + λ(m) G_{NN}(n, m)
+ λ(n − m) H_{NN}(n, m) + λ(m) λ(n − m) J_{NN}(n, m)].  (34)

Equation (33) shows that the presence of the fundamental mode n forces growth of the harmonic mode 2n. The function T (2n, n) acts like a driving force and its sign dictates if finger tip-splitting is favored or not by the dynamics. If T (2n, n) < 0, a_{2n} is driven negative, precisely the sign that leads to finger tip-widening and finger tip-splitting. If T (2n, n) > 0 growth of a_{2n} > 0 would be favored, leading to outwards-pointing finger tip-narrowing.

To investigate the influence of the non-Newtonian parameter δ on tip-splitting behavior at second order, we plot in Fig. 3 the behavior of T (2n, n) as a function of δ, for a few Fourier modes (n = 4, 5, 6). To simplify our analysis we adopt an instantaneous approach: we consider a particular Ṙ and R combination, using the identity ṘR = Q/2π, at the onset of growth of mode 2n [using the condition λ(2n) = 0] in the Newtonian limit δ = 0, where we know T (2n, n) is negative [15]. We see from Fig. 3 that the curve for mode n = 4 lies below the curves associated to n = 5, 6. This is an expected behavior since smaller values of n would mean more room for the existing fingers to split. We also observe that the curves associated to smaller n show a stronger angular inclination with respect to the horizontal δ-axis. Therefore, lower Fourier modes would be more sensitive to variations in δ.

Further inspection of Fig. 3 reveals that, for δ > 0, T (2n, n) becomes less negative as δ increases, meaning that the interface has less tendency toward tip-splitting. In contrast, for δ < 0 we observe that T (2n, n) becomes more negative as the magnitude of δ < 0 increases, indicating an enhanced tendency of the fingers to split at their tips. We conclude that tip-splitting is suppressed for shear thinning fluids, and enhanced for shear thickening ones. This fact is more clearly illustrated in Fig. 4. In Fig. 4 we plot the fluid-fluid interface for a certain time (t = 30 s), considering the interaction of two cosine modes (a fundamental n = 4 and its harmonic 2n = 8), for three diferrent values of δ: (a) δ = 0 (solid curve); (b) δ > 0 (black dashed curve); and (c) δ < 0 (gray dashed curve). From Fig. 4 it is evident that there is a stronger splitting in the shear thickening case.

These second order effects regarding tip-splitting are consistent with the first order effects described in section III A. With respect to shear thinning behavior, our analytical results agree with numerical simulations [12,13] and experiments [4,6] of fully nonlinear stages of interface evolution. In addition, we detect favored tip-splitting in the shear thickening case, a relevant nonlinear behavior not previously reported in the literature. We have also verified that, for a given mode n, smaller values of the surface tension parameter γ lead to enhancement (suppression) of splitting events for δ < 0 (δ > 0).

A physical mechanism that seems to be at work in shear thinning case has been proposed in references [11–13]: viscosity is lowered in high-shear regions, which are in front of the fingers. Hence the tip of a finger experiences less resistance than the part of the interface near the tip. This suggests that tip-splitting would be suppressed if the fluid shear-thins. Our second order mode-coupling results indicate that a similar mechanism is at work in the shear thickening case: the resistance would be increased at the finger tips, suggesting that tip-splitting would be increased.

2. Action of the parameter δ on side branching

Another relevant non-Newtonian effect that can be studied at second order refers to the side branching phenomenon. In the framework of a mode-coupling theory, side branching requires the presence of mode 3n. If the harmonic mode a_{3n} is positive and sufficiently large, it can produce interfacial lobes branching out sidewards which we interpret as side branching.

Consider the influence of the fundamental mode n, and its harmonic 2n, on the growth of mode 3n. The equation of motion for the cosine 3n mode is

ȧ_{3n} = λ(3n) a_{3n} + \frac{1}{2} S(3n) a_{n}a_{2n},  (35)

where the side branching function S(3n) = [T (3n, n) + T (3n, 2n)] can be easily obtained from Eq. (34). By analyzing Eq. (35) we observe that mode 3n can be spontaneously generated due to the driving term proportional to a_{n}a_{2n}, such that it enters through the dynamics even when it is missing from the initial conditions. The existence and phase of mode 3n depends on the interplay of the modes n and 2n. Side branching would be favored if a_{3n} > 0.

To study the growth of mode 3n as the non-Newtonian parameter δ is varied, we plot S(3n) as a function of δ in Fig. 5. We consider a particular Ṙ and R combination that corresponds to the onset of growth of mode 3n (i.e. obeying the condition λ(3n) = 0) in the Newtonian limit δ = 0. From Fig. 5 one can verify that in the weak shear limit we consider in this work, S(3n) is negative for all values of δ. As shown in section III B 1, starting with a fundamental mode a_{n}, the harmonic mode a_{2n} is driven negative. Hence the product S(3n)a_{n}a_{2n} in eq. (35) is positive, driving a_{3n} > 0.

Whether side branching actually occurs depends on the magnitude of a_{3n}. Indeed, our discussion shows the presence of a_{3n} > 0 even in the case of Newtonian fluids, where side branching does not occur. Why does side branching occur for shear thinning fluids? We speculate its origin may lie in the cubic form of our modified Darcy’s law, eq. (5), which will add a term of the form δa_{n}^{3} to ȧ_{3n}. Unfortunately, we have not calculated this term since we stopped our derivation at second order.

Furthermore, there is a delicate interplay of mode 3n and mode 2n, described by the evolution equation (35) and a similar expression for the growth of mode 2n,

ȧ_{2n} = λ(2n) a_{2n} + \frac{1}{2} [T (2n, n) a_{n}^{2} + S(2n) a_{n}a_{3n}],  (36)

where S(2n) = [T (2n, −n) + T (2n, 3n)] > 0. Hence, side branching via a positive a_{3n} will tend to drive a_{2n} positive (or at least less negative), inhibiting tip splitting but also reducing the growth rate of a_{3n} itself.

C. Rectangular geometry limit

It is interesting to study how the main interfacial features examined in the radial geometry behave if we consider flow of non-Newtonian fluids in rectangular Hele-Shaw cells. A mode-coupling equation describing the system can be obtained if we take the “rectangular geometry limit”: R → ∞ and Q → ∞, such that Q/(2πR) ≡ v_{∞} and n/R ≡ k remain constant, where v_{∞} is the flow velocity at infinity and k denotes the wave number of the disturbance. In this limit the interface evolution reverts to the evolution of the rectangular flow geometry, with linear growth rate given by

λ(k) = |k| [v_{∞}(1 + δv_{∞}^{2}) − γ|k|^{2}(1 + 2δv_{∞}^{2})] .  (37)

In addition, the only nonzero second order mode-coupling terms are

G_{N}(k, k′) = |k| [1 − sgn(k k′)] ,  (38)
G_{NN}(k, k′) = −|k|v_{∞}^{2} [1 − sgn(k k′)] ,  (39)

and

J_{NN}(k, k′) = \frac{v_{∞}}{2} \{3 − sgn[k′(k − k′)] + 2 sgn[k(k − k′)]\} .  (40)

Note that all rectangular limit expressions (37)-(40) are time-independent, and preserve the reflection symmetries shown in Eq. (31). Based on these findings, and using Eq. (27) we obtain a more compact mode-coupling equation for non-Newtonian flow in rectangular Hele-Shaw cells

ζ̇_{k} = λ(k) ζ_{k} + ∑_{k′≠k,0} [G_{N}(k, k′) + δ G_{NN}(k, k′)] ζ̇_{k′}ζ_{k−k′}
+ δ ∑_{k′≠k,0} J_{NN}(k, k′) ζ̇_{k′} ζ̇_{k−k′}.  (41)

The rectangular mode-coupling Eq. (41) is useful to study the influence of weak shear effects on the development of the fluid-fluid interface in rectangular cells. First, we examine finger tip-splitting related issues (finger narrowing/widening). As discussed in section III B 1, we analyze the influence of the fundamental k on the growth of its harmonic 2k. The equation of motion for the cosine mode 2k is

ȧ_{2k} = λ(2k) a_{2k} + \frac{1}{2} T (2k, k) a_{k}^{2},  (42)

where

T (2k, k) = δ [2k^{2}v_{∞}(v_{∞} − γ k^{2})^{2}] .  (43)

Note that the driving force term T (2k, k) vanishes in the Newtonian limit δ = 0. This agrees with the results obtained in Ref. [14] for Newtonian rectangular flow. Notice further, that for the non-Newtonian case (δ ≠ 0) the harmonic mode 2k grows spontaneously even if missing from the initial conditions. The selected sign of a_{2k} is given by the sign of T (2k, k) and thus is dictated by δ: if δ > 0 (δ < 0), T (2k, k) is positive (negative), and hence a_{2k} > 0 (a_{2k} < 0).

Considering the case δ > 0, this means that the fingers become narrower as the shear thinning behavior becomes more pronounced. Moreover, by inspecting Eq. (43) we observe that the width of the finger tips decreases for increasingly larger values of the flow velocity v_{∞}. Of course, these effects are minimized for larger values of the surface tension parameter γ. All these results are in agreement with recent experimental [16] and theoretical investigations [17,18] for weak shear-thinning flows in rectangular geometry, in which shear-thinning narrows the width of steady state fingers.

If we take k = k^{⋆} =≈ \sqrt{v_{∞}/3γ} as the fastest growing mode, then λ(2k) < 0 so exponential growth of a_{2k} is prevented. Unlike the case of radial flow where any mode n eventually goes unstable for sufficiently large R, the growth rates λ(k) in the rectangular geometry are time independent. Instead of exponential growth of a_{2k}, we instead expect its amplitude to saturate at the value T (2k, k)a_{k}^{2}/(−2λ(2k) obtained by setting ȧ_{2k} = 0 in eq. (42). Alternatively, one could specially prepare an initial condition with an initial perturbation of wavevector k sufficiently small that 2k still lies within the band of unstable modes (up to about \sqrt{v_{∞}/γ}), in which case mode 2k will be spontaneously generated and able to grow to a considerable magnitude.

For completeness, we also discuss the shear thickening case: for δ < 0 the function T (2k, k) is negative, favoring finger tip-widening. Even though this shear thickening behavior has not been studied experimentally, the broadening of the fingers in shear thickening flow in rectangular cells has been theoretically predicted by numerical simulations [18]. Our analytical results reinforce the correctness of such numerical predictions.

Finally, we briefly discuss side branching behavior in rectangular Hele-Shaw cells. We consider the influence of modes k and 2k on the evolution of mode 3k. In this case, the relevant equation of motion has the form

ȧ_{3k} = λ(3k) a_{3k} + \frac{1}{2} S(3k) a_{k}a_{2k},  (44)

where, as it was in the radial flow case, the coefficient S(3k) is the sum of T (3k, k) and T (3k, 2k)

S(3k) = δ [8k^{2}v_{∞}(v_{∞} − γ k^{2}) (v_{∞} − 4γk^{2})] .  (45)

The favored sign of a_{3k} depends on the sign of the product a_{k}a_{2k} and also on the sign of S(3k), which is positive (if δ is positive) in the band of wavevectors \sqrt{v_{∞}/4γ} < k < \sqrt{v_{∞}/γ} and negative outside this band. As we saw previously in the case of mode 2k, if the fundamental mode is taken as the fastest growing mode k = k^{⋆}, then exponential growth of the harmonic 3k is inhibited because λ(3k) < 0. Instead it will saturate at a magnitude of S(3k)a_{k}a_{2k}/(−2λ(3k).

IV. CONCLUDING REMARKS

Visually striking patterns arise when a less viscous Newtonian fluid displaces a more viscous non-Newtonian fluid in the confined geometry of a Hele-Shaw cell. These complex patterned structures are the result of finger tip-splitting and dendritic interfacial instabilities. Due to the complicated dynamics of the system, and also from the limitations imposed by purely linear analysis, the majority of the theoretical studies in this area of research rely heavily on sophisticated numerical methods. In this work, we developed an alternative theoretical approach to the problem which allowed us to address important nonlinear issues analytically.

We started our investigation by analyzing the role of fluid viscosity anisotropy on the development of the Saffman-Taylor instability in radial Hele-Shaw cells. To approach the problem analytically we considered the weak shear limit, and focused on the onset of nonlinear effects. In order to examine the influence of shear thinning and shear thickening on the shape of the emerging patterns, we derived a mode-coupling equation which is ideally suited to describe the weakly nonlinear interface evolution. Our analytical results show that finger tip-splitting is enhanced (diminished) in the case of shear thickening (thinning) fluids. We applied the rectangular geometry limit to our radial mode coupling equations of motion, and also studied finger-tip and side branching behavior in rectangular cells. Our main findings predict finger tip-widening (narrowing) for shear thickening (thinning) fluids.

In neither the radial nor the rectangular geometry limit were we able to show the appearance of side branching. We speculate that extending the current approach to third order could shed more light on this problem.

In summary, our analytical mode-coupling approach detects and predicts several features of the patterns formed in non-Newtonian Hele-Shaw flows in both radial and rectangular geometries. Furthermore, it predicts some behaviors not yet investigated in the literature, such as the evolution mechanisms involving tendency towards advanced finger tip-splitting and reduction of side branching for shear thickening fluids. We hope the main results presented in this paper will prompt further theoretical and experimental work on non-Newtonian Hele-Shaw flows, especially in the area of early stage transients in rectangular geometry where our theory makes specific predictions that have not been subject to experimental or numerical check.

ACKNOWLEDGMENTS

J.A.M. thanks CNPq (through its PRONEX Program) for financial support. Work of M.W. was supported in part by the National Science Foundation grant No. DMR-9732567.

APPENDIX A: SECOND ORDER MODE-COUPLING TERMS

This appendix presents the expressions for the second order Newtonian (N), and non-Newtonian (NN) mode-coupling coefficients which appear in Eq. (30)

F_{N} (n, m) = \frac{|n|}{R} \{\frac{Ṙ}{R} [\frac{1}{2} − sgn(nm)] − \frac{γ}{R^{3}} (1 − \frac{nm}{2} − \frac{3m^{2}}{2} )\} ,  (A1)
F_{NN}(n, m) = \frac{|n|}{R}Ṙ^{2}\{\frac{Ṙ}{R} [\frac{2α(m)}{m} \frac{|n|}{n} − \frac{2β(n)}{|n|} + \frac{f(n, m)}{2(|n| + 1)}]
− \frac{γ}{R^{3}} \frac{2|n|}{(|n| + 1)} (1 − \frac{nm}{2} − \frac{3m^{2}}{2} ) \},  (A2)
G_{N} (n, m) = \frac{|n|}{R} [1 − sgn(nm) − \frac{1}{|n|}] ,  (A3)
G_{NN} (n, m) = \frac{|n|}{R}Ṙ^{2} [\frac{2α(m)}{m} \frac{|n|}{n} − \frac{2β(n)}{|n|} + \frac{f(n, m)}{2(|n| + 1)}] ,  (A4)
H_{NN} (n, m) = \frac{|n|}{2(|n| + 1)} \frac{Ṙ^{2}}{R} f(n, m),  (A5)

and

J_{NN} (n, m) = \frac{|n|}{2(|n| + 1)}Ṙ f(n, m),  (A6)

where

f(n, m) = [\frac{|n|}{n}h (n, m) + k (n, m)] .  (A7)

REFERENCES

[1] P. G. Saffman and G. I. Taylor, Proc. Roy. Soc. London Ser. A245, 312 (1958).

[2] L. Paterson, J. Fluid Mech. 113, 513 (1981).

[3] B. Bensimon, L. P. Kadanoff, S. Liang, B. I. Shraiman, and C. Tang, Rev. Mod. Phys. 58, 997 (1986).

[4] K. V. McCloud and J. V. Maher, Phys. Rep. 260, 139 (1995).

[5] H. Van Damme, E. Lemaire, O. M. Abdelhaye, A. Mourchid, and P. Levitz, Nonlinearity and Breakdown in Soft Condensed Matter, edited by K. K. Bardan et al., Lecture Notes in Physics Vol. 437, (Springer, New York, 1994).

[6] A. Buka, P. Palffy-Muhoray, and Z. Racz, Phys. Rev. A 36, 3984 (1987).

[7] E Lemaire, P. Levitz, G. Daccord, and H. Van Damme, Phys. Rev. Lett. 67, 2009 (1991).

[8] H. Zhao and J. V. Maher, Phys. Rev. E 47, 4278 (1993).

[9] J. E. Sader, D. Y. C. Chan, and B. D. Hughes, Phys. Rev. E 49, 420 (1994).

[10] D. Bonn, H. Kelley, M. Ben Amar and J. Meunier, Phys. Rev. Lett. 75, 2132 (1995); D. Bonn, H. Kelley, M. Braunlich, M. Ben Amar, and J. Meunier, Physica A 220, 60 (1995).

[11] L. Kondic, P. Palffy-Muhoray, and M. J. Shelley, Phys. Rev. E 54, R4536 (1996).

[12] L. Kondic, M. J. Shelley, and P. Palffy-Muhoray, Phys. Rev. Lett 80, 1433 (1998).

[13] P. Fast, L. Kondic, M. J. Shelley, and P. Palffy-Muhoray, Phys. Fluids 13, 1191 (2001).

[14] J. A. Miranda and M. Widom, Int. J. Mod. Phys. B 12, 931 (1998).

[15] J. A. Miranda and M. Widom, Physica D 120, 315 (1998).

[16] A. Lindner, D. Bonn, and J. Meunier, Phys. Fluids 12, 256 (2000).

[17] M. Ben Amar and E. Corvera Poiré, Phys. Rev. Lett 81, 2048 (1998).

[18] M. Ben Amar and E. Corvera Poiré, Phys. Fluids 11, 1757 (1999).

FIGURE CAPTIONS

FIG. 1: Schematic configuration of radial flow in a Hele-Shaw cell. The inner fluid is Newtonian and has negligible viscosity. The outer fluid is non-Newtonian. All physical parameters are defined in the text.

FIG. 2: Linear growth rate (29) as a funciton of mode number n for three different values of the non-Newtoniam parameter: δ = 0 (black solid curve), +0.05 (black dashed curve), -0.05 (gray dashed curve). The injection rate Q = 2π cm^{2}/s, R = 1 cm, and γ = 1/200 cm^{3}/s. The units of λ(n) and δ are s^{−1} and (cm/s)^{−2}, respectively.

FIG. 3: Variation of T (2n, n) as a function of the non-Newtonian parameter δ, for three different Fourier modes (n = 4, 5, 6). The injection rate Q = 2π cm^{2}/s and γ = 1/100 cm^{3}/s. The units of T (2n, n) and δ are (cm s)^{−1} and (cm/s)^{−2}, respectively.

FIG. 4: Snapshot of the fluid-fluid interface (t = 30 s), for the interaction of two cosine modes n = 4 and 2n = 8. Three values of the non-Newtonian parameter are considered δ = 0 (black solid curve), δ = +0.01 (cm/s)^{−2} (black dashed curve), and δ = −0.01 (cm/s)^{−2} (gray dashed curve). Other physical parameters are: a_{n}(0) = 0.001 cm, R_{0} = 0.3 cm, γ = 0.025 cm^{3}/s, and Q = 3π cm^{2}/s. Splitting is favored for shear-thickening fluids (gray dashed curve).

FIG. 5: Variation of the side-branching function S(3n) = T (3n, n) + T (3n, 2n) with the non-Newtonian parameter δ, for modes n = 4, 5, 6. The injection rate Q = 2π cm^{2}/s and γ = 1/100 cm^{3}/s. The units of S(3n) and δ are (cm s)^{−1} and (cm/s)^{−2}, respectively.