Plasma Sheath Model in Fipy

[LOrellanaG]

Hi everybody!
I want to learn how to use Fipy and see if I can use it to simulate the plasma sheath phenomena. The following is what I want to model as first verificable example:

Continuity equation: $\nabla \cdot (n_i*u_i)=0$
Momentum equation: $u_i *(\nabla \cdot u_i)=-\partial_x\phi$
Poisson: $\nabla^2 \phi=n_i-\exp(-\phi)$

One solves for the variables $n_i$ (density), $u_i$ (ion speed), and $\phi$ (electric potential)
For a 1-Dimensional interval [0,600]
Boundary conditions: $\phi(0)=\phi_W$; $\phi(600)=0$; $\partial_x\phi(600)=0$; $u_i(600)=u_B$; $n_i(600)=n_0$.

This is the simple Bohm Model which I can compare with other program such as Comsol. The idea is to then be able to make it more complex as what there is in the literature.

Here is one implementation that I´ve tried:

from fipy import CellVariable, Grid1D, Viewer, DiffusionTerm,ImplicitSourceTerm,ConvectionTerm
from numpy import sqrt
from fipy.tools.numerix import exp

# Params
Te=2 # Elektron Temperatur
ep=8.85e-14 # epsilon 0
n=1e11 # Ladungsdichte
M=131*1.67e-27 # Xenon Mass
qe=1.6e-19 # Elektronsladung
Vw=10. # Potential an der Wand
lD=sqrt(ep*Te/qe/n) # debysche Länge
nlD=600 # X Definitionsbereich (Adimensional)
L=nlD*lD # X Definitionsbereich in physicalische Einheiten
nx=2000 # Nummer der Punkte in Definitionsbereich


# Mesh
dx=nlD/nx
mesh=Grid1D(dx=dx,nx=nx)
# Variablen
phi=CellVariable(mesh=mesh,name='Potential',value=0.)
ui=CellVariable(mesh=mesh,name='Vel',value=0.)
ni=CellVariable(mesh=mesh,name='n',value=0.)

#Equations
cont_eq=ConvectionTerm(coeff=ui*[[1]],var=ni)==0
mom_eq=ConvectionTerm(coeff=ui/2*[[1]],var=ui)==-phi.grad
poiss_eq=DiffusionTerm(coeff=1,var=phi)==ni-exp(-phi)

system_eq=cont_eq & mom_eq & poiss_eq

# Randbedingungen / Boundary
phi.constrain(0., mesh.facesRight)
phi.constrain(Vw, mesh.facesLeft)
ui.constrain(0.,mesh.facesRight)
ni.constrain(1.,mesh.facesRight)

# Solving
for i in range(10):  # iterate for convergence
    system_eq.sweep()

# Plot
viewer=Viewer(vars=(phi))
viewer.plot()

---

I have error when I run this this code, specifically in the part of the solver, system_eq.sweep(). First I had error that said that the coefficient of ConvectionTerm must be a vector, so I read somewhere else that one can write the coefficient (ui,) or ui*[[1]], but then I get the error "The coefficient must be rank 0 for a rank 0 solution variable."

If I use the notation for the coefficient (ui,) then I get the assertion error assert len(id1) == len(id2) == len(vector).

I believe it has something to do with the ConvectionTerm, but so far I haven´t been able to solve it from what I´ve read from similar problemss/discussions.
For the moment I would like only to check if the equations make sense for Fipy and solve this error that doesn´t allow me to at least run the solver.

Thank you for your time.

[guyer]

The issue isn't the ConvectionTerm, but rather the SourceTerm in mom_eq. phi.grad is rank-1, but your equation is rank-0. Change

-mom_eq=ConvectionTerm(coeff=ui/2*[[1]],var=ui)==-phi.grad
+mom_eq=ConvectionTerm(coeff=ui/2*[[1]],var=ui)==-phi.grad.dot([[1]])

At that point, you'll likely get a singular matrix error (exactly what will depend on which solver suite you're using). This can be seen with just cont_eq; with $u_i = 0$, then any $n_i$ is a solution to $\nabla\cdot(n_i * u_i) = 0$.

[LOrellanaG]

Guyer, thanks for the clarification.
I will have to recheck the equations to see how to get rid of the singular matrix error, but that´s another topic.
Thank you for your time!

[guyer]

Introducing time relaxation (with a TransientTerm) on one or more of your equations (I'd say continuity and momentum, but not Poisson) may help with the singular matrices.