Solving the P-N-P equations
Gib Bogle
behaviour of two ionic species (+ and - charge) in solution under the
influence of an applied potential field. I have a student who needs to
solve these equations as part of a larger problem. The FD method he is
using, in which he iteratively solves for the potential field and the
two concentration fields until convergence at each time step, is
extremely slow. It is necessary to use a very fine grid and small time
step. I am keen to hear from anyone else who has solved these
equations, to learn of any useful tricks to speed things up. One
obvious idea is to use a variable grid spacing (needed in one direction
only) - I'd like to know if there are pitfalls with this approach
related to possible loss of accuracy resulting from the way the
equations are discretized. Thanks.
Greg Heath
Decades ago the simulation group at the Stanford Institute of
Plasma Physics was able to increase the speed of the Poisson
equation solution by using the FFT whenever boundary conditions
in one direction were either periodic, Dirichlet or Neuman and the
boundaries in the orthogonal direction were uniform. The region was
either a rectangle (x,y) or an annulus (r,theta).
The principle investigator was the late Prof. Oscar Buneman.
I don't remember the journals in which they published.
Hope this helps.
Greg
William R. Frensley
A quick Google check verified what I had suspected: the P-N-P equations
are the same as the "Shockley Semiconductor Equations" used in
semiconductor device simulation. There is a great deal of published
work on numerical techniques for the latter, the vast majority of which
will send you in the wrong directions. One work which gives a very
useful view of the subject is Selberherr, Analysis and Simulation of
Semiconductor Devices (Springer 1984).
You have not specified the dimensionality of the problem: it appears to
be greater than one, so I will assume it to be two. (Three spatial
dimensions is not suitable for a student project.) You also imply
that you seek the transient (time-dependent) response of the system.
If this is simply a way to approach the steady-state solution, that
solution can be effectively solved using multi-dimensional Newton
techniques.
Now, some hardly known general principles are: (1) If your discrete
model requires a fine mesh to converge or avoid instabilities, you are
certainly using the wrong discretization. More comments on that below.
(2) Do not worry about discretization errors until you have a working,
robust simulation. An absurdly coarse mesh will tell you a great deal
more about the system than what you know right now; unless there are
large non-monotonic effects it will be qualitatively true; whether it is
"accurate" can be determined later. The unavoidable numerical
difficulties (due to ill-conditioning) increase as the number of mesh
nodes, so fine meshes are never an adequate solution.
(1) Problems that seem to be solved by using a finer mesh usually result
from discretizing a "small-signal" approximation to the differential
equation. This in fact what the textbooks tell you to do, but it is
exactly the wrong approach. You want to model the behavior of your
concentration fields between the meshpoints >as accurately as possible<
not just as a low-order polynomial. The way to do this is known in
semiconductor device simulation as the Scharfetter-Gummel scheme. Take
two adjacent meshpoints, assume that the potential varies linearly
between them, and analytically solve the Nernst-Planck (or
drift-diffusion) equation within this interval. This is readily done
using an inverse Boltzmann factor [exp(+V/kT) where V is the potential]
as an integrating factor. This gives you an expression for the flux
that reduces to the small-signal case for small potential differences,
but which correctly gives a barrier-limited Boltzmann form for large
potential differences.
(2) Using the Scharfetter-Gummel expression in the continuity equation
will give you a "transport operator" which contains the Laplace operator
from the diffusion equation, plus a correction for the drift term.
Again, if you really want a transient solution, what you need to do
is to use this operator in a ->fully implicit<- time-step algorithm.
That means you need to be able to solve the linear operator on a
two-dimensional discrete system. Selberherr explains the options for
doing so. I recommend direct solution using sparse-matrix techniques,
keeping the mesh as coarse as possible. The reason for keeping the
mesh coarse is that the "stiffness" of the time-stepping problem
increases with the number of mesh nodes. A finer discretization allows
shorter-wavelength components, which have faster decay rates in the
diffusion equation. That is the reason for choosing a fully-implicit
time evolution. It is stable at any size time step, and if you simply
step at an interval that is small compared to the times of physical
interest, you don't have to worry about excessive errors. (The issue
involves components of the solution that are decaying to zero, and
whether one is accurately representing their rate of decay. But if
they disappear before anything interesting happens in the system, who
cares how accurately they were represented?)
Now, while the fully-implicit scheme will stabilize the diffusion
equation, there is still a potential instability that operates via
the feedback through Poisson's equation. In the physical system,
this feedback is always negative and produces all the well-known
electrostatic screening effects in electrolytes, semiconductors,
plasmas, etc. But, if we time-step over an interval longer than the
characteristic time for this screening to operate, the corrections to
the potential overshoot, and numerical oscillations result. This
characteristic time is known as the "dielectric relaxation time"
and is equal to the dielectric constant (presumably of the solvent
in your case) divided by the (ionic) conductivity of the medium.
The simple solution is to keep the time step well below the shortest
dielectric relaxation time in the system (presumably in the regions
of highest concentration). If this is not feasible, then you will
need to use a fully-coupled solution scheme wherein Poisson's equation
is solved simultaneously with the transport equations.
- Bill Frensley
Gib Bogle
Thanks Greg, that might be a bridge too far.
Greg Heath
Why? Are your boundaries not coordinate aligned?
Hope this helps.
Greg
Gib Bogle
It's not so much my boundaries, as my student's (g), although there is
also a limit to how much effort I am prepared to devote to his research
topic.