  This work considers charged systems described by the modified
Poisson--Nernst--Planck (PNP) equations, which incorporate ionic steric effects
and the Born solvation energy for dielectric inhomogeneity. Solving the
steady-state modified PNP equations poses numerical challenges due to the
emergence of sharp boundary layers caused by small Debye lengths, particularly
when local ionic concentrations reach saturation. To address this, we first
reformulate the steady-state problem as a constraint optimization, where the
ionic concentrations on unstructured Delaunay nodes are treated as fractional
particles moving along edges between nodes. The electric fields are then
updated to minimize the objective free energy while satisfying the discrete
Gauss's law. We develop a local relaxation method on unstructured meshes that
inherently respects the discrete Gauss's law, ensuring curl-free electric
fields. Numerical analysis demonstrates that the optimal mass of the moving
fractional particles guarantees the positivity of both ionic and solvent
concentrations. Additionally, the free energy of the charged system
consistently decreases during successive updates of ionic concentrations and
electric fields. We conduct numerical tests to validate the expected numerical
accuracy, positivity, free-energy dissipation, and robustness of our method in
simulating charged systems with sharp boundary layers.
