1. Introduction
Foam drainage is the process of liquid channelling through a foam or froth. Gravity-driven drainage can extract liquid from a foam, until extremely low liquid fractions are reached. Such foams are unstable and eventually break. However, when a continuous flow of liquid, 
$Q$, is added to the top of a foam-filled column with cross-section 
$A$, an equilibrium between added and extracted liquid is maintained. The balance is controlled by the permeability 
$\alpha$ of the foam,
\begin{equation} \alpha = \frac{Q}{A} \frac{\eta_f}{\rho_f g}, \end{equation}
which in turn depends on the liquid fraction. Parameters 
$\eta _f$, 
$\rho _f$ and 
$g$ in (1.1) denote the fluid dynamic viscosity, density and gravitational acceleration, respectively. If a small flow of liquid is added homogeneously to the top cross-section, a static and homogeneous liquid fraction is established in the foam. However, when the flow of liquid exceeds a critical value, an inhomogeneous liquid fraction results and convection rolls set in. These convection rolls bear similarities to Rayleigh–Bénard cells occurring in thermally stratified liquids when the density gradient, respectively the Rayleigh number, exceeds the critical value (Chandrasekhar Reference Chandrasekhar2013). But in foam drainage, no vertical density gradient is present when the convection rolls occur. Consequently, the mechanism behind these rolls is still a matter of debate.
When liquid is added to the top region of an aqueous foam column, it drains downward through the network of Plateau borders. This effect is described by the drainage equation (Cantat et al. Reference Cantat, Cohen-Addad, Elias, Graner, Höhler and Pitois2013):
\begin{equation} \frac{\partial \phi_l}{\partial t} = - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} = -\frac{\rho_f}{\eta_f} \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha \boldsymbol{g}) + \frac{\gamma}{2 \delta_b R_b \eta_f} \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \frac{\alpha}{\phi_l^{3/2}} \boldsymbol{\nabla} \phi_l \right).\end{equation}
It links the local liquid fraction 
$\phi _l$ and the superficial liquid flow rate 
$\boldsymbol {q}$ with the properties of the liquid (density 
$\rho _f$, viscosity 
$\eta _f$, surface tension 
$\gamma$) and those of the foam (permeability 
$\alpha$, bubble radius 
$R_b$). The geometric constant 
$\delta _b=1.74$ results from the relation between the Plateau border cross-section and the liquid fraction. The drainage equation (1.2) considers two effects: the vertical flow due to gravity 
$g$ and a diffusion of the liquid fraction due to capillary pressure gradients. The diffusion term balances inhomogeneities in the liquid fraction.
Experiments in steady drainage (Hutzler, Weaire & Crawford Reference Hutzler, Weaire and Crawford1998; Vera, Saint-Jalmes & Durian Reference Vera, Saint-Jalmes and Durian2000; Hutzler et al. Reference Hutzler, Cox, Janiaud and Weaire2007) have investigated whether the foam structure remains static. They have observed that if the liquid fraction exceeds a certain critical value 
$\phi _{crit}$, a stable convection roll is formed, known as convective instability (CI). This is linked to an inhomogeneous liquid distribution over the channel cross-section, which cannot be explained by the drainage equation. The formation of this convection roll biases experiments in foam at a high liquid fraction, because it disturbs the assumption of homogeneous liquid fraction.
Although detailed studies have not been carried out in industrial foam and froth applications, CI can be anticipated to have a severe impact. To illustrate this, we take deep-froth flotation as an example, which is an important step in the beneficiation of valuable mineral particles, e.g. copper-bearing ores. Frequently, wash water is added to the particle-laden foam (froth) to remove hydrophilic contaminants, such as undesired entrained gangue particles. It can be presumed that the CI reduces the efficiency of wash water addition. As a result, the gangue particles are not completely washed out. Thus, the grade of the particle concentrate in the froth is reduced which has a significant economic impact.
The critical liquid fraction 
$\phi _{crit}$ for the onset of CI has been determined to depend on the bubble size (Hutzler et al. Reference Hutzler, Weaire and Crawford1998) and can be reduced by tilting the channel (Cox et al. Reference Cox, Alonso, Weaire and Hutzler2006). Convective instability results in increased effective drainage and can cause bubble sorting (Hutzler, Weaire & Shah Reference Hutzler, Weaire and Shah2000). However, the actual onset mechanism is still a matter of debate. Embley & Grassia (Reference Embley and Grassia2006) assume that a sudden scarification of individual Plateau borders leads to a localization of the drainage flow. Weaire & Hutzler (Reference Weaire and Hutzler2003) demonstrate that dilatancy combined with elastic deformation of the foam can cause a growing imbalance in the liquid fraction. Neethling (Reference Neethling2006) predicts that anisotropic drainage in sheared foam causes an imbalance in the liquid fraction.
Using a combined experimental, analytical and numerical approach, we revisit the onset of CI in detail and provide a conclusive explanation of the mechanism underlying CI. To facilitate understanding, the stages of the instability which are passed in a typical experimental run are illustrated in figure 1. Starting from a dry state (I), the drainage flow rate is increased in small steps, leading to a stepwise increase in the liquid fraction, which causes the foam to expand by the volume of additional liquid content (II). At each drainage level, a static foam structure is reached. Above a critical liquid fraction of 
$\phi _{crit} \approx 0.65\,\%$ the deformation becomes inhomogeneous (III). However, the foam structure still reaches a static state for each drainage level. In this static state the shear rate is zero, but the shear angle with reference to the initial stage I is non-zero. The inhomogeneity and shear angle increase with increasing drainage flow rate. Above a second critical liquid fraction of 
$\phi _{CI} \approx 1.05\,\%$, the critical yield strain is exceeded near the bottom of the column and a steady convection roll is formed in the lower region. In the upper region, the foam maintains its static inhomogeneous deformation (IV). With further increase in drainage flow rate, the convection roll covers more and more of the channel height but remains steady for a constant drainage flow rate.
(a) Stages for the onset of CI and (b) corresponding liquid fraction and drainage flow rate for a typical experimental run with stepwise increase of drainage flow rate, covering (I) the initial foam, (II) a homogeneous expansion, (III) an inhomogeneous expansion above the first critical liquid fraction 
$\phi _{crit}$ and (IV) the formation of a steady convection roll, the CI, above the second critical liquid fraction 
$\phi _{CI}$.
We claim that the primary instability is the transition from homogeneous (II) to inhomogeneous (III) static deformations of the foam structure, which is linked to an inhomogeneous liquid distribution. In stage III, the inhomogeneity grows in the vertical direction. If the column is long enough and the growth rate of the inhomogeneity with respect to the vertical position is high enough, the yield stress is exceeded at some vertical position. Below that position the convection roll sets in, corresponding to stage IV. Thus, the CI appears as a secondary instability following the primary one, the inhomogeneous drainage.
In this paper we present experimental findings on the static deformation of the foam structure and the steady velocity field in the case of CI. Our findings are compared with a linear stability analysis for the growth of inhomogeneities of liquid fraction. This analysis includes the anisotropic drainage, as predicted by Neethling (Reference Neethling2006). The identified unstable mode marks the transition from stage II to III. Numerical simulations of the drainage combined with the elastic deformation of the foam and anisotropic drainage (Neethling Reference Neethling2006) reproduce our experimental findings. Finally, we demonstrate the dependency of the critical liquid fraction 
$\phi _{CI}$ for the onset of CI on the initial liquid inhomogeneity as well as on the channel length.
2. Materials and methods
Experiments are carried out in a vertical foam channel with an effective length of 990 mm and a cross-section 
$T \times B = 30\ \textrm {mm}\times 100\ \textrm {mm}$. Cox et al. (Reference Cox, Alonso, Weaire and Hutzler2006) have demonstrated that a small tilt of the column can lead to CI. Thus, the tilt angle was controlled to be less than 
$0.05^{\circ }$. Bubbles are generated at the bottom of the channel by a tube with 20 holes, each 0.5 mm in diameter, submerged in the surfactant solution. The tube is loaded with pressurized air, generating bubbles with radii of 
$R_b = 2 \pm 0.2\ \textrm {mm}$ at a flow rate of 
$2\ \textrm {l}\ \textrm {min}^{-1}$.
The surfactant solution is deionized water with 35 mM sodium dodecylsulphate (SDS), generating stable foam. When keeping the foam for 6 hours at 0.5 % liquid fraction, the measured change in bubble size was less than 20 %. Nevertheless, fresh foam is generated for each run, which take less than 1 hour each.
At the top of the foam column, a steady liquid flow is imposed through four porous hollow cylinders with outer diameters of 16 mm, yielding a forced drainage configuration similar to that of Leonard & Lemlich (Reference Leonard and Lemlich1965). The volumetric liquid flow rate 
$Q$ is linked to the drainage flow rate 
$q=Q/({B}\,{T})$ by the channel cross-section 
$({B}\,{T})$. The drainage flow rate is also known as the superficial drainage velocity or liquid flux, and equals the volume of liquid that passes through a certain cross-section of the foam channel within a certain time.
Each cylinder is connected to one channel of a four-channel peristaltic pump. This enables each cylinder to be controlled independently. That is, a defined inhomogeneous distribution can be imposed by pumping liquid only through one, two or three of the four cylinders. The notation (oxxx) means that only the leftmost cylinder is charged with liquid. The liquid is taken from the liquid reservoir below the foam. In that way, the total volume of liquid and gas remains constant and independent of the volumetric liquid flow rate.
The liquid fraction is observed with four pairs of electrodes (
$5\ \textrm {mm}\times 45\ \textrm {mm}$ electrode area) at a horizontal cross-section 
$20$–
$65$ mm below the porous cylinders (see figure 2a,b), yielding the horizontal profile of the liquid distribution 
$\phi (y)$. From the four measured values 
$\phi _1, \dots , \phi _4$, the average liquid fraction 
$\phi _m= 0.25 ( \phi _1 + \phi _2 + \phi _3 + \phi _4 )$ and the first moment of liquid fraction
\begin{equation} M_{\phi} = \frac{1}{B \phi_m}\int_0^B \left(y-\frac{B}{2} \right) \phi(y) {\textrm{d}y} = \frac{37.5\ \textrm{mm}\left(\phi_4-\phi_1 \right) + 12.5\ \textrm{mm} \left(\phi_3-\phi_2 \right)}{ \left( \phi_1 + \phi_2 + \phi_3 + \phi_4 \right) } \end{equation}
are derived. In our measurements, we never achieved a perfectly homogeneous distribution of liquid fraction. Even when all four porous cylinders are charged, we find 
$M_{\phi }$ of the order of 1 mm. This corresponds to an average horizontal shift of the liquid content by only 1 % of the channel width. We spent a lot of effort to further reduce the inhomogeneity. For example, we measured the flow rate for each channel from our four-channel peristaltic pump. Also, we switched the channels. But, even if we maintained everything constant and only created fresh foam, we found the inhomogeneity changing and sometimes switching to the other side. So we believe that small initial values of 
$M_{\phi }$ can be considered as a random distortion which are always present in a random foam structure.
(a) Set-up including foam channel (1) of variable length, bubble generator (2), four-headed peristaltic pump (3), four porous cylinders (4), four pairs of conductivity electrodes (5) and a camera observing the lower part of the channel (6). (b) Close-up of the electrodes (5). (c) Images of the foam and (d) close-up of the upper, middle and lower section, respectively, at 
$q=0\ \mathrm {\mu }$m s
$^{-1}$ (stage I), 
$q=101\ \mathrm {\mu }$m s
$^{-1}$ (stage III) and 
$q=179\ \mathrm {\mu }$m s
$^{-1}$ (stage IV). The displacement of distinct elements of the foam structure under increasing liquid flow rate is marked. The green line serves as guide to the eye. At 
$q=179 \mathrm {\mu }$m s
$^{-1}$, a steady movement sets in at the lower section and the element is lost.
The foam structure inside the transparent channel is observed with a camera taking images with 
$1920\ \textrm {px} \times 1200$ px at 30 frames per second, yielding a spatial resolution of 
$3.3\ \textrm {px}\ \textrm {mm}$. Backlight illumination with a light-emitting diode panel permits exposure times down to 20 ms. Due to the illumination, the captured foam structure is a cumulative image in the 
$z$ direction. Consecutive images are analysed with the particle image velocimetry (PIV) algorithm implemented in DaVis 8.0 software. Particle image velocimetry (Adrian, Adrian & Westerweel Reference Adrian, Adrian and Westerweel2011) compares consecutive images to measure the shift of the foam structure between them.
The PIV analysis is employed in this work in two different cases. In case 1 (stage IV) the network of Plateau borders and vertices moves, while in case 2 (stages I, II and III), the network is static.
In case 1 the shift of the network structure between consecutive images is multiplied with the frame rate, yielding the velocity distribution 
$\boldsymbol {v}$ of the foam in the 
$x$–
$y$ plane. This velocity refers to the average velocity of the foam structure, i.e. to the average velocity by which the network of Plateau borders and vertices moves. In this case convection rolls are recorded (stage IV). Typically, the network velocity reaches a steady, non-static state in our experiments. In the case of low velocities, increments of up to 20 frames were used for the PIV analysis. In addition, averaging over 10 image pairs without outlier reduction is employed. This algorithm yields the velocity distribution 
$\boldsymbol {v}(x,y)$ in the 
$x$–
$y$ plane with a 5 mm spatial resolution. The velocity of the liquid inside the Plateau borders is not accessible. Due to the backlight illumination, the velocity is automatically averaged in the 
$z$ direction. For spatially resolved three-dimensional measurements of the foam velocity, other techniques such as ultrasound Doppler velocimetry (Nauber et al. Reference Nauber, Büttner, Eckert, Fröhlich, Czarske and Heitkam2018) or radiographic particle tracking (Lappan et al. Reference Lappan, Franz, Schwab, Kühn, Eckert, Eckert and Heitkam2020) are required.
In case 2, the network is static everywhere in the column. In this case, the PIV analysis mentioned above yields vanishing velocity 
$\boldsymbol {v}(x,y)=0$. However, even when the network is static, it has deformed compared with the initial state (stage I) of zero drainage flow. We increase the liquid flow rate stepwise and find at each step a different static network. When we then switched off the drainage flow, the network goes back to its initial shape. It is like a beam balance. For each set of weights (drainage flow) a certain static angle (shear angle) occurs. In order to analyse this static elastic deformation 
$\boldsymbol {U}(x,y)$ of the foam (figure 3), one image is taken from each drainage flow rate level, and these are combined and fed into the PIV algorithm. By this means, the static elastic shift of the foam network between consecutive levels of drainage flow rate can be computed (see figure 2d). This only works if the foam structure does not undergo a flowing movement, because in the case of flowing foam, the shift between consecutive drainage levels is too large to be tracked.
Static vertical displacement 
$U_x$ in stages II (
$q=30$–
$47\ \mathrm {\mu }$m s
$^{-1}$) and III (
$q=65$–
$148\ \mathrm {\mu }$m s
$^{-1}$), which are below the critical limit 
$q_{CI}=160\ \mathrm {\mu }$m s
$^{-1}$ for the onset of CI. (a) Liquid fraction and (b) resulting static vertical displacement for five drainage flow rates. (c) Static vertical displacement in six horizontal layers (dotted lines) for eight stepwise increasing drainage flow rates (oooo).
3. Stages II and III: static elastic deformation
Each experimental run is prepared by creating fresh foam and letting it rest for 3 minutes. Then, the pump is started, feeding a constant volumetric flow rate 
$Q$ to the top of the column, resulting in a constant, downward drainage flow rate 
$q=Q/(B\,T)$, which is identical to the vertical superficial liquid velocity. In this section, only homogeneous inflow (oooo) is considered. The drainage flow rate 
$q$ is increased in small steps of 
$20$–
$50 \mathrm {\mu }$m s
$^{-1}$, giving the foam column 3 minutes after each step to achieve a steady state. This experimental procedure, which avoids undesired drainage fronts, is sketched in figure 1(b).
In stages II and III, the stepwise increase in the drainage flow rate causes a corresponding stepwise static elastic deformation 
$\boldsymbol {U}$ of the foam. Figure 3(b) shows the measured static vertical displacement 
$U_{x}$ of the foam at various drainage flow rates. The corresponding initial liquid fraction distributions are depicted in figure 3(a). Figure 3(c) shows 
$U_{x}$ at selected horizontal lines for different drainage flow rates.When the drainage flow rate is switched off again, the foam relaxes back into its initial shape at stage I. In stage II below 
$q=65\ \mathrm {\mu }$m s
$^{-1}$, the displacement 
$U_{x}$ is fairly constant in the 
$y$ direction but increases in the 
$x$ direction. The relative displacement 
$U_x/x$ is of the order of 0.5 %. This homogeneous displacement is due to the change in the liquid volume inside the foam while the total volume of gas and liquid is kept constant. For example, the column with a height of 1000 mm displaces the foam-liquid level by 5 mm when the liquid fraction is increased by 0.5 %.
At around 
$q=65\ \mathrm {\mu }$m s
$^{-1}$ and 
$\phi _l \approx 0.6\,\%$, the transition to stage III occurs. A horizontal gradient is seen in the vertical displacement, which is equivalent to a shear deformation 
$\varepsilon _{xy}$ of the foam. This shear 
$\varepsilon _{xy}$ is predicted by Neethling (Reference Neethling2006) to generate a horizontal deflection 
$q_{y,aniso}$ of the vertical drainage flow rate 
$q_x$, which, in turn, increases the liquid imbalance and shear. This mechanism is discussed in § 5. The liquid imbalance grows with increasing drainage flow rate until, above 
$160\ \mathrm {\mu }$m s
$^{-1}$, CI sets in (see figure 4, predicated on the same run as figure 3).
Convective instability (stage IV): occurrence of convection rolls in the case of homogeneous inflow of liquid. (a) Steady liquid fraction distribution and (b) corresponding vertical velocity distribution in the channel for five drainage flow rates 
$q$, showing the onset of movement and, thus, the transition from stage III to stage IV at 
$q=179\ \mathrm {\mu }$m s
$^{-1}$. Note the different velocity scales in the contour plots.
4. Stage IV: steady convection rolls
4.1. Homogeneous inflow
The onset of CI is recognized as when foam velocities larger than 
$0.5\ \mathrm {\mu }$m s
$^{-1}$ are detected in a steady state at any point of the foam column. Similarly to Hutzler et al. (Reference Hutzler, Weaire and Crawford1998), the critical liquid fraction and drainage flow rate at the onset of the CI are determined. Figure 4 documents the velocity distribution of the foam in a steady state at certain drainage flow rates. In this particular run, the critical drainage flow rate for the onset of CI is between 
$q=160$ and 
$179\ \mathrm {\mu }$m s
$^{-1}$, corresponding to 
$\phi _{CI} = 1.05\,\%$. At 
$q=160\ \mathrm {\mu }$m s
$^{-1}$, the foam structure remains static. At 
$q=179\ \mathrm {\mu }$m s
$^{-1}$, the first convection roll is detected at the bottom region, constituting the onset of CI. This roll is steady. It does not grow in size or magnitude as long as the drainage flow rate remains constant. As 
$q$ increases further, the foam velocity in the convection roll increases, as does the extension of the roll in a vertical direction. At drainage flow rates above approximately 
$q=500\ \mathrm {\mu }$m s
$^{-1}$, the roll covers the entire height of the channel. This relation between the drainage flow rate 
$q$ and the roll extension supports our mechanism in figure 1: the inhomogeneity of liquid fraction in stage III requires a certain vertical length to grow and reach a critical imbalance. With higher liquid fraction, the inhomogeneity grows faster with respect to the vertical distance. In the cross-section, where a critical level of inhomogeneity (i.e. the yield strain; see § 6) is exceeded, the steady movement of foam sets in. The resulting convection roll covers the channel height from this critical cross-section to the bottom. At higher liquid fraction the inhomogeneity grows faster and, thus, the convection rolls sets in more close to the top.
4.2. Inhomogeneous inflow
In order to further investigate the growth of the inhomogeneity in the liquid fraction with respect to the vertical position, an inhomogeneous liquid fraction is deliberately imposed by feeding only some of the four porous cylinders. The corresponding liquid moment 
$M_{\phi }$ according to (2.1) is derived from the conductivity sensor below the porous cylinders.
Figure 5(a) shows how the drainage flow rate 
$q_{CI}$ and liquid fraction 
$\phi _{CI}$ for the onset of CI depend on the liquid moment. In the case of high initial inhomogeneity, 
$q_{CI}$ and 
$\phi _{CI}$ can be as low as 
$70\ \mathrm {\mu }$m s
$^{-1}$ and 0.65 %, respectively (dash-dotted line in figure 5a). These values are close to the transition from stage II to stage III in homogeneous drainage. This again supports our explanation: if the initial inhomogeneity is very high already, there is only a small increase in inhomogeneity required to reach an imbalance that exceeds the yield stress. In the given channel length, this can be achieved by a small growth rate with respect to the vertical position. As demonstrated in § 5, a liquid fraction just above the critical liquid fraction 
$\phi _{crit}$ will cause such a slow growth of the inhomogeneity with respect to the vertical direction, because it is fed by the anisotropic drainage. However, in the case of a lower initial inhomogeneity (symbol ‘o’ in figure 5a), even drainage flow rates up to 
$250\ \mathrm {\mu }$m s
$^{-1}$ and 1.2 % liquid fraction can result in a static foam structure (stage III). In this case, the liquid inhomogeneity does not grow fast enough to reach a critical imbalance in the provided channel length.
Liquid fraction 
$\phi _{CI}$ and corresponding drainage flow rate 
$q_{CI}$ for the onset of CI obtained for experiments with homogeneous and inhomogeneous inflow of liquid. The horizontal dash-dotted lines are guides to the eye and mark the derived values for 
$\phi _{crit}$ from the transition from stage II to stage III. (a) Onset in a channel of 990 mm in length depending on the initial first moment of the liquid fraction 
$M_{\phi }$ caused by the specific combination of porous cylinders charged with liquid (o, on; x, off). Multiple data points for ‘oooo’ correspond to repeat experimental runs. (b) Onset depending on channel length for three different channel lengths and three different combinations of porous cylinders charged with liquid. Error bars denote the standard deviation in multiple experimental runs.
Now, the length of the channel is reduced. Figure 5(b) documents the critical drainage flow rate for different channel lengths and initial inhomogeneities. The error bars denote the variation between several runs, each with fresh foam. If the initial liquid fraction exhibits strong inhomogeneity (ooxx), the reduced length has only a small influence on the critical drainage flow rate. Similar to high initial inhomogeneities in figure 5(a) the drainage flow rate has to be just high enough to maintain the inhomogeneity. No large channel length is required to reach a critical imbalance. However, for an initially homogeneous liquid fraction (oooo), the critical drainage flow rate increases significantly the shorter is the channel. Similar to small initial inhomogeneities in figure 5(a) the inhomogeneity requires a sufficient channel length to reach a critical imbalance. If the channel is shorter, the inhomogeneity has to grow much faster to reach a critical imbalance. To accomplish that, the liquid drainage flow rate has to be substantially higher.
5. Stability analysis
The transition from stage II to stage III marks the primary instability in foam drainage. In stage II small inhomogeneities in the inflow of liquid, i.e. small 
$M_{\phi }$, are damped and homogenized. In stage III such inhomogeneities grow in the vertical direction. In order to analyse whether inhomogeneities of the liquid fraction grow with respect to the vertical direction, a stability analysis is performed. It is based on the drainage equation (1.2), which is repeated for the sake of a coherent representation:
\begin{equation} \frac{\partial \phi_l}{\partial t} = - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} = -\frac{\rho_f}{\eta_f} \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha \boldsymbol{g}) + \frac{\gamma}{2 \delta_b R_b \eta_f} \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\frac{\alpha}{\phi_l^{3/2}} \boldsymbol{\nabla} \phi_l \right). \end{equation}
The foam permeability 
$\alpha$ also depends on the liquid fraction and the interface mobility. The surfactant used here, SDS, yields a mobile interface. Thus, in the present investigations the dependency is 
$\alpha \propto \phi _l^{3/2}$ (Lorenceau et al. Reference Lorenceau, Louvet, Rouyer and Pitois2009). Consequently, the term 
${\alpha }/{\phi _l^{3/2}}$ in (5.2) becomes a constant and is replaced by 
$K_{\alpha }$:
\begin{equation} \frac{\partial \phi_l}{\partial t} = - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} = -\frac{\rho_f K_{\alpha} \boldsymbol{g}}{\eta_f} \boldsymbol{\cdot} \boldsymbol{\nabla} (\phi_l^{3/2} ) + \frac{\gamma K_{\alpha}}{2 \delta_b R_b \eta_f} \left(\boldsymbol{\nabla}^2 \right) \phi_l,\quad K_{\alpha} = \frac{\alpha}{ \phi_l^{3/2}}. \end{equation}Surfactants leading to an immobile interface are considered in appendix B.
To account for the above-mentioned effect of anisotropic drainage, Neethling (Reference Neethling2006) suggests using an additional term:
\begin{equation} \boldsymbol{q}_{aniso} = 0.5 q_x \varepsilon_{xy}\boldsymbol{e}_y \approx 0.5 \frac{\rho_f g}{\eta_f} \alpha \varepsilon_{xy} \boldsymbol{e}_y.\end{equation}
A vertical drainage flow rate 
$q_x$ results in a horizontal drainage flow rate 
$q_y$, which is proportional to the local shear deformation 
$\varepsilon _{xy}$. Of course, the same mechanism would also cause an additional vertical drainage flow when a horizontal drainage flow is present. However, since gravitationally driven vertical drainage typically exceeds capillary-driven horizontal drainage, the reverse effect is neglected here. The extended drainage equation reads
\begin{align} \frac{\partial \phi_l}{\partial t} = - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} = -\frac{\rho_f K_{\alpha} g}{\eta_f} [\boldsymbol{e}_x \boldsymbol{\cdot} \boldsymbol{\nabla} (\phi_l^{3/2} ) + 0.5 \boldsymbol{e}_y \boldsymbol{\cdot} \boldsymbol{\nabla} (\phi_l^{3/2} \varepsilon_{xy})] + \frac{\gamma K_{\alpha}}{2 \delta_b R_b \eta_f} \left( \boldsymbol{\nabla}^2 \right) \phi_l. \end{align} In order to derive the growth of an initial horizontal inhomogeneity of the liquid fraction, the normal mode approach is applied to (5.4). For that purpose, the liquid fraction 
$\phi _l(x,y)$ is decomposed into a constant value 
$\phi _0$ and small perturbations 
$\phi _{k}(x)$ with wavenumbers 
$k$ which evolve in the vertical direction (cf. figure 6):
\begin{equation} \phi_l(x,y) = \phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky}. \end{equation}It is important to note that no temporal dependency is incorporated. This is motivated by our experimental observations, where the elastic deformation and the liquid fractions reached a steady state. As already mentioned, even in the case of highly symmetric liquid inflow, we never succeeded in generating a perfectly homogeneous liquid content in the conductivity sensor below the porous cylinders. Randomly, we found that the liquid content was higher at the left or right side of the channel. We assume that the asymmetry of the static foam structure which is in contact with the cylinders always generates small inhomogeneities. These slight inhomogeneities then remain steady over time. Consequently, in our experiments we always had steady but slightly asymmetric inflow conditions, which is reflected by the ansatz function.
Sketch of the considered system for linear stability analysis.
While the average foam weight is compensated for by a vertical pressure gradient, the perturbations lead to an imbalance in the gravitational force:
\begin{equation} \boldsymbol{f}_{g}(x,y) = \rho_f g [\phi_l(x,y)-\phi_0] \boldsymbol{e}_x = \rho_f g \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} \boldsymbol{e}_x . \end{equation}
For the stability analysis, we assume an infinitely long (high) channel. Consequently, this imbalance must not accumulate over the channel height. In each horizontal slice, the gravitational imbalance is compensated for by a local shear stress. The vertical force balance on a foam volume stretching between the left wall and the position 
$y=\tilde {y}$ yields
\begin{equation} \tau_{xy}(x,y=\tilde{y}) - \tau_{xy}(x,y=0) = - \int_0^{\tilde{y}} \boldsymbol{e}_x \boldsymbol{\cdot} \boldsymbol{f}_{g} \,{\textrm{d}y} = - \rho_f g \phi_{k}(x) \frac{1}{\textrm{i}k}\textrm{e}^{\textrm{i}k\tilde{y}}. \end{equation}At the vertical boundaries of the channel, a liquid film results in negligible static friction between the foam and wall. Consequently, the shear stress has to fulfil the boundary conditions
\begin{equation} \tau_{xy}(x,y=0)=0,\quad \tau_{xy}(x,y=B)=0. \end{equation}Hence, the shear stress equals
\begin{equation} \tau_{xy}(y) = - \rho_f g \phi_{k}(x) \frac{1}{\textrm{i}k}\textrm{e}^{\textrm{i}ky}. \end{equation} The shear modulus 
$G$ of foam depends on the liquid fraction. It can be estimated (Mason, Bibette & Weitz Reference Mason, Bibette and Weitz1995) by (5.10). In the present case, critical liquid fractions of 
$\phi _0 \leq 1\,\%$ are considered. In this small range the influence can be neglected, yielding
\begin{equation} G= \frac{\tau_{xy}}{\varepsilon_{xy}} = 1.4 (1-\phi_0)(0.36-\phi_0) \frac{\gamma}{R_b} \approx 0.5 \frac{\gamma}{R_b}.\end{equation}Combining (5.10) with (5.9) leads to the shear strain:
\begin{equation} \varepsilon_{xy}= \frac{\tau_{xy}}{G} = - K_1 \phi_{k}(x) \frac{1}{\textrm{i}k}\textrm{e}^{\textrm{i}ky},\quad \mbox{with}\ K_1=2.0 \frac{R_b}{\gamma} \rho_f g. \end{equation} Considering a steady liquid distribution, the term 
$\partial \phi _l / \partial t$ in (5.4) equals zero. Feeding the ansatz (5.5) into (5.4), neglecting the influence of capillary pressure on the vertical liquid transport and assuming that only small perturbations 
$\phi _{k}\ll \phi _0$ occur leads to
\begin{equation} 0 = \frac{3}{2}\frac{\rho_f g}{\eta_f}\sqrt{\phi_0} \phi_{k}'(x) + \left[0.5 \frac{\rho_f g}{\eta_f} K_1 \phi_0^{3/2} - \frac{\gamma k^2 }{2 \delta_b R_b \eta_f} \right] \phi_{k}(x). \end{equation}For a detailed derivation we refer the reader to appendix A. This ordinary differential equation is solved by
\begin{equation} \phi_{k}(x)=\phi_{k}(x=0)\,\textrm{e}^{Ax}. \end{equation}
The exponent 
$A$ represents the growth rate with respect to the vertical position 
$x$ and is given by
\begin{equation} A = \frac{2.0}{3.0 } \frac{R_b}{\gamma} \rho_f g \phi_0 - \frac{ \gamma k^2 }{3 \delta_b R_b \rho_f g \phi_0^{1/2}}. \end{equation}
Equations (5.13) and (5.14) show that small imbalances in 
$\phi _l$ grow exponentially with respect to the vertical distance 
$x$ provided that 
$A > 0$. Inspecting (5.14), this is the case for sufficiently small horizontal wavenumbers 
$k$. The smallest possible wavenumber in a channel that fulfils the boundary conditions (5.8) is 
$k_{min}={\rm \pi} /B$. At the critical point, 
$A=0$ and the liquid fraction 
$\phi _{0}=\phi _{crit}$. This yields the stability criterion for drainage in an infinitely long vertical foam channel of width 
$B$:
\begin{equation} \phi_{crit}^{3/2} = \frac{\gamma^2}{ \delta_b R_b^2 } \frac{0.5}{\rho_f^2 g^2} \frac{{\rm \pi}^2}{B^2}, \end{equation}
above which initial inhomogeneities grow exponentially. In the present case, we obtain 
$\phi _{crit} \approx 0.74\,\%$. This value corresponds well with the experimental data for the transition from stage II to stage III. The critical liquid fraction for the onset of the shear deformation of the foam in figure 3(c) is 
$\phi _{crit} \approx 0.6\,\%$.
6. Simulations
In order to investigate the interaction of stress, deformation and liquid distribution at the transition from stage II to stage III, a two-dimensional, phase-averaging numerical simulation of the experimental set-up is employed. The code computes the unsteady drainage equation, (5.4), in a finite-volume discretization with Euler-explicit time integration. In each time step, the simulation solves iteratively for the linear-elastic deformation 
$U$ and the stress tensor 
$\hat {\tau }$ of the foam with zero strain in the third dimension:
\begin{gather} \varepsilon_{xy} = \frac{1}{2} \left( \frac{\partial U_x}{\partial y} + \frac{\partial U_y}{\partial x} \right),\quad \varepsilon_{xx}=\frac{\partial U_x}{\partial x},\quad \varepsilon_{yy}=\frac{\partial U_y}{\partial y}, \end{gather}
\begin{gather}\tau_{xy}= \tau_{yx} = 2 G \varepsilon_{xy}, \end{gather}
\begin{gather}\tau_{xx}= 3 G K_{\nu} [(1-\nu) \varepsilon_{xx} + \nu \varepsilon_{yy}],\quad K_{\nu}=\frac{1}{\left(1+\nu \right) \left( 1-2\nu \right) }, \end{gather}
\begin{gather}\tau_{yy}= 3 G K_{\nu} [ ( 1-\nu) \varepsilon_{yy} + \nu \varepsilon_{xx}], \end{gather}
\begin{gather}\boldsymbol{\nabla} \hat{\tau} + \boldsymbol{f}_{g}=0, \end{gather}
balancing the gravitational load 
$\boldsymbol {f}_{g}$ from the liquid fraction contained. The shear modulus 
$G$ depends on the local liquid fraction according to (5.10). The Poisson ratio 
$\nu$ for incompressible foam would equal 0.5 and the compression module 
$K_{\nu }$ would be infinite. To avoid numerical problems resulting from an infinite compression module 
$K_{\nu }$, 
$\nu$ is set to the artificial value of 0.49. It has been found that small variations of 
$\nu$ have negligible effects on the results. The resulting local shear strain 
$\varepsilon _{xy}$ is fed back into the drainage equation. The pertinent boundary conditions are given in table 1.
Boundary conditions for the elastic simulations.
At the sidewalls, zero vertical stress, zero horizontal liquid flux and zero horizontal strain are imposed. At the top wall, zero strain and a prescribed vertical drainage flux 
$q_x(y)=q_0(1+q' \cos (y {\rm \pi}/B))$ with inhomogeneity 
$0 < q' <1$ are imposed. At the bottom, 20 % liquid fraction, zero horizontal stress and a hydrostatic vertical stress as a function of the displacement 
$U_x$ are imposed. Figure 7(a) shows the numerical results for a case of unstable drainage marked by a green star in figure 7(b). The initial distortion of the liquid fraction 
$\phi _l$ grows in the downward direction (figure 7a). The corresponding weight imbalance causes an inhomogeneous downward displacement 
$U_x$. Due to the hydrostatic boundary condition, the inhomogeneity of displacement is reduced close to the bottom (figure 7a). The horizontal gradient of 
$U_x$ corresponds to a shear strain 
$\varepsilon _{xy}$ (figure 7a), which is highest in the centre of the channel. The shear drives a horizontal drainage flow 
$q_y=q_{aniso}-q_{cap}$ (figure 7a), which is the anisotropic flow 
$q_{aniso}$ according to (5.3) minus the capillary-driven drainage flow 
$q_{cap}$ according to the second term on the right-hand side of (1.2). The horizontal drainage flow 
$q_y$ is positive in regions of high shear, feeding the imbalance. Close to the top and bottom, the boundary conditions inhibit the shear deformation, and a negative horizontal drainage flow occurs due to capillary forces. The ratio of the local strain 
$\varepsilon _{xy}$ to the local critical yield strain 
$\varepsilon _{y}=0.3 (0.36-\phi _l)$ (Saint-Jalmes & Durian Reference Saint-Jalmes and Durian1999) is highest in the lower part of the channel (figure 7a), where experimental investigations have found the onset of CI (figure 4).
Numerical simulation. (a) Liquid fraction, vertical displacement, shear strain, horizontal drainage flow rate and ratio between shear strain and local yield strain for an unstable case, marked in (b), at which the yield strain is exceeded. (b) Dependence of the liquid fraction 
$\phi _{CI}$ and drainage flow rate 
$q_{CI}$ for the onset of CI on liquid moment at the introduction for different initial distributions of 
$\phi _l$ and channel lengths 
$L$. The horizontal dash-dotted line is a guide to the eye and marks 
$\phi _{crit}$ for the transition from stage II to stage III.
Now, the initial inhomogeneity and total drainage flow rate are varied. If the local strain 
$\varepsilon _{xy}$ in the converged simulation exceeds the theoretical limit 
$\varepsilon _{y}$ at any place, the case is considered to show yielding and CI (stage IV). Figure 7(b) depicts the limit for the onset of CI for different channel lengths and initial moments of liquid distribution. It shows very good agreement with the experimental findings in figure 5. For high moments of the liquid fraction 
$M_{\phi }$, CI occurs at 
$\phi _{CI}=0.65\,\%$ liquid fraction and 
$q_{CI}=70\ \mathrm {\mu }$m s
$^{-1}$ drainage flow rate, regardless of the channel length. A small moment of 1.5 mm in a long channel causes CI at approximately 
$\phi _{CI}=1.7\,\%$ liquid fraction and 
$q_{CI}=250\ \mathrm {\mu }$m s
$^{-1}$ drainage flow rate. For shorter channels, the critical values are increased to 
$\phi _{CI}=2.7\,\%$ liquid fraction and 
$q_{CI}=500\ \mathrm {\mu }$m s
$^{-1}$ drainage flow rate.
7. Discussion
In this work, we have proven the mechanism for the onset of CI. Our experimental and numerical findings support the concept of an instability in forced drainage due to anisotropic drainage which was initially predicted by Neethling (Reference Neethling2006). The dependencies of the critical liquid fraction 
$\phi _{crit}$ on bubble radius, surface tension, gravity and channel size are identical to the findings of Neethling (Reference Neethling2006) for the assumption that the critical liquid fraction is negligible compared with the jamming point liquid fraction. However, the stability criterion that Neethling derived differs significantly from ours. In our case it would produce a critical liquid fraction of approximately 0.2 %, which is considerably lower than our experimental findings. We extend his theory by adding the concept of a growing inhomogeneity with respect to the vertical direction and by describing the CI as a secondary instability superimposing on that.
Initial inhomogeneities due to the liquid addition cause a horizontal gradient in the local gravitational force. This gradient is compensated for by a shear deformation of the foam structure. The shear deformation deflects the vertical drainage flow in the horizontal direction, increasing the imbalance. At the same time, capillary forces tend to reduce the liquid imbalance. If a critical liquid fraction is exceeded, the horizontal deflection exceeds the capillary effect and the initial inhomogeneity grows with respect to the vertical distance. This proposed mechanism is supported by the agreement of the critical liquid fraction from experiments, simulations and stability analysis, by the measurement and simulation of the distribution of inhomogeneous vertical displacement and by the dependency of the occurrence of CI on channel length and initial inhomogeneity.
In the experimental results (figure 3) and in the simulation (figure 7), we have observed the occurrence and vertical growth of such a shear deformation above a critical liquid fraction. This process does not necessarily involve a steady flow of the foam. Only if the channel is long and wide enough does the growing imbalance at some point exceed the yield stress of the foam and CI occurs below that point. However, the onset of CI is only a secondary instability of the primary instability, consisting of the shear deformation at 
$\phi _l > \phi _{crit}$. Our experimental and numerical observations of the CI have demonstrated that the onset of CI takes place close to the bottom of the channel (figure 4), where the imbalance is highest. Moreover, the liquid fraction 
$\phi _{CI}$ for the onset of CI in the case of an inhomogeneous introduction of liquid in figure 5(a) is 
$\phi _{CI} \approx 0.65\,\%$, which is very close to 
$\phi _{crit} \approx 0.74\,\%$. In line with the discussion above, a reduced channel height increases the critical liquid fraction for the onset of CI (figure 5b). Then again, a strong inhomogeneity in the liquid introduction reduces the channel height required to reach the onset of CI (figure 5b).
From our point of view, the presented findings on elastic deformation also rule out other existing explanations of CI. The idea that dilatancy might cause CI is disproved by the fact that the strongest shear is found in the centre of the channel (figure 3b). Therefore, dilatancy would cause the accumulation of liquid in the centre and not close to the channel side. Also, the idea of a sudden scarification of individual Plateau borders does not explain why the onset of CI is close to the bottom of the channel and why the size of the CI grows with increasing drainage flow rate.
The presented experiments on both the CI with strongly inhomogeneous liquid inflow and on the inhomogeneous elastic deformation together with the elastic simulation all consistently produce a critical liquid fraction of approximately 
$\phi _{crit}=0.65\,\%$. This critical liquid fraction marks the transition from stage II to stage III. In stage II, the static foam stretches homogeneously and initial inhomogeneities in liquid fraction are damped with respect to the vertical direction. In stage III, the static foam stretches inhomogeneously and initial inhomogeneities in liquid fraction grow with respect to the vertical direction. Due to the high scattering between experimental runs there is a certain uncertainty in the critical value resulting in an interval of 
$\phi _{crit} \approx 0.6\,\%$ to 0.65 %. Also, there is an uncertainty in distinguishing the different stages. In particular, it is difficult to assign the transition from stage III to stage IV to a fixed value of 
$\phi _{CI}$, because in some cases the foam network starts to flow and comes to a halt again. In figure 5 in some cases of high initial inhomogeneity, the CI sets in at liquid fractions 
$\phi _{CI}$ even below the theoretical critical liquid fraction 
$\phi _{crit}$. This might be due to such uncertainties or due to nonlinear effects under high imbalance and liquid fraction gradients.
Carrying out a linear stability analysis, we have established a stability criterion above which initial inhomogeneities in the liquid fraction grow with respect to the vertical position. The resulting critical liquid fraction of 
$\phi _{crit}=0.74\,\%$ is close to the experimental result of 
$\phi _{crit} \approx 0.6\,\%$ to 0.65 %, supporting the proposed mechanism. Our critical liquid fraction is a factor of five below the critical liquid fraction reported in Hutzler et al. (Reference Hutzler, Weaire and Crawford1998) for the onset of CI in cylindrical channels. One reason for the discrepancy is our significantly larger channel width 
$B$, which reduces the critical liquid fraction according to (5.15). Also, Hutzler et al. (Reference Hutzler, Weaire and Crawford1998) investigated the critical liquid fraction 
$\phi _{CI}$ for the onset of CI (transition from stage III to stage IV). By comparing figures 3 and 4 for homogeneous water inflow 
$\phi _{CI} \approx 1.05\,\%$ is demonstrated to be larger than 
$\phi _{crit} \approx 0.65\,\%$ for the transition from stage II to stage III. And finally, our channel is one order of magnitude longer, giving the inhomogeneity more distance to grow. Thus, smaller liquid fractions can lead to the onset of CI. The influence of the channel length on 
$\phi _{CI}$ is documented in figures 5(b) and 7(b).
We do not believe that the deviations between our critical liquid fraction and previous findings result from the small initial inhomogeneity of the liquid introduction. Despite considerable efforts, a small inhomogeneity of 
$M_{\phi } \approx 1\ \textrm {mm}$ has always been present. Presumably, this is a feature of the water addition to a random foam structure and cannot be avoided. Possibly, a similar inhomogeneity was also present in former studies. But since former studies did not measure the horizontal liquid distribution, this aspect may have gone unnoticed. Despite the different geometry and critical liquid fractions, the dependency of our stability criterion on the bubble size 
$\phi _{crit} \propto R_b^{-4/3}$ shows the same trend as earlier findings
This work takes only monodisperse foam into account. In industrial applications, foam is usually rather polydisperse. Polydispersity modifies the ingredients of the instability mechanism, namely the foam permeability, the anisotropic drainage term, the shear modulus and the capillary pressure term in the drainage equation. But none of these ingredients would be cancelled or negated in the case of polydispersity. Thus, we assume that our mechanism does apply to polydisperse foam as well, but the critical liquid fraction might change considerably.
For the stability analysis, 
$\alpha \propto \phi _l^{3/2}$ is assumed. This is only valid for some surfactants. However, as demonstrated in appendix B, the stability criterion is identical for 
$\alpha \propto \phi _l^{2}$. Since permeability acts on all terms on the right-hand side of the drainage equation, its effect is neutralized in steady state.
The stability analysis shows that any mode 
$k$ becomes unstable at a certain critical liquid fraction, which scales with 
$\phi _{crit} \propto k^{4/3}$ (see (5.15)). But, according to the linear stability analysis, the smallest modes show the highest growth rate in space and presumably suppress higher modes in many cases. Vera et al. (Reference Vera, Saint-Jalmes and Durian2000) have observed convective structures with higher modes in systems with bubble diameters below 
$100\ \mathrm {\mu }$m. In the present study, the case with the liquid being added by the two outermost porous cylinders (oxxo) corresponds to a mode 
$k=2 {\rm \pi}/B$ and shows superior stability (see figure 5a). However, the resulting convection roll then was not a double roll but a single roll filling the full cross-section (see figure 1a), corresponding to 
$k={\rm \pi} /B$.
While the one-dimensional drainage equation without anisotropic drainage does not show any signs of instabilities (Verbist, Weaire & Kraynik Reference Verbist, Weaire and Kraynik1996), we found an instability in two-dimensional steady drainage which presumably also occurs in three-dimensional steady drainage. This has strong implications for any experiment with columns of liquid foam under forced drainage. In such experiments, one should either measure the horizontal liquid fraction distribution or observe the elastic deformation of the foam.
The proposed mechanism of unstable drainage relies on the anisotropic deformation of the Plateau border network. In the case of very small bubbles or even solid grains and particles, (5.3) presumably loses its validity. This should be further investigated by drainage experiments under prescribed shear, and by simulations. It is possible that neutron imaging (Heitkam et al. Reference Heitkam, Rudolph, Lappan, Sarma, Eckert, Trtik, Lehmann, Vontobel and Eckert2018) or a combination of the Surface Evolver (Brakke Reference Brakke1992) with drainage equation could yield valuable insights into that phenomenon.
Acknowledgement
P. Jähnigen and B. Legrady supported the experiments.
Funding
The support of the Deutsche Forschungsgesellschaft (HE 7529/1-1) is gratefully acknowledged. Funding by the European Union's Horizon 2020 research and innovation programme under grant agreement no. 821265 is gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Solving the partial differential equation
This appendix demonstrates how to solve the partial differential equation given in (5.4). Note that (5.4) considers anisotropic drainage only in the horizontal direction, not in the vertical direction. Considering a steady liquid distribution, the term 
$\partial \phi _l / \partial t$ equals zero. Dividing by 
$K_{\alpha }$ yields
\begin{equation} 0 = - \frac{\rho_f g}{\eta_f} \left[ \frac{\partial \phi_l^{3/2}}{\partial x} + 0.5 \frac{\partial \varepsilon_{xy} \phi_l^{3/2}}{\partial y} \right] + \frac{\gamma}{2 \delta_b R_b \eta_f} ( \boldsymbol{\nabla}^2 ) \phi_l . \end{equation}
Feeding the ansatz (5.5) for 
$\phi _l$ and (5.11) for 
$\varepsilon _{xy}$ yields
\begin{align} 0 &=- \frac{\rho_f g}{\eta_f} \frac{\partial [\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky}]^{3/2}}{\partial x} \nonumber\\ &\quad + \frac{0.5 \rho_f g}{\eta_f} \frac{\partial K_1 \phi_{k}(x) \frac{1}{\textrm{i}k}\textrm{e}^{\textrm{i}ky} [\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky}]^{3/2}}{\partial y} \nonumber\\ &\quad + \frac{\gamma}{2 \delta_b R_b \eta_f} ( \boldsymbol{\nabla}^2) [\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky}]. \end{align}Carrying out the derivatives yields
\begin{align} 0 &= - \frac{3}{2}\frac{\rho_f g}{\eta_f} [ \phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} ]^{1/2} \phi_{k}'(x) \,\textrm{e}^{\textrm{i}ky} \nonumber\\ &\quad + \frac{0.5 \rho_f g}{\eta_f} \left[\vphantom{\left.\quad +\, \frac{3}{2} K_1 [\phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} ]^2 [ \phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} ]^{1/2}\right]} K_1 \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} [ \phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} ]^{3/2}\right.\nonumber\\ &\quad + \left.\frac{3}{2} K_1 [\phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} ]^2 [ \phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} ]^{1/2}\right] \nonumber\\ &\quad + \frac{\gamma}{2 \delta_b R_b \eta_f} \textrm{e}^{\textrm{i}ky} (\phi_{k}''(x) - k^2 \phi_{k}(x) ). \end{align} When the stability is analysed regarding small disturbances, 
$\phi _k$ is assumed to be small compared with 
$\phi _0$, yielding
\begin{align} 0 &= - \frac{3}{2}\frac{\rho_f g}{\eta_f} \left[ \phi_0 \right]^{1/2} \phi_{k}'(x) \,\textrm{e}^{\textrm{i}ky} \nonumber\\ &\quad + \frac{0.5 \rho_f g}{\eta_f} [ K_1 \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} [\phi_0 ]^{3/2}] \nonumber\\ &\quad + \frac{\gamma}{2 \delta_b R_b \eta_f} \textrm{e}^{\textrm{i}ky} (\phi_{k}''(x) - k^2 \phi_{k}(x) ). \end{align}
Now, one needs to compare 
$\phi _{k}''(x)$ and 
$k^2 \phi _{k}(x)$ from the third term of (A4). We assume an exponential growth of 
$\phi _{k}(x)$ in the 
$x$ direction (which we will find later):
\begin{equation} \phi_{k}(x) = \phi_{k}(x=0) \,\textrm{e}^{A x} \rightarrow \phi_{k}''(x) = \phi_{k}(x=0)A^2 \,\textrm{e}^{A x} = A^2 \phi_{k}(x). \end{equation}
In a linear stability analysis, one investigates the case where the perturbation just starts to grow. Thus, the growth rate 
$A$ is close to zero, i.e. much smaller than the wavenumber 
$k$. Consequently it holds that
\begin{equation} \phi_{k}''(x) = A^2 \phi_{k}(x) \ll k^2 \phi_{k}(x) \end{equation}
and we can neglect 
$\phi _{k}''(x)$ in (A4).
This yields
\begin{align} 0 &= - \frac{3}{2}\frac{\rho_f g}{\eta_f} \left[ \phi_0 \right]^{1/2} \phi_{k}'(x) \,\textrm{e}^{\textrm{i}ky} \nonumber\\ &\quad + \frac{0.5 \rho_f g}{\eta_f} [K_1 \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} [\phi_0]^{3/2}] \nonumber\\ &\quad - \frac{\gamma}{2 \delta_b R_b \eta_f} \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} (k^2 \phi_{k}(x)). \end{align} Dividing by 
$\textrm {e}^{\textrm {i}ky}$ and sorting the terms yields
\begin{equation} 0 = - \frac{3}{2}\frac{\rho_f g}{\eta_f} \phi_0^{1/2} \phi_{k}'(x) + \left(\frac{0.5 \rho_f g}{\eta_f} K_1 \phi_0^{3/2} - \frac{k^2 \gamma}{2 \delta_b R_b \eta_f} \right) \phi_{k}(x) \end{equation}and
\begin{equation} \phi_{k}'(x) = \left( \frac{1}{3.0} K_1 \phi_0 - \frac{k^2 \gamma }{3 \delta_b R_b \rho_f g } \phi_0^{-1/2} \right) \phi_{k}(x) . \end{equation}This ordinary differential equation is solved by
\begin{equation} \phi_{k}(x)=\phi_{k}(x=0)\,\textrm{e}^{Ax},\quad A = \frac{2.0}{3.0 } \frac{R_b}{\gamma} \rho_f g \phi_0 - \frac{ k^2 \gamma}{3 \delta_b R_b \rho_f g \phi_0^{1/2}}. \end{equation} Consequently, small imbalances of 
$\phi _l$ grow exponentially over the vertical distance 
$x$, if 
$A > 0$, i.e. if the horizontal wavenumber 
$k$ is sufficiently small. The smallest possible wavenumber in a channel that fulfils the boundary conditions (5.8) is 
$k_{min}={\rm \pi} /B$. This yields the stability criterion for drainage in an infinitely long vertical foam channel (see (5.15)):
\begin{equation} \phi_{crit}^{3/2} = \frac{\gamma^2}{ \delta_b R_b^2 } \frac{0.5}{\rho_f^2 g^2} \frac{{\rm \pi}^2}{B^2}, \end{equation}
above which initial inhomogeneities grow exponentially. In the present case, this yields 
$\phi _{crit} \approx 0.74\,\%$.
Appendix B. Other surfactants
In the above derivation, the relation 
$\alpha =K_{\alpha } \phi _l^{3/2}$ has been used, which is valid for the employed surfactant, SDS. However, for many other surfactants the liquid permeability 
$\alpha$ scales with 
$\alpha \propto \phi _l^{2}$. This changes the derivation, as the term 
$\alpha / \phi _l^{3/2}$ in (5.2) does not vanish. Thus, the second derivative is slightly more complex:
\begin{equation} (\boldsymbol{\nabla}^2) \phi_l \rightarrow \boldsymbol{\nabla} \phi_l^{1/2} \boldsymbol{\nabla} \phi_l , \end{equation}and also the exponents in (A1) change, yielding
\begin{equation} 0 = \frac{\rho_f g}{\eta_f} \left[ -\frac{\partial \phi_l^{\boldsymbol{2}}}{\partial x} - 0.5 \frac{\partial \varepsilon_{xy} \phi_l^{\boldsymbol{2}}}{\partial y} \right] + \frac{\gamma}{2 \delta_b R_b \eta_f} \boldsymbol{\nabla} \phi_l^{1/2} \boldsymbol{\nabla} \phi_l . \end{equation}Carrying out the derivations yields for the third term in (A3)
\begin{align} & \frac{\gamma}{2 \delta_b R_b \eta_f} \textrm{e}^{\textrm{i}ky} (\phi_{k}''(x) - k^2 \phi_{k}(x) ) \nonumber\\ &\quad \rightarrow (\phi_{k}'(x) \,\textrm{e}^{\textrm{i}ky})^{2} \tfrac{1}{2} (\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky})^{-1/2} + \textrm{i}k \phi_{k}''(x) \,\textrm{e}^{\textrm{i}ky} (\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky})^{1/2} \nonumber\\ &\qquad + (\textrm{i}k \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky})^{2} \tfrac{1}{2} (\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky})^{-1/2} \nonumber\\ &\qquad + \left(\textrm{i}k \right)^2 \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky} (\phi_0 + \phi_{k}(x) \,\textrm{e}^{\textrm{i}ky})^{1/2}. \end{align}
Under the assumption of small distortions 
$\phi _k \ll \phi _0$ and neglecting terms with 
$\phi _k^2$, this gives
\begin{equation} \frac{\gamma}{2 \delta_b R_b \eta_f} \left[ \phi_0 \right]^{1/2} \textrm{e}^{\textrm{i}ky} (\phi_{k}''(x) - k^2 \phi_{k}(x)), \end{equation}
the same term as in (A4), only with an additional factor 
$[\phi _0 ]^{1/2}$. Consequently, all the following steps can be carried out similarly, eventually yielding the growth rate:
\begin{equation} A = \frac{1}{2.0 } \frac{R_b}{\gamma} \rho_f g \phi_0 - \frac{ \gamma k^2 \phi_0^{1/2} }{4 \delta_b R_b \rho_f g \phi_0^{1}} = \frac{1}{2.0 } \frac{R_b}{\gamma} \rho_f g \phi_0 - \frac{ \gamma k^2 }{4 \delta_b R_b \rho_f g \phi_0^{1/2}}. \end{equation}The corresponding stability criterion equals
\begin{equation} \phi_{crit}^{3/2} = \frac{\gamma^2}{ \delta_b R_b^2 } \frac{0.5}{\rho_f^2 g^2} \frac{{\rm \pi}^2}{B^2}, \end{equation}which is identical to the stability criterion (A11) for SDS above. Consequently, the type of surfactant has no influence on the stability analysis.
Open access
