A higher order pressure-stabilized virtual element formulation for the Stokes-Poisson-Boltzmann equations

Electrokinetic phenomena in nanopore sensors and microfluidic devices require accurate simulation of coupled fluid-electrostatic interactions in geometrically complex domains with irregular boundaries and adaptive mesh refinement. We develop an equal-order virtual element method for the Stokes–Poisson–Boltzmann equations that naturally handles general polygonal meshes, including meshes with hanging nodes, without requiring special treatment or remeshing. The key innovation is a residual-based pressure stabilization scheme derived by reformulating the Laplacian drag force in the momentum equation as a weighted advection term involving the nonlinear Poisson–Boltzmann equation, thereby eliminating second-order derivative terms while maintaining theoretical rigor. Well-posedness of the coupled stabilized problem is established using the Banach and Brouwer fixed-point theorems under sufficiently small data assumptions, and optimal a priori error estimates are derived in the energy norm with convergence rates of order $\mathcal{O}(h^k)$ for approximation degree $k \geq 1$. Numerical experiments on diverse polygonal meshes – including distorted elements, non-convex polygons, Voronoi tessellations, and configurations with hanging nodes – confirm optimal convergence rates, validating theoretical predictions. Applications to electro-osmotic flows in nanopore sensors with complex obstacle geometries illustrate the method’s practical utility for engineering simulations. Compared to Taylor–Hood finite element formulations, the equal-order approach simplifies implementation through uniform polynomial treatment of all fields and offers native support for general polygonal elements.

💡 Research Summary

The paper addresses the computational challenge of simulating coupled fluid–electrostatic interactions described by the Stokes–Poisson–Boltzmann (SPB) system in geometrically complex domains such as nanopore sensors and micro‑fluidic devices. Traditional mixed finite‑element approaches require different polynomial orders for velocity and pressure (e.g., Taylor‑Hood elements) and often struggle with non‑convex, highly irregular meshes or hanging nodes, leading to cumbersome mesh generation and stabilization procedures.

To overcome these limitations, the authors develop a higher‑order equal‑order virtual element method (VEM) that works naturally on general polygonal meshes, including those with hanging nodes, without any special treatment. The central innovation is a residual‑based pressure‑stabilization scheme that rewrites the Laplacian drag term in the momentum equation as a weighted advection term involving the nonlinear Poisson–Boltzmann potential. By introducing a weight function β(φ) that depends on the electrostatic potential φ, the second‑order derivative is eliminated, and the resulting term ∇·(β(φ) u) can be discretized using only the boundary degrees of freedom inherent to VEM. This reformulation preserves the physical coupling between fluid viscosity and electrostatic forces while enabling a uniform polynomial approximation of degree k (k ≥ 1) for velocity, pressure, and potential.

Mathematical analysis proceeds in two stages. First, the linearized coupled system is shown to satisfy an inf‑sup condition thanks to the residual‑based stabilization, guaranteeing pressure stability. Second, for the full nonlinear problem, the authors apply the Banach fixed‑point theorem to prove that the nonlinear operator is a contraction under sufficiently small data (body forces, boundary data). Existence and uniqueness of the discrete solution are then secured via Brouwer’s fixed‑point theorem on the finite‑dimensional VEM space.

A priori error estimates are derived in the energy norm. Using standard VEM interpolation bounds and the stability of the pressure term, the authors obtain optimal convergence rates of order O(h^k) for all fields when the mesh size h tends to zero and the polynomial degree is k. The analysis holds for any shape‑regular polygonal mesh, regardless of convexity or the presence of hanging nodes.

Numerical experiments validate the theory on four representative mesh families: (i) highly distorted quadrilaterals, (ii) non‑convex polygons (e.g., L‑shaped and star‑shaped elements), (iii) Voronoi tessellations producing irregular polygons, and (iv) adaptive meshes with hanging nodes. In each case, L² and H¹ errors for velocity, pressure, and electrostatic potential exhibit the predicted O(h^k) decay for k = 1, 2, 3. The pressure field remains free of spurious oscillations, confirming the effectiveness of the residual‑based stabilization.

Finally, the method is applied to a realistic electro‑osmotic flow problem in a nanopore sensor featuring complex obstacle geometries. The equal‑order VEM accurately captures the electric double‑layer distribution and the resulting fluid motion without requiring mesh remeshing or special treatment of the irregular boundaries. Compared with a conventional Taylor‑Hood finite‑element implementation, the proposed VEM offers three decisive advantages: (1) a single polynomial order for all unknowns simplifies code development, (2) native support for arbitrary polygonal elements eliminates the need for mesh conformity, and (3) the pressure‑stabilization eliminates the inf‑sup restriction, allowing equal‑order approximations while preserving optimal accuracy.

In summary, the paper delivers a mathematically rigorous, high‑order, pressure‑stabilized virtual element formulation for the Stokes–Poisson–Boltzmann equations. It combines theoretical robustness (well‑posedness, optimal error bounds) with practical flexibility (general polygonal meshes, hanging nodes) and demonstrates its utility in engineering‑scale electro‑kinetic simulations. The approach opens the door to efficient, high‑fidelity modeling of coupled fluid‑electrostatic phenomena in increasingly complex micro‑ and nano‑scale devices.