读《Ideal MHD》（1）-磁流体力学方程组推导

0.前言-重新出发

   Freidberg的这本《Ideal Magnetohydrodynamics》是我进入等体圈接触到的第一批书里面的一本，当时（8年前）读的时候觉得很是吃力，现在决定重新读一遍，并做好笔记。 之前开会的时候会经常遇见Freidberg， 他非常的平和，绅士，向他请教问题时他都是娓娓道来。但是早些年读这本书真的是给我留下了非常惨痛的回忆。

文中如果有表述的不对的地方，欢迎大家指出改正。

1. 动理学方程

  等离子体是大量处于电离状态的粒子组成的体系。在描述大量粒子组成的体系时，引入粒子在相空间$(\mathbf{r},\mathbf{v})$的分布函数$f(\mathbf{r},\mathbf{v},t)$。$f(\mathbf{r},\mathbf{v},t)$的含义为：$f(\mathbf{r},\mathbf{v},t)d\mathbf{r}d\mathbf{v}$代表粒子空间位置在$\mathbf{r}\sim \mathbf{r}+d\mathbf{r}$，速度在$\mathbf{v}\sim \mathbf{v}+d\mathbf{v}$的粒子数目。等离子体中一般包括多种带正点的粒子，例如氢离子，氦粒子，还有带负电的电子。假设$\alpha$类粒子的分布函数为$f_{\alpha }(\mathbf{r},\mathbf{v},t)$，则$f_{\alpha }(\mathbf{r},\mathbf{v},t)$满足以下动理学方程，即Boltzman方程
$$\frac{\partial f_{\alpha }}{\partial t}+\mathbf{v}\cdot \frac{\partial f_{\alpha }}{\partial \mathbf{r}}+\frac{{\mathbf{F}}_{\alpha }}{m_{\alpha }}\cdot \frac{\partial f_{\alpha }}{\partial \mathbf{v}}={\left(\frac{\partial f_{\alpha }}{\partial t}\right)}_c\text{\hspace{1em}(1.1)}$$
这里$m_{\alpha }$是粒子的质量，等式右边${\left(\frac{\partial f_{\alpha }}{\partial t}\right)}_c$表示与$\alpha$粒子所发生的碰撞项。

2. 有些枯燥的速度矩方程推导

  这部分内容是从微观的Boltzman方程到宏观的MHD方程的关键，在推导的过程中可以了解MHD方适用的条件，所以还是非常重要的。这部分的推导过程，在《Ideal MHD》书中没有非常详细的过程，大概是因为学物理的都知道该怎么推导吧。不过对于我这个做数学出身的小菜鸟来说，着实是费了一番功夫的, 主要参考了资料[2]。在下面的推导过程中，会省去下标$\alpha$，并且对于等离子体，Boltzman方程中的受力即为洛伦兹力
$$\mathbf{F}=\frac{q}{m}(\mathbf{E}+\mathbf{v}\times \mathbf{B})$$
式中的$\mathbf{E}，\mathbf{B}$是自洽的电场强度和磁通量密度，$q$是粒子所带的电荷强度。
  对函数$\psi (\mathbf{v})$，定义速度矩
$$⟨\psi (\mathbf{v})⟩=\int \psi (\mathbf{v})f(\mathbf{r},\mathbf{v},t)d\mathbf{v}/n(\mathbf{r},t)$$
其中$n(\mathbf{r},t)=\int f(\mathbf{r},\mathbf{v},t)d\mathbf{v}$，表示在$t$时刻，位置$\mathbf{r}$处的粒子数密度。速度矩其实就是函数在速度空间的一个平均。下面先来推导零阶，一阶，二阶，三阶速度矩方程。

• 零阶矩
  $$n(\mathbf{r},t)=\int f(\mathbf{r},\mathbf{v},t)d\mathbf{v}$$
  为粒子密度函数。质量密度$\rho (\mathbf{r},t)=nm$。

• 一阶矩
  取$\psi (\mathbf{v})=\mathbf{v}$，则得到流体的平均速度
  $$\mathbf{u}(\mathbf{r},t)=\int \mathbf{v}f(\mathbf{r},\mathbf{v},t)d\mathbf{v}/n(\mathbf{r},t)=<\mathbf{v}>$$
  定义$\mathbf{w}=\mathbf{v}-\mathbf{u}(\mathbf{r},t)$,容易知道有$<\mathbf{w}>=0$，所以$\mathbf{w}$表示无规则热运动。

