[

# [

[
10 March 2010
###### Abstract

We discuss the creeping motion of plugs of negligible viscosity in rough capillary tubes filled with carrier fluids. This extends Bretherton’s research work on the infinite-length bubble motion in a cylindrical or smooth tube for small capillary numbers Bretherton_1961 (Bretherton 1961). We first derive the asymptotic dependence of the plug speed on the finite length in the smooth tube case. This dependence on length is exponentially small, with a decay length much shorter than the tube radius . Then we discuss the effect of azimuthal roughness of the tube on the plug speed. The tube roughness leads to an unbalanced capillary pressure and a carrier fluid flux in the azimuthal plane. This flux controls the relaxation of the plug shape to its infinite-length limit. For long-wavelength roughness, we find that the above decay length is much longer in the rough tube, and even becomes comparable to the tube radius in some cases. This implies a much-enhanced dependence of the plug speed on the plug length. This mechanism may explain the catch-up effect seen experimentally Ying_2008 (Ismagilov & Ying).

Plugs in rough capillary tubes]Plugs in rough capillary tubes: enhanced dependence of motion on plug length Q. Zhang, K. S. Turitsyn, and T. A. Witten]Q. ZHANG,K. S. TURITSYN,ANDT.A.WITTEN

Key Words: Micro-fluid dynamics, multiphase flows, drops and bubbles.

## 1 Introduction

The field of multiphase microfluidics has expanded greatly in the last decade because of the numerous applications of the technique in chemistry and biology Gunther_2006 (Gunther & Jensen 2006). Transporting reactants through microchannels in the form of bubbles or droplets has many advantages over the traditional single-phase systems. These include enhanced mixing rate, reduced dispersion and higher interface areas Gunther_2006 (Gunther & Jensen 2006). However there are still many technological problems that remain unsolved. A particular problem that aroused our attention is the control of spacing within a train of long droplets or plugs carried along the channel by a wetting fluid. Such plugs are separated from the tube by a thin lubricating film, whose thickness is controlled by the speed of motion Bretherton_1961 (Bretherton 1961). It is commonly observed Ying_2008 (Ismagilov & Ying) that the separation between the plugs changes as they move. Eventually one plug can coalesce with its neighbor, a process known as the catch-up effect. To our knowledge there is no widely accepted explanation of this effect.

A change in the distance between neighboring plugs is possible only if the fluid flux in the lubricating layer is different in the two plugs. Thus in order to understand the reason for the catch-up effect one needs to analyze the shape and dynamics of this lubricating layer. There are many physical mechanisms that could be responsible for the structure of the lubricating layer. In this work we analyze only one of them: irregularities of the channel shape. We argue that these irregularities, i.e. roughness of the tube wall, could be responsible for the fluctuations of the distance between the plugs. The influence of the substrate geometry or roughness on the lubricating layer thickness has been recognized for several decades and has been studied quite extensively (Stillwagon & Larson (1988); Schwartz & Weidner (1995); Kalliadasis, Bielarz & Homsy (2000); Mazouchi & Homsy (2001); Howell (2003)). In the present work we consider an effect beyond the direct affect on speeds. We find that the roughness enhances dramatically the sensitivity of the thin film width to the plug length, and thus produces speed variations between individual plugs that might eventually lead to catch-up coalescence events.

The earliest discussion of the roughness effect on the plug motion inside capillary tubes dates back to the seminal paper of Bretherton Bretherton_1961 (Bretherton 1961). He suggested that the roughness can be important if its amplitude is comparable to or larger than the lubricating layer thickness. Therefore the roughness effect is usually negligible for large enough tubes where other physical mechanisms are more important (see e.g. Gunther_2006 (Gunther & Jensen 2006), Ajaev_2006 (Ajaev & Homsy 2006) for review). However as modern microfluidic channels become smaller, surface roughness becomes more important relative to the channel width. Probably the first experimental observation that is attributed to the effect of roughness is reported in the work of Chen Chen_1986 (Chen 1986), where deviations from the Bretherton’s predictions were observed for small capillary numbers .

On the theoretical side the most widely studied effect of the surface roughness corresponds to the limit where the roughness amplitude is significantly smaller than the lubricating film thickness. In this case the roughness effect can be accounted for via effective boundary conditions: the rough surface can be modeled as a smooth one with non-zero slip length determined by the roughness properties (Einzel, Panzer & Liu (1990); Miksis & Davis (1994); DeGennes (2002)). Krechetnikov and Homsy Krechetnikov_2005 (Krechetnikov & Homsy 2005) have analyzed the effect of the slip length on the lubricating film profile, and showed that the effect is negligible for the slip length smaller than the Bretherton thin film thickness, but leads to a dramatic increase of the lubricating film thickness in the opposite limit: in this case the lubricating film thickness becomes comparable to the slip length and is independent of . In a broader context there have been a numerous studies on how the channel geometry affects the dynamics and shape of plugs. Wong et al. (Wong (1995); Wong (1995)) considered the plug motion in polygonal channels and showed that the shape of the plug strongly depends on the capillary number and that in the limit , the fluid flux in the lubricating film does not depend on . Hazel and Heil Heil_2002 (Hazel & Heil 2002) confirmed this effect by numerical simulations, and showed that a similar behavior is observed in tubes with elliptical cross-sections as well. It should also be noted that the surface roughness can have a non-hydrodynamic effect on the plug motion by changing wetting properties of the surface (Wenzel (1936); Herminghaus (2000); Bico, Tordeux & Quere (2001)) or dynamical destabilization of the lubricating layer due to fluctuations.

