A computational method to simulate mono- and poly-disperse two-dimensional foams flowing in obstructed channel

A modified phase-field model is presented to numerically study the dynamics of flowing foam in an obstructed channel. The bubbles are described as smooth deformable fields interacting with one another through a repulsive potential. A strength of the model lies in its ability to simulate foams with wide range of gas fraction. The foam motion, composed of about hundred two-dimensional gas elements, was analyzed for gas fractions ranging from 0.4 to 0.99, that is below and beyond the jamming transition. Simulations are preformed near the quasi-static limit, indicating that the bubble rearrangement in the obstructed channel is primarily driven by the soft collisions and not by the hydrodynamics. Foam compression and relaxation upstream and downstream of the obstacle are reproduced and qualitatively match previous experimental and numerical observations. Striking dynamics, such as bubbles being squeezed by their neighbors in negative flow direction, are also revealed at intermediate gas fractions.


Introduction
The flow of foam and froth plays a major role in many industrial applications such as mineral recovery (Lecrivain et al. 2015), food processing (Indrawati et al. 2008), fabric dyeing (Hussain and Wahab 2018), production of light-weight materials (Masi et al. 2014) and insulation components (Koniorczyk et al. 2017). Despite its industrial significance, many of the mechanisms, such as soft collisions, governing the foam dynamics remain to be elucidated (Cohen-Addad et al. 2013). One reason for this insufficient knowledge is Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00397-021-01288-y) contains supplementary material, which is available to authorized users. Professur für Bildgebende Messverfahren für die Energie-und Verfahrenstechnik, Technische Universität Dresden, 01062 Dresden, Germany that most established flow measurement techniques are not directly applicable to foam flow. Consequently, most experiments are conducted in two-dimensional configurations employing optical observation (Chen et al. 2019). Novel measurement techniques applied to foam, such as ultrasound Doppler velocimetry (Nauber et al. 2018), neutron imaging , positron emission particle tracking (Waters et al. 2008), nuclear magnetic resonance (Stevenson et al. 2010), or X-raying (Lappan et al. 2020) are not yet commonly established and suffer from limitations.
Numerical simulations are a valid alternative to investigate flowing foams in three-dimensional domains. Large advances in this field have been achieved with the Surface Evolver (Brakke 1992), which discretizes the air-liquid interfaces of the foam with a triangle mesh. Under a bubble volume conservation constraint, it then iteratively redistributes the triangle nodes to find the energy minima in terms of foam shape. Quasi-static dry foam motion around a spherical obstacle was simulated (Cox et al. 2006). The authors derived from two-dimensional simulations pressure fields, plastic events and drag and lift forces acting on the obstacle. Two-dimensional wet foam shearing in unobstructed channel was also simulated with the Surface Evolver (Jing et al. 2015). Bubble displacements in the transverse shear direction and elastic-plastic deformation of the foam were analyzed. Further numerical methods can also be found in the literature. The vertex model, where the dry foam is represented by a set of polygons, has been used to study foam shearing (Kern et al. 2004;Cantat 2011). In the latter work, two-dimensional simulations of foam shearing in an unobstructed channel were performed. The authors investigated the foam stress in the non-quasi static regime and found that the number of topological rearrangements, so-called T1 events, correlated with bubble elongation. The two-dimensional soft-disk model has also been applied to simulate shear foam flows (Durian 1995(Durian , 1997Sexton et al. 2011). In this model, the disk-like bubbles experience a repulsive spring force whenever an overlap occurs. In the work of Sexton et al. (2011), two-dimensional simulations were performed for a relatively wide range of shear rates and two distinct regimes with qualitatively different bubble dynamics reported. Lattice Boltzmann simulations (Benzi et al. 2009;Dollet et al. 2015) have been performed to determine the rate of plastic events in a Poiseuille foam flow. More recently, another approach was adopted by Kähärä (2017), who discretized each two-dimensional bubble into a series of segments subjected to surface tension, pressure gradient and contact visco-elastic forces. In this work, two dimensional foam shearing was also studied. These type of simulations can compute the structural response of the foam, stress, strain and velocity distributions. However, these methods do not behave properly for wide range of gas fraction within the foam. Thus, these models are not applicable if bubbles are not in contact with one-another. This typically occurs below the jamming point, which for randomly packed, poly-disperse bubbles corresponds to gas fractions ε 2D ≈ 0.85 (Katgert and van Hecke 2010;Desmond et al. 2013) and ε 3D = 0.64 (Cohen-Addad and Höhler 2014; Heitkam and Fröhlich 2019) in two and three dimensions, respectively.
Close to the jamming point, Heitkam and Fröhlich (2019) carried out numerical simulations of foam dynamics with the immersed-boundary method (Mittal and Laccarino 2005;Uhlmann 2008;Schwarz et al. 2015). These simulations also solved the Navier-Stokes equation for the interstitial liquid flow and thus, allowed to compute the formation of foam from rising bubbles. By suggesting a bubbles collision model, Heitkam and Fröhlich (2019) also investigated the shearing of bubble clusters. These above simulations are however restricted to spherical bubbles and, thus, loose validity if large bubble deformation occurs. To alleviate this, a modified phase-field model is here suggested as a numerical alternative to the immersed-boundary method. This model hence bridges the gap between the existing Surface Evolver and the available immersedboundary methods. On the one hand, it deals with interacting bubbles subject to large deformations and; on the other hand, with the movement of individual bubbles carried away in regions of low gas fraction. The method is here applied to mono-and poly-disperse foam flowing in an obstructed channel at gas fractions below, near and above the jamming point. Phase-field models constitute a class of efficient computational methods that allows one to investigate three-dimensional foam dynamics. In the following, we restrict however ourselves to two-dimensional systems and demonstrate the ability of the method to handle collision-dominated cases. Three-dimensional systems will be included in future works.

Foam model
The general idea behind the phase field modeling is to introduce a dimensionless identification function, also called order parameter, that continuously varies over a thin interfacial region. The sharp fluid interface between the gas and liquid phases is hence smeared out with a thin but non-zero transitional region. We refer the reader to the review article by Kim (2012) for further information on the phase-field modeling of multi-component flows. The model presently employed is a modified version of the phase field model by Nestler et al. (2008) that has been used in the field of biology, where cell-cell interactions are key in understanding the migrating dynamics of cell colonies (Alert and Trepat 2020). Each bubble i is here associated with an identification function φ i (x, t), where x = (x, y) is the spacial coordinate and t the time. As one moves from the bubble inner region to the outer region, φ i undergoes a smooth but rapid transition from unity to zero over the interfacial width ξ (Borcia et al. 2017). Figure 1 describes a typical smoothly transitioning identification function along the two coordinate axes crossing the bubble center of mass. The equilibrium foam state was obtained in a square box with a periodic boundary condition in the two spatial directions. Each bubble "sharp" interface, shown in blue in the figure, is hereafter defined by the isoline φ i = 0.5.
To compute the dynamic foam flow, each identification function is updating in time and space by solving where u is a velocity field, τ a response time associated with the foam dynamics and f b , respectively, f r the bulk and repulsive terms discussed below. The advection term u ·∇φ i on the left-hand side of Eq. 1 transports the foam, while the right-hand side term ensures an energetically stable foam conformation. The bulk term, which ensures φ i = 1 inside the bubble and zero elsewhere, is given by The area of each two-dimensional bubble is maintained through the correction φ c = 1/2 + α(Ã i /(Ã i ) 0 − 1)  (Nonomura 2012), where α > 0 is a growth rate,Ã i (t) = φ 2 i dr the modified area at time t (Mueller et al. 2019), dr an infinitesimally small area element and (Ã i ) 0 a constant value initially set. This correction causes the bubble area to expand, respectively shrink, until the conditionÃ i = (Ã i ) 0 is satisfied. The tilde is hereafter introduced to distinguish between the modified areaÃ i and its standard counterpart A i = φ i dr mostly used in the result section. The repulsive term in Eq. 1, which prevents bubble overlapping (Mueller et al. 2019), is given by where β > 0 controls the strength of the repulsion. The last term −ξ 2 ∇ 2 φ i is an interfacial contribution resulting in an energy excess across the smoothly spreading interface . A link between the tuning parameter β and experiments involving foam films stabilized by surface active agents (Andersson et al. 2010) can well be established. As the separating distance h between two bubble surfaces diminishes, surface active agents adsorbed at those interfaces rearrange resulting in the formation of thick adsorption layers (Peng et al. 2020). This surface interaction, expressed by the disjoining pressure in experimental tests, becomes significant when the film thickness drops below h < 100 nm (Stubenrauch and von Klitzing 2003). In a similar fashion, the parameter β controls the separating distance h, typically a few mesh elements, over which two neighboring bubbles experience strong repulsion. Linking the parameter α to any physical foam-associated quantity is however difficult. The term α is a pure numerical parameter that turns the Allen-Cahn Equation 1 into a conservative form (Aihara et al. 2019). It has been shown that, for a single bubble, the term α can safely be removed by using a Cahn-Hilliard formulation (Jeong and Kim 2017). This comes however at a cost that would be the arising complex biharmonic operator (∇ 4 ) in the interfacial term.
To ease the implementation, the immersed wall is treated as an additional immobile identification field φ w (x), which can be regarded experimentally as a free-slip hydrophilic wall boundary. That means, upon contact with the wall, there is no friction and the two normal vectors n w = ∇φ w and n i = ∇φ i tend to remain parallel. The general form of a smooth wall is given by a smooth hyperbolic tangent profile defined as where δ(x) is the field describing the shortest signed distance to the ideal sharp wall interface, with positive values in the fluid region (Luo et al. 2009). The choice for this wall representation is justified. In the simplified onedimensional scenario, Eq. 4 happens to be the equilibrium solution of 4φ w (φ w − 1/2)(φ w − 1) = ξ 2 (∂ 2 φ w /∂x 2 ), which corresponds to a free energy minima (Shinto 2012), hence the coefficient 4 preceding the bulk term in Eq. 2. The advection and Laplacian terms in Eq. 1 are discretized in space using a second-order finite difference and a compact fourth-order scheme, respectively. These two terms are integrated in time using an explicit first-order Euler and an implicit second-order Crank-Nicolson method, resulting in the following algebraic system where φ t i is the i-th identification field known at the beginning of each time step, φ t+δt i the field to be found and δt the time step.
The present model takes ingredients from the Lattice Boltzmann Method originally suggested by Benzi et al. (2009) for foam simulations running on Graphics-Processing-Units. The dynamics of the system are, here too, driven by the minimization of a free energy. The above transport model is implemented using the Portable Extensible Toolkit for Scientific Computation, a widely used open-source library for the scalable solution of partial differential equations (Balay et al. 2019). Parallel simulations with Message-Passing-Interface are presented in the following. With a two-dimensional flow domain discretized into about 40,000 nodes, about 100 gas bubbles and a 16-Central-Processing-Unit core available on most high performance computing centers, up to a month is the typical time needed for a complete dynamic simulation. Despite the excellent performance of the code, the number of gas bubbles (N b ) does set a limit on the computational gain. Solving the transport Eq. (1) N b times at each iteration is simply prohibitive. Compared to other foam simulation tools cited in the introduction, the present method has one major advantage. It is namely the capability to consider bubbly flow, wet and even dry foams coexisting in one single simulation. This allows to investigate in particular the flowing behavior of wet foam near and across the jamming limit.

Set-up
A periodic channel with height H and length L = 1.714H is obstructed by two semi-disks with radii R = 0.193H . A schematic of the channel and its coordinate system is given in Fig. 2a. The channel is filled with N b = 96 bubbles with specific size distributions and pre-defined gas fractions ε.
To achieve the initial equilibrium states, from which the dynamic simulations can start, several workflow steps are required. As illustrated in Fig. 2a, the identification functions φ i are first initialized to squares of equal sizes arranged in a chessboard-like structure with ε ≈ 0.5. At this stage, the exact value ε is not critical. The desired gas fraction is achieved by prescribing the initial areas (Ã i ) 0 to each individual bubble. These are randomly sampled using a standard uniform distribution, whose arithmetic meanĀ results from Equation 1 is then numerically solved with u = 0. Throughout this preliminary relaxation, the bubbles grow, respectively shrink in size. After being squeezed in the vertical channel direction to occupy the empty spaces ( Fig. 2bc), the poly-disperse bubbles eventually become spherical and rearrange into an energetically favorable state (Fig. 2d). The convergence is here achieved when |φ i (t + δt) − φ i (t)| < , with being an arbitrarily small value. In the presented simulations, it is set to = 0.1. A smaller value did not lead to convergence, because the bubble area is never constant but oscillates with a very small amplitude over time. All mono-and poly-disperse equilibrium states are obtained by restarting the state in Fig. 2d with appropriate areasÃ i . Table 1 lists the arithmetic mean of the equivalent diameterd and the corresponding standard deviation σ for each mono-disperse and poly-disperse configuration. Note that the standard deviation does not equal zero in the mono-disperse cases, because the area A i of a bubble never reaches the preset area (Ã i ) 0 perfectly. The equivalent diameter is given in a two-dimensional system by with A i = φ i dr. In the resulting equilibrium foam states illustrated in Fig. 3, different foam structures are well distinguishable. Away from the obstacle, a hexagonal closepacked two-dimensional crystalline structure is observed for the mono-disperse foam. As expected for the hexagonal pattern, each bubble is in contact with exactly six neighbors. With vanishing liquid fraction, Plateau borders eventually form (Forel et al. 2019). In the limiting case ε → 1, the bubbles take polygonal shapes, in-line with previous experimental observations Jones et al. 2013).
With the obtained equilibrium states as initial foam structures, the dynamic simulations are performed. To this end, the foam is advanced in time using a constant horizontal velocity u x . The advective term in Eq. 1 is hence set as The periodicity makes it difficult to drive the foam flow with a realistic pressure gradient. Numerical workarounds with the Surface Evolver have been suggested in the past. See for instance references (Raufaste et al. 2007;Boulogne and Cox 2011), where the foam is transported by small and targeted bubble area increments at each iteration. A constant flow rate could also here be achieved by setting a velocity field u that changes in time and probably space. Enforcing this is however difficult, especially at high gas fractions. We here decide to drive the flow by shifting all bubbles with a small incremental distance (u x δt) over the duration of one time step (δt). This is here achieved by prescribing the constant advection velocity u x in Eq. 8. The velocity u x is our attempt to reproduce the quasi-twodimensional experimental system in references (Dollet et al. , 2015. In such systems, a horizontal mono-layer bed of bubbles is confined between either a column of water and an upper plate or between two plates. Experimentally, the foam velocity is set at the inlet by imposing a constant flow rate. In reality, the Eulerian velocity u should vary in space. That is however not the case in the presented simulations. In the channel minima for example, that is at x = 0, u should increase because of mass conservation and at the wall, it should be zero. A constant u x might rather be reminiscent of the flotation process (Lecrivain et al. 2015), where swarms of gas bubbles rise in a liquid tank with a constant terminal velocity because of buoyancy. Nonetheless, this first work is an alternative and promising method for future foam simulations. Extension of this foam model to account for better flow conditions is currently on-going (Lecrivain et al. 2018), yet it requires additional efforts for gas fractions approaching unity, such as a dynamic mesh refinement algorithm to compute the liquid velocity in the interstitial spaces.

Results
Equation 1 is implemented in its non-dimensional form using the Peclet number given by P e = τ u x / , where u x is the constant advection velocity present in Eq. 8 and the size of one two-dimensional grid element. The Peclet number relates the convective to the diffuse transport. The limiting case P e → 0 corresponds to the scenario, in which the foam reaches an equilibrium in terms of minimum elastic energy at each time step. With a Peclet number set to P e = 0.01 in all presented simulations, the foam flow approaches the quasi-static limit, indicating that the bubble rearrangement is primarily driven by the soft collisions and not by the hydrodynamics, in turn warranting the absence of a Navier-Stokes solver at this stage. Besides, foams in flotation and fundamental experiments are performed well beyond the critical micelle concentration to ensure fully covered interfaces and stable lamellas. Thus, interfacial rheology and Marangoni stress might be neglected in this first attempt. The results are analyzed over one flow-through time, defined as T = L/u x , starting after a transient time period, that is 0.5 < t/T < 1.5. The growth rate and the repulsion parameter are respectively set to α = 1 and β = 1. Unless otherwise specified, the interfacial width equates the size of one grid element, that is ξ = . The time step δt is calculated from the Courant number here set to Co = (u x δt)/ = 10 −4 . The entire channel is discretized into 240×150 nodes in the axial and transverse flow directions, resulting in a grid-independent solution (see Appendix A).

Bubble trajectories and foam velocity
Using a periodicity correction (Bai and Breen 2008), the position of each center of mass is first calculated as Figure 4 shows the bubble trajectories X i (t) for selected simulated cases. At very low gas fraction (ε = 0.44), the bubbles tend to remain in the channel core. With increasing gas fraction, the foam expands vertically to gradually fill the stagnant lower and upper cavity zones between the obstacles. At gas fraction ε ≈ 0.7, the stagnant zone persists only downstream of the obstacles, in-line with previous observations (Deshpande and Barigou 2001). At same gas fraction, the bubble paths in the channel center do not seem to exhibit any ordered patterns. Crossing bubble trajectories are frequent and even more frequent in and further downstream of the channel minima, that is at position x/L = 0 where the vertical distance between the two obstacles is the shortest (See Fig. 4 Sample of bubble trajectories X(t), calculated at the bubble center of mass. The upper and lower rows illustrate the mono-and poly-disperse foams, respectively. ε is the gas fraction Fig. 2a). Such a chaotic behavior is reminiscent of twodimensional gas bubbles rising at low-Reynolds number (Esmaeeli and Tryggvason 1996). With increasing gas fraction (ε ≈ 0.84 and ε = 0.99), stagnant zones vanish and the bubbles eventually travel along relatively well-ordered non-intersecting paths, suggesting a "laminar" foam flow away from the obstacle (van der Net et al. 2009). In the channel minima, crossing trajectories occur irrespective of the gas fraction.
The instantaneous bubble velocity is conveniently calculated from the bubble path as V i (t) = (V x , V y ) i = dX i /dt. Figure 5 shows an instantaneous snapshot of the flowing foam overlaid with the velocity vector placed at each bubble center of mass. An animation can be found in the Supplementary Material. Figure 5 suggests an accelerating foam flow in the channel minima for high gas fractions. For comparison purposes, the time-averaged foam velocity field has been compared to that obtained with the potential flow theory. Given an incompressible and invicid single-phase fluid flowing over a infinite cylinder (Acheson 1990;Fang 2019), the laminar velocity is given in the polar system (e r , e θ ) by where r is the distance from the obstacle center and θ the polar angle. In the potential flow theory, the boundary condition for the velocity is set to u r (R, θ ) = 0. To calculate an average foam velocity at any given point x, the velocity of all bubbles moving in the vicinity of that point, that is |X i − x| < with being an arbitrary distance of some bubble diameters, is averaged in time. Figure 6 shows the comparison between the time-average foam velocity field with that obtained with the potential flow theory. The simulated foam velocity field somewhat agrees with the potential flow theory. The flow directions are similar both in the near-obstacle and the core channel regions. Near the obstacle, the foam flows in the tangential direction e θ , while in the core region, the foam flows in the channel streamwise direction. Differences to the potential flow are indeed excepted and results from the viscoelastic nature of the foam. These have been extensively discussed in previous experimental (Raufaste et al. 2007;Dollet and Graner 2007;Marmottant et al. 2008) and numerical investigations with the Surface Evolver (Cox et al. 2006). Figure 6 also shows, that with increasing gas fraction, the velocity generally decreases in magnitude. This is attributed to the increasing number of soft collisions, leading to kinetic energy dissipation. The model hence reflects well the pressure loss increasing with the overall gas fraction (Deshpande and Barigou 2000). A comparison to experimental data  in the vicinity of the spherical obstacle is provided in Appendix B.
To further assess the effect of the soft collisions, the mean foam velocityV x is plotted against the gas fraction in Fig. 7.
To calculateV x , we average the velocities V x of the bubbles, for which the center of mass is located near the periodic boundary, that is L/2 − λ < ∀X i < L/2 + λ, with λ being an arbitrarily small distance, typically in the order of some grid elements.
Below ε ≤ 0.5, the mean foam velocity equates the prescribed advection velocity u x . This is expected because  Foam velocity vectors (black) overlaid with those obtained with the potential flow theory (orange). The upper and lower rows illustrate the monoand poly-disperse foams, respectively. ε is the gas fraction of the very few collisions. With higher gas fractions, the mean foam velocity decreases linearly asV x /u x = −0.79ε, as illustrated by the straight line fitted with the least-squares method. One may also note that the normalized bubble velocity (V x ) i /u x fluctuates by about ±0.05. This occurs because the post-processed center of mass X(t) does not smoothly evolves with time due to the collisions.
The time-average velocity along the channel center line V c (x) provides further insight into the foam flow. For comparison purposes, we perform additional flow simulations with the open-source flow solver OpenFOAM (Weller et al. 1998), hereafter termed computational fluid dynamics (CFD) simulations. The available power-law model is used (Greenshields 2020) because of its relatively good accuracy in predicting foam flows in pipes (Rooki et al. 2015). This rheological model relates the apparent kinematic viscosity of the foam (ν) to the strain rate (γ ) as ν = kγ n−1 for ν min ≤ ν ≤ ν max , where k is the consistency coefficient and n the power-law index. We set the following properties to the foam, namely k = 0.001 m 2 s n−2 , n = 0.8 (Rooki et al. 2015), ν min = 10 −6 m 2 /s (magnitude order for water) and ν max = 10 −5 m 2 /s Fig. 7 Decrease of the mean foam velocityV x with increasing gas fraction ε. The foam velocity is normalized with the prescribed advection velocity u x . The linear curve is fitted to both mono-and poly-disperse data points (magnitude order for air). Fine changes in the minimum and maximum values of the viscosity did not show any significant differences in the resulting velocity field. On each wall a slip boundary condition is imposed, in-line with the present implementation. Equation 9 along with the numerical solutions obtained with two OpenFOAM solvers, namely "potentialFoam" and "nonNewtonianIcoFoam" are shown in Fig. 8.
At low gas fraction, that is ε < 0.52, the streamwise profile is significantly different to that obtained with the CFD simulations. The velocity V c (x) first decreases until it reaches the position x ≈ −0.25L to then follow a bell-shaped curve peaking downstream of the channel minima at about x ≈ R/L. The results obtained at low gas fraction should however be interpreted with caution, because the hydrodynamics associated with the viscous and inertial effects of the interstitial liquid are here left out. At intermediate gas fraction, that is ε ≥ 0.68, the streamwise velocity qualitatively matches that obtained with the power-law model. The velocity peaks at the channel minima, that is at position x/L = 0. Besides, the peak magnitude generally increases with the gas fraction. An important message emerging from Fig. 8 is that, with the exception of ε ≈ 0.84, the poly-dispersity has a minor effect on the foam flow. Near this intermediate gas fraction, that is ε = 0.84 for the mono-disperse foam, the streamwise velocity curve V c (x) exhibits a relatively wide double-peak across the obstacle length. Such a double peak is not seen in the corresponding poly-disperse scenario.
The bubble rearrangement near the obstacle is complex and can not be precisely analyzed with only time-averaged data. Figure 9 shows for instance a bubble being squeezed in the backward flow direction. The bubble ejection mechanisms are solely attributed to the soft collisions occurring at intermediate gas fraction.

T1 events
The plastic behavior of foams is caused by irreversible topological rearrangements, known as T1 events. According to Cantat et al. (2013), a T1 event occurs when four Plateau borders meet at one node. The configuration becomes unstable and leads to a spontaneous bubble rearrangement reducing the surface energy of the system. We visually searched in the image sequence of the driest configurations (ε = 0.99) for the locations and times, at which T1 events occurred (Dennin 2004). The results are depicted in the Fig. 10 for the mono-and poly-disperse foams. In both cases, the frequency of T1 events increase near the obstacles, with a further increased occurrence in the upstream region. This finding matches those previously reported in previous numerical (Cox et al. 2006) and experimental studies Jing et al. 2015). In terms of T1 events, the effect of poly-dispersity is here noticeable. On the one hand, T1 events observed in the flowing mono-disperse foam occur mostly at constant locations. This is because the same-sized bubbles tend to follow somewhat organized trajectories with only few crossing paths. On the other hand, small bubbles evolving in the poly-disperse foam are more likely to cross trajectories (see Fig. 4). This in-turn introduces sporadicity as well as a wider spacial distribution of the topological changes near the obstacles. A total of 241 and 285 T1 events are here reported for the mono-and poly-disperse scenarios, respectively. This converts to a 18% increase in the number of T1 events occurring as a result of foam poly-dispersity.
The obstacle can be regarded as an artificial wall roughness. Consequently, the bubbles at the wall would remain immobile while those near the channel centerline would move faster. However, as can be seen in Figs. 6 and 7, this did not occur. We had hoped to observe an accumulation of T1 events along some horizontal lines. The spatial distribution of the T1 events does however not point to any shear localization or shear banding (Divoux et al. 2016). Figure 10 shows that the T1 events occur mostly up-and downstream of the obstacle, not supporting the prediction of shear banding. Three-dimensional simulations and an appropriate subgrid wall boundary condition would probably be necessary to further investigate shear banding.