• 二阶矩
  取$\psi (\mathbf{v})=\mathbf{v}\mathbf{v}$,这是一个二阶张量，一共有9个分量。
  定义
  $$\begin{array}{rl}\mathbf{P}& =nm<\mathbf{v}\mathbf{v}>=nm<(\mathbf{u}+\mathbf{w})(\mathbf{u}+\mathbf{w})>\\ & =nm<\mathbf{u}\mathbf{u}>+nm<\mathbf{w}\mathbf{w}>\\ & :=nm\mathbf{u}\mathbf{u}+\mathbf{p}\end{array}$$
  这里$\mathbf{p}$是热压强张量。根据对称性，$\mathbf{p}$只有6个独立的分量。如果体系处于局域热平衡状态，那么其分布函数为局域Maxwell分布
  $$f(\mathbf{r},\mathbf{v},t)=n(\mathbf{r},t){\left(\frac{m}{2\pi T(\mathbf{r},t)}\right)}^{3/2}\exp \left(-\frac{m{w}^2}{2T(\mathbf{r},t)}\right)$$
  利用分布函数可以计算得到
  $$p_{kk}=<w_{kk}>=n(\mathbf{r},t)T(\mathbf{r},t)=p(\mathbf{r},t)$$
  $$nm<w^2>=\sum _{k=1}^3{p}_{kk}=3nT=3p$$
  即$\mathbf{p}$的对角项是热压强。
  这里还可以得到一个有物理意义的量：体系的总动能密度
  $$K=<\frac{1}{2}nm{v}^2>=\frac{1}{2}\sum _{k=1}^3{P}_{kk}=\frac{1}{2}nm{u}^2+\frac{3}{2}p$$
  上式中的第一项为单位体积流体平均运动动能，第二项为热运动动能。
  现在我们定义粘滞应力张量
  $$\mathbf{\Pi }=\mathbf{p}-p\mathbf{I}$$
  那么有关系式
  $$\mathbf{P}=nm<\mathbf{v}\mathbf{v}>=nm<(\mathbf{u}+\mathbf{w})(\mathbf{u}+\mathbf{w})>=nm\mathbf{u}\mathbf{u}+p\mathbf{I}+\mathbf{\Pi }$$

• 三阶矩
  $\psi (\mathbf{v})=m\mathbf{v}\mathbf{v}\mathbf{v}$, 和分布函数对速度空间做积分可以得到一个三阶张量，一共有27个分量。 这里我们只关心有明确物理意义的分量
  $$\begin{array}{rl}\mathbf{Q}& =\frac{1}{2}nm<v^2\mathbf{v}>=\frac{1}{2}nm<v^2(\mathbf{u}+\mathbf{w})>\\ & =K\mathbf{u}+\frac{1}{2}nm<v^2\mathbf{w}>\\ & =K\mathbf{u}+nm\mathbf{u}\cdot <\mathbf{w}\mathbf{w}>+\frac{1}{2}nm<w^2\mathbf{w}>\end{array}$$
  如果我们定义
  $$\mathbf{q}=\frac{1}{2}nm<w^2\mathbf{w}>$$
  那么有
  $$\mathbf{Q}=K\mathbf{u}+\mathbf{u}\cdot \mathbf{p}+\mathbf{q}$$
  上式中，$K\mathbf{u}$的含义为流体宏观流动带走的总动能，$\mathbf{u}\cdot \mathbf{p}$的含义为流体宏观流动时压力张量做的功率, $\mathbf{q}$为热流矢量。 当流体宏观速度$\mathbf{u}$为零时，$K\mathbf{u}$和$\mathbf{u}\cdot \mathbf{p}$都为零，但是$\mathbf{q}$不为零，它是因为碰撞产生的热量从高温流体元到低温流体元的流动。