In this work we consider capillary tubes with roughness smaller than the lubricating film thickness. In this regime the roughness cannot produce any significant effect on the width of the lubricating layer. However, as we will show, in the rough tube the decay length of the plug shape to its infinite-length limit is increased dramatically. In some cases, the decay length becomes comparable to or even larger than the tube radius. The ripples of the plug shape propagate from the front region to the rear region of the plug, leading to a much-enhanced dependence of the fluid flux on the plug length. On the level of the plug train this effect translates the variations of plug speeds to the fluctuations of the distance between them. Thus it can be responsible for the catch-up events. We also note that this is essentially a non-local effect. To our knowledge this effect has not been discussed perviously.

We will illustrate this effect by using a very simple system, a single air plug of finite length moving through an almost cylindrical tube. For the sake of simplicity, throughout the paper we neglect the effects of gravity, inertia and surfactant concentrations, assuming that the only relevant stresses are due to surface tension and viscosity. We also limit our analysis to the small capillary number regime: . We assume that the tube roughness can be modeled as weakly non-circular cross sections of the tube, which does not change in the longitudinal direction and can be described by a roughness function with being the azimuthal angle. We will consider the situation where the deformation is smooth enough for the lubrication approximation to be valid, that is, , i.e. the deformation is small compared to the Bretherton’s thin film thickness. However, we expect that our scaling results can be safely extrapolated to the regime: .

The structure of this paper is as follows. We start our analysis by revisiting the classical Bretherton solution for the semi-infinite air plug moving in the perfectly cylindrical tube. Then we follow the finite-plug analysis of Teletzke Teletzke_1983 (Teletzke 1983) and extend it to a wider range of plug length. The main finding of this section is the exponential suppression of the speed corrections due to the plug length variations: , where measures the length of the plug, and is the characteristic decay length of the plug shape to its infinite-length limit for the smooth tube case. In Section 3 we explain that the presence of small tube roughness produces a whole spectrum of relaxation modes with different decay lengths . The relaxation mode with the decay length provides the largest contribution to the dependence on the plug length. In Section 4 we estimate of the amplitude of the plug shape distortion due to the roughness. Then we estimate the effect of tube roughness on the finite plug speed in Section 5. We conclude in Section 6 by discussing the limitations of the model and proposing future research directions.

## 2 Plug speed dependence on the plug length in smooth tubes

Suppose we have a smooth cylindrical capillary tube with radius . It is initially filled with a liquid of viscosity , which is called the carrier fluid in the following discussion. Now an air plug with negligible viscosity is forced into the tube and moves slowly together with the carrier fluid due to a pressure gradient along the tube. As discussed in detail later, the plug speed is different from the speed of the carrier fluid far away from the plug Jensen_2004 (Hodges, Jensen & Rallison 2004). In this section, we will investigate the effect of the plug length on the plug speed . Both gravitational and inertial effects are assumed to be negligible.

### 2.1 Bretherton-Teletzke mechanism

In Bretherton’s seminal work Bretherton_1961 (Bretherton 1961), he studied the limit case of semi-infinite air plug motion. He showed that a dimensionless number , called capillary number, controls the plug motion:

 Ca=μUγ, (2.0)

where is an interfacial tension between the plug and the carrier fluid. Bretherton’s steady state solution of a semi-infinite plug is illustrated in Fig.1. For convenience, we have chosen the plug as the frame of reference. In this frame, the tube wall moves with velocity . The shape of a semi-infinite plug consists of three regions: a spherical cap (static region) is connected to a thin film region of almost constant liquid film thickness via a transition remobilization region. Since the internal viscosity of the plug is zero, there can be no shear stress at the plug surface. Thus the carrier fluid in the thin film region moves without shear at the wall speed . Bretherton also found that the lubricating film thickness decays exponentially with a decay length , denoted as Bretherton_1961 (Bretherton 1961). Thus is a measure of the length of the remobilization region. Both capillary force and viscous force are important in the remobilization region, where the plug shape is determined by an equilibrium between these two forces. In the limit of small capillary numbers, , Bretherton derived the dependence of on and Bretherton_1961 (Bretherton 1961):

 h∞=0.643R⋅(3Ca)2/3. (2.0)

Obviously, both and are much less than for small capillary numbers. Relative motion between the plug and the carrier fluid is governed by an incompressibility constraint. This constraint states that, in the steady state, the total volume of the plug and the carrier fluid between two arbitrary cross sections along the tube doesn’t change over time. Mathematically, it can be translated into the equation: . By keeping only the lowest order approximation, we get:

 U/V ≃ 1+2h∞/R (2.0) = 1+1.286(3Ca)2/3.

Strictly speaking, Eqn.(2.0) is only valid in the limit of the semi-infinite plug. For a finite plug, we expect Eqn.(2.0) to be modified as Teletzke_1983 (Teletzke 1983): , with different from its limiting value . Physically speaking, represents the average carrier fluid flux per unit circumference. As shown in detail later, depends on the plug length and approaches as we increase the length.

The motion of a finite-length air plug in a capillary tube was first studied by Teletzke Teletzke_1983 (Teletzke 1983). In the following discussion length, velocity, and pressure are made dimensionless by using tube radius , plug speed , and capillary pressure . Capillary number was kept small in Teletzke’s analysis: . By imposing a lubrication approximation and appropriate boundary conditions (non-slip condition at the tube wall and stress-free condition at the plug surface), the balance of viscous and capillary pressure gradients yields Teletzke_1983 (Teletzke 1983):

 d3hdz3=3Ca⋅(h−h∗)h3, (2.0)

