A finite element solver for a thermodynamically consistent electrolyte model

TL;DR

This work develops a finite element solver based on a thermodynamically consistent electrolyte model that captures multicomponent ion transport with steric, solvation, and pressure coupling. Grounded in non-equilibrium thermodynamics, the model enforces mass conservation, charge neutrality, and entropy production, and it extends beyond classical Nernst--Planck frameworks by incorporating a variational formulation and a momentum balance coupled to electrostatics. Implemented in FEniCSx, the solver demonstrates robust convergence in 1D and 2D settings, validates against benchmark problems, and reveals key electrolyte phenomena such as electric double layer formation, ion saturation due to solvation, and the impact of Debye length and compressibility on transport. The framework, validated against benchmarks and compared to NP models, provides a physically faithful tool for simulating complex electrochemical systems and is released with reproducibility materials for broader adoption.

Abstract

In this study, we present a finite element solver for a thermodynamically consistent electrolyte model that accurately captures multicomponent ionic transport by incorporating key physical phenomena such as steric effects, solvation, and pressure coupling. The model is rooted in the principles of non-equilibrium thermodynamics and strictly enforces mass conservation, charge neutrality, and entropy production. It extends beyond classical frameworks like the Nernst-Planck system by employing modified partial mass balances, the electrostatic Poisson equation, and a momentum balance expressed in terms of electrostatic potential, atomic fractions, and pressure, thereby enhancing numerical stability and physical consistency. Implemented using the FEniCSx platform, the solver efficiently handles one- and two-dimensional problems with varied boundary conditions and demonstrates excellent convergence behavior and robustness. Validation against benchmark problems confirms its improved physical fidelity, particularly in regimes characterized by high ionic concentrations and strong electrochemical gradients. Simulation results reveal critical electrolyte phenomena, including electric double layer formation, rectification behavior, and the effects of solvation number, Debye length, and compressibility. The solver's modular variational formulation facilitates its extension to complex electrochemical systems involving multiple ionic species with asymmetric valences. We publicly provide the documented and validated solver framework.

Figures (15)

• Figure 1: Schematic of the electrical double layer (EDL) near a positively charged electrode. Negatively charged ions accumulate near the electrode surface, forming the compact layer, while the bulk electrolyte remains electrically neutral with a uniform distribution of positive and negative ions.
• Figure 2: Convergence behavior of the finite element solver for a one-dimensional incompressible ternary electrolyte system. The plot shows the $e_{2}-$ and $e_{\infty}-$error norms as functions of mesh resolution (number of elements, $n_x$) on a logarithmic scale.
• Figure 3: Validation study comparing the thermodynamically consistent model with the classical NP model in a one-dimensional ternary electrolyte system. (a) Number density profiles of anions ($n_A$) and cations ($n_C$) for three model variants: ideal mixture ($\kappa = 0$), solvated ions ($\kappa = 8$), and the classical NP model in Poisson--Boltzmann formulation. (b) Electrostatic potential ($\varphi$) and pressure ($p$) distributions for the same configurations. Dots indicate the reference solution from dreyer2018bulk. The results demonstrate close agreement with the reference and highlight the divergence behavior in the NP model, emphasizing the improved physical consistency of the extended model.
• Figure 4: Comparison between (a) classical NP model, and (b) DGM model: solution profiles of anion atomic fraction at different potential differences.
• Figure 5: Comparison between DGM and the classical NP model: profiles of electric potential ($\varphi$), pressure ($p$) and atomic fraction ($y_{\alpha}$) at different electric potentials (a) $\delta \varphi=1$, and (b) $\delta \varphi = 10$.