至此为止，我们刚刚做好了准备工作。下面要将在(1.1)式两端乘以$\psi (\mathbf{v})$，并在速度空间做积分
$$\begin{array}{r}\int \psi (\mathbf{v})\frac{\partial f}{\partial t}d\mathbf{v}+\int \psi (\mathbf{v})\mathbf{v}\cdot \frac{\partial f}{\partial \mathbf{r}}d\mathbf{v}+\int \psi (\mathbf{v})\frac{q}{m}(\mathbf{E}+\mathbf{v}\times \mathbf{B})\cdot \frac{\partial f}{\partial \mathbf{v}}d\mathbf{v}=\int \psi (\mathbf{v}){\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}\end{array}$$
这里我们一项一项来分析，第一项第二项比较简单，这里需要注意的是我们是在广义坐标$(\mathbf{r},\mathbf{v},t)$下来考虑的，所以
$$\int \psi (\mathbf{v})\frac{\partial f}{\partial t}d\mathbf{v}=\frac{\partial }{\partial t}\int \psi (\mathbf{v})fd\mathbf{v}=\frac{\partial }{\partial t}(n<\psi (\mathbf{v})>)$$
在分析第二项的时候，这里将$\frac{\partial }{\partial \mathbf{r}}$换成$\nabla$，两个符号是等价的
$$\int \psi (\mathbf{v})\mathbf{v}\cdot \nabla fd\mathbf{v}=\nabla \cdot \int \psi (\mathbf{v})\mathbf{v}fd\mathbf{v}=\nabla \cdot (n<\psi (\mathbf{v})\mathbf{v}>)$$
第三项因为是对$\mathbf{v}$求微分，所以处理起来稍微要仔细一些，这里需要用到张量的运算公式$\nabla \cdot (\mathbf{A}\mathbf{B})=(\nabla \cdot \mathbf{A})\mathbf{B}+(\mathbf{A}\cdot \nabla )\mathbf{B}$，$\nabla \cdot (f\mathbf{A})=\nabla f\cdot \mathbf{A}+f\nabla \cdot \mathbf{A}$。为了记号简便，下面的公式中$\frac{\partial }{\partial \mathbf{v}}$会用${\nabla }_{\mathbf{v}}$来替代。所以第三项为(注意这里推导的时候没有因子$\frac{q}{m}$)
∫ ψ ( v ) ( E + v × B ) ⋅ ∇ v f d v = ∫ ψ ( v ) { ∇ v ⋅ [ ( E + v × B ) f ] − ∇ v ⋅ ( E + v × B ) f } = ∫ ψ ( v ) { ∇ v ⋅ [ ( E + v × B ) f ] } = ∫ ∇ v ⋅ [ ψ ( v ) ( E + v × B ) f ] − ( E + v × B ) f ⋅ ∇ v ψ ( v ) d v = [ ψ ( v ) ( E + v × B ) f ] 边 界 − ∫ ( E + v × B ) f ⋅ ∇ v ψ ( v ) d v \begin{aligned} &\int\psi(\mathbf{v})(\mathbf{E}+\mathbf{v}\times\mathbf{B})\cdot\nabla_{\mathbf{v}}fd\mathbf{v} \\&=\int\psi(\mathbf{v}) \left\{\nabla_{\mathbf{v}}\cdot\left[(\mathbf{E+v\times B})f\right]- \nabla_{\mathbf{v}}\cdot(\mathbf{E+v\times B})f\right\} \\&=\int\psi(\mathbf{v}) \left\{\nabla_{\mathbf{v}}\cdot\left[(\mathbf{E+v\times B})f \right] \right\} \\&=\int \nabla_{\mathbf{v}}\cdot\left[\psi(\mathbf{v})(\mathbf{E+v\times B})f\right]-(\mathbf{E+v\times B})f\cdot\nabla_{\mathbf{v}}\psi(\mathbf{v}) d\mathbf{v} \\&=\left[\psi(\mathbf{v})(\mathbf{E+v\times B})f\right]_{边界}-\int (\mathbf{E+v\times B})f\cdot\nabla_{\mathbf{v}}\psi(\mathbf{v}) d\mathbf{v} \end{aligned} ​∫ψ(v)(E+v×B)⋅∇v​fdv=∫ψ(v){∇v​⋅[(E+v×B)f]−∇v​⋅(E+v×B)f}=∫ψ(v){∇v​⋅[(E+v×B)f]}=∫∇v​⋅[ψ(v)(E+v×B)f]−(E+v×B)f⋅∇v​ψ(v)dv=[ψ(v)(E+v×B)f]边界​−∫(E+v×B)f⋅∇v​ψ(v)dv​
上的推导中，由于$\mathbf{E}=\mathbf{E}(\mathbf{r},t)$, 所以${\nabla }_{\mathbf{v}}\cdot \mathbf{E}=0$, 并且${\nabla }_{\mathbf{v}}\cdot (\mathbf{v}\times \mathbf{B})=0$。因为当$\mathbf{v}\to \infty$时，$f(\mathbf{v})\to 0$，所以上式中的边界项为零。综合以上的推导，可以得到最后的速度矩方程