wherever . For small , this includes both the thin film region and the remobilization region. For a plug of fixed length, is a self-consistent parameter whose value is determined by the fact that remobilization regions must match with spherical caps at the plug ends. Suppose the front and rear remobilization regions are separated by a distance for a finite plug. Teletzke proposed a method to solve the above differential equation numerically in the limit of short plugs: Teletzke_1983 (Teletzke 1983). Numerical results obtained by using Teletzke’s method are plotted in Fig.3. As expected, approaches as the plug length increases. However, the trend is not monotonic. The origin of this oscillating behavior will be clear in the next section.

Teletzke’s method did not address the asymptotic behavior for finite plugs with length: . So the knowledge of the dependence of the plug speed on the plug length is incomplete at the current stage. We only understand two limit cases: Bretherton’s semi-infinite plugs and Teletzke’s short plugs with . In the next section, we complete this previously unstudied part of the plug speed problem, by determining the asymptotic behavior of for plugs with large .

### 2.2 Asymptotic corrections to the semi-infinite plug speed

Due to rotational symmetry, the problem of determining plug shape is reduced to calculating a lubricating film thickness . The governing equation of is given by Eqn.(2.0). For small capillary numbers , this differential equation is valid in both the thin film and the remobilization region, as noted above. Its solution has to match spherical caps at the ends of remobilization regions. By using the rescalings: and , we may rewrite Eqn.(2.0) as:

 d3ηdξ3=η−1η3. (2.0)

In the above substitutions, is the remobilization length of the finite plug defined as: . Comparing with Eqn.(2.1), we observe that measures thickness in units comparable to the asymptotic thin film thickness . Likewise evidently measures lengths in units comparable to the remobilization length . We denote the plug length in these units as . So in this section we will study long plugs with .

As we approach the thin film region, becomes arbitrarily close to 1: . So we can simplify Eqn.(2.0) as follows: . This autonomous equation has solutions of the form: , where . It has three roots and , which correspond to three solution modes given below:

 ψ1(ξ) = exp(ξ), ψ2(ξ) = exp(−ξ/2)cos(√3ξ/2), ψ3(ξ) = exp(−ξ/2)sin(√3ξ/2). (2.0)

So general solutions are linear combinations of these modes, and we have:

 η=1+3∑i=1Ciψi(ξ), (2.0)

where coefficients are determined from boundary conditions away from the thin film region. We note for future reference that for large positive , is dominated by the growing function.

Far from the thin film region treated above, the following relation is valid: . For very small , lubrication approximation () is still applicable in this regime Teletzke_1983 (Teletzke 1983). So Eqn.(2.0) can be approximated as follows: . Integration of this equation yields parabolic profiles:

 ηF,R=12PF,Rξ2+QF,Rξ+SF,R, (2.0)

where superscripts and stand for front and rear solutions, respectively. As defined by Teletzke Teletzke_1983 (Teletzke 1983), axial (-direction) positions of and in Fig.2, correspond to axial positions of the lowest points of front and rear parabolas. And he defined the distance between and as the plug length . We will keep Teletzke’s definition in the following discussion.

Since both ends of the plug are spherical caps for small , mean curvatures far from the thin film region must match those of spherical caps. As a result, the following constraint needs to be satisfied near the spherical cap: . We may rewrite this constraint in terms of reduced variables: . By using asymptotic solutions of given in Eqn.(2.0), we get:

 (3Ca)2/3h∗⋅PF,R=1. (2.0)

Since is fixed for a given plug, the above equation implies that . So we can drop the superscript from now on. For a semi-infinite plug, we shall use Bretherton’s result for . This leads to the limit value of as the plug length goes to infinity: . As derived in detail in the Appendix of this paper, the dependence of the curvature on the plug length for finite long plugs is given by: , where functions and are linear combinations of modes and , defined in Eqn.(A 0). By using the constant curvature constraint derived above, we have:

 h∗h∞ = P(ξp)P∞ (2.0) = 1−2.15δ2(ξp)+0.75δ3(ξp).

This is the dependence of on the dimensionless plug length . It can be easily transformed to represent the dependence on via the scaling relation: . Thus the thickness oscillates around with an amplitude that decays as increases, and with a decay length of order . Strictly speaking, the above derivation is only valid in the regime: , where and are arbitrarily small. For illustrative purpose, we extrapolate our analytic results to Teletzke’s short-plug regime and make a comparison between data obtained via two methods in Fig.4. Deviation between two sets of data is obvious for short plugs, where Teletzke’s method is more accurate. For short plugs in the Teletzke’s regime, contributions from all modes are important. Different modes interact nonlinearly, which leads to a much stronger finite-length effect, as shown in Fig.3. Around the data obtained via Teletzke’s method begin to merge with our curve. This is the regime where the prediction in Eqn.(2.0) becomes valid.

In this section, we completed the previously unstudied part of the finite plug speed problem. Together with Teletzke’s method, we can now predict how fast plugs of different length will move in smooth cylindrical tubes. Evidently the speed dependence on the plug length is governed by the solution modes in the thin film region. For long plugs studied in this section, their speeds are the most sensitive to the mode with the longest wavelength: . In smooth tubes, this is much smaller than the tube radius for small capillary numbers of interest. The asymptotic corrections of the finite plug speed is exponentially weak, with a decay length . For very long plugs with , the speed is almost the same as that of the semi-infinite one. As discussed in the next section, tube roughness can modify the spectrum of solution modes in the thin film region dramatically so that may become comparable to or larger than . This is crucial in understanding some interesting experimental phenomena involving the relative motion between finite plugs of different length in non-smooth capillary tubes.

## 3 Relaxation modes in the thin film region with tube roughness

