Chapter 04: Conservation laws¶

Fluid theory¶

Once collisions are numerous enough to average particle motion to that of a fluid, we can use a simpler model. This model is now using macroscopic quantities like the mass density $\rho$ $(kg/m^3)$, the velocity $\vec v$ $(m/s)$ or the kinetic energy density $\rho v^2$ $(J/m^3)$.

Kinetic theory¶

The distribution function¶

As we saw in the previous chapter, the best way to define collisions is to use particles which are labeled using invariant quantities of motion (rest mass $m$ and charge $q$) as well as quantities that varies with motion (the position $\vec r$ and the velocity $\vec v$).

If we are at a fix position is space $\vec r=(x,y,z)$ we can count all the particles inside an infinitely small small volume $dV=dxdydz$. Let's say that this number is $dN$. Once all these particles have been counted we can sort them as a function of their velocities. We end up with a distribution of particles across all available velocities at the location $\vec r$. Now, this distribution $f$ is a function of both $\vec r$ and $\vec v$ and we can write it as $f(\vec r,\vec v,t)$. $f$ is simply the number or particles per $m^3$ and $(m/s)^3$.

One can use the momentum $\vec p=m\vec v$ instead and then $f$ is now function of both $\vec r$ and $\vec p$. $$f=f(x,y,z,p_x,p_y,p_z,t)$$ So $f$ has a total of six coordinates plus time. To differentiate between the two different spaces (space and momentum space) the distribution function is said to be a function of space $(x,y,z)$ and phase space $(p_x,p_y,p_z)$.

Boltzmann equation¶

For a closed system with no source or sink of:

1. particles (i.e. mass conservation)
2. momentum (i.e. momentum conservation)
3. energy (i.e. energy conservation)

we have $$df(\vec r,\vec p,t)=0\tag{4.1}$$

Note that Eq. $(4.1)$ is a total derivative since both $\vec r$ and $\vec p$ are both function of time. Evidently, we can use the chain rule to get $$df(\vec r,\vec p,t)=\frac{\partial f}{\partial t}dt+\frac{\partial f}{\partial \vec r}d\vec r+\frac{\partial f}{\partial \vec p}d\vec p$$

Clearly, we can rewrite this equation using the fact that $$\frac{\partial f}{\partial \vec r}d\vec r=\frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz=\vec \nabla_{\vec r}\,f.d\vec r\tag{4.2}$$ and $$\frac{\partial f}{\partial \vec p}d\vec p=\frac{\partial f}{\partial p_x}dp_x+\frac{\partial f}{\partial p_y}dp_y+\frac{\partial f}{\partial p_z}dp_z=\vec \nabla_{\vec p}\,f.d\vec p\tag{4.3}$$

So we are left with $$\frac{\partial f}{\partial t}dt+\vec \nabla_{\vec r}\,f.d\vec r+\vec \nabla_{\vec p}\,f.d\vec p=0$$ Clearly, we can divide by $dt$ and use the fact that $$\frac{d\vec p}{dt}=\vec F$$ to get $$\frac{\partial f}{\partial t}+\frac{\vec p}{m}. \vec \nabla_{\vec r}\,f+\frac{\vec F}{m}.\vec \nabla_{\vec v}\,f=0\tag{4.4}$$

Note that, besides using $d\vec p/dt=\vec F$, there is no physics added to Eq. $(4.4)$. It's all math! From now on, we are going to use velocity phase space coordinate $\vec v$ instead of momentum phase space coordinate $\vec p$, so eq. $(4.4)$ becomes $$\frac{\partial f}{\partial t}+\vec v. \vec \nabla_{\vec r}\,f+\frac{\vec F}{m}.\vec \nabla_{\vec v}\,f=0\tag{4.4}$$

So when does the physics comes into play? During interaction between particles, a.k.a. collisions. These collisions are source or sink of particles. For instance, consider collisions that push particles from a velocity $\vec v_1$ to a velocity $\vec v_2$, then $f(\vec r,\vec v_1,t)$ decreases with time while $f(\vec r,\vec v_2,t)$ increases with time. For now, we encapsulate collisions inside a symbolic function $$\Big(\frac{\partial f}{\partial t}\Big)_{coll}$$ and we get, $$\frac{\partial f}{\partial t}+\vec v. \vec \nabla_{\vec r}\,f+\frac{\vec F}{m}.\vec \nabla_{\vec v}\,f=\Big(\frac{\partial f}{\partial t}\Big)_{coll}$$ Since $\vec F/m=\vec a$ we get the final equation, known as Boltzmann equation, $$\frac{\partial f}{\partial t}+\vec v. \vec \nabla_{\vec r}\,f+\vec a.\vec \nabla_{\vec v}\,f=\Big(\frac{\partial f}{\partial t}\Big)_{coll}\tag{4.5}.$$

