## Abstract

The steady-state solutions of the viscous froth model for foam dynamics are analysed and shown to be of finite extent or to asymptote to straight lines. In the high-velocity limit, the solutions consist of straight lines with isolated points of infinite curvature. This analysis is helpful in the interpretation of observations of anomalous features of mobile two-dimensional foams in channels. Further physical effects need to be adduced in order to fully account for these.

## 1. Introduction

Foams are widely used in traditional applications such as fire-fighting and froth flotation (Kraynik 1988; Prud’homme & Khan 1996). More recently, they have begun to find new applications in the emerging field of discrete microfluidics (Garstecki *et al.* 2004; Drenckhan *et al.* 2005), in which small volumes of gas (or, equivalently, liquid) can be manipulated within narrow channels.

When monodisperse bubbles are confined in narrow channels with a low liquid fraction, they generally form ordered structures (Garstecki *et al.* 2004; Drenckhan *et al.* 2005; Garstecki & Whitesides 2006) (figure 1). They may be readily manipulated in networks of such channels. Such a system is suggestive of a practical system of microfluidics (Weaire & Drenckhan 2007), so the motion of the bubbles, when driven by a pressure, is of practical as well as of basic interest. The desire to predict the structure and dynamics of such a foam arises from the requirement to design both chemical formulation and container geometry with maximum efficiency, to perform specific functions.

In a rectangular channel whose width is much greater than its depth, the structure in question may be regarded as essentially two-dimensional. That is, each bubble touches both of the bounding plates and the soap films that span the gap are perpendicular to those plates. Accordingly, a two-dimensional model, previously developed in the physics of foams (Kern *et al.* 2004) has been used in this area. This ‘viscous froth’ model is only a provisional or skeletal one (including in particular linear relationships when power laws may be more realistic), but it has succeeded in shedding light on some observed phenomena that are seen when the flow velocity is large enough to depart from the quasi-static condition (Kern *et al.* 2004; Cox 2005; Green *et al.* 2006; Grassia *et al.* 2008). However, it then appeared that significant qualitative details of observed structures were apparently not compatible with the model and it was not clear what needed to be added to it. The present study, which explores the high-velocity limit, was stimulated in part by these anomalous structures.

We review the model in §2 and find its steady-state solutions for single lines in §3. Even this apparently straightforward problem presents quite a wealth of mathematical detail. Simplicity is restored in the high-velocity limit (Grassia *et al.* 2008), which is described in §4. Any of these steady-state solutions may be used to provide pieces of a more general solution in which the lines meet at 120^{°}, allowing comparison with the experiment, as discussed in §5.

## 2. The viscous froth model

A static dry foam can be represented by the ideal two-dimensional soap froth (Weaire & Hutzler 1999), in which each film is a circular arc, and three arcs meet at a vertex at equal angles of 120^{°}. The elastic force on each film, which is the product of curvature κ and surface (or line) tension γ, balances the pressure difference Δ*p* across it according to the Laplace law, Δ*p*=γκ. The viscous froth model (VFM) extends the ideal soap froth to predict the dynamics of such structures at a finite deformation rate (see figure 2). In the VFM, *three* forces act on each element and sum to zero, since inertia is neglected. They are due to surface tension, pressure difference and a drag force that resists motion normal to the film.

Consider a point on a soap film: a force balance in the direction of the normal leads to the following evolution equation:
2.1
where λ is the drag coefficient and *v*_{n} the normal velocity. The exponent β is usually set to unity (the linear VFM, as considered here) to speed up the numerical algorithm and to simplify the analytical theory (Kern *et al.* 2004).

Note that the 120^{°} rule for internal vertex angles is an essential feature of the model, since three equal surface tension forces, and no others, act on a point.

## 3. Steady-state motion

This analysis is concerned with the steady-state motion with a constant velocity. It generalizes earlier treatments of steady-state solutions of curvature-driven growth (the VFM with zero pressure difference), reviewed in Weaire & McMurry (1996), and soap film motion in channels (Green *et al.* 2006). Here we state the key results, relegating derivations and details to appendix A.

We look for solutions of the VFM for motion at a constant velocity of magnitude *V* . We use Cartesian axes and seek curves that translate in the positive *x*-direction by writing equation (2.1) in the following form (Cox & Mishuris in press):
3.1
where *x*_{y} is the derivative of *x*, etc. We normalize the *x* and *y* coordinates by the length scale
3.2
and write
3.3
(cf. the mobility parameter *a* of (Grassia *et al.* 2008, equation (18))) to give
3.4

