PaperlayerView exact revision on arXiv

Source-linked AI summary

A dynamic pore-network model with discrete microbubble transport

Kristoffer Skjelanger, Anna Mareike Kostelecky, Maziar Veyskarami, Timo Koch

arXiv:2610.09888v1cs.CEphysics.flu-dyn

TL;DR

The paper asks how apparent gravitational forces shift gas transport between connected pathways and discrete bubbles in a PEM-electrolyser porous transport layer. It couples a dynamic pore-network model to a bubble tracker, using pore-scale Stokes simulations for closures. In network simulations, the bubble-carried gas fraction increases with Γ, reaching approximately 40% of produced gas for Γ > 100.

  • Problem

    It remains unclear how strongly apparent gravitational forces shift gas transport between connected-phase pathways and discrete bubbles in a PEM-electrolyser porous transport layer.

    Source: §1 · p. 1

  • Method

    The paper couples a dynamic two-phase pore-network model to a discrete-bubble tracker and derives bubble-motion and throat-impact closures from Stokes simulations.

    Source: §1 · p. 2, §4 · p. 7

  • Results

    Approximately 40% of produced gas leaves the network as bubbles when Γ > 100, and the bubble-carried fraction increases as Γ increases.

    Source: §8 · p. 23

  • Takeaways & Limitations

    In the simulated PEM-electrolyser PTL networks, increasing Γ shifts more produced gas toward discrete-bubble transport, largely as the detachment radius decreases and bubbles move more freely.

    Source: §8 · p. 23

    The detachment process is simplified: growth location is unresolved and the effective acceleration is assumed constant and perpendicular to the growth surface, making detachment radius constant throughout the network.

    Source: §7.8 · p. 22

Abstract

from arXiv · show

We develop a hybrid pore-network model (PNM) for discrete microbubble transport coupled to connected two-phase flow in porous media. We apply the model to simulate flow in the porous transport layer of a rotating proton exchange membrane (PEM) electrolyser. At sufficiently large apparent gravitational forces, bubbles growing at the catalyst layer can detach at sizes smaller than the characteristic pore size and move through liquid-filled pores via buoyancy and advection. The model couples a dynamic two-phase PNM to a Lagrangian bubble tracker. Pore-scale Stokes simulations in square throats provide closures for the bubble buoyancy velocity, the hydraulic resistance of bubble-occupied throats, and the reciprocal coupling between liquid flow and bubble motion. Verification tests assess conservation, pressure-flow coupling, bubble growth, advection, merging, trapping, and transfer to the connected gas phase. The network simulations show that increasing apparent gravitational forces shifts a larger fraction of the produced gas toward discrete-bubble transport. The increase mostly comes from a reduced detachment radius, which decreases bubble size and lets them move more freely through the network.

1. Introduction

The paper addresses how apparent gravitational forces shift gas transport between connected pathways and discrete bubbles in PEM-electrolyser porous transport layers. It develops a hybrid dynamic pore-network model to track both representations and their interactions.

1. Introduction

  • The balance between connected-phase and discrete-bubble transport under strong apparent gravity remains unclear for the PEM-electrolyser porous transport layer.Connected gas pathways commonly form by capillary invasion under stationary conditions, while experiments report improved performance under centrifugation, attributed presumably to faster detachment and gas removal.Source: §1 · p. 1
  • Interface-resolving simulations capture individual interfaces but are computationally infeasible for network-scale simulations involving many bubbles.Pore-network models offer a lower-cost reduced description, assigning storage and primary variables to pore bodies and fluxes to throats.Source: §1 · p. 1
  • Prior network models address individual bubble mechanisms but do not couple discrete tracking to dynamic two-phase flow with conservative gas transfer and liquid-flow feedback.The identified gap is the combination of a discrete-bubble tracker with a dynamic two-phase pore-network model.Source: §1 · p. 2
  • The study derives bubble-motion rules from Stokes simulations in square throats, assuming small capillary number and Bond number slightly below unity.The model treats bubble shape as determined by volume and throat geometry, and derives local rules for buoyancy and advection.Source: §1 · p. 2
  • The hybrid model couples a dynamic two-phase pore-network model to a discrete bubble tracker and uses three Stokes-derived closure functions.Discrete bubbles affect pore volume, throat conductance, liquid pressure-flow coupling, and exchange of mass with connected gas; the formulation is described as conservative.Source: §1 · p. 2

2. Problem setup and assumptions

The section sets assumptions for a rotating electrolyser model and motivates discrete-bubble transport when apparent gravity is strong. It establishes bubble-motion and shape simplifications used for the pore-throat closures.

