Nitsche method for the Stokes-Poisson-Boltzmann equation with Navier slip boundary condition

Abstract: We study the Stokes--Poisson--Boltzmann equations with Dirichlet and Navier boundary conditions. The system consists of the incompressible Stokes equations coupled with a nonlinear Poisson--Boltzmann equation through electrostatic forcing and convective transport effects. To handle the Navier boundary conditions in a unified framework, we employ Nitsche's method for their weak imposition within a conforming finite element setting. We derive a consistent and stable discrete formulation and establish the well-posedness of the resulting problem. By carefully choosing the penalty parameters, the bilinear form is shown to be coercive and continuous. A priori error estimates are proved in the natural energy norms, yielding optimal-order convergence under suitable regularity assumptions. Furthermore, we develop residual-based a posteriori error estimators that incorporate element residuals, inter-element jump residuals, and boundary residuals arising from the Nitsche formulation. The estimators are shown to be reliable and locally efficient. Numerical experiments confirm the theoretical results and demonstrate the robustness and accuracy of the proposed method for the Stokes--Poisson--Boltzmann system.

• The paper develops a symmetric Nitsche finite element method that weakly enforces Navier slip conditions for the coupled Stokes–Poisson–Boltzmann system and proves discrete stability under a sufficiently large penalty parameter.
• The analysis establishes existence, uniqueness, quasi-optimal a priori error bounds, and optimal convergence rates, including second-order convergence for pressure and potential in P2–P1–P2 experiments under small-data assumptions.
• The paper introduces a residual-based a posteriori estimator that is reliable and locally efficient, with adaptive tests resolving re-entrant corners and boundary layers, although its effectivity shows some dependence on viscosity.

Model problem and continuous formulation

The paper analyzes the stationary Stokes–Poisson–Boltzmann (SPB) system on a bounded Lipschitz domain $\Omega \subset \mathbb{R}^n$, $n=2,3$, describing electro-osmotic flow of an incompressible electrolyte. The momentum equation $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$ is coupled to a nonlinear transport–diffusion equation for the double-layer potential $\psi$ through the charge density $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$ and convective terms $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$ and $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$. A notable modeling choice is that the $\varepsilon\Delta\psi$ forcing is rewritten using the potential equation itself, so that no integration by parts of $\Delta\psi$ is required in the momentum balance; this simplifies the subsequent analysis.

The boundary splits into a Dirichlet part $\Gamma_D$ (no-slip, homogeneous potential) and a Navier slip part $n=2,3$0, where $n=2,3$1 together with tangential friction conditions $n=2,3$2. Three structural assumptions are imposed: uniform boundedness of the potential (A1), global Lipschitz continuity of $n=2,3$3 with both upper and lower Lipschitz constants $n=2,3$4 on the bounded range of $n=2,3$5 (A2), and $n=2,3$6 (A3). Assumption (A2) is strong — it requires a priori knowledge that $n=2,3$7 stays within a bounded interval where $n=2,3$8 is bi-Lipschitz — and it underpins all well-posedness results.

The weak formulation seeks $n=2,3$9 with $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$0 incorporating the normal constraint on $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$1. Continuity of all forms follows from Sobolev embeddings ($$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$2 for $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$3 in 2D, $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$4 in 3D); ellipticity of $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$5 on the divergence-free kernel $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$6 and the inf-sup condition for $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$7 are inherited from prior work by Bansal et al., yielding a global inf-sup condition for the linearized form $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$8 via standard Babuška–Brezzi arguments.

Discrete scheme and well-posedness

The discretization uses Taylor–Hood-type spaces: velocity and potential in continuous piecewise polynomials of degree $$-\mu\Delta\boldsymbol{u} + \nabla p = \boldsymbol{f} - \varepsilon\Delta\psi\,\boldsymbol{E}$$9, pressure of degree $\psi$0. Since $\psi$1 does not enforce $\psi$2 on $\psi$3, the method is non-conforming in this sense, and the symmetric Nitsche variant weakly imposes the normal condition through consistency, adjoint-consistency, and penalty terms scaled by $\psi$4. The discrete energy norm augments the $\psi$5 norm with the facet term $\psi$6.

The stability analysis proceeds along familiar lines: inverse and trace inequalities give boundedness of $\psi$7, $\psi$8, and the trilinear forms in the mesh-dependent norm; ellipticity of $\psi$9 on the discrete kernel $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$0 holds provided the penalty parameter satisfies $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$1 with $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$2; and a discrete inf-sup condition holds for $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$3 uniformly in $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$4. These combine into a global inf-sup bound for $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$5 with constant $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$6.

Well-posedness of the nonlinear coupled system is established by decoupling: for fixed $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$7 in a ball $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$8 of radius $\mathcal{K}(\psi) = k_0\sinh(k_1\psi)$9, the flow subproblem is solved via Babuška–Brezzi theory; for fixed $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$0 in a ball $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$1, the potential subproblem is solved via Minty–Browder, exploiting strong monotonicity inherited from assumption (A2). Banach's fixed-point theorem then closes the loop. Crucially, uniqueness requires two small-data assumptions: a smallness condition on $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$2 relative to $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$3 and the model constants, and a contraction condition $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$4. This is a genuine restriction: existence and uniqueness are guaranteed only for sufficiently small data and permittivity not too small relative to the coupling strength, and the analysis does not address large-data regimes or multiple solutions.