There are three classes of solutions, the details of which are given in appendix A. Note that, given the form of equation (3.1), we may choose each curve to originate at (0,0). They are classified according to their shape at the origin, illustrated schematically in figure 3; for a given α, there may be more than one solution, distinguished by the tangent at the origin.

The most straightforward steady states are simply straight lines (in which the middle term of equation (3.4) is zero) with slope
3.5
as in figure 3*a*. These solutions are found only for |α|≤1; for α=0 (i.e. no pressure difference), they are parallel to the *x*-axis, whereas for α=±1 they are parallel to the *y*-axis.

In addition, there are curves that evolve from the origin with , as *x*→0. For 0<α≤1, they asymptote to the above straight lines (equation (3.5)). Otherwise they reach a region of zero slope at which we terminate our calculation to leave curves of finite extent with accessible endpoints (see figure 4*a*). The details of the envelopes, which are the locus of the endpoints, are presented in appendix A.

Finally, there are solutions with *y*∼±*x*^{2} at the origin, shown in figure 4*b*. For positive α, these lie above the *x*-axis, and changing the sign of α results in a reflection in this axis. For |α|≤1, they again asymptote to the above straight lines (equation (3.5)). As for the previous case, we follow the curves up to the point at which the respective derivatives tend to infinity, and the corresponding points lie on essentially the same envelope.

The coincidence of the envelopes can be seen by appealing to figure 3: each curve of finite extent in figure 3*b*, with square-root behaviour at the origin, approaches the envelope as a parabola. Thus, it is also represented in figure 3*c*, i.e. it is one of the curves with quadratic behaviour at the origin, albeit with a shift and possible change of orientation. For example, the dashed curve for α=−0.5 in the quadrant *x*>0,*y*>0 in figure 3*b* can be translated to its (dashed) equivalent in the quadrant *x*<0,*y*<0 in figure 3*c*.

Each curve of finite extent can also be extended with any other solution for the same value of α. For example, the curve for α=1.5 in the quadrant *x*>0,*y*>0 in figure 3*c* can be continued by the curve for α=1.5 in the quadrant *x*<0,*y*>0 in figure 3*b*. This is elaborated in figure 1.

As the velocity increases, the extent of the finite, curved, solutions tends to zero, whereas the (normalized) straight-line solutions remain invariant. The other solutions of infinite extent tend towards piece-wise straight-line forms, as discussed below.

Before proceeding, we first illustrate the use of the steady-state solutions with the example of a single bubble confined within a channel and attached to a moving side wall, shown in figure 1. The ends of the interface enclosing the bubble are fixed to the wall at *y*=0 a distance 2*L* apart, and the wall is moved in the positive *x*-direction with speed *V* . We use a Surface Evolver-based implementation (Brakke 1992; Cox 2005) to solve the VFM to find the bubble shape, with parameters γ=λ=1, , time step Δ*t*=1×10^{−5} and circa 60 line segments. The pressure *p* in the bubble is varied to enforce a constraint of fixed area, *A*_{b}=π*L*^{2}/2, although for the values of *V* considered here *p* never differs greatly from γ/*L*. Figure 1*a* shows bubble shapes for different *V* , which cease to be physical when *V* is sufficiently large that the angle between the trailing edge and the wall goes to zero. We find this critical velocity to be *V* ≈2.3 (figure 1*b*): the angles at the leading and trailing edges are equal only at low velocity, and as the velocity increases, the angle at the trailing edge decreases more quickly.

We show in figure 1*c* how the solutions of figure 4 can be used to predict the bubble shape. In the case , the parameter α is equal to one. Thus, the shape of the leading edge is given by part of the curve for α=1 in figure 4*b*. Since the pressure difference acts in the opposite direction on the trailing edge, the latter’s shape is given by part of the curve for α=−1 in figure 4*a*. The two are joined at the point where the slope is vertical and truncated to give two endpoints the required distance apart. The simulation result is clearly consistent with the derived steady-state solutions (figure 5*c*).

## 4. High-velocity limit

Here we examine more carefully the forms that emerge as *V* tends to infinity, by reference to figures 3 and 4. It is easiest to consider α to be constant which, in an experiment, would imply that Δ*p* is increased in proportion to *V* , according to equation (3.3) (i.e. fixing α while varying the pressure difference (Green *et al.* 2006)) while λ is kept constant. This analysis builds upon the results of Grassia *et al.* (2008), who showed that the curvature is pushed towards the side walls and gave approximate expressions for the film shape there.

The only effect of increasing *V* to a high value is then to uniformly shrink the curves progressively wherever they have finite curvature. The normalized infinite (or semi-infinite) solutions have significant curvature only close to the origin. The curvature here is of order 1/*r*≈α/*L*_{0} (=Δ*p*/γ), and the curvature is concentrated within a distance of the origin that is of order *r*. Now consider the evolution of the curve as the velocity goes to infinity. The curved parts become ever sharper, as *r* goes to zero.