Fig.1 Change of f as a function of time inside the volume dV due to a spacial variation of f along the x axis and a velocity parallel to that gradient.

The figure above shows that $f$ changes in time due to a change in f along the x-direction and a velocity $v_x$ that moves this change insdie the volume dV.

Vlasov equation¶

Vlasov realized that to properly capture Coulomb collisions, one should not use the collision operator $(\frac{\partial f}{\partial t})_{coll}$ but split the charge particle distribution function into two distribution functions:

1. one for the electrons, $f_e$
2. one for the ions, $f_i$

We suppose here that there is only one ions species with ionization $Z_i$. However, if other collisions (e.g. charged particles-neutrals) are not negligible then they have to be taken into account.

The total system used by Vlasov let the two distribution functions interact via the electric and magnetic fields that these particles generate. $$\begin{split} \frac{\partial f_e}{\partial t} + \vec v.\vec \nabla_{\vec r}\, f_e - \frac{e}{m_e}\left(\vec {E}+\vec v\times\vec {B}\right).\vec \nabla_{\vec v}\,f_e = 0 \\ \frac{\partial f_i}{\partial t} + \vec v.\vec \nabla_{\vec r}\,f_i +\frac{ Z_i e}{m_i}\left(\vec E+\vec v\times\vec B\right).\vec \nabla_{\vec v}\,f_i = 0 \end{split} \tag{4.6} $$ So how do we find these global fields?

We used Maxwell's equations $$ \begin{split} \vec \nabla\times\vec {B} =\mu_0\vec {j}+\mu_0\varepsilon_0 \frac{\partial\vec {E}}{\partial t} \\ \vec \nabla\times\vec {E} =-\frac{\partial\vec {B}}{\partial t} \\ \vec \nabla\cdot\vec {E} =\frac{\rho_E}{\varepsilon_0} \\ \vec \nabla\cdot\vec {B} =0 \end{split} \tag{4.7}$$ By global, Vlasov meant that these fields are the fields created by all the particles at the location $(x,y,z)$ but average over all the particle momenta at this location.

The global fields are defined by the charge density $\rho_E(x,y,z)$ and by the electrical current $\vec j(x,y,z)$. Note that these quantities are now independent of the particle momenta. Indeed, this feat is only possible by using $$\rho_E=e\int(Z_i f_i-f_e)d^3v \\ \vec j=e\int(Z_i f_i \vec v -f_e \vec v)d^3v $$ Basically, we integrated the distribution function over all velocities at each location $(x,y,z)$ to get local quantities where the velocity information has been encapsulated. This approach reduces the amount of information we have to deal with.

Moment equations¶

What do we obtain if we integrate the distribution function $f_i$ or $f_e$ over all velocities?

We obtain the number density in $particles/m^3$ $$n_{i,e}=\int f_{i,e}(x,y,z,v_x,v_y,v_z,t)d^3v\tag{4.8}$$ Note that we used the subscript $i$ and $e$ in the same equations rather than writing one equation for the electrons and one for the ions.

What happens if the integrate the flux of the distribution function $\vec \Gamma_{f_{i,e}}=f_{i,e}\vec v$ over all velocities ?

We obtain the total flux of particles $\vec \Gamma_{i,e}$ $$\vec \Gamma_{i,e}=\int f_{i,e}(x,y,z,v_x,v_y,v_z,t)\vec v d^3v$$ Since we have defined the flux as a number density times a velocity, then we can get a flow velocity $\vec u_{i,e}$ from dividing the flux by the particle density $n_{i,e}$ $$\vec u_{i,e}=\frac{\vec \Gamma_{i,e}}{n_{i,e}}$$ which we can rewrite as $$n_{i,e}\vec u_{i,e}=\int f_{i,e}(x,y,z,v_x,v_y,v_z,t)\vec v d^3v\tag{4.9}$$ Note that the flow velocity $\vec u$ is different for ions and electrons, while the velocity $\vec v$ is a phase space coordinate of the distribution function, like $\vec x$ is a spacial coordinate. So if $\vec x$ has no subscripts, then $\vec v$ should not either!

What is the physical meaning of $m_{i,e}n_{i,e}v^2$?

This is the kinetic energy density. What happens if we multiply the distribution function by $\vec v^2/2$?