In this section, we investigate the shape of the thin film region of semi-infinite plugs in tubes with roughness. Here tube roughness is defined as surface deformation of the tube wall, whose shape is not cylindrical. For simplicity, the only type of roughness that we will consider in the following discussion is azimuthal roughness. A tube with the azimuthal roughness preserves its translation symmetry in the axial direction, but cross sections are no longer circular. A sketch of such a tube is plotted in Fig.5(A). The cross section of the rough tube is parameterized in the polar coordinate system as: , where is the smooth tube radius and is the polar angle. The roughness function measures the deformation. In order to isolate the regime of small roughness we may express as: . Here is some small parameter much less than the Bretherton’s thin film thickness , and is a function of order unity. The tube roughness breaks rotational symmetry in the azimuthal plane. As a result, the plug shape now depends on both the and the coordinate: . Local lubricating film thickness is defined by the distance between the tube wall and the plug surface: , as shown in Fig.5(B). As in the previous section, we nondimensionalize the system by using the smooth tube radius , the plug speed and the capillary pressure . As derived later in this section, the semi-infinite plug shape in the thin film region can be expressed as:

 h(z,θ)≃H∞+∑kCkψk(θ)exp(νkz) (3.0)

where and are eigenfunctions and eigenvalues governing the relaxation of the plug shape to the cylindrical asymptote . Relaxation mode coefficients are determined by boundary conditions away from the thin film region. In Section 4, we will estimate the order of magnitude of near the transition between the thin film and the remobilization region.

### 3.1 The plug shape equation in rough tubes

As in the case of the smooth tube, the plug shape is dictated by the local carrier fluid flux :

 →q(z,θ)=∫h(z,θ)−e(θ)→v(z,θ,r)dr, (3.0)

where is the carrier fluid velocity. As shown in Fig.5(B), has both an axial component and an azimuthal component . In the limit of lubrication approximation, and with the same boundary conditions as in the smooth case, these fluxes are given by deGennes_2003 (DeGennes, Brochard & Quere 2003):

 qz = −(h+e)−(h+e)33Ca∂zP, qθ = −(h+e)33Ca∂θP, (3.0)

where is the capillary pressure on the plug surface. In our reduced units (Howell (2003)): . Eqn.(3.0) resembles the usual Poiseuille’s law for thin layer fluid fluxes deGennes_2003 (DeGennes, Brochard & Quere 2003). Rearranging terms in the equation of , we get:

 ∂3zh=3Ca(h+e)3⋅(qz+h+e)−∂z(1+∂2θ)h. (3.0)

This is the counterpart of the Teletzke’s plug shape equation, Eqn.(2.0), in rough tubes. Without roughness, is a constant and this equation is sufficient to determine the plug shape . However, now depends on both and coordinates, the flux is no longer constant and another equation describing the variation of is needed to determine the plug shape. This additional equation is the incompressibility constraint: , which leads to the differential equation:

 ∂zqz = ∂θ(h+e)33Ca∂θP (3.0) = −∂θ(h+e)33Ca∂θ[(1+∂2θ)h+∂2zh].

Combining Eqn.(3.0) and (3.0), we get a complete set of differential equations dictating both the plug shape and the flux :

 (3.0)

The first two equations are simply definition equations of variables and . And the parameter is given by: . Evidently, and derivatives are separated to two sides of the equation.

In the smooth tube the plug surface has rotational symmetry, so the last equation of Eqn.(3.0) is simplified as: . Obviously, it implies a constant flux . Defining this flux as and using it in the third equation of Eqn.(3.0), we get:

 ∂3zh=3Ca⋅(h−h∗)h3−∂zh. (3.0)

Finally, we saw in the previous section that in the region of interest varies over a length scale of about . So the slope term in Eqn.(3.0) is negligible compared to the other two terms for small capillary numbers. And the Teletzke’s smooth plug shape equation, Eqn.(2.0), is obtained after dropping .

Returning to the case of nonzero , for convenience, we define a plug-shape vector: , where the superscript indicates the matrix transpose operation. Then we can rewrite the above set of differential equations in a concise matrix form:

 ∂z→ψ = ⎛⎜ ⎜ ⎜ ⎜⎝01000010D−1−(1+∂2θ)0D−1−∂θD∂θ(1+∂2θ)0−∂θD∂θ0⎞⎟ ⎟ ⎟ ⎟⎠→ψ (3.0) ≡ M(h)→ψ.

The nonlinearity of this equation arises because of the dependence of on . In the following discussion, the roughness function is assumed to be of a simple sinusoidal form:

 e(θ)=e0cos(mθ), (3.0)

where are positive integers indicating different roughness modes. The roughness amplitude is assumed to be small in the regime of interest: . As discussed in detail later, Eqn.(3.0) can be reduced to a linear form in the thin film region. More general roughness functions can be expanded in terms of the Fourier series and analyzed correspondingly.

### 3.2 Asymptotic shape of semi-infinite plugs

In this section, we consider the asymptotic shape of the front region of a semi-infinite plug of arbitrarily long length . A similar analysis can easily be applied to the rear region of the semi-infinite plug. One immediate concern with these long plugs is that they may be unstable against fission via the Rayleigh instability deGennes_2003 (DeGennes, Brochard & Quere 2003). This instability is independent of the steady state effect that we discuss here. Indeed, for small capillary numbers, the growth rate of the Rayleigh instability becomes negligible deGennes_2003 (DeGennes, Brochard & Quere 2003).