Gas pressure
A strength of the phase-field model is its ability to conveniently calculate a surface tension force F i (x, t) that is non-zero across the diffuse bubble interface and from which a pressure field p i (x, t) is derived. The surface tension force, made here non-dimensional with the interfacial width ξ and the surface tension γ , is calculated as Kim (2007) and Jeong and Kim (2017) where n i = ∇φ i is the vector normal to the smooth bubble boundary, κ i = ∇ · (n i /|n i |) the local curvature and an arbitrarily small value. An advantage of the above formulation is that the field F i is solely function of the identification function φ i . In the quasi-static limit, as is the case here because of the low Peclet number, it is fair to assume that the pressure gradient cancels out the surface tension force as By applying the divergence operator on each side of Eq. 11, that is ∇ 2 p i = ∇ · F i , a pressure field for each bubble is determined. Equation 10 is here validated against the solution of the Young-Laplace equation. For an infinite cylinder, the theoretical pressure difference is given by |p in − p out | = γ /r (Kim 2005; Joshi and Jaiman 2018), where r = d/2 and p in and p out are the pressures inside and outside the bubble, respectively. In Fig. 11, a set of equilibrium solutions are compared to their theoretical counterparts for increasing bubble diameter d. With an interfacial width set to ξ = 2 and a minimum of 15 grid points discretizing the bubble diameter, the results are in excellent agreement with the theory. In all other simulations, the interfacial width is set to ξ = . With such, relatively good results are still achieved. The blue curve is best-fit and is hereafter used as the reference pressure p 0 , so that for an uncompressed disk-like bubble, the normalized pressure p i /p 0 equates unity in the gas and zero outside.
To calculate an instantaneous pressure field p(x, t), we proceed as follows. The algebraic system given by Eq. 11 is first solved using a Newton-like non-linear solver (Balay et al. 2019) for each bubble, that is N b times. All bubbleassociated pressure fields are then summed as p = p i . The resulting overlay is shown in Fig. 12. At low gas fraction, that is here ε = 0.44, the pressure inside each bubble equates the reference pressure p 0 . With increasing gas fraction, that is here ε = 0.68, the pressure equates unity mostly downstream of the obstacle where expansion occurs. In that channel area, the bubbles have few, if any, neighbors and hence remain spherical. Right upstream of the obstacles, the pressure increases by 10 to 20 %. At intermediate gas fraction, that is ε = 0.83, the pressure inside the bubble equates 1.2p 0 and increases up to 1.5p 0 upstream of the obstacle. For the last scenario, that is ε = 0.99, the gas pressure increases even more and its distribution inside each gas bubble is highly heterogeneous. The pressure contours in Fig. 12 do not reveal any noticeable differences between the mono-and poly-disperse scenarios. Even though the simulations approach here the quasi-static limit, the surface tension force in Eq. 10 could be used as source term in the Navier-Stokes momentum equation.