We get the density of kinetic energy $\epsilon_{i,e}$ (not to confuse with $\varepsilon$ which is the permitivity of materials). $$\epsilon_{i,e}=\int f_{i,e}(x,y,z,v_x,v_y,v_z,t)\frac{\vec v^2}{2} d^3v\tag{4.10}$$ This energy density of the fluid is made up of energy at the macroscopic scale (i.e. kinetic energy) and energy at the miscroscopic scale, i.e. magnetic energy and internal energy. In our case, internal energy is simply pressure. $$m_{i,e}\epsilon_{i,e}=\frac{1}{2}m_{i,e}n_{i,e}\vec u_{i,e}^2+\frac{3}{2}p_{i,e}$$ Since $\epsilon_{i,e}$, $n_{i,e}$ and $\vec u_{i,e}$ are know from the previous moment equations, we can define the pressure as $$\frac{3}{2}p_{i,e}=m_{i,e}\epsilon_{i,e}-\frac{1}{2}m_{i,e}n_{i,e}\vec u_{i,e}^2$$ Here $p_{i,e}(\vec r,t)$ is the isotropic pressure. Note that $\frac{3}{2}p_{i,e}$ corresponds to the internal (i.e. microscopic) energy density, while $\frac{1}{2}m_{i,e}n_{i,e}\vec u_{i,e}^2$ correspond to the flow kinetic energy density.

Let's compute numerically the particle density for the distribution function given by $$f(\vec r,\vec v,t)=A_0\exp\Big(-\sqrt{\frac{x^2+y^2+z^2}{L_0^2}}\Big)\exp\Big(-\frac{v_x^2+v_y^2+v_z^2}{v_0^2}\Big)$$ where $A_0=1000\,particles/m^6s^3$, $L_0=1\mu m$ and $v_0=0.1\,m/s$.

So the answer is $n(x,y,z)=5.568327996831708\exp\Big(-\sqrt{\frac{x^2+y^2+z^2}{L_0^2}}\Big)\,particles/m^3$. Can we compute the previous integral analytically?

Yes. The integral of $\exp(-x^2)$ from $-\infty$ to $+\infty$ is

$$\int_{-\infty}^{+\infty}\exp(-x^2)dx=\sqrt{\pi}$$

Since we know that $$\int\int\int_{-\infty}^{+\infty}\exp\Big(-\frac{v_x^2+v_y^2+v_z^2}{v_0^2}\Big)dv_xdv_ydv_z=\Big(\int_{-\infty}^{+\infty}\exp(-\frac{v_x^2}{v_0^2})dv_x\Big)^3$$

we get $$\int\int\int_{-\infty}^{+\infty}\exp\Big(-\frac{v_x^2+v_y^2+v_z^2}{v_0^2}\Big)=(v_0\sqrt{\pi})^3$$

So the answer is $n(x,y,z)=A_0(v_0\sqrt{\pi})^3\exp\Big(-\sqrt{\frac{x^2+y^2+z^2}{L_0^2}}\Big)\,particles/m^3$. The particles have a Gaussian distribution.

Two fluid magnetohydynamics¶

Now that we have integrated $f$, $f\vec v$ and $f\vec v^2$ over all possible velocities, let's redo the same thing with $df$, $df\vec v$ and $df \vec v^2$.

Conservation of mass¶

The conservation of mass in obtained by integrating Vlasov's equation $$\int\int\int_{-\infty}^{+\infty}\frac{\partial f_e}{\partial t} + \vec v.\vec \nabla_{\vec r}\, f_e +\vec a.\vec \nabla_{\vec v}\,f_ed^3v=\int\int\int_{-\infty}^{+\infty}\Big(\frac{\partial f}{\partial t}\Big)_{e\,coll}d^3v\tag{4.11}$$ Since the left hand side of Eq. $(4.11)$ is integrated for all possible velocities, we will get the same results for electrons and ions. However, the right hand side of Eq. $(4.11)$ is dependent on cross-sections and has to be treated separately for ions and electrons.

When we integrate eq. $(4.11)$ over all velocities, we find the first conservation equation, called the conservation of mass or continuity equation $$\frac{\partial n_e}{\partial t}+\vec \nabla . \big(n_e\vec u_e\big)=S_e\tag{4.12}$$ $S_e$ is the sources of particles in this equations and comes from $$S_e=\int\int\int_{-\infty}^{+\infty}\Big(\frac{\partial f}{\partial t}\Big)_{e\,coll}d^3v$$ Note that these collisions are all but Coulomb collisions in the context set by Vlasov. If $S_e$ is positive then the electron density $n_e$ at the location $(x,y,z)$ is increasing or the flux of particle leaving this location is non-zero and points away from this volume. This is the meaning of the divergence operator $\vec\nabla.\bullet$

Physically, what can be considered a source of electrons at a given location?

Ionizations.