$$\begin{array}{rl}\frac{\partial }{\partial t}(n<\psi (\mathbf{v})>)& +\nabla \cdot (n<\psi (\mathbf{v})\mathbf{v}>)-\frac{nq}{m}\mathbf{E}\cdot <{\nabla }_{\mathbf{v}}\psi (\mathbf{v})>-\frac{nq}{m}<(\mathbf{v}\times \mathbf{B})\cdot {\nabla }_{\mathbf{v}}\psi (\mathbf{v})>\\ & =\int \psi (\mathbf{v}){\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}\end{array}\text{\hspace{1em}(2.1)}$$

3. 终于要看到磁流体方程组了

   下面我们要从(2.1)式出发，导出流体方程组了。 在这之前，需要假设：等离子体不存在粒子（电子，离子）之间的电离，复合等情况。也就是说（2.1）式中的碰撞项只包含弹性碰撞。等离子是由处于电离状态的电子、离子组成的粒子体系，如果其中的电子和离子没有达到平衡，那么电子和离子应作为两种不同的粒子体系，这样将得到双流体方程。如果电子和离子达到平衡，那么可以将等离子看做是电中性的流体，此时将得到单流体方程组。

3.1 双流体方程组

   接下来还是慢慢地来推导, 我们先来推导某种粒子的方程，然后再考虑不同的粒子。 首先在(2.1)式中，取$\psi (\mathbf{v})=1$, 结合之前准备工作推导的各阶矩，有
$$\frac{\partial n}{\partial t}+\nabla \cdot (n\mathbf{u})=\int {\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}$$
分析下右端的碰撞项，${\left(\frac{\partial f}{\partial t}\right)}_c$表示粒子由于碰撞（根据假设，只考虑弹性碰撞）导致分布函数的变化。$\int {\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}$的含义就是$t$时刻，在空间位置$\mathbf{r}$，粒子数密度由于弹性碰撞的净增量。因为只考虑弹性碰撞（没有电离，复合等过程），所以这一项为零。从而得到了连续性方程
$$\frac{\partial n}{\partial t}+\nabla \cdot (n\mathbf{u})=0$$
上式两端乘以粒子的质量$m$，则可以得到质量守恒方程
$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0$$
这里用到了关系式$\rho =nm$。
   接下来，在(2.1)式中，取$\psi (\mathbf{v})=m\mathbf{v}$, 那么可以得到

$$\begin{array}{rl}\frac{\partial }{\partial t}(nm\mathbf{u})& +\nabla \cdot (nm<\mathbf{v}\mathbf{v}>)-nq(\mathbf{E}+\mathbf{u}\times \mathbf{B})\\ & =\int m\mathbf{v}{\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}\end{array}$$
上式中，记$\mathbf{F}=q\mathbf{E}+\mathbf{u}\times \mathbf{B}$，为洛伦兹力。$nm<\mathbf{v}\mathbf{v}>=nm\mathbf{u}\mathbf{u}+p\mathbf{I}+\mathbf{\Pi }$。 这里还需要分析一下右端的碰撞项
$$\begin{array}{rl}\int m\mathbf{v}{\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}& =\int m(\mathbf{u}+\mathbf{w}){\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}\\ & =\int m\mathbf{w}{\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}:=\mathbf{R}\end{array}$$
这一项表示由粒子做无规则热运动时碰撞所导致的动量的改变，在宏观体现为摩擦力。综上，可以得到动量方程
$$\begin{array}{r}nm\frac{\partial \mathbf{u}}{\partial t}+nm\nabla \cdot (\mathbf{u}\mathbf{u})=n\mathbf{F}-\nabla p-\nabla \cdot \mathbf{\Pi }+\mathbf{R}\end{array}$$
引入随体导数
$$\frac{d}{dt}\equiv \frac{\partial }{\partial t}+\mathbf{u}\cdot \nabla$$
动量方程可以表达为
$$\begin{array}{r}nm\frac{d\mathbf{u}}{dt}=n\mathbf{F}-\nabla p-\nabla \cdot \mathbf{\Pi }+\mathbf{R}\end{array}$$
   快接近尾声了，再接下来，在(2.1)式中，取$\psi (\mathbf{v})=\frac{1}{2}m{v}^2$, 那么可以得到