2. Problem setup and assumptions

  • The rotating-electrolyser setup is treated as isothermal, with fluid properties fixed to representative PEM titanium-PTL values at 60°C.The material-property values are used throughout the paper for all developed models.Source: §2 · p. 2, Table 1 · p. 3

2.1. Apparent gravitational forces

  • The model treats centrifugal acceleration as uniform across the PTL and neglects Coriolis effects on bubble closures and network-scale motion.The effective acceleration is expressed through dimensionless G-level Γ and assumed normal to the catalyst surface.Source: §2.1 · pp. 2–3

2.2. A buoyancy-driven bubble flow regime

  • Strong apparent gravity can detach bubbles smaller than typical pores, allowing buoyant motion through liquid and entry into smaller throats.The detachment model balances acceleration-induced buoyancy against surface tension; with a representative contact radius of 0.1 μm, predicted detachment radii are 10–100 μm.Source: §2.2 · pp. 3–4

2.3. Bubble motion and deformation

  • The assumed regime supports fixed-shape, throat-local closures, while the study derives them from Stokes equations despite possible inertial effects at high Reynolds number.At Γ = 300, the estimated Bond number is approximately 2 ⋅ 10^-1 and is not large enough for gravity to dominate bubble shape; inertia-induced deformation is neglected.Source: §2.3 · pp. 4–5

2.4. Bubble shapes in square capillaries

  • Square-throat bubbles are spherical below the prescribed wall-clearance limit and elongated above it, with a uniform dimensionless clearance of 0.01.The closure coefficients are calibrated for this clearance, and the study does not establish independence from it; resistance near the shape transition may be sensitive to the parameter.Source: §2.4 · p. 5

2.5. Gas-liquid interface conditions for bubbles

  • The model assumes immobile, no-slip bubble interfaces to represent contamination expected in an electrolyser PTL.The resulting closures apply to strongly immobilised interfaces rather than clean, mobile bubble interfaces.Source: §2.5 · p. 5

3. Hybrid pore-network model with discrete bubble tracking

The hybrid model combines a dynamic two-phase pore-network description with a discrete-bubble tracker. Its balances account for bubble motion and for the effects of discrete bubbles on connected-phase flow and mass transfer.

3. Hybrid pore-network model with discrete bubble tracking

  • The hybrid model combines an isothermal, non-compositional two-phase pore-network model with a discrete oxygen-bubble tracker.Its primary variables include wetting-phase pressure, non-wetting saturation, bubble position, and bubble volume.Source: §3 · p. 5
  • Bubble velocity is decomposed into buoyancy and advection components, and each bubble has its own mass and momentum balance.The balance formulation includes bubble mass sources and throat-scale motion.Source: §3 · p. 6
  • The model couples bubbles to continuum flow through pore-volume occupancy, liquid-conductance changes, buoyancy-induced pressure, and liquid displacement between pore bodies.When a bubble transfers to connected gas, its mass is conserved and continuum saturations are updated accordingly.Source: §3 · p. 6
  • Throat conductances distinguish fully water-saturated flow, wetting-phase corner flow after invasion, and non-wetting gas flow in invaded throats.The conductance laws draw on square-throat models, including the cited proposals of Patzek and Silin, Ransohoff and Radke, Zhou et al., Bakke and Øren, and Tørå et al..Source: §3 · pp. 6–7

4. Bubbles in pore throats: closure relations from Stokes simulations

Pore-scale Stokes simulations provide three geometry-dependent closures for bubble motion, bubble-occupied throat conductance, and reciprocal pressure–advection coupling. The fitted closures remain within 5% relative error of the simulation results.

4. Bubbles in pore throats: closure relations from Stokes simulations

  • The three fitted closures lie within 5% relative error of the Stokes simulation results.Together they describe buoyancy velocity, force-free conductance in bubble-occupied throats, and reciprocal pressure–advection coupling.Source: §4 · p. 7

4.1. Simulation setup

  • Stokes simulations use a stationary-bubble frame to obtain buoyancy-driven and pressure-driven flow closures under distinct force and flux constraints.The buoyancy case sets net liquid flux to zero; the advection case makes the bubble force-free.Source: §4.1 · p. 7

4.2. Buoyancy velocity model

  • The buoyancy-velocity fit combines a square-channel wall correction with additional resistance from corner and film flow as bubble length increases.The fitted model is shown with relative errors in Figure 3.Source: §4.2 · p. 8, Figure 3 · p. 9