In the limit , the solutions therefore consist entirely of straight lines. Consider, for example, the solution that we have referred to as ‘square-root, α=0.5’ in figure 3*b*. In this case, the limiting high-velocity solution consists of a pair of straight lines inclined at ±30^{°} to the *x*-axis and joined at the origin, where the curvature is infinite. The case that is referred to as ‘quadratic, α=0.5’ in figure 3*c* consists of just a single line, also inclined at 30^{°}. All other solutions or parts of solutions vanish in the limit , so we are left with essentially a single possibility consisting of straight lines, arising from the three cases of figure 3.

### (a) Applying the high-velocity solutions

In appropriate cases, we would hope to use these simple solutions directly. Note that we may cut any piece of the infinite curves of figure 3 to fit certain boundary conditions. We may also use such pieces to construct high-velocity solutions with vertices, for structures such as the one in figure 1.

The main subtlety in doing so relates to the infinite curvature at the kink of the solution in figure 3*b,c* in the high-velocity limit. If a solution is excised that has this as an endpoint, the tangent at the rear end of this curve is not well defined. By referring to figures 3 and 4, we see that the direction of the tangent may be taken to lie anywhere in a broad range of angles.

To demonstrate these results, we first perform representative VFM calculations on a single film.

### (b) Single film with moving endpoints

Consider a single film (or line) that experiences a pressure difference Δ*p* and whose endpoints, located a distance 2*L* apart in the *y*-direction, move at a constant velocity *V* in the *x*-direction. We solve the VFM numerically, as described above, to determine the steady-state shape of the film in this case, for fixed α and γ=λ=*L*=1. Figure 6*a* shows the convergence towards the predicted shape as the velocity (and pressure difference) increases for α=0.5. In the special case α=0, the films consist of segments of Mullins’ solution.

### (c) Single film between parallel walls

Now consider a single film that experiences a pressure difference Δ*p* that pushes it along a straight channel of width 2. In addition to the drag on the bounding plates, with drag coefficient λ, we introduce a further drag on the side walls of the channel proportional to velocity, with coefficient λ_{sw}. There is then a new evolution equation for the points of the film at the side walls
4.1
where θ is the angle between the film and the wall.

The numerical implementation of the VFM described above is used to find the shape and velocity of the film for a given pressure difference Δ*p* and drag coefficient λ_{sw}=0.01. For these parameter values, we find that the velocity increases in direct proportion to Δ*p*, giving α≈1. Figure 6*b* shows the convergence with increasing velocity towards the film shape expected from figure 4 for this value of α. That is, the film becomes vertical away from the wall and the angle θ decreases to zero.

### (d) Staircase structure

We next consider the extended array of bubbles depicted in figure 1; in order to mimic experiments (Drenckhan *et al.* 2005), we choose to find the steady-state shape of this 211 or staircase structure at finite flow rate, with sidewall drag as in equation (4.1).

To avoid end-effects, we use a sample that is periodic in the *x*-direction, in a channel of width 2*L*=7/4 with a bubble area . A pressure gradient is applied by invoking an additional contribution of Δ*p** to the pressure difference across each film. (This results in a drop of 2Δ*p** between adjacent bubbles on the same side of the channel.) For sufficiently high Δ*p**, the angle θ goes to zero and the system becomes unstable. In what follows, we shall concentrate on the configuration close to this unstable limit, denoting the critical values of pressure and velocity by and *V*_{crit}, respectively.

Applying the arguments of the previous sections leads us to expect the limiting structure at a high velocity to take a surprising form: it returns to the original static or low-velocity form of figure 1. This is consistent with what we find (figure 7): note that curvature is here concentrated at the points of contact with the side walls, so that a contact angle consistent with the drag at the side walls is maintained.