As with the smooth tube case discussed above, we may fix the origin of coordinate at the transition between the thin film and the remobilization region. Far inside the thin film region where , the front semi-infinite plug restores the translation symmetry in the direction. It implies vanishing derivatives in this limit, i.e. . Evidently, and are zero. The third equation of Eqn.(3.0) leads to the condition: . And the asymptotic plug shape is determined by the last equation of Eqn.(3.0): . One obvious solution is that , where is a constant. This solution corresponds to a cylindrical plug centered inside the tube. There are two other solutions: and , which correspond to cylindrical plugs shifted laterally inside the tube. The translation symmetry far inside the thin film region leads to the cylindrical plug shape with constant capillary pressure on the plug surface and vanishing flux . For convenience, we focus on the centered asymptote. The implication of non-centered solutions will be addressed in the Discussion section.

The centered solution measures the magnitude of average per unit circumference:

 ¯q ≡ 12π∫2π0qz(∞)dθ (3.0) = 12π∫2π0(−H∞−e)dθ = −H∞.

So has the same physical meaning as in the smooth tube case, and dictates the relative motion between the plug and the carrier fluid: . For tubes with small roughness, the difference between and is also very small: . The leading order correction in due to the tube roughness appears in the second order of .

### 3.3 Eigenmodes of relaxation in the thin film region

At finite values, the semi-infinite plug shape deviates from its cylindrical asymptote. The plug-shape vector is rewritten in the following form:

 →ψ = →ψ(∞)+δ→ψ (3.0) ≡ {H∞+δh,δh′,δh′′,−H∞+δqz}T,

where measures the amount of shape correction. In the thin film region, satisfies the relation: , where the superscript indicates different vector components. So the plug shape equation (3.0) can be linearized by using the approximation: , and we get:

 ∂zδ→ψ = ⎛⎜ ⎜ ⎜ ⎜⎝01000010D−11−(1+∂2θ)0D−11−∂θD1∂θ(1+∂2θ)0−∂θD1∂θ0⎞⎟ ⎟ ⎟ ⎟⎠δ→ψ (3.0) ≡ M1δ→ψ,

Evidently, we can solve the above equation via the technique of separation of variables. The general solution of is given by:

 δ→ψ=∑kCk→ψk(θ)exp(νkz), (3.0)

where and are eigenvectors and eigenvalues of the matrix operator , with . Thus in the thin film region, are the eigenmodes dictating the relaxation of the plug surface to its cylindrical asymptote. The speed of relaxation is determined by the decay length of each mode , defined as . The larger the value of , the slower of the relaxation process. Relaxation modes with positive are dominant for the front semi-infinite plug. Likewise, negative ones are dominant for the rear semi-infinite plug.

In the regime of vanishingly small roughness , Eqn.(3.0) can be simplified further as follows:

 ∂zδ→ψ = ⎛⎜ ⎜ ⎜ ⎜⎝01000010D−10−(1+∂2θ)0D−10−D0(∂2θ+∂4θ)0−D0∂2θ0⎞⎟ ⎟ ⎟ ⎟⎠δ→ψ (3.0) = M0δ→ψ,

where the coefficient is now given by: . Evidently the matrix operator has no explicit dependence on the coordinate, so its eigenvectors correspond to different harmonics in the Fourier series of . Due to the form of the roughness function chosen, we may focus on cosine series solutions with the same symmetry: . Here is the harmonic number. Up to the first order in the roughness , we need to analyze the zeroth and first harmonics, with each harmonic leading to four eigenmodes:

 δ→ψ ≃ δ→ψ(0)+δ→ψ(1) (3.0) = 3∑k=0C(0)k→ψ(0)kexp(ν(0)kz)+4∑k=1C(1)k→ψ(1)kcos(mθ)exp(ν(1)kz),

where superscripts (0) and (1) stand for the zeroth and the first harmonics respectively. The zeroth harmonic solution preserves the rotational symmetry and satisfies the equation:

 ∂zδ→ψ(0)=⎛⎜ ⎜ ⎜ ⎜⎝01000010D−10−10D−100000⎞⎟ ⎟ ⎟ ⎟⎠δ→ψ(0). (3.0)

The above matrix operator is a four-by-four matrix and does not depend on the form of the roughness function. So the zeroth harmonic spectra of different roughness functions are the same. For small capillary numbers, four eigenvalues of this matrix are given by: , , , and . There is no plug surface relaxation associated with the mode . It corresponds to the degree of freedom to fix the radius of the cylindrical asymptote as . The other three modes , and correspond to the ones found in the previous section for the smooth tube, up to a difference in the choice of units. For illustrative purpose, the zeroth harmonic spectrum is plotted in Fig.6. Obviously, and are the relaxation modes with the longest decay length in this spectrum.

A major difference between the rough tube and the smooth one comes from the first harmonic solution . It satisfies the equation:

 ∂zδ→ψ(1)=⎛⎜ ⎜ ⎜ ⎜⎝01000010D−10−1+m20D−10D0(m2−m4)0D0m20⎞⎟ ⎟ ⎟ ⎟⎠δ→ψ(1). (3.0)

Four eigenvalues are solutions of the characteristic equation:

 ν4+(1−2m2)ν2−D−10ν+m4−m2=0. (3.0)

The first harmonic spectra with the roughness mode number 5, 6, and 7 are plotted in Fig.7. Compared to the zeroth harmonic spectrum shown in Fig.6, the mode with the longest decay length now lies near the origin and has much longer decay length. For small values of , i.e. long wavelength roughness, becomes comparable to or even larger than the tube radius . Obviously, in the regime of interest. So is the new slowest decay mode in rough tubes after analyzing both the zeroth and the first harmonic spectra. So we will drop the superscript and refer to it as from now on.