4.3. Throat conductance model

  • Bubble effects on throat conductance are represented as an excess length, with excess lengths for multiple bubbles assumed to add in series.The fitted excess-length model is shown in Figure 4.Source: §4.3 · pp. 8–9, Figure 4 · p. 9

4.4. Advection velocity and viscous pressure model

  • A reciprocal coupling function links buoyancy-induced pressure drop and advection-driven bubble velocity, allowing both to be obtained from one fitted coefficient.The closure fit and residuals are shown in Figure 5.Source: §4.4 · pp. 9–10, Figure 5 · p. 11

5. Numerical algorithm

The hybrid time step couples an implicit pore-network solve with bubble tracking, iterating until bubble positions, network state, events, and invasions are consistent. Bubbles are injected and grown at source pores, then move through throats under buoyancy and advection, with rules for entry, merging, and transfer to connected gas.

5. Numerical algorithm

  • The coupled time step iterates the network solve and bubble tracker, updating events and invasions until the coupled state is consistent.The inner iteration holds events and invasions fixed; the outer iteration updates them, and the time step is reduced if convergence fails.Source: §5 · p. 12

5.1. Spawning and bubble growth

  • At source pores, gas starts or enlarges a discrete bubble until detachment, unless connected gas is present, in which case the source feeds that phase directly.Bubble gas is assumed to have fixed density.Source: §5.1 · p. 12

5.2. Bubble movement

  • Bubbles move through one throat at a time under buoyancy and liquid advection, selecting the next throat by minimum resistance.If a bubble exceeds the throat’s inscribed radius, it enters only when pressure forcing exceeds the capillary-pressure increase associated with deformation.Source: §5.2 · pp. 13–14
  • Bubbles merge on contact or overtaking, while a bubble too large for its pore or entering a continuously gas-invaded region transfers irreversibly to connected gas.Transfer also applies when the entry condition fails; the resulting non-wetting saturation can invade neighboring pores.Source: §5.2 · p. 14

6. Model verification

The model is verified with analytical checks of conservation and coupled pressure-flow behavior, plus tests of discrete bubble events. The setup represents oxygen generation at the bottom of a water-saturated network and gas removal at the top.

6.1. Initial and boundary conditions

  • The network starts fully water-saturated, with oxygen generated at bottom source pores and removed at top outlets; side boundaries impose no flow.The bottom boundary can represent an open system with prescribed pressure difference or a closed reactive boundary with no phase flux.Source: §6.1 · pp. 14–15, §6 · p. 14

6.2. Verification

  • Tests 1–4 check bubble-induced flux, storage, growth, advection, and pressure coupling against analytical relations; reported fluxes and pressure effects match to machine precision.The tests include a moving bubble in a chain, growth and detachment, advection, and a buoyancy-driven bubble in a loop network.Source: §6.2 · pp. 15–16
  • The loop-network test reproduces analytical pressure-drop magnitudes and bubble-velocity changes to machine precision under the closure laws.The test probes reciprocal pressure-flow coupling as a buoyancy-driven bubble moves through a bypassed throat.Source: §6.2 · p. 16
  • Tests 5–8 verify merging, continuous creation, trapping, and throat narrowing, reporting conserved bubble volumes and, where checked, velocities matching expected solutions to machine precision.The trapping test transfers bubbles to the non-wetting phase while conserving mass.Source: §6.2 · p. 16

7. Gas transport at high apparent gravitational forces

Across the network simulations, increasing apparent gravity raises the fraction of gas transported as bubbles, while pore geometry and trapping shape how much reaches the outlet. The section also reports time-step checks and modeling assumptions that constrain these predictions.

7.1. Simulation setup

  • The 5 × 5 × 5 network generates gas at a constant rate in bottom-boundary source pores and is run to a defined steady state.Its height is comparable to that of a real electrolyser PTL; steady state requires no recent invasion events and outlet flow within 2% of the source rate.Source: §7.1 · pp. 16–17

7.2. Pore network geometry

  • Networks are built from a structured pore grid with randomly deleted throats and pore radii sampled from a log-normal distribution.Throat radii are calculated from adjacent pore radii and their separation using an averaging relation.Source: §7.2 · pp. 17–18

7.3. Time-step independence

  • For buoyancy-dominated transport, the estimated maximum time step keeps a typical bubble from crossing more than half a throat, and refinement beyond the estimate changed exit ratio by less than 1%.The estimate may need adjustment when advection dominates; the reported test found little impact from further refinement.Source: §7.3 · p. 18

