Numerical Methods for Solving Nonlinearly Coupled Poisson Equations in Dual-Continuum Modeled Porous Electrodes

Porous electrodes are widely used in electrochemical systems, where accurately determining electric potentials, particularly overpotentials, is essential for understanding electrode behavior. At the macroscopic scale, porous electrodes are typically modeled using a dual-continuum approach, treating the porous solid phase and the liquid electrolyte as spatially superimposed domains. Determining potential distributions requires solving two Poisson equations that are nonlinearly coupled through Butler-Volmer kinetics under galvanostatic and potentiostatic operating modes. Under galvanostatic operation, these equations form an underconstrained singular system due to all-Neumann boundary conditions, posing numerical challenges. This paper systematically presents numerical methods for solving nonlinearly coupled Poisson equations in dual-continuum porous electrodes, with a particular focus on galvanostatic solutions. We mathematically establish solution uniqueness in terms of the potential difference between the electrode and electrolyte (or overpotential), as well as the individual potentials up to a shared constant shift. To resolve the nonuniqueness of the solution, we introduce three numerical approaches: (1) Lagrange Constrained Method (LCM), (2) Dirichlet Substitution Method (DSM), and (3) Global Constraining Method (GCM), where GCM enables solving the overpotential without imposing an explicit system reference potential. Additionally, we develop both decoupled and fully coupled nonlinear solution strategies and evaluate their computational performance in both homogeneous and heterogeneous conductivity cases. The presented numerical methods are general for addressing similar underconstrained nonlinear systems. A Python implementation is provided at https://github.com/harrywang1129/porous_electrode_solver.