基于Fokker-Planck方程的等离子体模拟ALI SHAJII，DANIEL SMITHMKS Instruments, 90 Industrial Way, Wilmington, MA, USA, 01887 对于各个计算组织，等离子体的模拟一直是个极大的挑战，有很多不同近似程度的模拟计算方法。

包括完整的动力学计算方法，流体近似方法和关于漂移扩散方程的方法。

近几年来，有人用Fokker-Planck方程处理等离子体中的电子，同时把离子当作流体进行耦合计算，获得了很好的计算结果。

本章我们将介绍基于通用Fokker-Planck方程的计算求解过程，并通过一个具体实例得到电容放电过程的电子密度分布。

希望通过该简单模型使读者对等离子放电建模过程有个初步的了解。

1．引言各种工业等离子体应用“过程”中都存在一个关键步骤[2][3]。

历史上曾采用各种不同方法对等离子体进行简化建模，分别对应于不同层面问题所需准确性[3][5][7]。

这些层面包括：●完整的动力学模型（多组分Boltzmann方程）[4]；●使用Monte-Carlo方法的颗粒模拟[3]；●Fokker-Planck近似[1][7]；●多尺度流动模型（也被称作漂移扩散模型）[3]。

出于种种原因，使得等离子体的建模和模拟非常困难。

首先，最直接的使用多流体方程的模型不能反应相关的等离子体物理过程。

其次，“水动力学”系数完全取决于研究的特定问题，不能作为纯气体或液体的常数简单测量。

最重要的一点是，完整的动力学模型包括Boltzmann方程，计算求解非常困难。

对于完整动力学模型和流动模型之间的需求空白，通常采用Fokker-Planck (FP)近似或者Monte Carlo (MC)颗粒模拟。

这两种方法可以在所需计算复杂度和捕获等离子体重要物理细节之间找到一个很好的平衡。

本章的主要目的是展现用COMSOL Multiphysics求解FP方程的功能。

为了对该问题给出一个整体认识，我们把侧重点集中在一个简单的例子上。

特别是在第二节我们通过一个简单的例子对FP 方程给出一个直观描述，将它用于布朗运动的颗粒模拟。

本章最主要的贡献在于介绍了如何在外部电场情况下对电子动力学过程进行建模。

最后，在第四节中对如何使用COMSOL Multiphysics 实现该模型给出了详细的讨论。

2．一维FP 和Langevin 方程对FP 方程直接求导可以得到Langevin 方程[1]。

考虑浸没在流体中的“布朗”颗粒，如果颗粒足够小，会同时受到两种力。

一是颗粒和流体介质间的粘性力，它会降低平均颗粒速度。

二是颗粒与流体“分子”间随机碰撞的力。

该布朗颗粒的运动控制方程如下[1]：()t υγυ=-+Γ (1)其中γ是阻尼系数，随机项()t Γ表示颗粒与背景流体的连续碰撞。

对于这个简单的例子，通常我们假设粘性力线性依赖于颗粒速度。

同时，根据随机近似，Langevin 力必须满足以下方程[1]：()0()(')(')t t t q t t δΓ=ΓΓ=-其中2/q kT m γ= [1]，表示整体平均，T 是流体温度，k 是Boltzmann 常数，m 是布朗颗粒的质量。

同时注意到力Γ可以很容易的通过MATLAB 函数randn （见第四节）实现。

给定适当的颗粒数和初始条件后，使用MC 方法可以很容易的求解方程(1)。

为了对颗粒在时间t 时的速度进行精确静态测量，颗粒数通常需要超过一百万。

最简单的初始条件为0(0)t υυ==。

总之，Langevin 方程描述了背景流体介质中布朗颗粒的运动。

当然，如果没有随机力，颗粒的路径为0t e γυυ-=。

但是由于Γ的存在，方程(1)对于很多颗粒算出的解通常是条分布曲线，分布的宽度由q 决定[1]。

FP 方程为MC 方法提供了另一种思路。

包含随机力的Langevin 方程可以等价与以下偏微分方程：2w q w w w v t v v v γ∂∂∂∂⎛⎫⎛⎫=++ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭(2)初始条件和边界条件为：0(0)()()0w t v v w v δ==-→±∞=(3)其中w 是颗粒在速度空间中的分布概率，也就是颗粒出现的概率，在无限大系统中，速度通常表示为wdv 。

初始条件(3)如图1所示。

求解该偏微分方程等价于对无限个颗粒进行Monte-Carlo 模拟，以获得t >0时的概率分布函数。

图2给出了两种方法的比较。

MC 方法针对初始条件为00υ=的情况，求解了20000个常微分方程。

这还不是强非线性阻尼力和三维空间情况下的颗粒运动，即使如此，对于这个简单的例子，求解偏微分方程的时间和常微分方程一样。

当颗粒痕迹的常微分方程是非线性时，或者需要从FEA模拟中查询外部力的值时，MC方法的计算时间就会变得非常大。

图1 初始条件w图2 各时间步长下Fokker-Planck 和Monte Carlo 方法的对比作为例子，我们考虑修改后的Langevin 方程：3()t υγυ=-+Γ(4)和相应的FP 方程[1]：121D w w w D w D t v v vv γ∂∂∂∂∂⎛⎫⎛⎫=++ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭ (5)其中31D v γ=-，2/2D q =。

这种情况下的分布函数如图3所示。

即使计算一百万个颗粒，常微分方程的分布函数也存在很大的数值噪声，而FP 结果就很平滑。