Note that injecting electrons at this location does not count because the represent an electron flux and this would be accounted for in $\vec \nabla . \big(n_e\vec u_e\big)$. When $S_e$ is negative then there is a sink of electrons. Physically, what can be considered a sink of electrons at a given location?

Recombinations.

It is usual to break $S_e$ down into two positive quantities, $G_e$, the term corresponding to electron sources or gains, and $L_e$, the term corresponding to electrons sinks or losses. So, we get $$\frac{\partial n_e}{\partial t}+\vec \nabla . \big(n_e\vec u_e\big)=G_e-L_e.\tag{4.13}$$

One way to understand conservation of mass (with no subscripts) $\frac{\partial n}{\partial t}+\vec \nabla . \big(n\vec u\big)=0$ is to look at Fig. 1 and replace $f$ with $n$. If $n$ changes along one spacial direction and there is also a velocity along this direction then $n$ changes with time. In other words, $n\vec v$ represents a flux of particles. If this flux of particles is forced inside a volume, then the number of particles inside that volume increases with time, i.e. $\frac{\partial n}{\partial t}$ is positive.

Clearly, we can repeat the whole process for ions $$\frac{\partial n_i}{\partial t}+\vec \nabla . \big(n_i\vec u_i\big)=G_i-L_i.\tag{4.14}$$

And we can also apply this method when particles are neutrals (i.e. gas) and we get $$\frac{\partial n_g}{\partial t}+\vec \nabla . \big(n_g\vec u_g\big)=G_g-L_g\tag{4.15}$$

If we look at ionizations only, we can write $$G_e=\nu_{iz}n_e$$ or $$G_e=n_en_g\sigma<v>_{ionization}$$ Note that $G_e=L_g$, $G_e=G_i$, together with $G_g=0$ (no recombinations) and $L_g=G_i$. This effectively couples all three equations together.

If ions can be ionized multiple times, then we will have a conservation equation for each charged ion species (e.g. single ionized, double ionized, ...). But all the conservation equations are still coupled via their source terms.

Conservation of momentum¶

We now look at what happens when Vlasov equations for ions and electrons are multiplied by $m_{i,e}\vec v$. Integrating over all velocities we get conservation of momentum. $$\frac{\partial n_e\vec u_e}{\partial t}+\vec \nabla \cdot (n_e\vec u_e \vec u_e)=\frac{1}{m_e}\big(-en_e[\vec E +\vec u_e \times \vec B]-\vec \nabla p_e+\vec F|_{e\,coll}\big)$$ and $$\frac{\partial n_i\vec u_i}{\partial t}+\vec \nabla \cdot (n_i\vec u_i \vec u_i)=\frac{1}{m_i}\big(Z_ien_i[\vec E +\vec u_i \times \vec B]-\vec \nabla p_i+\vec F|_{i\,coll}\big)$$

Here by $\vec u\vec v$ we mean the outer product (i.e. vector product not dot product) which is a matrix $$\begin{align}\vec u\vec v = \begin{bmatrix}u_x \\ u_y \\ u_z \end{bmatrix} \begin{bmatrix}v_x & v_y & v_z\end{bmatrix} = \begin{bmatrix}u_xv_x & u_xv_y & u_xv_z \\ u_yv_x & u_yv_y & u_yv_z \\ u_zv_x & u_zv_y & u_zv_z \end{bmatrix}.\end{align}$$ Note that $\vec \nabla \cdot \big(\vec u \vec v\big )$ is now a vector!

While the momentum equation seems more complex than mass conservation, it is in fact the same equation but where $n$ as been replaced by $n\vec u$. So one can see $mn\vec u \vec u$ is the flux of momentum $mn\vec u$ traversing a surface with velocity $\vec u$.

The operator vector $\vec F|_{coll}$ represents the force density of collisions with respect to other ions or neutrals. We will ignore such effects here. Further, we also dropped in this treatment secondary effects, such as viscosity. In this context, the pressure is given by our integral calculus and has the following form. $$p_{i,e}= \frac{1}{3}m_{i,e}n_{i,e}\Big<\big(v_{i,e}-u_{i,e}\big)^2\Big>_{\vec v}$$ where $\Big<\big(\bullet\big)^2\Big>_{\vec v}$ is the average operator over all velocities.

If we now suppose that the electron and ion plasmas behave like an ideal gas then this operator yields for both pressures, $p_e=n_eeT_e$ and $p_i=Z_in_ieT_i$.

In summary, we get conservation laws for momentum $$\frac{\partial n_e\vec u_e}{\partial t}+\vec \nabla \cdot (n_e\vec u_e \vec u_e)=\frac{1}{m_e}\big(-en_e\vec E +\vec j_e \times \vec B-\vec \nabla p_e\big)\tag{4.16}$$ and $$\frac{\partial n_i\vec u_i}{\partial t}+\vec \nabla \cdot (n_i\vec u_i \vec u_i)=\frac{1}{m_i}\big(Z_ien_i\vec E +\vec j_i \times \vec B-\vec \nabla p_i\big)\tag{4.17}$$ after ignoring collisions and using the fact that $qn_q\vec u_q=\vec j_q$, the current density for the particle with charge q.

Conservation of energy¶

We now look at what happens when Vlasov equations for ions and electrons are multiplied by $m_{i,e}v^2$ and then integrated over all velocities. This gives conservation of energy for both species. $$\frac{\partial \epsilon_e}{\partial t}+\vec \nabla \cdot (\epsilon_e\vec u_e)=\frac{1}{m_e}\big(-\vec \nabla \cdot Q_e+\vec j_e.\vec E\big)\tag{4.18}$$ Clearly, we have here a similar structure in the equation, mainly a change in time od a certain quantity $\sigma$ is related to its flux $\sigma \vec u$ via the divergence operator $\vec \nabla$. And as before we have introduced a new quantity $\vec Q_e$ that is the flux of energy at the microscopic scale, i.e. $\frac{1}{2}m_en_e<u_e^2\vec u_e>_e$, also called heat flux. The energy density source term $\vec j_e.\vec E$ is the work of the electric field on the moving electrons.

We see here that we were not able to define $\vec Q_e$ from any previous quantities. So we will need to come up with a model for this effect, purely caused by collisions. An popular approximation for the heat flux is $$\vec Q_e=\kappa_e\vec \nabla T_e\tag{4.19}$$ where $\kappa_e$ is the electron heat conductivity.

Of course, we get a similar equations for the ions $$\frac{\partial \epsilon_i}{\partial t}+\vec \nabla \cdot (\epsilon_i\vec u_i)=\frac{1}{m_i}\big(-\vec \nabla \cdot Q_i+\vec j_i.\vec E\big)\tag{4.19}$$

Summary¶

For any fluid, ions or electrons, we have shown that an explicit change in time of any moment ($n$, $n\vec u$ or $\epsilon$) can only be cause by an explicit change of flux in space ($n\vec u$, $n\vec u \vec u$ or $\epsilon \vec u$) or by explicit cources and sinks. Sources and sinks comes from collisions or from electromagnetic sources (according to Vlasov theory).

The physical meaning of the moments on the left hand side of every equation is tightly bound to the integral operator, but the moments and their flux can alsways be defined as long as an integral operator can be defined. This is not true for most terms on the right hand side. If they come from collision operators, the collision operators have to be defined physically (non trivial task when dealing with plasmas) and some quantities have to be purely and simply defined: the closure quantities (like $p$ or $\vec Q$).

$$\begin{align} &\frac{\partial n_e}{\partial t}+\vec \nabla . \big(n_e\vec u_e\big)=0\\ &\frac{\partial n_i}{\partial t}+\vec \nabla . \big(n_i\vec u_i\big)=0\\ &\frac{\partial n_e\vec u_e}{\partial t}+\vec \nabla \cdot (n_e\vec u_e \vec u_e)=\frac{1}{m_e}\big(-en_e\vec E +\vec j_e \times \vec B-\vec \nabla p_e\big)\\ &\frac{\partial n_i\vec u_i}{\partial t}+\vec \nabla \cdot (n_i\vec u_i \vec u_i)=\frac{1}{m_i}\big(Z_ien_i\vec E +\vec j_i \times \vec B-\vec \nabla p_i\big)\\ &\frac{\partial \epsilon_e}{\partial t}+\vec \nabla \cdot (\epsilon_e\vec u_e)=\frac{1}{m_e}\big(-\vec \nabla \cdot Q_e+\vec j_e.\vec E\big)\\ &\frac{\partial \epsilon_i}{\partial t}+\vec \nabla \cdot (\epsilon_i\vec u_i)=\frac{1}{m_i}\big(-\vec \nabla \cdot Q_i+\vec j_i.\vec E\big)\\ &\vec \nabla\times\vec {B} =\mu_0\vec {j}+\mu_0\varepsilon_0 \frac{\partial\vec {E}}{\partial t} \\ &\vec \nabla\times\vec {E} =-\frac{\partial\vec {B}}{\partial t} \\ &\vec \nabla\cdot\vec {E} =\frac{\rho_E}{\varepsilon_0} \\ &\vec \nabla\cdot\vec {B} =0\\ \end{align}$$

Conservation equations for a single fluid¶

The major difficulty in going from two fluid equations to equations for one single fluid is the disparity of both temporal and spacial scales between the ions and the electrons. Because they have the same charge, share similar momentum source or sink, but because the electrons are so much lighter, they react more quickly, 2000 times more quickly if we compare protons and electrons. As a results, two-fluid equations have to be solved on the electron time scale, requiring a numerical treatment that calls for very small time steps.

Because the bulk of the mass lies with the ion fluid, macroscopic effects are most visible at the ion scale. Yet the equations have to be solved very slowly, to account for electron physics, and their interactions with the ions.

To simplify the problem, we assume macroscopic ambipolarity ($\vec u_i\approx\vec u_e$) and quasi-neutrality ($n_e\approx Z_in_i$). We use here the symbol $\approx$ (i.e. almost equal) allows for the existance of an electrical current density $\vec j=Z_ie\vec u_i-en_e\vec u_e$, which can be seen as a difference in flow velocity between ions and electrons. Perfect neutrality and ambipolarity forces the electrical current to be $\vec j=\vec 0$.

Because of these two assumptions, we can redefine our electron distribution function so it takes into account many time cycles on the electron scale inside volumes compatible with the ion scale length. By doing so, we virtually have not change the left hand side of Eqs. $(4.13)$, $(4.16)$ and $(4.18)$. We have put the mathematical burden on the right hand side, i.e. on the closure quantities (e.g. $p$ or $\vec Q$ )and the collision operators. Assumptions will be made on these quantities latter. In this case, we can rewrite all our moments as $$\begin{align} &n=\sum_{i,e}n_{i,e}\\ &\rho=\sum_{i,e}n_{i,e}m_{i,e}\\ &m=\sum_{i,e}m_{i,e}\\ &n\vec u=\frac{1}{\rho}\sum_{i,e}\rho_{i,e}u_{i,e}\\ &\vec j=Z_ien_i\vec u_i-en_e\vec u_e\end{align}$$

One way to keep as much electron physics as possible without breaking the bank is to look at the impact of the electrons on the electric and magnetic fields felt by the ions. This is called Ohm's law.

Ohm's law¶

Ohm's law (or the generalized Ohm's law as we are going to study it) is simply that. It is a law that averages over electron motion at the microscopic (fast) scale to give the effective electric field felt by the ion at the macroscopic (ion) scale.

One might think that quasi-neutrality means that there is no electric field inside a plasma. In fact, it states that there is no source of electrostatic field inside a plasma, i.e. $\rho_E=Z_ien_i-en_e=0$. However, this is not because $\vec \nabla\cdot\vec E=0$ that there is no electrostatic field. It just means there exist a vector $\vec A_E$ such that $\vec E=\vec \nabla\times\vec A_E.$

The simplest case is coming from electron collisions with the ions inside materials. When averaged over many collisions, the force experience collectively by the electron is some kind of drag (fluid view). This is called resistivity in physics. If there is no external electric field, thermal motion averages out to zero. If there is an external electric field, then there is net motion in the direction of the field. But this flow is not proportional to the average force $<-e\vec E>$. This is because the ion are deflecting more eletrons than the electric field can entrain the electrons. As a result, the flow is reduced and depends on the lectric field and the type of material via the factor of proportionality $\eta$, called resistivity $$\vec E=\eta \vec J$$

Ohm was inspired by Fourier's work on heat conduction (thermodynamics or statistical approach) and use experimental measurements to deduce the linear correspondance between field and current. While $\vec j$ is a consequence of the electric field in this framework, the reverse is also true. A flowing current (generated by gravity for instance) can generate an eelctric field.

More generally, electric fields can be generated by different phenemona tightely connected to electron physics. We will look at the following in the next chapter:

1. Dynamo
2. Hall effect
3. Electron pressure (also known as the Biermann battery term)
4. Electron inertia

Conservation laws¶

While the meaning of electron closure quantities and collisional operators can be recast on the ion time scale, their mathematical definitions become very restrictive. We lose a sizable portion of electron physics by doing so. If the physical system indeed allows it, then the equations will closely match the physical behavior of this system. If not, then it will be a poor approximation. Using the one-fluid quantities as defined above, we get $$\begin{align} &\frac{\partial n}{\partial t}+\vec \nabla . \big(n\vec u\big)=0\tag{4.20}\\ &\frac{\partial n\vec u}{\partial t}+\vec \nabla . \big(n\vec u\vec u\big)=\frac{1}{m}(\vec j\times\vec B-\vec\nabla\cdot p)\tag{4.21}\\ &\frac{\partial \epsilon}{\partial t}+\vec \nabla \cdot (\epsilon\vec u)=\frac{1}{m}\big(-\vec \nabla \cdot Q+\eta j^2\big)\tag{4.22} \end{align}$$

Note that the electric field has been completely eliminated from these equations, by using Ohm's law ($\vec E=\eta \vec j$). Further, when ignoring electromagnetic radiations ($\varepsilon_0\mu_0\frac{\partial \vec E}{\partial t}\approx \vec 0$), it can be dropped from Maxwell's equations $(4.7)$ yielding

$$\begin{align} &\frac{\partial\vec {B}}{\partial t}=-\vec \nabla\times\big( {\frac{\eta}{\mu_0}}\vec \nabla\times\vec {B}\big) \tag{4.23}\\ &\vec \nabla\cdot\vec {B}=0 \tag{4.24} \end{align} $$

Hydrodynamics¶

Before solving flows where electromagnetism plays an important role, let's first look at what happens when flows are solely driven by kinetic energy, i.e. $B=0$ and $E=0$. In this case, we need to solve Eqs. $(4.20)$ to $(4.22)$ with $\vec j\times \vec B=\vec 0$. We will suppose also that all the heat is transferred via convection rather than conduction so we also ignore $\vec Q$.

Finite volume¶

To do this we recast our conservation laws using a more general framework. Regardless of their sources (even if fields are involved), Eqs. $(4.20)$ to $(4.22)$ can be written as $$\frac{\partial U}{\partial t}+\vec \nabla \cdot \vec F(U)=S_U\tag{4.25}$$ Here $U$ can be a scalar or a vector. If we discritize our domain into an ensemble of adjacent domains $\Omega_j$ with boundary $\partial\Omega_j$ we can rewrite Eq. $(4.25)$ as $$\iiint_{\Omega_j}\frac{\partial U}{\partial t}dV+\iiint_{\Omega_j}\vec \nabla \cdot \vec F(U)dV=\iiint_{\Omega_j}S_U\,dV$$

Using Gauss's theorem $$\iiint_V\left(\mathbf{\nabla}\cdot\mathbf{F}\right)\,dV=\iint_S\vec F\cdot\vec n\,dS $$ where $\vec n$ is the normal to $\partial\Omega_j$ we get

$$\iiint_{\Omega_j}\frac{\partial U}{\partial t}dV+\iint_{\partial\Omega_j}\vec F(U)\cdot \vec n\,dS=\iiint_{\Omega_j}S_U\,dV$$

By using the theorem we actually turned the spatial derivative inside a volume into a surface integral of the function itself. To first order we also suppose that $U$ and $S_U$ are constant inside the volume (this is actually why we discretized our equations) so we get

$$\frac{\partial U_j}{\partial t}+\frac{1}{V_j}\iint_{\partial\Omega_j}\vec F(U)\cdot \vec n\,dS=S_{U_j}\tag{4.26}$$

In our case, $U$ is $n$, $n\vec u$ and $\epsilon$. $\vec F(U)$ is $n\vec u$, $n\vec u\vec u$ or $\epsilon\vec u$.

Hydrodynamics code¶

We are now ready to design a code to handle hydrodynamics flows. Ignoring electric fields and magnetic fields allows us to not only simplify the code but also run the code faster, as we discuss later.

We also use numba to improve code performance (it is python after all) and some python widget to monitor the remaining time of the computation.

Keep in mind that such computation are very intensive on large computers. Using python with neither accelerator (i.e. graphics card) nor parallelization will make any computation inherently slow. So be wise when you choose your resolution.

The code presented here is two-dimensional. It can be used in slab mode, i.e. $z$ is the direction of symmetry, or in cylindrical mode, i.e. $\phi$ is the direction of symmetry. In other words, in slab geometry any quantity is such that $$U(x,y,z)=U(x,y).$$ In cylindrical geometry we have $$U(r,\phi,z)=U(r,z).$$

Python housekeeping¶

As usual we import our python libraries. Note the new addition of ipywidget that may have to be installed using the procedure discussed in Ch01.

Definitions and constants¶

For practical reason we will define some constants that will be use to remember which component in the array correspond to which quantity, e.g. the temperature can be accessed using tp at the location i,j as Q[tp,i,j]. rh is the number density, mx is the momentum density in the x-direction, my and mz in the y and z directions. en is for the energy density, tp the temperature. nQ is the total number of variables, from rh to tp.

mu is the atomic number of the fluid, in this case aluminum. aindex is the adiabatic index.

Here we define our dimensional parameters. L0 is the geometrical scale length, t0 the time scale, n0 the density scale, v0 is the velocity scale, p0 the pressure scale, te0 the temperature scale. We also define a couple of floor quantities under which no density, temperature and pressure can go.