There is an apparent approximate scaling for the critical velocity and pressure difference in the high-velocity limit
4.2
which, with the exception of the pre-factors, can be derived from equations (2.1) and (4.1) in the limit in which . This scaling for *V*_{crit} is in agreement with the results in Grassia *et al.* (2008) for a single film (cf. §4*c*).

## 5. Discussion

We are therefore able to distinguish two important cases: the case |α|≤1 consists of solutions that asymptote to a straight-line solution, while for |α|>1 they do not. In the former case, we find the ‘inverted wing’ and the ‘hockey stick’ shown in figure 3. In the latter, all solutions are of finite extent. For any given system, with a fixed drag coefficient λ, *V* increases with Δ*p* (not necessarily in proportion) so that α does not vary much (Grassia *et al.* 2008). If α is sufficiently small in magnitude, then as *V* or Δ*p* tends to infinity any curved region shrinks as its curvature too becomes infinite.

The analysis of the previous sections has proved useful in understanding the way in which the results of numerical simulations vary with *V* . But can they, in the end, shed light on the existing experimental results that provoked this study? We show these in figure 8; they may be summarized as follows:

the lines are almost straight,

internal vertex angles clearly deviate from 120

^{°}.

The first result is entirely consistent with the kind of high-velocity solutions that we have described. At first sight, it might seem that the existence of undefined tangents at endpoints is the source of the anomalous internal angles, but a careful application of the rules that follow from §4 shows this to be impossible. That is, these deviations from the 120^{°} angle rule *cannot* be explained by the present model. We may ask what *additional* effects might be responsible for them. Although we have no definite evidence as yet, a prime candidate is the effect of a drag force acting at the vertex, which is in reality a junction of finite extent. One may envisage a vertex correction in the form of a force acting at the vertex point. Such forces have been posited before, but without any experimental support.

If such a connection is considered, we find that the angles observed in figure 8 (notwithstanding the difficulty of measuring such an angle from an image of a foam) are such that the required vertex force must act in the direction of the velocity, rather than opposing it. This *negative* drag on a vertex seems unphysical at first. However, it may be taken to represent the *change* in the local drag as a *finite* vertex (attached to a finite transverse Plateau border) is replaced by the point vertex of the present model. There is no reason why the finite vertex should not suffer a lesser drag force. Further experiments are required to test this suggestion and interpretation.

Experiments on single films in micro-channels (Raufaste 2007; Raven 2007) show shapes similar to those found in figure 6*a*, and further analysis of these experiments may be useful in trying to observe the curvature of the films close to the wall. In addition, the solutions derived here may be combined in different ways to shed light on, for example, the motion of single bubbles in channels.

The present appreciation of the properties of the model is a step forward in modelling microfluidic networks, but we would not wish to underestimate the further factors that may come into play in due course, particularly when shearing states are considered. They include at least film stretching (Marangoni, etc.) effects, finite Plateau borders, flow in Plateau borders, surface viscosity and nonlinear drag laws. For example, the VFM as originally formulated included a power-law form for the velocity, in order to accommodate the Bretherton-type analysis of bubble motion in tubes (Bretherton 1961), and Grassia *et al.* (2008) found that for β<1 there is an increase in the width of the region in which the curvature is concentrated at steady state. Nevertheless, the present model provides an adequate platform to begin to investigate designs of networks and their components for experimental testing.

## Acknowledgements

We thank W. Drenckhan for stimulating discussion and the provision of experimental results. Financial support from (D.W.) the European Space Agency (MAP AO-99-108: C14914/02/NL/SH, MAP AO-99-075: C14308/00/NL/SH) (G.M.), the Wales Institute of Mathematical and Computational Sciences and (S.J.C.) EPSRC (EP/D048397/1, EP/D071127/1) is gratefully acknowledged.

## Appendix A. Derivation of steady-state solutions

It is apparent that the VFM is closely related to motion by local curvature (Glazier & Weaire 1992). Indeed, it was designed to be so in the limit in which pressure differences between cells go to zero. Motion by mean curvature (Brakke 1977; Weaire & McMurry 1996), also known as the curve-shortening problem and described by the eikonal-curvature equation, is often invoked to describe the motion of grain boundaries in crystals. Mullins (1956) gave steady-state solutions in this limit (i.e. equation (2.1) with Δ*p* set to zero) for translation, rotation and expansion. Here, we generalize the solutions for translation to the case of the VFM.