To better understand the origin of these slow relaxation modes in tubes with roughness, we propose a following physical explanation of . The remobilization region ends at about distance from the spherical cap. At this transition between the remobilization region and the thin film region, the plug shape has not reached its cylindrical asymptote yet. We assume that the plug surface has a typical height variation of about . A sketch of such a plug cross section is given in Fig.5(B). Non-cylindrical plug shape leads to unbalanced capillary pressure in the azimuthal plane and a non-zero lateral flow . In the thin film region, the plug surface relaxes toward the cylindrical asymptote via this flow . Suppose that the relaxation length scale is given by . It can be estimated as the product between the carrier fluid speed in this region and a relaxation time scale over which smoothes the plug surface. In the following analysis, we will estimate the length scale and obtain its scalings with other system parameters. For this part of discussion only, we will work in the lab units.

The expression of the lateral flux is given by Eqn.(3.0). In the lab units, we have:

 qθ=−(h+e)33μR∂θP. (3.0)

As mentioned above, the capillary pressure is no longer a constant in the azimuthal plane for the non-cylindrical plug shape. The length scale in the direction is defined by the roughness wavelength . So we may estimate derivatives as: . Since the -direction plug surface adjustment mostly happens in the remobilization region, we expect that derivatives are typically much larger than derivatives in the thin film region. The limit of this assumption will be derived later. So the capillary pressure is approximated as follows:

 P = −γ[∂2zh−(1+∂2θ)hR2] (3.0) ∼ γΔh(R/m)2.

By definition, the surface relaxation time scale is given by:

 T≃ΔhR−1∂θqθ∼μ(R/m)4γh3∞. (3.0)

The relaxation length scale is estimated as the product between and : . We may simplify this estimation of further by using the scaling of and get:

 L∗∼RCa⋅m4. (3.0)

For a self-consistent derivation, ’s value is constrained by the approximation used in Eqn.(3.0). This is equivalent to the condition: , which yields the relation: . By using the estimation of obtained above, we get the valid regime of the roughness mode number: . It corresponds to capillary tubes with long wavelength roughness. Obviously is much longer than the smooth remobilization length in this regime of roughness.

The length scale is actually the longest decay length in the eigenvalue spectrum derived previously. It corresponds to the smallest solution of Eqn.(3.0), which goes to zero as approaches zero. For this root, the term dominates the higher power terms of , and Eqn.(3.0) simplifies as follows:

 −D−10⋅λ−1max+m4−m2=0. (3.0)

For an order of magnitude estimation, we may drop the term relative to and get: . This is the same scaling relation as that of , up to a difference in the choice of units. Thus and measure the same length scale in rough tubes. The scaling argument presented above physically explains the origin of large values in the rough tube. This new length scale reflects the slow relaxation of the plug shape via a vanishingly small lateral fluid flux . Historically, people has observed some evidence of this new length scale in non-cylindrical tubes Wong (1995).

## 4 Amplitude of the plug surface bumpiness in the remobilization region

In the previous section, we analyzed the semi-infinite plug shape in the thin film region and found that the plug surface decays to its cylindrical asymptote via different relaxation modes. To proceed further, we must understand how the rough tube wall imposes a non-cylindrical plug shape in the first place. For small roughness, the plug shape is dominated by the zeroth and the first harmonic solutions:

 h(z,θ)≃H∞+3∑k=1C(0)kexp(ν(0)kz)+4∑k=1C(1)kcos(mθ)exp(ν(1)kz), (4.0)

where the coefficients and are determined from the boundary conditions away from the thin film region. In this section, we will treat the front and the rear semi-infinite plugs simultaneously, with the understanding that only the modes with the correct sign of are excited for a single semi-infinite plug. The coordinate origin is fixed at the transition between the thin film and the remobilization region as in Section 2. In that section we have shown that the zeroth harmonic coefficients are of the order of . We will see later in this section that the first harmonic coefficients are of the order of . The exact numerical solutions will be presented in our future work.

In the remobilization region, the governing equation of the plug shape is derived by combining the last two equations in Eqn.(3.0):

 ∂z(h+e)−13Ca[∂z(h+e)3∂z+∂θ(h+e)3∂θ](1+∂2z+∂2θ)h=0. (4.0)

This equation can be simplified further by noticing the fact that in the remobilization region, varies over about in the direction and about in the direction. For long wavelength roughness with , the following estimation of derivatives is valid : . In other words, in the remobilization region the plug shape varies much more rapidly in the direction than in the direction. So we may keep only derivatives in Eqn.(4.0), and get:

 (4.0)

Bretherton’s smooth tube solution satisfies the above differential equation with . In the regime of small roughness, Eqn.(4.0) can be solved perturbatively by expanding around the smooth tube solution : . Here and are the leading order plug shape correction attributed to the tube roughness. For vanishingly small roughness, the relation holds for arbitrary . Since we are interested in the coefficients of the first harmonic solution, we will focus on the shape correction for the rest of this section. Keeping only the lowest order approximation, satisfies the following equation:

 (∂zh30∂3z)u(1)+(3h20∂3zh0−3Ca)∂zu(1)+[∂z(3h20∂3zh0)](u(1)+e0)=0. (4.0)

This is an inhomogeneous differential equation. The correct boundary conditions of are that it goes to zero as approaches the spherical cap, and matches the relaxation modes in Eqn.(4.0) as approaches the thin film region. The general solution of Eqn.(4.0) takes the form: . Here is a dimensionless function satisfying the homogeneous equation and goes to 1 as approaches the spherical cap. By using the rescalings as in the smooth tube: , and , we get the following equation of :

 (∂ξη30∂3ξ)f+(3η20∂3ξη0−3)∂ξf+[∂ξ(3η20∂3ξη0)]f=0. (4.0)

