EPJ manuscript No.
(will be inserted by the editor)
Simulations of two-dimensional foam
rheology: localization in linear Couette flow
and the interaction of settling discs
A. Wyn1, I.T. Davies1, and S.J. Cox1,2 a
1 Institute of Mathematical and Physical Sciences, Aberystwyth University, Ceredigion
SY23 3BZ, Wales, UK
2 UMR 5588 ? Laboratoire Spectrom?etrie Physique, B.P. 87, F-38402 St. Martin d?H`eres
Cedex, France.
January 7, 2008
Abstract. Surface Evolver simulations of flowing two-dimensional
foams are described. These are used for two purposes. Firstly, to
extract the location of the T1s, the changes in bubble topology that
occur during plastic flow. It is shown that the T1s are localized in
space, becoming more so as the polydispersity of the foam decreases.
Secondly, the sedimentation of two circular discs through a foam un-
der gravity is studied. If the discs are sufficiently close, they begin to
interact and one moves behind the other during their descent.
PACS. 47.57.Bc Foams and Emulsions ? 83.80.Iz Emulsions and
Foams
1 Introduction
Liquid foams are familiar from domes-
tic use and important in industrial ap-
plications including ore-separation and
enhanced oil recovery [1, 2]. They are
elasto-visco-plastic complex fluids with
a highly nonlinear response to applied
forces: at low strain they deform elas-
tically, like a solid, while above a yield
stress they flow like a viscous liquid [3].
Of all complex fluids, liquid foams pro-
vide one of the most experimentally ac-
cessible systems for study, since bubbles
are objects that can have millimetric di-
mensions. Moreover, Plateau?s laws [4]
meanthat the internalstructureofa foam
is well understood, at least at the level of
the network of films. Foams thus provide
a prototypical complex fluid.
a email: foams@aber.ac.uk
However, given the degree of disorder
within the foam structure and the com-
plex response, it makes sense to first con-
sider two-dimensional (2D) foams, such
as can be made by squeezing a foam be-
tween parallel glass plates until it con-
sistsof a singlelayerofbubbles [5].Other
realizationsofa 2D foamincludethe bub-
ble raft of Bragg and Nye [6], promoted
recently by Dennin and co-workers [7, 8,
9], and the hybrid method of Cyril Stan-
ley Smith [10] and Fortes and co-workers
[11, 12]. A theme of current research is
exploring the different responses of each
of these experimental setups [13, 14], re-
quiring an understanding in particular of
the effects of liquid content.
The mathematical idealization of a
two-dimensional foam is, however, clear:
a dry 2D foam at equilibrium consists
of bubbles with fixed areas surrounded
2 A. Wyn et al.: 2D Foam Rheology
by films that are circular arcs meeting
threefoldatanglesof120? (figure1).These
rules are consequences of minimization
of energy [15], which is in this case the
total film length multiplied by surface
tension. This model, and various approx-
imations to it, have long been used for
simulation [16, 17, 18, 19, 20, 21]. Here,
we use the Surface Evolver software [22]
to simulatewith highaccuracyfoamscon-
sisting of many hundreds of bubbles, to
predict the plastic response of 2D foams.
At low strain a foam responds as an
elastic medium. That is, the shear stress,
given as a sum of surface tension contri-
butions in the films [23, 24], increaseslin-
early with strain. As the strain increases,
the foam begins to yield and bubbles be-
gin to slide past each other in plastic
events known as T1 topological changes
(figure 1) [25]. These occur when a film
shrinks to zero length and a fourfold ver-
tex is formed. Such a vertex is unsta-
ble, andimmediately dissociatesinto two
threefoldverticeswiththe connecting film
now perpendicular to the vanishing one.
In the notation of Wang et al. [9], two
bubbles thatwerenearest-neighboursbe-
come next-nearest neighbours, and vice
versa. Each T1 event contributes a drop
in both total film length and stress.
LocalizationofT1 events,also referred
to as shear banding, has been described
in experiments in an annular wide-gap
Couette viscometer [26]. After an initial
transient, the majority of T1 events oc-
cur close to the inner moving wall. Sim-
ilar results have been found in simula-
tions [27]. In linear Couette shear be-
tween parallel side-walls, there is also lo-
calization of T1s. In this geometry, since
the shear stress should be homogeneous
there is no preferred location for the lo-
calized region based upon the boundary
conditions,confirmedbyPottsmodel [28]
and Surface Evolver [19] simulations.
The presence of phenomena such as
shear localization presents a non-trivial
obstacle to the development of contin-
uum models for foam rheology [29, 30].
Approaches such as the theory of shear
transformationzones[31] alsorequirethat
the local dynamics is first understood.
In ?2 we predict the width of the lo-
calized region in linear Couette shear,
and its dependence on the area disorder
of the foam. This is, in effect, a predic-
tion of the degree to which the foam is
fluidized under shear. In the limit of zero
area disorder - a monodisperse foam -
T1 events tend to occur in a very nar-
row band and shear-induced crystalliza-
tion is evident. We show here that mak-
ing a foam more polydisperse widens the
localized region and can thus reduce the
amount of static foam present.
We characterizethepolydispersity,or
volumetricdisorder,ofa foamby the sec-
ond moment of the distribution of bub-
ble areas A:
?2(A) =
angbracketleftbigg(A??A?)2
?A?2
angbracketrightbigg
(1)
where??denotesanaverageoverthe whole
foam. In contrast to the disorder in the
number of sidesnofeachbubble, ?2(n) =
?(n ? 6)2?, which varies in time due to
T1s, the area disorder is fixed in each
of our simulations. That is, we exclude
inter-bubblegasdiffusion(coarsening)and
film collapse. The elastic response of a
foam is characterized by the shear mod-
ulus [32], which decreases by up to 10%
at both high topological and high volu-
metric disorder [23].
In a further effort to understand the
response of a foam, we consider a ge-
ometry in which we know approximately
where the T1s will occur, and ask what
is the interaction between the foam flow
and an embedded object (?3).
A number of authors have studied
the flow of a 2D foam past a fixed object,
with both experiments [12, 33, 34, 35]
and simulations [36, 35]. The drag and
lift forces on the object are due to a
number of contributions. At low veloc-
ity the dominant ones are the force from
the tensions in the films attached to the
object and the pressures in the bubbles
that touch it. For a circular object in
the centre of a channel the drag force
increases with object diameter [36, 37]
and decreaseswith increasing liquid frac-
tion [35]. For asymmetric objects such as
an aerofoils [12], and for circular objects
A. Wyn et al.: 2D Foam Rheology 3
Fig. 1. An ideal 2D foam consists of films that are circular arcs, meeting at 120?. The
sequence of images shows a T1 topological change. The intermediate step with the fourfold
vertex is energetically unstable. The orientation of a T1 is characterized by the lines joining
bubble centres adjacent to the deleted and created films, marked by dashed lines.
close to one of the walls of the channel,
there is in addition a lift force. For an
aerofoil this is a negative lift [12], while
for the circularobjectit pointsawayfrom
the wall [36].
Related work in three dimensions is
concerned with single spheres, with a di-
ameter larger than the bubble size, be-
ing pulled [38, 39, 40] or dropped [41, 42]
through a foam.
In ?3 we examine the the interaction
betweentwo circulardiscsfalling through
a foam under their own weight. We aim
to answer questions such as the condi-
tionsunder whichtwo objectsfalling thr-
ough a foam are mutually attracted or
repelled, as has been done for a num-
ber of viscoelastic fluids [43, 44, 45]. The
answers will give guidance in determin-
ing the effect of a wake in a discretized
elasto-plastic fluid.
To guide our intuition we recall work
on the flow of a foam past an ellipse: Dol-
let et al. [46] found that the only sta-
ble orientation of an ellipse was with its
long axis parallel to the direction of flow.
This is a feature of elastic fluids [47]. Is it
therefore the case that the plastic events
arenot significant in determining this as-
pect of the foam response, and we can
treat it as an elastic liquid?
2 The localization of
topological changes in linear
Couette shear
We describe simulations of the slow lin-
ear Couette shear of a 2D foam confined
betweenparallelwalls(figure2).The area
dispersity ?2(A) is varied to determine
how the width of the localized region de-
pends upon this parameter.
2.1 Method
We usethe SurfaceEvolver[22]in a mode
in which each film is represented as a cir-
cular arc. We use foams of N = 1120
bubbles, in a channel of width W = 3.2
and length L = 8.0, giving an average
bubble size of?A? = 0.0229and about 21
bubbles between the walls. The value of
surface tension, which should be thought
of as a line tension with units of energy
per unit length, is taken equal to one
throughout. A realistic foam structure is
found byminimizing the totalfilmlength
subject to the prescribed bubble areas.
For some parameter values we doubled
the number of bubbles in the x? direc-
tion to ensure that the results were not
affected by the possibility of system-wide
avalanchesofT1 events(data not shown).
The simulation procedure is as fol-
lows. A Voronoi construction [48] is first
used to generate a fully periodic tessel-
lation of the plane. Bubbles at the top
and bottom are sequentially deleted un-
til the required number of bubbles re-
mains. This structure is imported into
the Surface Evolver and peripheral films
constrained to one of the two side-walls,
a distance W apart. New bubble areas
are determined randomly from a Weibull
distribution:
f(A;?,?) = ??
parenleftbiggA
?
parenrightbigg??1
e?(A/?)?. (2)
4 A. Wyn et al.: 2D Foam Rheology
L
Wy
x
Fig. 2. An example of the foams used to simulate linear Couette shear in a channel, with
?2(A) = 0.175. The channel is periodic in the x direction. To shear the foam the side-wall
at the top of the image is moved to the right (positive x-direction) in small increments d?.
Those films that meet the side-walls have their ends pinned to the wall (no-slip condition).
The parameter ? > 1 determines the
area dispersity and the parameter ? is
chosen as ? = 1.115?A?(so that the peak
of the distribution is close to A = ?A?).
The second moment of this distribution
is
?2(A) =
?
parenleftbigg
1 + 2?
parenrightbigg
?
parenleftbigg
1 + 1?
parenrightbigg2 ?1 (3)
where ? is the Gamma function. The
limit ? ? ? corresponds to a monodis-
perse foam (?2(A) = 0); decreasing ?
leads to increasingly polydisperse foams.
Note that since our foam sample is finite,
the value of ? chosen for each simula-
tion can lead to slightly different values
of ?2(A).
The initial structure for each simula-
tion is found by reducing the total film
length to a local minimum. During this
minimization T1s are triggered by delet-
ing each film that shrinks below a cer-
tain length lc and allowing a new films to
form to complete the process. The criti-
cal length lc is a measure of liquid frac-
tion ? [36], but we keep it small enough
here (lc = 0.005throughout, correspond-
ing to ? = 2.6?10?4) that it should not
affect the results [27].
To shear the foam, a small step in
strain is applied by moving one of the
confining walls a distance d?, moving all
vertices affinely, and then reducing the
film length to a minimum. In this way,
the foam passes through a sequence of
equilibrium states, appropriate to an ap-
plied strain with strain rate much lower
than the rate of equilibration after T1s.
The value d? = 0.0078was used through-
out, and the foam sheared up to a total
strain of at least ? = 5.
2.2 Position of T1s
Figure 3 shows the T1 positions in foams
at three representative values of ?2(A).
Plotting the y position of a T1 against
strain, or the number of iterations, indi-
cates that for each value of ?2(A) there
is an initial transient that lasts up to ap-
proximately unit strain. In the monodis-
perse case shown in figure 3(a), the T1s
mostly occur close to the moving wall.
Although thisis notthe casefor allof our
simulations in the monodisperse limit,
we found that monodisperse foams usu-
ally localize near one of the walls (see fig-
ure4(a)). Asthe polydispersityincreases,
the width ofthe localizedregionincreases
and it often occurs further from the walls
(figure 4(a)). At large values of polydis-
persity (small ?) the T1s occur almost
throughout the channel. For the inter-
mediate value of ?2(A) shown in figure
3(b), we note that plotting the data on y
vs x axes illustrates a slight undulation
in the localized region.
We measure the width lw of the lo-
calized region as follows. After the tran-
A. Wyn et al.: 2D Foam Rheology 5
(a)
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  1  2  3  4  5  6  7  8x position
yp
os
iti
on
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  100  200  300  400  500  600
 0  1  2  3  4  5
yp
os
iti
on
Strain
Iterations
(b)
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  1  2  3  4  5  6  7  8x position
yp
os
iti
on
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  100  200  300  400  500  600
 0  1  2  3  4  5
yp
os
iti
on
Strain
Iterations
(c)
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  1  2  3  4  5  6  7  8x position
yp
os
iti
on
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  100  200  300  400  500  600
 0  1  2  3  4  5
yp
os
iti
on
Strain
Iterations
Fig. 3. Locations of topological changes in simulations at three values of ?2(A), shown on
y vs x axes (left-hand column) and on y vs ? (strain or number of iterations) axes (right
hand column). (a) Monodisperse (?2(A) = 0,? ? ?). (b) Moderately polydisperse (the
foam in figure 2, ?2(A) = 0.175). (c) Highly polydisperse (?2(A) = 0.561). The vertical
bars show the position and width of the localized region after the transient; they are
centred at the average y position of the T1s, averaged in groups of Ns = 100 iterations,
and their total height encompasses the foam width within which 90% of T1s occur.
sient, taken to be the first 200 iterations,
we find the mean y position of the T1s
in bins of Ns iterations. The localisa-
tion width wl is the interval in y position
within which 90% of T1s are found. We
find that Ns = 100 gives the best mea-
sure of wl, that is, it balances the need
to have many points in each bin with the
desire to accurately reflect the width of
the evolving localization.
Figure 4(b) shows the increase of lo-
calisation width wl/W with disorder. It
is clear that at high disorder, and for
narrow (low W foams), the localized re-
gion may encompass the whole foam. At
low disorder, crystallization is more fre-
quent, and localization usually occurs in
a narrow band. We find the following
rule of thumb:
wl
W ?
radicalbig
?2(A). (4)
On a few occasions (data not shown) we
found that two narrow localized regions
persisted up to strains of about 3.
6 A. Wyn et al.: 2D Foam Rheology
(a)
 0
 0.5
 1
 1.5
 2
 2.5
 3
 3.5
 0  0.1  0.2  0.3  0.4  0.5  0.6
wradicalbig
?2(A)
yp
os
iti
on
(b)
 0.001
 0.01
 0.1
 1
 1e-005  0.0001  0.001  0.01  0.1  1
 0.1
 1
 10
wl
/W
w/
radicalbig
?A
?
?2(A)
Fig. 4. (a) The position of the centre of the localized region, given as the average of
the mean y position of the T1s in each group of Ns = 100 iterations after the transient.
Each point corresponds to one full simulation. As the foam becomes less disordered, the
T1s localize closer to the walls, but without showing a preference for either the stationary
or moving wall. (b) The width of the localized region wl/W, given as the width of foam
within which 90% of T1s occur after the transient, increases with area disorder ?2(A).
Note the log axes. The right-hand axis indicates the localization width in terms of bubble
diameters. The straight line has slope one-half: wl/W = 0.8
radicalbig
?2(A).
2.3 Angular dependence of T1s
Wang et al. [9] recorded the orientation
of the disappearing film and the newly-
createdfilmduring eachT1 eventin their
experiments on a sheared bubble raft.
This is done by drawing a line between
the centres of the bubbles neighbouring
each film, in consecutive images brack-
eting the T1 event, as in figure 1. The
distributions of angles show peaks at an-
gles of about 45? and 135? respectively
to the walls of the channel. The height of
the peak increases with shear-rate, and
at low shear rate a small ?knee? appears
at about 90?, i.e. parallel to the direction
of shear.
From our simulations we can extract
the same data, in the limit of low shear-
rate. Figure 5 shows these distributions,
with data obtained after the transient in
each simulation. We find the same peak
for the filmsthat disappear,but the most
probable orientation for new films is at
30?. Moreover, no detail is seen around
90?. We therefore believe that this dis-
crepancy is due to the high liquid con-
tent in the bubble raft experiments, un-
attainable with the methods described
here.
Although we do not probe the effect
of shear-rate on the height of the peak
(as could be done with the Viscous Froth
Model [49]), we see that it is strongly
affected by polydispersity: monodisperse
foams exhibit a much higher peak. This
is perhaps suggestiveof hexagonal order-
ing, although we did not test this explic-
itly (e.g. by measuring ?2(n), with n the
number of sides of a bubble).
2.4 Other predictions
Our ultimate goal is to predict the re-
gion of a foam where localization will oc-
cur from the initial structure. Given that
there are many T1s throughout the foam
during the transient, this is difficult. To
begin, we seek ways of characterizing the
foam structure, and then following these
characterizations through each simula-
tion to detect robust changes when the
foam localizes. In particular, these meth-
ods should be able extractable from a
single image of a foam at a given time,
rather thanrequiringthe trackingofbub-
bles from video analysis.
Our structural measures extract the
distribution ofvariousquantitiesasa func-
tion of y. In particular, we seek to char-
acterize the area distribution of a foam.
Five methods were investigated:
A. Wyn et al.: 2D Foam Rheology 7
(a)
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0.12
 0.14
 0  20  40  60  80  100  120  140  160  180
?2(?) = 0.43
?2(?) = 0.13
?2(?) = 0.076
?2(?) = 0.051
?2(?) = 0.0065
?2(?) = 0.0017
?2(?) = 0.00016
Orientation (degrees)
PD
F
(b)
 0
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0  20  40  60  80  100  120  140  160  180
?2(?) = 0.43
?2(?) = 0.13
?2(?) = 0.076
?2(?) = 0.051
?2(?) = 0.0065
?2(?) = 0.0017
?2(?) = 0.00016
Orientation (degrees)
PD
F
Fig. 5. The distribution of the angle made with the x axis by a line joining the centres
of the two bubbles neighbouring a film that (a) disappears and (b) is created during a T1
event. The former peaks at around 30? and the latter at 135?.
1. wefindthe centre(xc,yc) ofeachbub-
ble by averaging the coordinates of
itsvertices,and calculatethe histogram
of yc in 20 bins;
2. for each value of y0 ? [0,W] we cal-
culate the average area Ay of those
bubbles that intersect the line y =
y0;
3. wefindthe centre(xc,yc) ofeachbub-
ble by averaging the coordinates of
its vertices, and calculate the local
foam disorder ?2(A), from (1), in 20
bins based upon yc;
4. for each value of y0 ? [0,W] we cal-
culate the average length Ly of the
line y = y0 that is covered by each
bubble, sometimes referred to as the
linear intercept method [50];
5. we calculate the texture tensor [51]
based upon bubble centres in 20 bins.
Figure 6 shows these structural mea-
sures for the foam in figures 2 and 3(b)
at the beginning and end of the simu-
lation, i.e. before and after localization
has occurred. All measures are uniform
at the beginning of the simulation, sug-
gesting that it is not possible to predict
where a foam will localize.
Method 4, Ly, is the only 1D struc-
tural measure to give a clear indication
that the foam has localized. Note that
it does have the disadvantage that large
fluctuationsareobservedfor orderedstruc-
tures.
The texture tensor, method 5, is the
tensorial equivalent ofLy. It is more sen-
sitive than the latter but more difficult
to extract from the data. We shall return
to it in future work, in addition to other
measures such as the local stress in the
foam.
3 Movement of discs through a
foam
We describe simulations that probe the
interaction between macroscopic objects
falling through a foam. The system un-
der study consists of two circular discs,
whosediametersareequalandlargerthan
the bubble size. Recall that an elliptical
object tries to align itself with the foam
flow [46]. Here, we show that when the
discs are sufficiently close, one of them
moves behind the other.
Thediscs? motioniscommencedfrom
a positionnear the top of a monodisperse
foam. They descend under the action of
three forces, defined in figure 7: (i) grav-
ity; (ii) the resultant tension force Fn
due to the network of films that contact
each obstacle; (iii) the resultant pressure
force Fp from the bubbles that touch
eachobstacle.The filmstouching the discs
are not uniformly distributed around the
circumference: as figure 8(a) shows, they
bunch up behind the obstacle. It is this
inhomogeneity that leads to the resul-
tant network and pressure forces.
For each disc the network force is a
sum over those films j that touch the
disc. Each film meets the disc perpen-
dicularly, and makes an angle ? with the
8 A. Wyn et al.: 2D Foam Rheology
 0
 0.5
 1
 1.5
 2
 2.5
 3
 3.5 1 2 3 4 5
yp
os
iti
on
Measure
Fig. 6. The five measures of foam structure are shown from left to right for the foam
in figure 3(b), both at the beginning and the end of the simulation. It is clear that the
foam is initially fairly isotropic, and that none of the methods 1 to 3 indicate that any
structural change has occurred during the evolution. In contrast, the peak around y ? 1.7
in Ly (method 4) corresponds to the region where the T1s have localized. The tensorial
texture tensor (method 5) also shows this feature, as well as indicating that the bubbles
are more closely aligned (cf. figure 5).
y direction [36]. Thus:
Fn = ?
summationdisplay
filmsj
(sin?j,cos?j). (5)
The pressure force is a sum over all bub-
bles touching the obstacle:
Fp =
summationdisplay
bubblesk
pklk(sin?k,cos?k). (6)
where pk is the pressure inside the bub-
ble, lk is the length of contact, and ?k
is the angle that the inward normal at
the midpoint of the line of contact makes
with the y-direction.
3.1 Method
We perform quasi-static simulations as
described in?2.1, using a foam with N =
727bubbles, width W = 0.792andlength
L = 1. The bubble size is therefore A ?
1 ? 10?3 (it shrinks slightly in propor-
tional to the disc size, since the total
area of the foam and two disc system
is constant). The cut-off length for T1
events is lc = 0.002, corresponding to a
dry foam with ? ? 1?10?3. The channel
is periodic in the y direction, parallel to
the direction of gravity, and we stop the
simulations before either of the discs re-
turns to the top of the foam. Films that
meet the side-walls have that end fixed
throughout each simulation (no-slip con-
dition). The ends of films that touch the
discs are free to slide, so as to be able
to make an equilibrium 90? angle (slip
condition).
Our dimensionless units are chosen
so that the line tension ? has the value
1. We choose the discs to have equal ar-
eas in the range 2A to 7A and equal
weights of w = 10 irrespective of their
size. We first ensured that this value of
weight is large enough that the discs are
not brought to a halt by the opposing
forces due to film tensions and bubble
pressures.
Two starting configurations are cho-
sen, as shown in figure 8. The disc cen-
tres are initially separated by a distance
di, either horizontally (i = 1) or verti-
cally (i = 2). In the first configuration,
there is a possibility of a small lift force
due to, and perpendicular to, the walls
[36], acting to push the discs together;
we quantify this in ?3.2 and show that
A. Wyn et al.: 2D Foam Rheology 9
Fig. 7. The position of each disc evolves under the gravitational, tension and pressure
forces shown.
the discs are far enough from the wall
that it is negligible. The second case has
the advantage that the lift force on each
disc should be zero on average.
The simulations proceed as follows.
A foam containing the two discs in their
starting positions is relaxed to equilib-
rium, using the method described in [35].
At each iteration, the resultant forces on
each disc in the x and y directions are
calculated and the disc centres moved
according to
?x = ?(Fnx +Fpx),
?y = ?parenleftbigFny +Fpy +wparenrightbig, (7)
where the subscripts denote the x and y
components of the forces. The parame-
ter ? = 5 ? 10?4 measures how far the
centres are moved at each iteration. The
foam perimeter (energy) is then brought
back to a local minimum with the discs
fixed.This comprisesoneiteration,which
is repeated until a disc reaches the bot-
tom of the simulation cell. Representa-
tive paths of two discs are shown in fig-
ure 9.
3.2 Wall effects
To estimate the effect to which the lat-
eral movement observed is due to the
wall, we ran simulations for each con-
figuration with and without one of the
discs present. For a simulation of con-
figuration 1 in which the discs are far
enough apart that they do not interact,
in contrast to figure 9(a), figure 10(a)
shows that the motion of the left-hand
disc changes little when the right-hand
disc is removed. We therefore surmise
that the wall has no influence on the mo-
tion of the discs here.
In configuration 2, figure 10(b) shows
that the lower disc perturbs the foam in
such a way that the upper disc moves
more quickly than if it were not present.
Neither disc moves sideways to a great
extent.
3.3 Varying disc size
We now fix the initial centre-to-centre
separation of the discs and measure how
the separationevolvesduring thedescent
for discs of different size. For configura-
tion 1, figure 11 shows that, in general,
the area of the discs makes little differ-
ence. In all but one case one of the discs
falls behind the other one. That our sys-
tem always chooses the left-hand disc is
probablyan artefact due to the foam cre-
ation step.
Fixing the initialcentre-to-centresep-
aration of the discs in configuration 2
leads to the result, shown in figure 12,
that the distance between the discs is re-
duced more quickly when the discs are
smaller. This can be attributed to the
fact thatsmaller discsexperiencea smaller
drag force [36], and are therefore more
affected by small changes in the foam
structure in the wake of another disc.
10 A. Wyn et al.: 2D Foam Rheology
Fig. 8. Two discs in a monodisperse 2D foam confined in a channel of width W. (a)
Configuration 1, in which the discs start side-by-side, with a distance d1 between their
centres and an angle ? made by the line joining their centres. (b) Configuration 2, in
which the discs start one above the other a distance d2 apart.
 0
 0.2
 0.4
 0.6
 0.8
 1
 0  0.2  0.4  0.6  0.8
 0
 0.2
 0.4
 0.6
 0.8
 1
 0  0.2  0.4  0.6  0.8
x positionx position
yp
os
iti
on
yp
os
iti
on
Fig. 9. The motion of the disc centres in two typical simulations. Disc area is 4.5A. (a)
Configuration 1, with d1 = 0.08. Here, both discs move a small distance to the right, and
the disc that began on the left of the foam advances more slowly and moves behind the
right-hand disc. (b) Configuration 2, with d2 = 0.2. The discs barely deviate to the sides,
but the higher disc moves slightly faster in the wake of the lower one.
A. Wyn et al.: 2D Foam Rheology 11
(a)
 0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0  50  100  150  200
Second Disc Removed
Second Disc Included
rep
lac
em
en
x
po
sit
ion
Iterations
 0
 0.2
 0.4
 0.6
 0.8
 1
 0  50  100  150  200
Second Disc Removed
Second Disc Included
yp
os
iti
on
Iterations
(b)
 0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0  50  100  150  200
Second Disc Removed
Second Disc Included
x
po
sit
ion
Iterations
 0
 0.2
 0.4
 0.6
 0.8
 1
 0  50  100  150  200
Second Disc Removed
Second Disc Included
yp
os
iti
on
Iterations
Fig. 10. The position of the disc centres with disc area 3.5A. (a) In configuration 1, the
presence of the right-hand disc (with d1 = 0.24) does not affect the descent of the left-hand
one. (b) In configuration 2, the upper disc descends more quickly when the lower disc (a
distance d2 = 0.2 away) is present, but they do not deviate sideways.
(a)
 0
 1
 2
 3
 4
 5
 0  50  100  150  200  250  300
Se
pa
ra
tio
n
d1
/?
A
Iterations
2A
2.5A
3A
3.5A
4.5A
7A
(b)
pi/2
pi/4
0
-pi/4
 0  50  100  150  200  250  300
An
gle
,?
Iterations
2A
2.5A
3A
3.5A
4.5A
7A
Fig. 11. (a) The separation, measured in bubblewidths, of the disc centres in configuration
1 for a range of disc areas, expressed in multiples of the bubble area. The initial separation
is close to 3 bubble widths. On average there is a slight increase in the separation between
the discs, and no clear trend with disc size. (b) The angle ? made by the line joining the
disc centres. In all but one case the left-hand disc moves behind the right-hand one.
3.4 Varying disc separation
The results for configuration 1 in figure
11 may be difficult to interpret because
the distance between the edges of the
discs varies as well as their areas. This
means that there are a different number
ofbubbles betweenthe discsin eachcase.
We now fix the disc size and vary the ini-
tial separation between the edges of the
discs. Figure 13 emphasizes that discs
which startclosertogether aremorelikely
to interact, and for one to move behind
the other.
12 A. Wyn et al.: 2D Foam Rheology
 3
 3.5
 4
 4.5
 5
 5.5
 6
 6.5
 7
 0  50  100  150  200  250  300
Se
pa
ra
tio
n
d2
/?
A
Iterations
2A
3.5A
4A
5A
5.5A
Fig. 12. The separation, measured in bubble widths, of the two discs in configuration 2.
The upper disc ?catches up? the lower disc, and it does so more quickly if the discs are
smaller.
 0
 0.5
 1
 1.5
 2
 0  50  100  150  200  250  300
0.1W
0.2W
0.3W
0.4W
No
rm
ali
ze
d
Se
pa
ra
tio
n
Iterations
pi/2
pi/4
0
-pi/4
 0  50  100  150  200  250  300  350
0.1W
0.2W
0.3W
0.4W
An
gle
,?
Iterations
Fig. 13. The separation, normalized by the initial separation, and the angle made by the
line joining the centres of the two discs in Configuration 1. The initial separations in the
key are given in terms of the width W of the foam. For the closest discs, the left-hand disc
moves behind the right-hand one, while for the next closest pair the reverse occurs.
4 Conclusion
We have shown that localization in the
linear Couetteshearof a 2Dfoamishighly
dependent on the polydispersity in bub-
ble areas. At high dispersity, the local-
ized region can extend throughout the
foam, with the width of this region de-
pendent upon the square-root of the sec-
ond moment of area disorder, ?2(A). A
signature of the localization is given by
a 1D measure of the polydispersity, Ly,
enabling us to probe how bubbles read-
just to the shear and allow T1s to collect
in specific regions of the foam.
The results begin to explain why the
simulations of Kabla and Debregeas [52]
indicatethat T1soccur closeto the walls:
the area dispersity is low in their case, so
that we expect the localised region to be
near a wall.
Our simulationsof two discs descend-
ing through a foam probe the interac-
tions between falling objects. We have
shown that small particles move faster
in the wake of another, and that for suf-
ficiently close pairs of particles one is at-
tracted into the wake of the other. Given
the elastic nature of the interaction, and
the discrete nature of the foam which al-
lows single bubbles to detach from the
disc and possibly move upwards, we are
currentlyworkingtowardsextracting the
bubble displacement field to test for the
phenomenon of a ?negative wake? [53].
It remainsto investigatethe effects of
a polydisperse foam, although we expect
these to be small, and the interaction be-
tween more particles/discs. We are cur-
rently pursuing an experimental realiza-
tion of this system, to ascertain the ap-
plicability of the results given here.
For both systems studied here, and
indeed for foam research in general, the
influence of liquid is important. Extend-
ing this work to wet foams presents a
A. Wyn et al.: 2D Foam Rheology 13
challenge to simulation, although this is
itself overshadowed by the demands of
three-dimensional calculations on wet
foams.
Acknowledgements
We thank K. Brakke for developing, dis-
tributing andsupporting theSurfaceEvolver
and for providing the script to gener-
ate Voronoi input. This work benefited
from discussions with D. Weaire, K. Kr-
ishan and F. Graner. Financial support
is gratefully acknowledged from EPSRC
(EP/D048397/1,EP/D071127/1)andthe
British Council Alliance programme.
References
1. J.J. Bikerman. Foams: Theory and
Industrial Applications. Reinhold,
New York, 1953.
2. R.K. Prud?homme and S.A. Khan,
editors. Foams: Theory, Measure-
ments and Applications, volume 57
of Surfactant Science Series. Marcel
Dekker, New York, 1996.
3. D. Weaire and S. Hutzler. The
Physics of Foams. Clarendon Press,
Oxford, 1999.
4. J.A.F. Plateau. Statique
Exp?erimentale et Th?eorique des
Liquides Soumis aux Seules Forces
Mol?eculaires. Gauthier-Villars,
Paris, 1873.
5. C.S. Smith. The Shape of Things.
Scientific American, 190:58?64,
1954.
6. L. Bragg and J.F. Nye. A dynamical
model of a crystal structure. Proc.
R. Soc. Lond., A190:474?481, 1947.
7. J. Lauridsen, M. Twardos, and
M. Dennin. Shear-induced stress re-
laxation in a two-dimensional wet
foam. Phys. Rev. Lett., 89:098303,
2002.
8. M. Twardos and M. Dennin. Com-
parison between step strains and
slow steady shear in a bubble raft.
Phys. Rev. E, 71:061401, 2005.
9. Y. Wang, K. Krishan, and M. Den-
nin. Statistics of microscopic yield-
ing in sheared aqueous foams. Phil.
Mag. Letts., 87:125?133, 2007.
10. C.S. Smith. Grain shapes and other
metallurgical applications of topol-
ogy. In Metal Interfaces, pages 65?
108. American Society for Metals,
Cleveland, OH, 1952.
11. M.E. Rosa and M.A. Fortes. De-
velopment of bamboo structure in a
2D liquid foam. Europhys. Lett., 41:
577?582, 1998.
12. B. Dollet, M. Aubouy, and
F. Graner. Inverse Lift: a sig-
nature of the elasticity of complex
fluids. Phys. Rev. Lett., 95:168303,
2005.
13. M.F. Vaz and S.J. Cox. Two-
bubble instabilities in quasi-two-
dimensional foams. Phil. Mag.
Letts., 85:415?425, 2005.
14. Y. Wang, K. Krishan, and M. Den-
nin. Impact of boundaries on veloc-
ity profiles in bubble rafts. Phys.
Rev. E, 73:031401, 2006.
15. J.E. Taylor. The structure of singu-
larities in soap-bubble-likeand soap-
film-like minimal surfaces. Ann.
Math., 103:489?539, 1976.
16. D. Weaire and J.P. Kermode.
Computer simulation of a two-
dimensional soap froth I. Method
and motivation. Phil. Mag. B, 48:
245?249, 1983.
17. F. Bolton and D. Weaire. The ef-
fects of Plateau borders in the two-
dimensional soap froth. II. General
simulation and analysis of rigidity
loss transition. Phil. Mag. B, 65:
473?487, 1992.
18. T. Okuzono, K. Kawasaki, and
T. Nagai. Rheology of Random
Foams. J. Rheol., 37:571?586, 1993.
19. S.J. Cox, D. Weaire, and J.A.
Glazier. The rheology of two-
dimensional foams. Rheol. Acta, 43:
442?448, 2004.
20. I. Cantat and R. Delannay. Dynam-
ical transition induced by large bub-
bles in two-dimensional foam flows.
Phys. Rev. E, 67:031501, 2003.
21. S.J. Cox. The mixing of bubbles in
two-dimensional foams under exten-
14 A. Wyn et al.: 2D Foam Rheology
sional shear. J. Non-Newtonian Fl.
Mech., 137:39?45, 2006.
22. K. Brakke. The Surface Evolver.
Exp. Math., 1:141?165, 1992.
23. S.J. Cox and E.L. Whittick. Shear
modulus of two-dimensional foams:
The effect of area dispersity and dis-
order. Eur. Phys. J. E, 21:49?56,
2006.
24. A.M. Kraynik, D.A. Reinelt, and
F. van Swol. The structure of ran-
dom monodisperse foam. Phys. Rev.
E, 67:031403, 2003.
25. D. Weaire and N. Rivier. Soap,
cells and statistics?random patterns
in two dimensions. Contemp. Phys.,
25:59?99, 1984.
26. G. Debr?egeas, H. Tabuteau, and
J.M. di Meglio. Deformation and
flow of a two-dimensional foam un-
der continuous shear. Phys. Rev.
Lett., 87:178305, 2001.
27. S.J. Cox. Simulations of Two-
Dimensional Foam under Cou-
ette Shear (preprint), 2007.
http://hdl.handle.net/2160/323.
28. Y. Jiang, P.J. Swart, A. Saxena,
M. Asipauskas, and J.A. Glazier.
Hysteresis and Avalanches in two-
dimensional foam rheology simula-
tions. Phys. Rev. E, 59:5819?5832,
1999.
29. E. Janiaud, D. Weaire, and S. Hut-
zler. Two dimensionalfoam rheology
with viscous drag. Phys. Rev. Lett.,
97:038302, 2006.
30. P. Marmottant and F. Graner. An
elastic, plastic, viscous model for
slow shear of a liquid foam. Euro.
Phys. J. E, 23:337?347, 2007.
31. M.L. Falk and J.S. Langer. Dy-
namics of viscoplastic deformation
in amorphous solids. Phys. Rev. E,
57:7192?7205, 1998.
32. N.P.Kruyt. Onthe shearmodulus of
two-dimensional liquid foams: a the-
oretical study of the effect of geo-
metrical disorder. J. Appl. Mech.,
74:560?567, 2007.
33. B. Dollet, F. Elias, C. Quilliet,
C. Raufaste, M. Aubouy, and
F. Graner. Two-dimensional flow of
foam around an obstacle: Force mea-
surements. Phys. Rev. E, 71:031403,
2005.
34. B. Dollet and F. Graner. Two-
dimensional flow of foam around
a circular obstacle: local measure-
ments of elasticity, plasticity and
flow. J. Fl. Mech., 585:181?211,
2006.
35. C. Raufaste, B. Dollet, S. Cox,
Y. Jiang, and F. Graner. Yield
drag in a two-dimensional foam flow
around a circular obstacle: Effect of
liquid fraction. Euro. Phys. J. E, 23:
217?228, 2007.
36. S.J. Cox, B. Dollet, and F. Graner.
Foam flow around an obstacle: sim-
ulations of obstacle-wall interaction.
Rheol. Acta., 45:403?410, 2006.
37. B. Dollet, F. Elias, C. Quil-
liet, A. Huillier, M. Aubouy, and
F. Graner. Two-dimensional flows of
foam: dragexerted oncircular obsta-
cles and dissipation. Coll. Surf. A,
263:101?110, 2005.
38. J.R. de Bruyn. Transient and
steady-state drag in foam. Rheol.
Acta, 44:150?159, 2004.
39. I. Cantat and O. Pitois. Mechani-
cal probing of liquid foam ageing. J.
Phys.: Condens. Matter, 17:S3455?
S3461, 2005.
40. I. Cantat and O. Pitois. Stokes ex-
perimentin a liquid foam. Phys. Flu-
ids, 18:083302, 2006.
41. S.J. Cox, M.D. Alonso, S. Hutzler,
and D. Weaire. The stokes ex-
periment in a foam. In P. Zitha,
J. Banhart and G. Verbist, editor,
Foams, emulsions and their appli-
cations, pages 282?289. MIT-Verlag,
Bremen, 2000.
42. H. Tabuteau, F.K. Oppong, J.R. de
Bruyn, and P. Coussot. Drag on
a sphere moving through an aging
system. Europhys. Lett., 78:68007,
2007.
43. J. Feng, P.Y. Huang, and D.D.
Joseph. Dynamical simulation of
sidementation of solid particles in an
oldroyd-b fluid. J. non-Newt. Fl.
Mech., 63:63?88, 1996.
44. S. Daugan, L. Talini, B. Herzhaft,
and C. Allain. Aggregation of par-
ticles settling in shear-thinning flu-
ids. Part 1. Two-particle aggrega-
A. Wyn et al.: 2D Foam Rheology 15
tion. Euro. Phys. J. E, 7:73?81,
2002.
45. B. Gueslin, L. Talini, B. Herzhaft,
Y. Peysson, and C. Allain. Aggrega-
tion behavior of two spheres falling
through an aging fluid. Phys. Rev.
E, 74:042501, 2006.
46. B. Dollet, M. Durth, and F. Graner.
Flowof foampast anelliptical obsta-
cle. Phys. Rev. E, 73:061404, 2006.
47. J. Wang and D.D. Joseph. Potential
flow of a second-order fluid over a
sphere or an ellipse. J. Fl. Mech.,
511:201?215, 2004.
48. K. Brakke. 200,000,000
Random Voronoi Polygons.
www.susqu.edu/brakke/papers/voronoi.htm,
1986. Unpublished.
49. N. Kern, D. Weaire, A. Martin,
S. Hutzler, and S.J. Cox. Two-
dimensional viscous froth model for
foam dynamics. Phys. Rev. E, 70:
041411, 2004.
50. L. Arnaud, J. Weiss, M. Gay, and
P. Duval. Shallow-ice microstruc-
tureatDome Concordia,Antarctica.
Ann. Glaciol., 30:8?12, 2000.
51. M. Aubouy, Y. Jiang, J.A. Glazier,
and F. Graner. A texture tensor
to quantify deformations. Granular
Matter, 5:67?70, 2003.
52. A. Kabla and G. Debregeas. Quasi-
static rheology of foams. part 1. os-
cillating strain. J. Fluid Mech., 587:
23?44, 2007.
53. O. Hassager. Negative wake behind
bubbles in non-Newtonian liquids.
Nature, 279:402?403, 1979.