# [

###### 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:

(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):

(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:

(2.0) | |||||

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):

(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:

(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:

(2.0) |

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

(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:

(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:

(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:

(2.0) | |||||

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:

(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 :

(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):

(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:

(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:

(3.0) | |||||

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:

(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:

(3.0) | |||||

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:

(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:

(3.0) | |||||

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) | |||||

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:

(3.0) | |||||

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

(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:

(3.0) | |||||

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:

(3.0) | |||||

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:

(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:

(3.0) |

Four eigenvalues are solutions of the characteristic equation:

(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:

(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:

(3.0) | |||||

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

(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:

(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:

(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:

(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):

(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:

(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 :

(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:

(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:

(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 :

(5.0) | |||||

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:

(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:

(5.0) | |||||

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:

(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:

(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: