[笔记] 分布函数, Vlasov方程和PIC

[latexpage]

Basic notions

$df/dt$: 全微分,只随时间变化,没有空间坐标依赖关系 $\Rightarrow$ 描述的是随体运动坐标系中的物理量的性质(Lagrangian坐标系)。

转换到时间$t$-固定空间$x$坐标系(Eulerian坐标系)中描述,流体元$f(x,t)$的时空演化: $df/dt = \partial f/\partial t + \partial f/\partial x \cdot dx/dt = \partial f/\partial t + u_x \cdot \partial f/\partial x $. 其中,$\partial f / \partial t$是固定坐标系中的时间变化。

分布函数$f(x,v,t) \Rightarrow df/dt = \partial f/\partial t + u_x \cdot \partial f/\partial x + F/m \cdot \partial f/ \partial v $.

如果:1)$df/dt=0$,即无碰撞;2)$F=q/m(E+v\times B)$,即作用力为电磁力, 则$df/dt=0$这个方程叫做Vlasov Eq.

这里的记号$x,v$均应看作三维空间的矢量,每个矢量都有沿三个正交坐标方向的分量,如果采用笛卡尔坐标($e_x,e_y,e_z$),则$x=(x,y,z)$和$v=(v_x,v_y,v_z)$. Vlasov codes就是在Eulerian坐标系上的($x,v$)6D相空间中用有限差分方法来数值求解Vlasov Eq.这个偏微分方程,需要计算和存储6D网格上每个格点的数据及其随时间的演化,计算量很大。

The meaning of distribution function

$f(x,v,t)$代表在时间$t$、位置$x$、速度为$v$的粒子数目。更严谨的说法是$fdxdv$代表$t$时刻在$x$和$x+dx$之间以及$v
$和$v+dv$之间的相空间体积元中的粒子数。它是“相空间”坐标($x,v$)和时间$t$的函数。如果默认考虑单位体积,即$f \equiv f/Vol$,$Vol$是整个体系的体积,则$f$代表粒子密度。如果$f$是归一化的,即$\int fdv=1$,则$f$代表几率。

对应通常密度$n=n(x,t)$来理解:如果对速度积分,即考虑在时间$t$和位置$x$包含了所有速度的(单位体积)粒子数目,则$\int_{-\infty}^{\infty}f(x,v,t)dv=n(x,t)$,这就是通常所说的密度$n(x,t)$,它是四个标量变量即位置$x$和时间$t$的函数,因为这个密度是通常坐标空间单位体积的粒子数,只随位置$x$和时间$t$发生变化。

本来完全描述一个N个粒子-体系的分布函数的形式应该是$F(x_1,v_1,\cdots,x_N,v_N,t)$,相空间是6N维。但是根据统计物理得出的结论,对于每个种类(电子或离子)的粒子来说,用一个单粒子分布函数就可以很好地描述这类粒子体系的性质,相空间只有6维,前提是粒子间的关联(碰撞)效应很弱,可以用微扰处理。

Liouville theorem

如前所述,全微分$df/dt$就解释为随体(相空间流体元)运动的坐标系中$f$随时间的变化率。那么为啥除了有碰撞外,相空间密度$f$不随时间变化,即有Liouville定理:$df/dt=0$ ?

我的理解是:因为在随体坐标系看来,$f$只是$t$的函数了,但是$f(t)$表示在时间$t$的粒子数目,如果没有碰撞或热运动,流体元中的所有粒子一起运动,流体元中的粒子数目不随时间变化;如果考虑单位体积,就是说相空间密度不随时间变化。如果有碰撞或热运动,这个随体的“体”就会逐渐瓦解,不能再代表一团流体元整体一起运动,流体元中的粒子数目$f$会随时间$t$发生变化,即此时$df/dt\neq 0$.

固定($x,v$)坐标系来看,相空间分布(即$f$)会随时间变化,这只是说$\partial f / \partial t \neq 0$. 而说无碰撞时分布函数不随时间变化,说的是$f$的全微分为零。