$$\begin{array}{rl}\frac{\partial }{\partial t}(n<\frac{1}{2}m{v}^2>)& +\nabla \cdot (n<\frac{1}{2}m{v}^2\mathbf{v}>)-\frac{nq}{m}\mathbf{E}\cdot <m\mathbf{v}>-\frac{nq}{m}<(\mathbf{v}\times \mathbf{B})\cdot m\mathbf{v}>\\ & =\int \frac{1}{2}m{v}^2{\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}\end{array}$$
之前我们在计算二阶矩的时候，定义了$K=<\frac{1}{2}nm{v}^2>$，计算三阶矩的时候有$\mathbf{Q}=<\frac{1}{2}nm{v}^2\mathbf{v}>$。同样这里还需要分析下右端的碰撞项
$$\begin{array}{rl}\int \frac{1}{2}m{v}^2{\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}& =\int \frac{1}{2}m(u^2+2\mathbf{u}\cdot \mathbf{w}+w^2){\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}\\ & =\mathbf{u}\cdot \mathbf{R}+Q\end{array}$$
上式中，定义了$Q=\int \frac{1}{2}m{w}^2{\left(\frac{\partial f}{\partial t}\right)}_{c}d\mathbf{v}$，表示粒子的弹性碰撞引起的热能交换。综上得到了能量方程
$$\begin{array}{r}\frac{\partial K}{\partial t}=-\nabla \cdot \mathbf{Q}+nq\mathbf{E}\cdot \mathbf{u}+\mathbf{u}\cdot \mathbf{R}+Q\end{array}$$
结合质量守恒方程，动量方程，在等离子体处于局部热平衡假设的前提下，上式可以简化为关于$T$的方程，
$$\frac{3}{2}n\frac{\partial T}{\partial t}+\mathbf{p}\cdot \nabla \cdot \mathbf{u}+\nabla \cdot \mathbf{q}=Q$$
上面推导的式子中加上粒子种类的标识$\alpha$，再结合Maxwell方程组，至此，我们终于得到了双流体方程组:
$$\begin{array}{rl}& \frac{\partial n_{\alpha }}{\partial t}+\nabla \cdot (n_{\alpha }{\mathbf{u}}_{\alpha })=0\\ & n_{\alpha }{m}_{\alpha }\frac{d{\mathbf{u}}_{\alpha }}{dt}=n{\mathbf{F}}_{\alpha }-\nabla p_{\alpha }-\nabla \cdot {\mathbf{\Pi }}_{\alpha }+{\mathbf{R}}_{\alpha }\\ & \frac{3}{2}n\frac{d{T}_{\alpha }}{dt}+{\mathbf{p}}_{\alpha }\cdot \nabla \cdot {\mathbf{u}}_{\alpha }+\nabla \cdot {\mathbf{q}}_{\alpha }=Q_{\alpha }\\ & \nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t}\\ & \nabla \times \mathbf{B}={\mu }_0\sum n_{\alpha }{q}_{\alpha }{\mathbf{u}}_{\alpha }+\frac{\partial \mathbf{E}}{\partial t}\\ & \nabla \cdot \mathbf{E}=\frac{1}{{ϵ}_0}\sum n_{\alpha }{q}_{\alpha }\\ & \nabla \cdot \mathbf{B}=0\end{array}$$
不知道是不是篇幅有点长，编辑文章的时候实在是很卡。所以理想磁流体方程的推导就放到下篇。

===================================================================
参考资料：
[1] (最重要的当然是读的这本书啦)Jeffrey P. Freidberg, Ideal Magnetohydrodynamics
[2]https://max.book118.com/html/2018/0512/166044967.shtm