7.4. Model predictions for varying Γ

  • As Γ rises and r_d/r̄_p falls, bubble transport increases; smaller bubbles travel farther before trapping, while invasion dominates as r_d/r̄_p approaches 1.The invasion patterns also show fewer invaded pores at higher Γ, with invasion reaching higher pores.Source: §7.4 · p. 18, Figure 9 · p. 19, Figure 10 · p. 20

7.5. Model predictions for varying mean pore radiī 𝑟𝑝

  • At fixed Γ, varying mean pore radius qualitatively reproduces the trend from varying Γ, but the difference indicates that network heterogeneity also influences bubble transport.Bubble movement depends on whether pores accommodate the bubbles and whether the throat-entry condition is met.Source: Figure 11 · p. 20, §7.5 · p. 20

7.6. Model predictions for varying Γ with fixed 𝑟𝑑∕̄𝑟𝑝

  • With detachment radius fixed at 50 μm, the bubble-phase fraction increases with Γ through the throat-entry condition, while the exit rate is barely affected.This enlarged detachment radius highlights cases where bubbles are close to throat size and entry conditions determine mobility.Source: Figure 12 · p. 21, §7.6 · p. 21

7.7. Larger pore networks and simulation times

  • In 10 × 10 × 10 networks, inlet bubble-transport ratios resemble the smaller networks, but outlet ratios are lower because more bubbles can become trapped.The larger-network invasion patterns otherwise resemble the smaller-network results, with higher Γ associated with invasion farther up the network and increased bubble transport.Source: §7.7 · p. 21, Figure 13 · p. 22

7.8. Model limitations

  • The model fixes gas injection and detachment assumptions, transfers stuck bubbles to connected gas, and uses Stokes closures limited to small Reynolds numbers.The transition to connected gas is unidirectional, and closure assumptions include uniform minimum film thickness, negligible Coriolis effects, and negligible dynamic deformation.Source: §7.8 · pp. 22–23

8. Conclusion

The paper introduces a hybrid pore-network model that couples dynamic two-phase flow with discrete bubble tracking and verifies its component behaviors. It predicts greater bubble-carried gas at higher effective gravity, while dynamic bubble shapes affected by flow remain outside its scope.

8. Conclusion

  • The hybrid model couples a dynamic two-phase pore-network model with a discrete bubble tracker, using Stokes-based closures for bubble movement and impact in square throats.Verification tests isolate model features; the model is intended for regimes where flow has negligible impact on bubble shape.Source: §8 · p. 23
  • Flow-driven dynamic bubble shapes are beyond the study’s scope but matter in spatially larger systems.Source: §8 · p. 23
  • At Γ > 100, approximately 40% of produced gas leaves the network as bubbles, with bubble transport increasing as effective gravitational strength rises.The increase is associated with a reduced detachment radius and with stronger gravity forcing bubbles into throats they otherwise could not enter.Source: §8 · p. 23

Data and code availability

The paper identifies a repository for its source code, scripts, and simulation data. Its appendix describes how square-channel bubble geometries and meshes are prepared for Stokes simulations.

Data and code availability

  • Source code, reproduction scripts, and simulation CSV data are available in the same repository, with additional output available from the authors upon reasonable request.Source: Data and code availability · p. 24

A. Bubble rising in a square pore throat

  • The Stokes-closure dataset uses precomputed bubble shapes and Stokes solutions in a narrow square channel, with auxiliary calculations for small-bubble and long-bubble limits.The simulations are reported as mesh-converged to a tolerance considered acceptable for the application.Source: A. Bubble rising in a square pore throat · p. 27

A.1. Bubble shape and fluid domain

  • Bubble surfaces are generated with Surface Evolver and subtracted from a square channel to create tetrahedral fluid meshes, using a minimum film thickness of 0.01 r_ij.Only one quarter of the domain is meshed and simulated because the problem is axisymmetric along the throat.Source: §A.1 · p. 27

A.2. Stokes system and boundary conditions

The appendix specifies Stokes-flow boundary conditions for the bubble-channel problem and numerical procedures for finding terminal bubble velocity and viscous pressure drop. The velocity calculation balances vertical force; the pressure-drop calculation also enforces zero wetting-phase flux.

A.2. Stokes system and boundary conditions

  • The Stokes system uses a no-slip bubble surface, prescribed inlet and outlet tractions, and fixed outer walls that move oppositely to the bubble in its co-moving frame.Source: §A.2 · pp. 27–29