After the rescaling both and its derivatives are of the order of unity in the remobilization region, so all coefficients in Eqn.(4.0) are of the order of unity. Together with the boundary condition of near the spherical cap, the solution must be of the order of unity in the remobilization region. By matching with relaxation modes in the thin film region, we get an order of magnitude estimation of the first harmonic coefficients: . Thus the semi-infinite plug shape in the thin film region can be rewritten as follows:

 h(z,θ)≃H∞(1+3∑k=1~C(0)kexp(ν(0)kz))+4∑k=1~C(1)ke(θ)exp(ν(1)kz), (4.0)

where all coefficients and are of the order of unity. The above estimation is valid for both the front and the rear semi-infinite plug.

## 5 Finite plug motion in rough tubes

In this section, we qualitatively study the dependence of the finite plug motion on the plug length in the rough tube, following the logic of Section 2. For the finite plugs, relaxation modes of both signs of are excited in the thin film region. Far inside the thin film region where the plug shape is almost cylindrical, all the other relaxation modes are vanishingly small compared to the slowest decay mode . Thus we may focus on the effect of . By defining a coordinate whose origin is at some position far inside the thin film region, the plug shape is approximated as follows:

 h(z,θ)≃H∞+~Ce(θ)exp[(~z−ℓ)/λmax], (5.0)

where . The coefficient is of the order of unity as derived in the previous section. Here the variable measures the half length of the finite plug. For plugs used in microfluidic experiments, is usually comparable to or larger than the tube radius Ying_2008 (Ismagilov & Ying). As mentioned previously, the relative motion between the plug and the carrier fluid is governed by the azimuthal average of :

 U/V = 1−1π∫2π0qzdθ (5.0) = 1−2¯q.

Far inside the thin film region, the plug shape is almost cylindrical, which leads to the following estimation of the carrier fluid fluxes: , and . By using the approximate plug shape in Eqn.(5.0), we have: . Evidently, the average flux vanishes in this linear approximation. Thus any effect on appears in the second or higher order of the roughness: . The asymptotic dependence of the plug speed on the plug length is exponentially small for long plugs, and we have:

 U/V≃1+2H∞[1+~Be20H∞exp(−2ℓ/λmax)]. (5.0)

The above form of length dependence implies an interesting resonance effect in the rough tube. To illustrate this effect, we may estimate the speed difference between two plugs with length and . We also assume that the length difference between these two plugs, defined as , is small such that . By using Eqn.(5.0), we derive their speed difference as follows:

 (U2−U1)/V ≃ 2~Be20[exp(−2ℓ2/λmax)−exp(−2ℓ1/λmax)] (5.0) ≃ 2~Be20(−Δℓℓ1)⋅[exp(−2ℓ1/λmax)⋅(2ℓ1/λmax)],

The function of in the last square bracket above is not monotonic. It has a maximum at . For plugs used in microfluidics measurements, we have: . Thus the speed difference of these plugs is maximized in rough tubes with . By using the scaling of derived in Section 3, we get the roughness mode of resonance: .

In the rest of this section, we will use an example to illustrate how the tube roughness enhances the sensitivity of the plug motion to the plug length. Suppose that we have two finite plugs of length: and . Strictly speaking, Eqn.(5.0) is valid for the regime of vanishingly small roughness. However, we expect this estimation to be at least qualitatively correct for tubes with large roughness Thus for computation simplicity, we may set the coefficient to be 1. In the smooth tube case, and . And we get the speed difference between these two plugs less than with . Thus the finite plug length effect is negligible in the smooth tube. On the other hand for tubes with long wavelength roughness, as shown in Fig.6, we have parameter values: and . The plug speed difference can be up to a few percent, which is significantly larger than that in the smooth tube.

## 6 Discussion

In the previous sections we have argued that roughness in a microchannel can qualitatively alter capillary motion in the channel. We demonstrated the new effects in the context of an air plug being pushed through the channel. Roughness creates new modes of relaxation of the fluid interface. These modes make major changes in the Bretherton-Teletzke mechanism that sets the film thickness and governs the plug’s speed. In particular, they alter the nature of the the remobilization region that dictates the thickness of the liquid film between the plug and the wall. These modes can thus alter the effect of other variations in the system, such as the addition of surfactants or the replacement of the air plug by a viscous fluid.

Our focus above has been to consider the effect of plug length on its speed. This basic effect has practical consequences as noted in the Introduction. Microfluidic devices for chemical processing create long trains of droplets that need to remain separated over long distances. However in practice these droplets can have slight differences of speed. This can make droplets merge, thus disrupting the desired behavior. Since the length of these droplets is not precisely controlled, the length differences could in principle be responsible for the undesired speed differences. However, such a length effect is incompatible with the classic Bretherton explanation for the speed. We showed above that roughness changes this conclusion qualitatively.

This effect complements previously-considered effects of channel shape. It differs from the roughness effect identified by Homsy et al. Krechetnikov_2005 (Krechetnikov & Homsy 2005), in which roughness alters the effective boundary condition at the wall. We showed here that roughness can also have long-range effects on the thin film behavior. These effects do not require the qualitative changes in channel shape found, e.g., for square channels (Wong (1995); Wong (1995)). Remarkably, even high-wavenumber roughness can alter the long-range effects if the capillary number is sufficiently small, as shown in Section 5.

Roughness may be viewed as a relevant perturbation to the plug viewed as a dynamical system. The fast relaxation to the thin film within a narrow remobilization region seen in a Bretherton plug depends crucially on its assumed azimuthal symmetry. When this symmetry is broken by roughness, the relaxation takes on a qualitatively new behavior. In this sense the symmetric plug is unstable against the generic un-symmetric behavior discussed above. For this reason it appears important to consider roughness when estimating any new aspect of plug behavior.