We now define our geometrical functions. They are used quite often so we will compile them inside the code to accelerate computations. To avoid writing two codes, one for slab symmetry and one for cylindrical symmetry, we wrote this code in for slab symmetry and in Cartesian coordinates, i.e. $(x,y,z)$. When using the code for problem with cylindrical symmetry then $$\begin{align} &x\rightarrow r\\ &y\rightarrow z\\ &z\rightarrow \phi\end{align}$$

The rc function computes the radius in cylindrical coordinates, which appears in source terms and in the conservation equations. Clearly, this is simply $x$ as mentioned above. That's why rc returns xc. To turn the code into slab symmetry, just have the function rc return 1.

We now add radius and angle function for axi-symmetric distribution in slab geometries, not ot be confused with rc above. Note the the letter c, for center, is now gone.

Finite volume method and time advance algorithms¶

We now compute the sources for all our conservation equations.

Here do the time advance of the conservation equations. The flux are used here as seen in Eq. $(4.26)$. The flux are computed later.

We now write the function that computes the flux flx and the freezing speeds cfx along the x-direction. These speeds correspond to the largest speed found in the hyperbolic equation at the given location $(x,y)$. Using freezing speed avoids using a Riemann solver.

We now write the function that computes the flux fly and the freezing speeds cfy along the y-direction.

We now smooth out the different fluxes using a total diminishing variation scheme (or TVD]. Basically, this schemes applies numerical viscosity to smooth out sharp transitions caused by shocks or sharp interfaces. Typically, sharp transitions cause numerical oscillations (called Gibbs phenomenon). Note that this scheme is conservative, in the sense that it does conserve all global quantities $n, n\vec u,\epsilon$.

We truly compute all the fluxes here.

We now limit flows to avoid minor instabilities caused by round errors.

Domain and variable definitions¶

Everything you find below should be the only lines of code you need to change. We now initialize the computational domain. ngu correspond to the number of points at the boundary. Leave it alone. nx and ny are the number of grid points along the x- and y-directions. lxu and lxd are the upper and lower corners of the mesh along the x-direction. lyu and lyd are the corners for the y-direction. dxi and dyi are the inverse of the grid resolution dx and dy. Since we use both 1/dx and /dy a lot and it is more expansive to compute inverse, we compute them once and that's it. Since python is interactive, we could see them as constants here.

Every time you want to restart a new computation you need to recompile what's below.

As discussed earlier, it is important to keep the time step dt small enough to keep the code stable. For that we use a combination of flow speed and sound speed. There is a stability parameter that is called the Courant–Friedrichs–Lewy condition (or CFL). The fastest speed $U$ in the computation (the maximum of the flow speed and sound speed) defines the largest time step $\Delta t$ that can be taken as a function of the resolution $\Delta x$. $$ UC_{CFL}=\frac{\Delta x}{\Delta t}$$ $C_{CFL}<1$. Here we take it 0.1 for second order time advance. Note that higher resolution (i.e. $\Delta x$ smaller) requires to take $\Delta t$ to be smaller. So every time the resolution is increased by a factor of 2 in both x and y, the computational time is $2^3=8$ times longer.

To find the right value for your problem, start with a large number, e.g. 0.5, and decrease that number until the code runs with not error. Then decrease it even more until the final solution of the problem look the same. Then you will have the optimum value for dt. When if cfl is too small, there are a lot of time steps and sharp features in the solution will start to be washed out. So it is important to find the right cfl number. But once it is found then similar problems can be solved using the same condition.

Below we defined the initial ti and final tf times. n_out corresponds to to the total number of outputs. The rest are variable initialization. Note the most important one, the time t.

And now all our arrays are defined here. Note that we use a mask MMask for parts of the simulation space were we want to preclude motion. That will come handy later in the course.

Initial conditions¶

That's for an empty domain.

Boundary conditions¶

Below we define our boundary conditions. In this example, we do not use the time t and focus only on steady state solutions. The first block define the default boundary conditions for the whole domain. Note that the conditions for the y-direction are pretty simple. These boundary conditions correspond to an open boundary. everything can go through it.

Boundary conditions are a bit more complex for the x-direction. We are using a cylindrical coordinate system with open boundary conditions on the right but reflecting boundary conditions on the left. That's because nothing can pass through the axis of symmetry $r=0$.

At the end we have the boundary conditions that make up our problem, e.g. material free falling toward the axis of our coordinate system. Note that we are using very large density because the author works in high energy density plasmas. the initial parameters L0, t0 and n0 can be changed to accommodate other regimes. But, chances are that the CFL scaler cflm may have to be adjusted.

Computation and outputs¶

Let's define some plotting parameters

And here we go!