Compression dynamics
Compression refers here to the local change in liquid fraction. The volume of each bubble is conserved by the correction term in Eq. 1. The foam dynamics can be further studied using the bubble circularity, which indicates not only the degree of deviation from a perfect disk, but is also a measure of the local shear stress within the foam . A number of indices have been proposed to measure the circularity of an object (Blott and Pye 2008). We here use the area-based parameter originally defined in two dimensions by Cox (1927) as where P i is the perimeter of the i-th bubble. At low gas fraction and away from the obstacle, the circularity should equate unity because bubbles are circular in shape. The circularity of squeezed bubbles will however vary between zero and unity. Figure 13 shows a snapshot of the foam colored by circularity and Fig. 14 the corresponding timeaveraged field. As expected, the circularity is mostly c = 1 away from the obstacle at low (ε = 0.44 and 0.68) and intermediate gas fractions (ε = 0.83). In these three scenarios, the bubble deformation occurs mostly upstream of the obstacle. At higher fractions (ε = 0.99), the circularity decreases in the entire domain because of the hexagonal shape of the bubbles. A pronounced decrease in circularity is also observed at higher fraction near the obstacle. In this area, the hexagons elongate in the flow direction in-line with previous simulations in obstructed channels with the Surface Evolver (Cox 2015).
An additional measure of foam disorder is the number of neighbors each bubble has. In the spirit of the phase-field model, a natural way to compute the number of neighbors is to look at the value φ i φ j | i =j > 0 from the smooth bubble contours (Mueller et al. 2019). Our tests showed however, that this neighbor detecting method does not work well at intermediate gas fraction because the shortest separation distance between two bubble contours is typically in the order 1 to 2 . For these reasons, the shortest distance was calculated directly from the sharp bubble surface defined Fig. 13 Instantaneous snapshot of foam colored by circularity. The upper and lower rows illustrate the mono-and poly-disperse foams, respectively. ε is the gas fraction  Figure 15 shows results of the neighbor-detecting algorithm. A slight overlap between the isolines only occurs at gas fractions approaching unity, that is ε → 1. In such cases, the overlap did not exceed one grid size . In the following, a bubble contact is counted whenever the shortest distance between two isolines is lower than one grid size element.
As seen in Fig. 15, at low gas fraction (ε = 0.44 and 0.68) the bubbles have no neighbors, which coincides with a gas pressure p = p 0 and a circularity c = 1 worked out earlier in Figs. 12 and 13, respectively. Since the gas fraction is below the jamming point, the bubbles are not in permanent contact with one another. At high gas fraction (ε = 0.99), each bubble of the mono-disperse foam has exactly 6 neighbors in the channel core, corresponding the crystalline structure. Figure 16, which illustrates the timeaveraged number of neighbors along the axial direction x, reveals striking features for each gas fraction.
At low and intermediate gas fractions, that is below ε < 0.8 (blue curves), the foam turns into a compact state upstream of the obstacle to then relax across the channel minima (−R < x/L < R). At high gas fraction, that is ε > 0.9 (green curves), the foam exhibits in the entire channel a highly ordered state with nearly six neighbors throughout, corresponding to the hexagonal packing. However, downstream of the obstacle the packing is disturbed and diverges from its hexagonal structure. Striking, yet different features are observed for the in-between gas fraction (ε = 0.84, orange curve). The foam is in a compact state away from the obstacle and relaxes to some extent as it moves Fig. 15 Neighbor-detection. Each blue marker indicates that the shortest separation distance falls below the grid size . The integer inside each bubble refers to the number of neighbors N. The upper and lower rows illustrate the mono-and poly-disperse foams, respectively. The term ε is the gas fraction through the channel minima. Foam poly-dispersity only a has minor effect in terms of average neighbor numbers, when compared to the reference mono-disperse case.