这里出现的静态噪声是有限颗粒MC 方法模拟的经典问题（也就是说数据的后处理是很精细、很关键的步骤）。

图3 各时间步长下基于非线性Langevin方程的Fokker-Planck和Monte Carlo方法的对比下面，根据本节介绍的背景知识，我们来看一个基于FP方程的更为复杂的例子。

3．外加场中的电子动力学现在我们考虑电子流在电容两极板间外加电场作用下的动力学过程。

一个极板保持常电压，另一个极板接地。

假设电子流从接地极板侧以一定速度引入，并在阳极侧积累。

并且两极板间的整个区域充满了温度为T的静止流体。

图4给出了该几何结构。

图4 模拟域示意图模型主要假设为：● 忽略由于电子和背景介质碰撞所产生的离子； ● 电子间的相互作用可以忽略； ● 在一维情况下求解该问题。

为了精确求解物理空间（y ）和速度空间（v y ）的电子分布，需要求解FP 方程。

在本例中我们还要计算电子密度、电子流和电容两极板间的静电场。

在通常情况下，FP 方程比静电场问题多需要一个空间维数，所以我们需要同时求解一个一维静电场和二维FP 方程。

电子分布的控制方程为：22yy y y y y F kT ww w w v w v m v y mv v γγ⎛⎫⎛⎫∂∂∂∂-=--++ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭ (6)对于该电场ρε∇⋅=E (7)用标量形式电势和电子密度表示为22()e en y V y ε∂=∂ (8)其中()e y n y wdv ∞-∞=⎰(9)并且()y VF y ey∂=∂ (10)注意这里的n e 和F y 都是y 的函数。

边界条件和初始条件为：()()00,y y w y v v v δ==-(11) (),0y w y v =±∞=(12) ()00V y == (13) ()10.3V y ==(14)其中m 是电子质量，k 是Boltzmann 常数，T 是背景流体温度，e 是电子电荷，0ε是自由空间介电常数，γ是阻尼系数。

同时，无限大系统中，在速度v y 和位置y 处找到一个电子的概率为wdvdy 。

方程(6)和(8)是两个耦合的不同空间维数的偏微分方程。

由于COMSOL Multiphysics 允许不同空间维数的任意数量几何形状进行耦合，所以我们现在同时求解这两个方程。

通过投影耦合变量方法将电子密度引入静电场问题中，通过导出耦合变量的方法将颗粒上的力引入FP 问题中。

概括起来讲，电子流以速度v 0进入“电容放电”区域，同时受到三种力：Γ，v γ-和V ∇。

这些力导致电子速度和密度分布沿y 方向变化，所以与之前简单线性分布例子相比，当修改极板电势V 时，不存在这种情况。

图5给出了函数w 的解。

注意到初始时刻，当电子进入计算域时速度降低。

这种非线性、非单调曲线是由如图6所示的阴极电子屏蔽所导致的。

查看域中间位置附近的最大电子密度。

图8表明电容的电子电势在最大n e 处有最小值。

图5 概率分布图6 电子密度图7 平均速度图8 平行极板间电子电势随位置的分布最后，给出电子流量：y y wv dv ∞-∞Φ=⎰由于质量守恒，流量在y 方向应该保持为常数。

图7沿y 方向的电子流量表明该过程中质量守恒。

在下一节中我们将详细给出在COMSOL 中如何实现该模型。

但是在开始之前应该注意到，虽然以上模型不是对等离子体的完整模拟，但是通过模拟它可以给出了很多重要的特性，即，一个更完整的模型仍然应该包括FP 方程，但是现在必须耦合离子动力学过程。

同时，在很多情况下，FP 方程都需要在二维空间中求解。

FP 方程的速度依赖性采用相同的处理方法[7]。

4．COMSOL Multiphysics 的实现该模型用到COMSOL Multiphysics 的一些高级功能。

这些问题包括：● 速度空间的拓扑域是无限的，而物理空间只有厘米量级。

这会导致几何结构具有极高的横纵比；● 两个偏微分方程的空间维数不同；● 必须对所有的y 在速度空间中进行一系列线性积分，然后把得到的值用于静电场问题；● 同样的，静电场问题产生的力必须施加到Fokker-Planck 域物理和速度空间的每一个点；● 两个物理场之间存在强耦合，即：电子密度是外部施加力的强函数，且外部施加力是电子密度的强函数。

以上每一点都会带来潜在的问题，但是COMSOL Multiphysics 可以通过一些方法来克服这些困难。

第一点通过调节偏微分方程(6)和(8)的尺度来克服。

这使得几何结构具有较低的横纵比，从而提高网格质量。

尺度调节后的偏微分方程为：22T y y y C y y v w w dV w wv w v v v y dy v v α⎛⎫⎛⎫∂∂∂∂-=--++ ⎪ ⎪∂∂∂∂-⎝⎭⎝⎭(15) 2y V wdv β∞-∞∇=-⎰(16)其中T v =C v d γ=T eV mv dαγ=2000ed n V βε=m 是电子质量，k 是Boltzmann 常数，T 是背景流体温度，e 是电子电荷，0ε是自由空间介电常数，d是平行极板间的距离，V0是下极板电压，γ是阻尼常数，n0是y＝0处的电子密度。

我们通过模型导航栏建立两个几何结构，几何结构1为二维结构，独立变量y，v y和z，偏微分方程通用形式的因变量为w。

然后再增加一个几何结构2，使用独立变量x，y和z，偏微分方程通用形式的因变量为V。