Though our work suggests that roughness is important, our quantitative exploration of this importance has been very limited. First, we considered only azimuthal roughness, ignoring any dependence of the wall position. We believe the nonlocal effects identified above are not qualitatively altered for generic roughness, we have not addressed this issue explicitly. We consider only perturbative effects of roughness in lowest order. Thus we have no systematic calculation for the speed, which requires the next order of approximation. Our analysis of the remobilization region was limited to a scaling discussion, with no systematic results. We didn’t carry out more detailed calculations because the perturbation approach itself has serious flaws. It becomes qualitatively inadequate for low capillary number, where the unperturbed film thickness becomes smaller than any roughness. In this regime the interface is mainly supported by asperities where the film thickness is much smaller than its average. In order to get quantitative information on how roughness and length affect speed, we plan a numerical integration of Eqn.(3.0).

Our analysis unearthed a peculiar behavior for the lowest-wavenumber perturbation, namely . This mode amounts to a lateral displacement of the plug in a circular channel. Remarkably this perturbation does not relax since it does not alter the circular shape of the plug. Centering effects of particles in tubes have been discussed Brenner_1963 (Brenner 1963), but we are not aware of an analysis suited to Bretherton plugs.

The roughness effects shown here suggest a new regime of interest in microfluidic capillary flow. In this rough-channel regime, the symmetry of the Bretherton plugs is lost, but the fluid coating the interface is still a thin film in lubrication flow. Thus these channels are distinct from, e.g., the square channels discussed by Wong (Wong (1995)), where the lubrication approximation is inadequate. In addition to capillary number and the plug aspect ratio, two dimensionless variables characterize a system in this regime. The first is the roughness wavelength relative to the Bretherton film thickness . The second is the amplitude relative to this thickness. We have noted above the special effects to be explored when the decay length is comparable to the plug length. We have also suggested an asperity-dominated regime to be expected when the dimensionless amplitude becomes large. Naturally the chief experimental impact of these effects comes when the plug is replaced by a droplet of viscous liquid. The effects anticipated above should all have counterparts for such liquid droplets.

Our predictions, though qualitative, are subject to experimental tests. An obvious test is to create tubes with controlled corrugated profiles and verify that when the roughness amplitude is comparable to the predicted Bretherton film thickness, the speed increment changes significantly relative to that of a smooth tube with the same cross section. One can also create corrugations with a wavelength where sensitivity to plug length is expected and look for significant variations in when the roughness amplitude is comparable to the Bretherton thickness. Such tests could be done at the microfluidic scale or at larger scales.

More primitive observations using microfluidic channels are also suggested by our analysis. One may readily prepare circular tubes with 200 micron cross section and controllable flows with capillary numbers of order Ying_2008 (Ismagilov & Ying). For such flows, the predicted Bretherton thickness is 60 nanometers, so that roughness amplitudes of this order should be significant. AFM measurements on these tubes show roughness with of order and amplitude of order 50 nanometers. This amplitude is evidently large on the scale of interest. Since the roughness varies randomly along the tube, one expects statistical fluctuations in the speed along the tube. Some of these fluctuations can be due to differences in cross sectional area, so that they cause the fluid speed to change with position. In addition one expects variations in owing to our mechanism. The indication that these variations are caused by roughness is that is determined by position. Thus e.g. large values occur at particular places along the tube.

## 7 Conclusion

The subtle interplay of nonlocal forces that controls a microfluidic droplet is a classic example of the power of simple capillary flow to produce nontrivial structures. It appears from our analysis that weak corrugation of the guiding channel can add unexpected richness to this picture. The addition of transient effects, binary interaction between droplets and surfactants seems likely to add even more to this richness. Exploring these effects appears promising both for microfluidic applications and for basic understanding of what capillary flow can do.

The authors are grateful to Prof. Rustem Ismagilov and Dr. Ying Liu, whose microfluidic experiments inspired this work. This work was supported in part by the National Science Foundation’s MRSEC Program under Award Number DMR 0820054. And by a grant from the US-Israel Binational Science Foundation.

## A

In this appendix, we derive the dependence of the rescaled curvature on the finite plug length . All quantities are in the dimensionless units defined by the rescalings in Section 2. To set the stage for our treatment of finite plugs, we analyze the behavior of for the semi-infinite plug treated by Bretherton. We proceed from a semi-infinite plug with the front cap. Far inside the thin-film region, and we may use the general solution in Eqn.(2.0). In the front region of the thin film, the contributions from and are arbitrarily small relative to that from . So the plug shape takes the form:

 η=1+C1ψ1(ξ). (A 0)

Consistency with this form implies: , where and stand for the second and the first derivative of with respect to respectively. We may thus find the profile by integrating forward from an arbitrary origin to , starting from an initial value of very close to . The asymptotic curvature is necessarily finite and independent of axial translation in the small-amplitude limit. Specifically, . It is convenient to define as the for which the contribution in Eqn.(A 0) extrapolates to 1, i.e. . Thus the thin film region of Eqn.(A 0) is where , and . Then from the profile thus determined numerically, one finds that Teletzke’s point has position: . Evidently, is independent of .

A similar analysis can be applied to a semi-infinite plug with the rear cap. Here, the contribution is negligible in the rear region of the thin film. So the plug shape takes the form:

 η=1+C2ψ2(ξ)+C3ψ3(ξ). (A 0)

As with the front region, we may find the rear profile by integrating backward from the thin-film region where . Far inside the thin-film region the and contributions are necessarily small, but their ratio varies sinusoidally with position, with a wavelength comparable to the remobilization length . As above, we define a such that . Then in the thin film region takes the form:

 η=1+e−(ξ−ξR)/2⋅[cos(√3(ξ−ξR)/2)+C3sin(√3(ξ