Conclusions
A modified phase-field method has been successfully introduced for the simulation of flowing mono-and polydisperse foams. In this method each bubble is described with a smooth transitioning identification function φ i . This enables investigation of a wide variety of scenarios, ranging from the strong deformation of highly packed bubbles down to loosely packed and even individual bubbles. The transition between both states, such as the formation of dry foam from rising bubbles, could also well be covered with the present model. This is a rather difficult task to overcome with other currently available numerical methods. The model was applied to investigate the compression and relaxation dynamics of foam as it passes through a channel with confinement. Many of the properties observed experimentally were reproduced. Foam states ranging from low (ε = 0.44), intermediate (ε = 0.84) to high gas fractions (ε = 0.99) matched the previously reported kinetics of two-dimensional foam flowing around a cylindrical obstacle. The model does however have some limitations. A disadvantage of the model is, that surface tension, interfacial rheology (for instance sorption of surface active agents) and interaction forces between neighboring bubbles do not appear explicitly in the formulation. Instead, they have to be accounted for by the repulsive parameter β in Eq. 3.
The next step will include the coupling with a Navier-Stokes solver. In the present case, the lamellas between the bubbles require special treatment, because the film lamellas in realistic foam are too thin to be spatially resolved. Alternatively, a subgrid model for the tangential stress between sliding bubbles could be incorporated. Also foam coarsening by gas diffusion could be implemented by adding a source term to account for the gas transfer based across the overlapping area formed by two adjacent identification functions φ i and φ j and by changing the growth rates α i and α j accordingly. Future work will also include an extension to three-dimensional simulations. This can be easily achieved with the phase-field model, because in contrast to other methods, e.g. Volume-of-Fluid methods, the transport equation for the scalar identification functions do not differ substantially between two and three dimensions. With a fully applicable phase-field model, the rheology of wet foam near and across the jamming point could well be approached to extend the experimental findings of other authors (Herzhaft et al. 2005;Heitkam and Fröhlich 2019). Also, the reversibility and the relaxation dynamics of foam deformation close to the jamming point could be further investigated and compared to space experiments (Langevin 2017).
The code implementation together with selected raw data have been made available (Lecrivain 2020).