更严格深刻的理解,可以参考B. V. Somov的Plasma Astrophysics, Part I: Fundamentals and Particle一书,开篇就是Liouville定理,并给出了数学推导及守恒条件:

在相空间运动的粒子应遵守连续性方程,
\begin{equation}
\partial f / \partial t + \nabla \cdot (f\dot{X}) = 0,
\end{equation}
其中$X={x,v}$. 可以类比电荷连续性方程$\partial \rho / \partial t + \nabla \cdot j = 0$(其中$j=\rho v$)来理解:$f \rightarrow \rho$. 用${x,v}$将连续性方程的$X$具体写出:
\begin{equation}
\partial f / \partial t + \nabla_x \cdot (f \dot{x}) + \nabla_v \cdot (f \dot{v}) = 0.
\end{equation}
其中,左边第二项:
\begin{equation}
\nabla_x \cdot (f \dot{x}) = \nabla_x \cdot (f v) = f \nabla_x \cdot v + (v\cdot \nabla_x)f = 0 + (v\cdot \nabla_x)f
\end{equation}
上式最后那个0是因为在相空间中$x$和$v$是独立变量。再看上上式左边第三项:
\begin{equation}
\nabla_v \cdot (f \dot{v}) = f \nabla_v \cdot \dot{v} + \dot{v} \cdot \nabla_v f.
\end{equation}
所以连续性方程现在可写成
\begin{equation}
\partial f / \partial t + (v\cdot \nabla_x)f + f \nabla_v \cdot \dot{v} + \dot{v} \cdot \nabla_v f =0.
\end{equation}
到目前为止都还是严格的,没有任何假设和近似。

假如,现在开始我们只考虑上式左边第三项为零的情形,即:
\begin{equation}\label{fcondition}
f \nabla_v \cdot \dot{v} = 0.
\end{equation}
这样连续性方程变成
\begin{equation}
\partial f / \partial t + (v\cdot \nabla_x)f + \dot{v} \cdot \nabla_v f =0.
\end{equation}
这就是Liouville方程,而且$\partial / \partial t + v\cdot \nabla_x + \dot{v} \cdot \nabla_v $ 刚好就是全微分$d / dt$.

所以,相空间密度守恒即Liouville定理成立的关键就是满足条件(\ref{fcondition}),也就是说要求\textbf{某个方向的加速度或者力,与对应方向的速度无关}。可以验证,Lorentz力和重力满足这个条件,对相空间密度守恒来说是保守力;而摩擦力与碰撞力不满足,因为它们与速度有关,换言之,对相空间密度守恒来说是耗散力。\footnote{Mar. 27, 2017 Mon}

连续性方程是从物理出发得到的最根本性的方程,它不含任何假设。由它可以推导出Liouville和Vlasov方程。

Distribution function and particle-in-cell simulation

现在有了$df/dt=0$,即$f$在无碰撞时为常数,这意味着粒子在相空间中是沿等$f$线运动的,也就是说$f$的等值线就是粒子在相空间($x,v$)运动的轨道。如果考虑2D相空间($x,v_x$),$f$的等值线是一条线,e.g.,所有粒子速度都为$v_0$,或产生一个沿$x$传播的波动,速度分布为一条正弦曲线$v_0\sin x$。考虑6D相空间,则$f$张成一个3D的超曲面,粒子轨道沿这个超曲面。随着系统演化,这个等值曲面可能发生拉伸、收缩、变形,但始终是3D的。用一定数量的Finite Phase-Fluid Elements (FPFE),即一定数量的在相空间中的有限大小流体元(或称为宏粒子),来近似(approximate)或取样(sample)这个等$f$曲面,再跟踪它们的分布演化,就可得到整个粒子系统的演化信息,所以在相空间中的FPFE方法与通常说的PIC方法等价,不同的是前者中的宏粒子是相空间中的有限流体元,而后者中的宏粒子代表一大堆德拜球内的真实粒子。