We classify the solutions of equation (3.4) with respect to their behaviour at origin (since we can set *x*(0)=0 without loss of generality). We conclude that
where the following three cases can be distinguished:
A 1
There are in fact no other cases. A further condition on the derivative at the origin is required to distinguish the different branches of the solution.

Solutions in the case ϖ=1 are straight lines given by equation (3.5), found only for |α|≤1.

#### (a) Case ϖ=2

To seek further, curved, steady-state solutions, we rewrite equation (3.1) with the substitution:
A 2
In the case ϖ=2, we adopt the following boundary conditions:
A 3
Depending on the value of the parameter α, three different families of solutions with finite gradient are possible (see figure 4*a*): finite length solutions for α<0 and α>1 and infinite length solutions for 0≤α≤1 that tend to the respective straight-line solutions discussed above. Note that the limiting cases are the well-known Mullins solution (α=0) or a line parallel to the *y*-axis when α=1.

This behaviour, as well as the shape of the envelope, can be confirmed by an analysis of the corresponding differential equation (A 2). Under conditions (A 3), this takes the form
A 4
with argument restricted to the interval . *G*(*t*) itself takes the following form (Green *et al.* 2006), depending on the value of α
A 5
where . For any α<0, *G*(*t*) is a positive increasing and bounded function in the interval [0,π/2]. Then the solution to equation (3.4) under conditions (A 3) takes the form
A 6
where *G*^{−1} is the inverse to the function *G*. In the case α>1, the function *G*(*t*) decreases in the interval [0,π/2], and the solution takes the same form (A 6) but now *x*(*y*)<0 for *y*>0. In both of these cases, the curves are of finite extent and we can compute their envelope, consisting of the points (*y*_{*},*x*(*y*_{*})), from equation (A 6). On the other hand, it follows from equation (A 4) that as *y*→*y*_{*} for any fixed α<0 or α>1. As the point (*y*_{*},*x*(*y*_{*})) can always be shifted to the origin of a new coordinate system, this means that the end of any curve for a given value of α (satisfying conditions (A 3)) can be continued by another curve from the remaining class of solutions, which satisfy slightly different conditions (A 7), described below. Note that the solutions obtained can be extended to negative *y*, as *x*∈*C*^{2}[0,*y*_{*}). In fact, owing to the symmetry in equation (3.4), *x*(−*y*)=*x*(*y*).

#### (b) Case ϖ=1/2

The third case of equation (A 1), ϖ=1/2, cannot satisfy conditions (A 3) since *x*_{y}(*y*) is in the interval (cf. equations (A 2) and (A 4)). This solution has infinite derivative at the origin, and we therefore apply the conditions
A 7
The corresponding equation for the derivative *y*_{x} can be evaluated in a way similar to equation (3.4), to give
A 8
Note that the case α=0 gives the trivial solution *y*(*x*)≡0 and is excluded from the following analysis. The solutions to equation (A 8) are invariant with respect to changing the signs of both *y* and α together; thus, solutions for negative α may be found by reflecting those for the corresponding positive value of α in the *x*-axis.

We now have , and hence
A 9
This can be re-written in the form
A 10
with the following solution for |α|<1 and *x* such that |e^{−x}−1|<|α|^{−1}:
A 11
Then, from equation (A 10), the following asymptotic formula holds for any |α|≤1:
A 12
That is, as , all curves again asymptote to straight lines, with the same slope as before (equation (3.5)).

#### (c) Envelope

Finally, we note that the non-zero points where the gradient *y*_{x} is zero (ϖ=2) or infinite (ϖ=1/2) coincide, since they both correspond to a limit of integration in equation (A 2), but from different sides (i.e. one can either integrate from the point at which *y*_{x}=0 to the point at which , or alternatively from to *y*_{x}=0). Thus, we can give an expression for the envelope. In parametric form, for the case ϖ=1/2 (figure 4*b*), it is
A 13

For negative α, one needs to change the sign of *y*_{*}. It is possible to eliminate the parameter α to have explicit representations in Cartesian coordinates, but we deliberately stay with the parametric form since it shows clearly the value of the physical parameter α for every point on the curve.

## Footnotes

- Received February 3, 2009.
- Accepted April 15, 2009.

- © 2009 The Royal Society