Appendix A: Grid dependency test
A grid dependency test is here presented. Five Cartesian grids discretizing the channel into N x = 112, 144, 176, 208, 240 and N y = 70, 90, 110, 130, 150 nodes are considered, where N x and N y are the number of nodes in the axial and transverse directions of the channel flow. In all simulations, the ratio of the sharp bubble diameter to the channel height is kept constant and corresponds to the mono-disperse scenario ε = 0.68 illustrated in Fig. 5, that isd/H = 0.12. The interfacial width is set to the size of one grid element, that is ξ = and is the only non-constant value. The time-average velocity along the channel center line V c (x) is determined for the fives grids. As seen in Fig. 17, all simulations deliver similar velocity profiles. In this work, simulation data obtained with the finest grid, that is N x = 150 and N y = 240, are discussed.

Appendix B: Comparison to reference literature data
We compare the time-averaged velocity field of the driest mono-disperse foam to experimental data. In the original experiment , a mono-layer bed of bubbles is confined between a lower liquid pool and an upper transparent wall. A cylinder is placed in the channel center and the bubbles flow around it. The experimental scenario is somewhat different to that presently simulated. Here, two half-disks are located at the wall. In the numerical Fig. 18 Foam velocity vectors (black) overlaid with those obtained experimentally (orange). The foams are mono-disperse and dry in both the simulation and the reference experiment . The ratio of the equivalent bubble to obstacle diameter equate d/(2R) = 0.37 and 0.15, respectively. ε is the gas fraction and experimental scenarios, the ratios of the equivalent bubble to obstacle diameter have the same magnitude order and equate d/(2R) = 0.37 and 0.15, respectively. Figure 18 compares the numerical and experimental velocity vectors, that have been digitized from the publication. Despite the differences between the two systems in terms of obstacle position, the magnitude and the direction of the foam velocity qualitatively match. The discrepancies tend to increase as one approaches the obstacle.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.