Minimal documentation

README.md

@@ -0,0 +1,14 @@
+## Overview
+
+`glowDischargeSolver` is a collection of python functions supporting
+the approximation of plasma systems in "glow discharge" regime.  The
+model takes the form of a fluid model that uses the drift-diffusion
+approximation so that the state variables consist of a set of species
+number densities and an electron temperature.  For more details on the
+physical model, see the LaTeX documentation in the `doc` directory.
+
+These equations are discretized in space using Chebyshev collation and
+in time using fully implicit schemes, either backward Euler or
+Crank-Nicolson.  Further, to support the time-periodic case,
+time-periodic-shooting is implemented to accelerate converged to the
+time-periodic solution.

doc/.gitignore

@@ -0,0 +1,11 @@
+vc.tex
+*~
+glow.aux
+glow.fdb_latexmk
+glow.fls
+glow.log
+glow.pdf
+glow.bbl
+glow.blg
+glow.out
+glow.toc

doc/Makefile

@@ -0,0 +1,23 @@
+PAPER    := glow.pdf
+TEX_SUFS := .aux .log .nav .out .snm .toc .vrb .fdb_latexmk .bbl .blg .fls
+GITID    := $(shell git describe --always 2> /dev/null)
+TEXSRC   := $(wildcard *.tex bolsig/*.tex fluid/*.tex hybrid/*.tex)
+FIGURES  := $(wildcard figures/*.pdf figures/*.jpg)
+BIBS     := $(wildcard bibs/*.bib)
+
+%.pdf:%.tex $(TEXSRC) $(BIBS) $(FIGURES) gitinfo2.pm
+	echo "hello"
+ifdef GITID
+	./gitinfo2.pm
+endif
+	latexmk -pdf $<
+
+all: $(PAPER)
+
+clean:
+	$(RM) $(foreach suf, ${TEX_SUFS}, $(PAPER:.pdf=${suf}))
+ifdef GITID
+	$(RM) vc.tex
+endif
+	$(RM) $(PAPER)
+

doc/appendix.tex

@@ -0,0 +1,77 @@
+\appendix
+\chapter{Non-dimensionalization}
+
+For the non-dimensionalization of the system of equations, the following dimensionless variables are used
+%
+\begin{align*}
+t' = \frac{t}{\tau}, \quad \phi' = \frac{\phi}{V_0}, \quad x' = \frac{x}{L},  \quad n_{i,e}' = \frac{n_{i,e}}{n_{p0}}, \quad n_{\alpha}' = \frac{n_{\alpha}}{n_{p0}}, \quad \left(n_{e}e_{e}\right)' = \frac{n_{e}e_{e}}{n_{p0} e_0}, \quad ne_t' = \frac{ne_t}{n_{Ar} e_0}.
+\end{align*}
+%
+The Poisson equation governing the
+electric potential $\phi$ takes the following form:
+%
+\begin{align*}
+-\pp{^2\phi}{x^2} = \frac{\rho_c}{\epsilon_0} = \frac{q_e}{\epsilon_0} \left(n_i - n_e\right) \quad &\Rightarrow \quad \pp{^2\phi}{x^2} = - \frac{q_e}{\epsilon_0} \left(n_i - n_e\right) \quad \Rightarrow \quad \frac{V_0}{L^2}\pp{^2\phi'}{x'^2} = - \frac{q_e}{\epsilon_0} \left(n_i' - n_e'\right) n_{p0} \\
+\quad &\Rightarrow \quad \pp{^2\phi'}{x'^2} = \color{red}\boxed{\color{black} - \frac{q_e n_{p0} L^2}{\epsilon_0 V_0}} \color{black}\left(n_i' - n_e'\right) .
+\end{align*}
+%
+The species transport equations can be non-dimensionalized as follows:
+%
+\begin{align*}
+\pp{n_{\alpha}}{t} + \pp{}{x} \left( - z_{\alpha} \mu_{\alpha} n_{\alpha} \pp{\phi}{x} - D_{\alpha} \pp{n_{\alpha}}{x}\right)= k n_{Ar} n_{e}  \ &\Rightarrow \ \frac{1}{\tau}\pp{n_{\alpha}'}{t'} +  \frac{1}{L}\pp{}{x'} \left( - z_{\alpha} \mu_{\alpha} n_{\alpha}' \frac{V_0}{L}\pp{\phi'}{x'} - D_{\alpha} \frac{1}{L} \pp{n_{\alpha}'}{x'}\right) = k n_{Ar} n_{e}' \\
+\ &\Rightarrow \ \pp{n_{\alpha}'}{t'} - \pp{}{x'} \left( z_{\alpha} \color{red}\boxed{\color{black} \mu_{\alpha} \frac{\tau V_0}{L^2} } \color{black} n_{\alpha}'\pp{\phi'}{x'} + \color{red}\boxed{\color{black} D_{\alpha} \frac{\tau}{L^2}} \color{black}  \pp{n_{\alpha}'}{x'}\right) = \color{red}\boxed{\color{black} k n_{Ar} \tau } \color{black} n_{e}' .
+\label{eqn:species}
+\end{align*}
+%
+The electron energy equation can be non-dimensionalized as below:
+%
+\begin{align*}
+\pp{}{t} \left( n_e e_e \right)
++
+\pp{}{x} \left( - \frac{5}{3} D_e \pp{\left(n_e e_e\right)}{x} \right. \\
+\left. - \frac{5}{3}  z_e \mu_e \pp{\phi}{x} \left(n_e e_e\right) \right)
+&=
+\underbrace{q_e \left(-z_{e} \mu_{e} n_{e} \pp{\phi}{x} - D_{e} \pp{n_{e}}{x} \right) \pp{\phi}{x}}_{\textrm{Joule heating}} \\
+&-
+\underbrace{\frac{3}{2} k_B n_e \frac{2 m_e}{m_b} \left( T_e - T_g \right) \bar{\nu}_{e,b}}_{\textrm{Elastic}} \\
+&-
+\underbrace{q_e \sum_{j=1}^{N_r} \Delta E^e_j G_j}_{\textrm{Inelastic}} \ \Rightarrow
+\end{align*}
+\begin{align*}
+{\color{blue}\frac{n_{p0} e_0}{\tau}} \pp{}{t'} \left( n_e e_e \right)' 
++
+{\color{blue}\frac{1}{L}} \pp{}{x} \left( - \frac{5}{3} D_e {\color{blue}\frac{n_{p0} e_0}{L}} \pp{\left(n_e e_e\right)'}{x'}  \right. \\
+\left. - \frac{5}{3}  z_e \mu_e {\color{blue}\frac{V_{0}}{L}} \pp{\phi'}{x'} {\color{blue}n_{p0} e_0} \left(n_e e_e\right)' \right)
+&=
+\underbrace{q_e \left(-z_{e} \mu_{e}  {\color{blue}\frac{V_{0} n_{p0}}{L}} n_{e}' \pp{\phi'}{x'} - D_{e}  {\color{blue}\frac{n_{p0}}{L}} \pp{n_{e}'}{x'} \right) {\color{blue}\frac{V_{0}}{L}} \pp{\phi'}{x'}}_{\textrm{Joule heating}} \\
+&-
+\underbrace{\frac{3}{2} k_B n_e \frac{2 m_e}{m_b} \left( T_e - T_g \right) \bar{\nu}_{e,b}}_{\textrm{Elastic}} \\
+&-
+\underbrace{q_e \sum_{j=1}^{N_r} \Delta E^e_j G_j}_{\textrm{Inelastic}} \ \Rightarrow
+\end{align*}
+\begin{align*}
+\pp{}{t'} \left( n_e e_e \right)' 
++
+\pp{}{x} \left( - \frac{5}{3} {\color{red}\boxed{\color{black} D_e \frac{\tau}{L^2}}} \color{black}\pp{\left(n_e e_e\right)'}{x'}  \right. \\
+\left. - \frac{5}{3}  z_e  {\color{red}\boxed{\color{black} \mu_e \frac{\tau V_{0}}{L^2}}} \color{black}\pp{\phi'}{x'} \left(n_e e_e\right)' \right)
+&=
+\underbrace{q_e \left(-z_{e} {\color{red}\boxed{\color{black} \mu_{e} \frac{V_{0} \tau}{L^2}}} \color{black}n_{e}' \pp{\phi'}{x'} - {\color{red}\boxed{\color{black} D_{e} \frac{\tau}{L^2}}} \color{black} \pp{n_{e}'}{x'} \right) {\color{red}\boxed{\color{black} V_{0}}} \color{black} \pp{\phi'}{x'}}_{\textrm{Joule heating}} \\
+&-
+\underbrace{\frac{3}{2} {\color{blue}\frac{\tau}{n_{p0} e_0}} k_B n_e \frac{2 m_e}{m_b} \left( T_e - T_g \right) \bar{\nu}_{e,b}}_{\textrm{Elastic}} \\
+&-
+\underbrace{q_e {\color{blue}\frac{\tau}{n_{p0} e_0}} \sum_{j=1}^{N_r} \Delta E^e_j G_j}_{\textrm{Inelastic}} \ \Rightarrow
+\end{align*}
+%
+The non-dimensionalization of the ideal gas law is presented below:
+%
+\begin{align*}
+p = n_e k_B T_e + k_B T_g \sum_{\alpha \in h} n_{\alpha}  \quad \Rightarrow \quad {\color{blue} n_{Ar} k_B T_0} \ p' = {\color{blue} n_{p0} e_0} \frac{2}{3} \left(n_e e_e\right)' + {\color{blue} n_{p0} e_0}  \frac{2}{3} \left(k_B T_g \sum_{\alpha \in h} n_{\alpha} \right)'.
+\end{align*}
+%
+\begin{align*}
+\pp{ne_t}{t} = &- \frac{5}{3} \left( D_{e}\pp{n_e e_e}{x} + z_{e}\mu_{e} \pp{V}{x} n_e e_e \right) \\
+ &- \frac{5}{3} \sum_{\alpha \in h} \left\{ \left( D_{a}\pp{\left( \frac{3}{2} k_B  n_a T_g\right)}{x} + z_{a} \mu_{a} \pp{V}{x} \left( \frac{3}{2} k_B n_a T_g\right) \right)  \right\} \\
+ &+ \sum_{\alpha \in h} \left\{ \varepsilon_{\alpha} \left( - z_{\alpha} \mu_{\alpha} n_{\alpha} \pp{V}{x} - D_{\alpha} \pp{n_{\alpha}}{x} \right) \right\} +
+\underbrace{\sum_{\alpha} q_e z_{\alpha} F_{\alpha} E}_{\textrm{Joule heating}}
+\end{align*}
+%

doc/bibs/glow.bib

@@ -0,0 +1,10 @@
+@article{panneer2015computational,
+  title={Computational modeling of the effect of external electron injection into a direct-current microdischarge},
+  author={Panneer Chelvam, Prem Kumar and Raja, Laxminarayan L},
+  journal={Journal of Applied Physics},
+  volume={118},
+  number={24},
+  pages={243301},
+  year={2015},
+  publisher={AIP Publishing LLC}
+}

doc/bolsig/eqns.tex

@@ -0,0 +1,112 @@
+\section{BOLSIG+ Governing Equations}
+The electron Boltzmann equation may be written as follows:
+%
+\begin{equation}
+\pp{f}{t}
++
+\mbf{v} \cdot \nabla f
+-
+\frac{e}{m} \mbf{E} \cdot \nabla_{\mbf{v}} f
+=
+C[f],
+\label{eqn:electron-boltzmann-full}
+\end{equation}
+%
+where $f(\mbf{x}, \mbf{v}, t)$ denotes the electron distribution
+function, $\mbf{v}$ is the velocity vector, $\nabla f$ is the spatial
+gradient of $f$ (i.e., wrt $\mbf{x}$), $\nabla_{\mbf{v}} f$ is the
+velocity gradient of $f$ (i.e., wrt $\mbf{v}$), $e$ is the absolute
+value of the charge of an electron, $m$ is the mass of an electron,
+$\mbf{E}$ is the electric field, and $C[f]$ denotes the collision
+term.
+
+The development here is based upon the description of Hagelaar and
+Pitchford~\cite{} (referred to as HP for short) and the BOLSIG+ User's
+Manual~\cite{}.
+
+To begin, it is assumed that the electric field and collision
+probabilities are spatially uniform.  In this situation, HP asserts
+that the distribution function $f$ is 1) ``symmetric in velocity space
+around the electric field direction'' and 2) ``may vary only along the
+field direction'' in physical space.  The meaning of 1) is somewhat
+unclear, but in the subsequent analysis, a spherical coordinate system
+is used and no possible dependence on the azimuthal angle is included.
+Thus, the combined import of these assumptions is to reduce the
+dependence of $f$ from 6D plus time to 3D plus time.  Specifically,
+letting $z$ denote the direction of the electric field, $v$ denote the
+velocity magnitude, and $\theta$ denote the polar angle between the
+electric field direction and the velocity direction, we have $f = f(z,
+v, \theta, t)$.  In this situation,~\ref{eqn:electron-boltzmann-full}
+simplifies to
+%
+\begin{equation}
+\pp{f}{t}
++
+v \cos \theta \pp{f}{z}
+-
+\frac{e}{m} E \left( \cos \theta \pp{f}{v} - \frac{\sin \theta}{v} \pp{f}{\theta} \right)
+=
+C[f].
+\label{eqn:electron-boltzmann-reduced}
+\end{equation}
+%
+
+In preparation for subsequent development, let $\xi = \cos \theta$.
+Then,~\ref{eqn:electron-boltzmann-reduced} becomes
+%
+\begin{equation}
+\pp{f}{t}
++
+v \xi \pp{f}{z}
+-
+\frac{e}{m} E \left( \xi \pp{f}{v} + \frac{(1-\xi^2)}{v} \pp{f}{\xi} \right)
+=
+C[f].
+\label{eqn:electron-boltzmann-xi}
+\end{equation}
+%
+To continue, the following two-term expansion is introduced:
+%
+\begin{equation}
+f(v, \xi, z, t) = f_0(v, z, t) + \xi f_1(v, z, t).
+\label{eqn:two-term}
+\end{equation}
+%
+To derive governing equations for $f_0$ and $f_1$, we apply the
+Galerkin method in the $\xi$ coordinate.  Specifically,
+substitute~\eqref{eqn:two-term} into
+~\eqref{eqn:electron-boltzmann-xi} and integrate against $1$ and
+$\xi$.  Rewriting the result in terms of a rescaled energy variable
+$\epsilon \equiv (v^2 / \gamma^2)$, where $\gamma^2 = 2e / m$, gives
+%
+\begin{gather}
+\pp{f_0}{t}
++
+\frac{\gamma}{3} \sqrt{\epsilon} \pp{f_1}{z}
+-
+\frac{\gamma}{3} \frac{E}{\sqrt{\epsilon}} \pp{}{\epsilon} \left( \epsilon f_1 \right)
+=
+C_0
+\label{eqn:f0}\\
+%
+\pp{f_1}{t}
++
+\gamma \sqrt{\epsilon} \pp{f_0}{z}
+-
+\gamma E \sqrt{\epsilon} \pp{f_0}{\epsilon}
+=
+C_1
+\label{eqn:f1}
+\end{gather}
+%
+where
+%
+\begin{gather*}
+C_0 = \frac{1}{2} \int_{-1}^{1} C[f_0 + \xi f_1] \, d\xi\\
+C_1 = \frac{3}{2} \int_{-1}^{1} \xi C[f_0 + \xi f_1] \, d\xi.
+\end{gather*}
+%
+Note that~\eqref{eqn:f0} and~\eqref{eqn:f1} are equivalent to HP (5)
+and (6), respectively.  Although, to match the forms in HP more
+precisely, we must examine the collision terms in more detail.
+\todo[inline]{Collisions}

doc/bolsig/overview.tex

@@ -0,0 +1,6 @@
+This chapter documents the models and data that support determination
+of transport coefficients and reaction rates.  For electron transport
+properties and electron-impact reactions, the coefficients and rates
+are determined using the BOLSIG+
+code\footnote{http://www.bolsig.laplace.univ-tlse.fr/}, which solves
+an approximate form of the electron Boltzmann equation.

doc/fluid/bcs.tex

@@ -0,0 +1,74 @@
+\section{Boundary Conditions} \label{sec:fluid-bcs}
+
+\subsection{Species Continuity}
+We impose conditions on the species fluxes at the walls.  The electron
+flux $F_e$ is given by
+%
+\begin{equation*}
+F_e \mbf{n}_x
+=
+\frac{1}{4} n_e \left( \frac{8 k_B T_e}{\pi m_e} \right)^{1/2}
+-
+\sum_{\alpha=1}^{N_s} \gamma_{\alpha} F_{\alpha} \mbf{n}_x
+\end{equation*}
+%
+where $\mbf{n}_x$ denotes the outward pointing unit normal (in 1D,
+$\mbf{n}_x =-1$ on the 'left' and $\mbf{n}_x = 1$ on the 'right') and
+$\gamma_{\alpha}$ is the secondary electron emission coefficient for
+species $\alpha$.  \todo[inline]{We need values for the
+  $\gamma_{\alpha}$ parameters!}
+
+The ion flux, $F_i$, is given by
+%
+\begin{equation*}
+F_i \mbf{n}_x
+=
+\frac{1}{4} n_i \left( \frac{8 k_B T_g}{\pi m_i} \right)^{1/2}
++
+n_i \max(0, \mu_i E \mbf{n}_x),
+\end{equation*}
+%
+\todo[inline]{Check inconsistency in refs for ion BC with Raja}
+
+Finally, the all neutral species fluxes ars given by
+%
+\begin{equation*}
+F_{\alpha} \mbf{n}_x
+=
+\frac{1}{4} n_{\alpha} \left( \frac{8 k_B T_g}{\pi m_{\alpha}} \right)^{1/2}
+\end{equation*}
+%
+
+\subsection{Electric Potential}
+Dirichlet conditions on the electric potential are imposed on both
+sides.  For instance, assuming $x=0$ is ground and $x=L$ has an
+imposed fluctuating voltage:
+%
+\begin{equation*}
+\phi(0,t) = 0, \quad \phi(L,t) = V_0 \sin(2 \pi f t),
+\end{equation*}
+%
+where $V_0$ is the voltage amplitude and $f$ is the frequency.
+
+\subsection{Energy}
+For the electron energy equation, the total surface energy flux at the
+wall is given by
+%
+\begin{equation*}
+Q_{e,wall} \mbf{n}_x
+\equiv
+\left. \left( \frac{5}{2} k_B T_e F_e - \kappa_e \pp{T_e}{x} \right)\right|_{wall} \mbf{n}_x
+= 
+\frac{5}{2} k_B T_e F_e \mbf{n}_x,
+\end{equation*}
+%
+where $F_e \mbf{n}_x$ is taken from the species continuity boundary
+condition.
+
+Isothermal wall conditions are imposed on the gas temperature:
+%
+\begin{equation*}
+T_g(0,t) = T_g(L,t) = T_{g,wall},
+\end{equation*}
+%
+where $T_{g,wall}$ is a specified value.