A priori error estimates

The convergence proof uses a Céa-type argument built on a solution-dependent bilinear form $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$5 obtained by freezing the trilinear couplings at admissible iterates. An inf-sup estimate for this form holds on balls $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$6 of radius $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$7, and Galerkin orthogonality of the nonlinear residual yields quasi-optimality:

$$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$8

under the additional regularity-size assumptions $$\boldsymbol{c}_1(\boldsymbol{w};\psi,\boldsymbol{v}) = \int_\Omega (\boldsymbol{w}\cdot\nabla\psi)\boldsymbol{E}\cdot\boldsymbol{v}$$9, $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$0 with $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$1. Combined with Lagrange interpolation estimates, this gives optimal-order convergence

$\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$2

for solutions in $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$3. The constant $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$4 depends only on $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$5, $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$6, and $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$7, so the estimate is robust in the viscosity, although the smallness assumptions on the solution norms remain embedded in the argument.

A posteriori estimation

The residual-based estimator combines element residuals $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$8 (momentum), $\boldsymbol{c}_2(\boldsymbol{w};\psi,\phi)$9 (mass), $\varepsilon\Delta\psi$0 (potential equation), interior flux jumps $\varepsilon\Delta\psi$1 and $\varepsilon\Delta\psi$2, and — specific to the Nitsche treatment of slip — boundary residuals $\varepsilon\Delta\psi$3 (tangential friction defect) and $\varepsilon\Delta\psi$4 (normal velocity defect), weighted as $\varepsilon\Delta\psi$5. Data oscillation terms $\varepsilon\Delta\psi$6 account for piecewise polynomial approximation of $\varepsilon\Delta\psi$7 and $\varepsilon\Delta\psi$8.

Reliability rests on a global inf-sup stability result for the linearized operator in the continuous space, valid for frozen states of sufficiently small norm, plus a decomposition $\varepsilon\Delta\psi$9 into conforming and non-conforming parts, where the non-conforming component is controlled solely by the normal-jump indicator via $\Delta\psi$0. Under smallness assumptions on $\Delta\psi$1 and $\Delta\psi$2, the reliability bound

$\Delta\psi$3

holds. Efficiency is proved with classical Verfürth bubble-function techniques element-wise and facet-wise, including the trace-residual bound $\Delta\psi$4, giving local efficiency up to oscillation. Together these establish that $\Delta\psi$5 is a computable, reliable, and locally efficient error control for adaptive refinement — the first such estimator, to the authors' knowledge, for the SPB system with weakly imposed slip conditions.

Numerical experiments

All computations use FEniCS with MUMPS and Newton iteration, with a maximum-strategy Dörfler marking adaptive cycle. In the manufactured-solution test on the unit square with the $\Delta\psi$6–$\Delta\psi$7–$\Delta\psi$8 pair, the observed orders match theory: pressure and potential converge at rate 2.00 throughout, while velocity converges at rate ~2.03 for $\Delta\psi$9 and at superconvergent rates (~2.4–2.8) for $\Gamma_D$0. The effectivity index remains bounded but behaves differently across viscosities: approximately constant near 5.7 for $\Gamma_D$1, whereas for $\Gamma_D$2 it grows from about 0.97 to about 2.11 over the refinement sequence — still bounded, but indicating that the estimator's quality degrades somewhat in the low-viscosity regime rather than being uniformly sharp.

Three further experiments exercise the adaptive framework qualitatively. On non-convex C-, L-, and T-shaped domains with re-entrant corners, adaptive refinement concentrates degrees of freedom near corners where uniform refinement yields suboptimal rates. On a triangular domain with manufactured exponential boundary layers ($\Gamma_D$3), adaptivity restores optimal convergence once layers are resolved, with a stable effectivity index. Finally, a pipe-flow problem with a circular obstruction carrying the Navier condition demonstrates applicability to non-homogeneous inflow/outflow configurations, which the homogeneous-boundary analysis does not directly cover.

Limitations and open questions

Several restrictions qualify the results. Well-posedness and contraction require small data and a bounded, bi-Lipschitz regime for $\Gamma_D$4 (assumption A2), so large forcing, strong coupling, or potentials outside the assumed range fall outside the theory. The a priori and a posteriori bounds additionally assume small solution norms ($\Gamma_D$5, small $\Gamma_D$6, $\Gamma_D$7), and the reliability proof uses $\Gamma_D$8 control of the velocity that is not established from the primitive data. Only the symmetric Nitsche variant is analyzed; non-symmetric variants, curved boundaries with the associated geometric consistency issues treated by Gjerde and Scott, and time-dependent extensions remain open. Whether the effectivity index can be made uniformly bounded independently of $\Gamma_D$9 is not resolved by the numerical evidence presented.

Conclusion

This work provides a complete numerical analysis pipeline — well-posedness, optimal a priori bounds, and reliable/efficient residual-based a posteriori estimation — for the Stokes–Poisson–Boltzmann system with Navier slip conditions imposed by a symmetric Nitsche method. The theoretical guarantees are conditional on small-data and small-solution assumptions typical of nonlinear coupled elliptic systems, and the numerical experiments confirm the predicted rates while exposing a mild viscosity dependence of the effectivity index. The main open questions concern relaxing the smallness assumptions, extending the analysis to non-symmetric Nitsche variants and curved boundaries, and establishing $n=2,3$00-robust efficiency of the estimator.