采用FPFE方法的PIC的好处:

1)和Vlasov code比,大大减小了计算量,因为只用sample和follow有粒子分布或有感兴趣的物理过程发生的3D等$f$超曲面区域,而不用计算和存储在6D相空间中即使为零或不感兴趣的每个格点,从这个意义来说,这样的PIC codes可看作“压缩”或“Lagrangian”的Vlasov codes,而且用较少的宏粒子FPFE就能准确地取样超曲面;

2)很容易推广到one (macroparticle) to one (real particle) simulation的情形,只需把粒子间的碰撞力加上就行。

但如何保证这个方程是沿着3D等$f$超曲面演化的呢?

等$f$意味着没有加热(热运动造成相空间的diffusion),因此保证系统总能量守恒就保证了等$f$的条件。系统总能量由粒子动能、静电能、电磁能、磁场能(Lorentz力自动保证能量守恒)组成。如果运动方程满足正则方程,则总能量即总Hamiltonian守恒。第$n$个FPFE(宏粒子)的正则运动方程为
\begin{equation}
dx_n/dt = p/m\gamma,
\end{equation}
\begin{equation}
dp_n/dt=F(+F_{collision}).
\end{equation}
再通过特定的电流分配方法,可保证电荷守恒,即保证满足连续性方程$\partial \rho / \partial t + \nabla \cdot j = 0$,那么电磁场变化的Maxwell方程组
\begin{equation}
\partial E / \partial t = c\nabla \times B ## 4\pi j,
\end{equation}
\begin{equation}
\partial B / \partial t = -c \nabla \times E
\end{equation}
也满足能量守恒。这样得到的系统总能量守恒。

current interpolation to the grid: 粒子有一定的中心位置,当粒子运动,沿运动路径积分得到电流密度,通过保证电荷守恒的方法,将电流密度分配到网格上。然后计算网格上的电场。(电荷密度是定义在网格中心的,通过形状因子将粒子位置分配到中心)

force interpolation to the actual particle position: 将网格上的电场分配到当前粒子位置,其分配方法和电流密度被分配到网格上的方法相同。通过正则运动方程由场推动粒子得到动量。\footnote{Mar. 21, 2017}

Why no collisions in PIC

描述N个粒子的分布函数$ F^N(\textbf{x}_1,\cdots,\textbf{x}_N,\textbf{p}_1,\cdots,\textbf{p}_N) $是6N维的。

\textbf{如果}满足粒子间碰撞(inter-particle correlations/collisions)效应很弱(可微扰处理),则:

1)用单种粒子的分布函数$ f(\textbf{x},\textbf{p}) $来描述系统就足够了,只有6维;

2)粒子的运动满足Boltzmann-Vlasov Eq.: $ \partial f /\partial t + (P/\gamma m) \nabla f + (F/m) \nabla_p f = s.t. $,其中s.t.是粒子间碰撞项。

而sample the distribution function by a set of Finite Phase-Fluid Elements (FPFEs): $ f(\textbf{x},\textbf{p}) = \sum_n WEIGHT_{n}^{ph} SHAPE^{ph}(\textbf{x}-\textbf{x}n,\textbf{p}-\textbf{p}_n); \; d \textbf{x}_n / d t = \textbf{p}/ m \gamma; \; d\textbf{p}_n /d t = \textbf{F}{ext} + \textbf{F}_{St} $,这是PIC中采用的相对论粒子运动方程,它比Boltzmann-Vlasov Eq.还来得基本,因为Boltzmann-Vlasov Eq.的前提是粒子间关联或碰撞足够小,而FPFEs可以推广到处理包含了碰撞在内的one(numerical particle)-to-one(real particle) simulation!

为啥通常PIC里面没有包含碰撞呢?因为只考虑了外加电磁场的力$ \textbf{F}{external} $这一项,而忽略了$ \textbf{F}{St} $这一项,它是带电的粒子与粒子之间的静电库仑力,即inter-particle correlations/collisions.

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注