A.3. Calculating the bubble velocity

  • Terminal bubble velocity is found by iteratively adjusting wall velocity until the vertical force residual reaches zero.The iteration accounts separately for the gas bubble’s weight because the gas phase is not explicitly solved.Source: §A.3 · p. 29

A.4. Calculating the viscous pressure drop

  • Viscous pressure drop is calculated by jointly adjusting inlet pressure and wall velocity until both vertical force and wetting-phase flux are zero.A two-variable Jacobian is used because hydrostatic boundary pressure alone produces a spurious flux in the bubble co-moving frame.Source: §A.4 · pp. 29–30

A.5. Mesh independence

The Stokes-flow closure results become more accurate as the mesh is refined, and changing throat length has little effect on the calculated quantities. The study also checks closure accuracy against an analytical limit.

A.5. Mesh independence

  • Finer meshes reduce RMS relative errors in buoyancy velocity, relative conductance, and closure 𝜙 against Richardson-extrapolated limits.For the finest mesh, the error in 𝜙 relative to the analytical limit is approximately 1%.Source: §A.5 · p. 30, Figure 18 · p. 30
  • Changing throat length from 100 to 300 μm did not substantially affect the calculated quantities.The comparison fixes throat radius at 25 μm and uses a mesh cell size of 0.8 μm.Source: §A.5 · p. 30, Figure 19 · p. 31

B. Corner flow factor 𝑘

The corner-flow factor is obtained from a geometric Laplace problem, while the supplementary passages describe the flow-resistance formulation, bubble-flow simulations, distributions, and verification material.

B. Corner flow factor 𝑘

  • For the long-bubble asymptote, the geometric corner-flow factor is k = 0.428, and 𝜙∞ depends only on geometry, not Stokes velocity or pressure data.The reported k range is 0.427–0.429 across a 3–12 cell refinement across the film.Source: B. Corner flow factor 𝑘 · p. 31

C. Flow resistance 𝛽

  • The wetting-phase resistance formulation allows an intermediate wetting–non-wetting boundary condition by setting f₂ between 0 and 1.The passage identifies f₂ = 1 as no-slip and f₂ = 0 as free-slip.Source: C. Flow resistance 𝛽 · p. 32

S1. Bubble velocity decomposition

  • Figure S1 illustrates the additive buoyancy and advection velocity components when Δp = −10 Pa and Γ = 10 drive motion in the same direction.Source: Figure S1 · p. 33

S2. Stokes flow around a bubble

S3. Bubble distributions

S4. Details of the model verification tests

Test 1 - local flux conservation

The verification tests check local flux conservation during bubble transfer and the modeled growth and detachment response to a gas source. Bubble growth tracks the imposed source rate until detachment.

Test 1 - local flux conservation

  • When a bubble crosses a pore throat midpoint, its volume transfers between pore bodies while water flows proportionally in the reverse direction to fill the vacated space.The test uses a fixed-radius 50 μm bubble moving at 0.05 m/s.Source: Test 1 - local flux conservation · p. 36, Figure S7 · p. 36

Test 2 - bubble growth and detachment

Test 3 - bubble advection

The verification cases show how bubbles alter network flow, move, merge, become trapped, and transfer into the connected gas phase. Supplementary studies also report time-step behavior and the computational cost of larger networks.

Test 3 - bubble advection

Test 4 - viscous pressure drop

  • In the looped network, a bubble’s viscous pressure drop induces wetting flux in the loop and contributes to bubble velocity.The viscous pressure drop is applied across the rest of the network when the bubble moves in the bottom or top throat.Source: Test 4 - viscous pressure drop · p. 38, Figure S10 · p. 38

Test 5 - bubble merging / Test 6 - continuous bubble creation

Test 7 - bubble trapping

Test 8 - throat narrowing

S5. Hybrid PNM invasion patterns

  • For Γ = 100 and the stated network parameters, invasion patterns are consistent when Δt_max remains below the estimated 1.2 × 10^-4 s threshold.Source: Figure S15 · p. 41, Figure S16 · p. 42

S6. Solver statistics for the time-step study

  • As Δt_max decreases, runtime increases; the fraction of accepted two-cycle states increases with Δt_max.The solver statistics cover the 10 s simulations and include gas-exit rates and bubble-fate ratios.Source: Figure S18 · p. 43, Figure S17 · p. 43

S7. Steady-state runtime

  • The 10 × 10 × 10 network costs approximately 30 times more computationally per simulated second than the 5 × 5 × 5 network; the implementation was not performance-tuned.Source: S7. Steady-state runtime · p. 43
Loading 2610.09888v1…