Numerical Solution of nonlinear P-B Equation in unbounded domain for determining the EDL potential distributions around a spherical particle

Question

For my project I am studying a paper, namely "Perturbation solutions for the nonlinear Poisson–Boltzmann equation with a higher order-accuracy Debye–Huckel approximation" by Cunlu Zhao, Qiuwang Wang and Min Zeng, Zeitschrift für angewandte Mathematik und Physik volume 71, 140 (2020), DOI: 10.1007/s00033-020-01367-9.
I am stuck at a problem: I tried to solve a nonlinear P-B equation in an unbounded domain, in order to determine the EDL potential distributions around a spherical particle, by using the finite difference scheme given in that paper. But I found I am unable to figure out how specify boundary conditions on an infinite domain. Precisely the equation is the following $$ \frac{1}{r^2}\frac{d}{dr}\bigg(r^2 \frac{d\psi}{dr}\bigg)=\sinh(\psi), $$ while the boundary conditions are $$ \begin{cases} \psi|_{r=k}=\zeta\\ \psi|_{r\rightarrow\infty}=0 \end{cases}. $$

I am not figuring out how to construct the grid points in an infinite domain and solve it. Please help me dealing with this equation. Thanks in advance!

Answer 1

For large $r$ the equation gets close to $$ ψ''=\sinhψ\implies ψ'^2=2\coshψ+C $$ and with the asymptotic convergence at infinity $C=-2$ so that $$ ψ'^2-(2\sinh(ψ/2))^2=0. $$ The desired stable behavior for the far-field solution corresponds to the factor $$ψ'+2\sinh(ψ/2)=0.$$ This can be taken as upper boundary condition at some largish but finite upper interval end.

Let's try an example with $ψ(0.5)=1$.

def sys(r,u):
    ψ,θ = u
    return θ, sinh(ψ)-2*θ/r

def bc(ua,ub):
    ψa,θa = ua
    ψb,θb = ub
    return ψa-1, θb+2*sinh(ψb/2)

a,b = 0.5, 5 
r = a+logspace(-3,0,10)*(b-a)
u = [1-r/b, -1/b+0*r]

res = solve_bvp(sys,bc,r,u)

This works and gives a simple falling and asymptotic curve

Answer 2

For the case of a single flat plate, the potential varies as $\psi\propto \exp(-\kappa y)$, both for small applied surface potential ("Debye-Hückel theory") and for the full nonlinear Poisson Boltzmann equation (Gouy-Chapman). In the paper you mentioned, the authors have nondimensionalized all lengths by the Debye length $\kappa^{-1}$. Hence, $\psi(\bar{y})\propto \exp(-\bar{y})$. This means that already at $\bar{y}=5$, say, the potential is $\exp{(5)}\approx 148$ times smaller than at the electrode at $\bar{y}=0$ where the potential is applied. So presumably you will see that the difference between applying your second boundary condition either at $\psi(\bar{y}=5)=0$ and $\psi(\bar{y}=6)=0$ are very small. As Kurt G. says:

Choose a large upper bound , choose finite grid points and experiment with those parameters until your solution doesn't change anymore . Done.

The case of a spherical particle will be very similar. Chose a certain dimensionless sphere radius $K$ and apply the second boundary condition at $r=K+10$. You will see that the potential at $K+5$ will already be close to zero. Now set the second boundary condition at $r=K+6$ and observe that the differences between these two solution will be small.