基于Fokker-Planck方程的等离子体模拟.doc

1、 基于Fokker-Planck方程的等离子体模拟ALI SHAJII，DANIEL SMITHMKS Instruments, 90 Industrial Way, Wilmington, MA, USA, 01887对于各个计算组织，等离子体的模拟一直是个极大的挑战，有很多不同近似程度的模拟计算方法。包括完整的动力学计算方法，流体近似方法和关于漂移扩散方程的方法。近几年来，有人用Fokker-Planck方程处理等离子体中的电子，同时把离子当作流体进行耦合计算，获得了很好的计算结果。本章我们将介绍基于通用Fokker-Planck方程的计算求解过程，并通过一个具体实例得到电容放电过程的电子

2、密度分布。希望通过该简单模型使读者对等离子放电建模过程有个初步的了解。1引言各种工业等离子体应用“过程”中都存在一个关键步骤23。历史上曾采用各种不同方法对等离子体进行简化建模，分别对应于不同层面问题所需准确性357。这些层面包括：l 完整的动力学模型（多组分Boltzmann方程）4；l 使用Monte-Carlo方法的颗粒模拟3；l Fokker-Planck近似17；l 多尺度流动模型（也被称作漂移扩散模型）3。出于种种原因，使得等离子体的建模和模拟非常困难。首先，最直接的使用多流体方程的模型不能反应相关的等离子体物理过程。其次，“水动力学”系数完全取决于研究的特定问题，不能作为纯气体或

3、液体的常数简单测量。最重要的一点是，完整的动力学模型包括Boltzmann方程，计算求解非常困难。对于完整动力学模型和流动模型之间的需求空白，通常采用Fokker-Planck (FP)近似或者Monte Carlo (MC)颗粒模拟。这两种方法可以在所需计算复杂度和捕获等离子体重要物理细节之间找到一个很好的平衡。本章的主要目的是展现用COMSOL Multiphysics求解FP方程的功能。为了对该问题给出一个整体认识，我们把侧重点集中在一个简单的例子上。特别是在第二节我们通过一个简单的例子对FP方程给出一个直观描述，将它用于布朗运动的颗粒模拟。本章最主要的贡献在于介绍了如何在外部电场情况下

4、对电子动力学过程进行建模。最后，在第四节中对如何使用COMSOL Multiphysics实现该模型给出了详细的讨论。2一维FP和Langevin方程对FP方程直接求导可以得到Langevin方程1。考虑浸没在流体中的“布朗”颗粒，如果颗粒足够小，会同时受到两种力。一是颗粒和流体介质间的粘性力，它会降低平均颗粒速度。二是颗粒与流体“分子”间随机碰撞的力。该布朗颗粒的运动控制方程如下1： (1)其中是阻尼系数，随机项表示颗粒与背景流体的连续碰撞。对于这个简单的例子，通常我们假设粘性力线性依赖于颗粒速度。同时，根据随机近似，Langevin力必须满足以下方程1：其中 1，表示整体平均，T是流体温度

5、，k是Boltzmann常数，m是布朗颗粒的质量。同时注意到力可以很容易的通过MATLAB函数randn（见第四节）实现。给定适当的颗粒数和初始条件后，使用MC方法可以很容易的求解方程(1)。为了对颗粒在时间t时的速度进行精确静态测量，颗粒数通常需要超过一百万。最简单的初始条件为。总之，Langevin方程描述了背景流体介质中布朗颗粒的运动。当然，如果没有随机力，颗粒的路径为。但是由于的存在，方程(1)对于很多颗粒算出的解通常是条分布曲线，分布的宽度由q决定1。FP方程为MC方法提供了另一种思路。包含随机力的Langevin方程可以等价与以下偏微分方程： (2)初始条件和边界条件为： (3)其

6、中w是颗粒在速度空间中的分布概率，也就是颗粒出现的概率，在无限大系统中，速度通常表示为wdv。初始条件(3)如图1所示。求解该偏微分方程等价于对无限个颗粒进行Monte-Carlo模拟，以获得t0时的概率分布函数。图2给出了两种方法的比较。MC方法针对初始条件为的情况，求解了20000个常微分方程。这还不是强非线性阻尼力和三维空间情况下的颗粒运动，即使如此，对于这个简单的例子，求解偏微分方程的时间和常微分方程一样。当颗粒痕迹的常微分方程是非线性时，或者需要从FEA模拟中查询外部力的值时，MC方法的计算时间就会变得非常大。图1 初始条件w图2 各时间步长下Fokker-Planck和Monte

7、Carlo方法的对比作为例子，我们考虑修改后的Langevin方程： (4)和相应的FP方程1： (5)其中，。这种情况下的分布函数如图3所示。即使计算一百万个颗粒，常微分方程的分布函数也存在很大的数值噪声，而FP结果就很平滑。这里出现的静态噪声是有限颗粒MC方法模拟的经典问题（也就是说数据的后处理是很精细、很关键的步骤）。图3 各时间步长下基于非线性Langevin方程的Fokker-Planck和Monte Carlo方法的对比下面，根据本节介绍的背景知识，我们来看一个基于FP方程的更为复杂的例子。3外加场中的电子动力学现在我们考虑电子流在电容两极板间外加电场作用下的动力学过程。一个极板保

8、持常电压，另一个极板接地。假设电子流从接地极板侧以一定速度引入，并在阳极侧积累。并且两极板间的整个区域充满了温度为T的静止流体。图4给出了该几何结构。图4 模拟域示意图模型主要假设为：l 忽略由于电子和背景介质碰撞所产生的离子；l 电子间的相互作用可以忽略；l 在一维情况下求解该问题。为了精确求解物理空间（y）和速度空间（vy）的电子分布，需要求解FP方程。在本例中我们还要计算电子密度、电子流和电容两极板间的静电场。在通常情况下，FP方程比静电场问题多需要一个空间维数，所以我们需要同时求解一个一维静电场和二维FP方程。电子分布的控制方程为： (6)对于该电场 (7)用标量形式电势和电子密度表示

9、为 (8)其中 (9)并且 (10)注意这里的ne和Fy都是y的函数。边界条件和初始条件为： (11) (12) (13) (14)其中m是电子质量，k是Boltzmann常数，T是背景流体温度，e是电子电荷，是自由空间介电常数，是阻尼系数。同时，无限大系统中，在速度vy和位置y处找到一个电子的概率为wdvdy。方程(6)和(8)是两个耦合的不同空间维数的偏微分方程。由于COMSOL Multiphysics允许不同空间维数的任意数量几何形状进行耦合，所以我们现在同时求解这两个方程。通过投影耦合变量方法将电子密度引入静电场问题中，通过导出耦合变量的方法将颗粒上的力引入FP问题中。概括起来讲，电

10、子流以速度v0进入“电容放电”区域，同时受到三种力：，和。这些力导致电子速度和密度分布沿y方向变化，所以与之前简单线性分布例子相比，当修改极板电势V时，不存在这种情况。图5给出了函数w的解。注意到初始时刻，当电子进入计算域时速度降低。这种非线性、非单调曲线是由如图6所示的阴极电子屏蔽所导致的。查看域中间位置附近的最大电子密度。图8表明电容的电子电势在最大ne处有最小值。图5 概率分布图6 电子密度图7 平均速度图8 平行极板间电子电势随位置的分布最后，给出电子流量：由于质量守恒，流量在y方向应该保持为常数。图7沿y方向的电子流量表明该过程中质量守恒。在下一节中我们将详细给出在COMSOL中如何

11、实现该模型。但是在开始之前应该注意到，虽然以上模型不是对等离子体的完整模拟，但是通过模拟它可以给出了很多重要的特性，即，一个更完整的模型仍然应该包括FP方程，但是现在必须耦合离子动力学过程。同时，在很多情况下，FP方程都需要在二维空间中求解。FP方程的速度依赖性采用相同的处理方法7。4COMSOL Multiphysics的实现该模型用到COMSOL Multiphysics的一些高级功能。这些问题包括：l 速度空间的拓扑域是无限的，而物理空间只有厘米量级。这会导致几何结构具有极高的横纵比；l 两个偏微分方程的空间维数不同；l 必须对所有的y在速度空间中进行一系列线性积分，然后把得到的值用于静

12、电场问题；l 同样的，静电场问题产生的力必须施加到Fokker-Planck域物理和速度空间的每一个点；l 两个物理场之间存在强耦合，即：电子密度是外部施加力的强函数，且外部施加力是电子密度的强函数。以上每一点都会带来潜在的问题，但是COMSOL Multiphysics可以通过一些方法来克服这些困难。第一点通过调节偏微分方程(6)和(8)的尺度来克服。这使得几何结构具有较低的横纵比，从而提高网格质量。尺度调节后的偏微分方程为： (15) (16)其中m是电子质量，k是Boltzmann常数，T是背景流体温度，e是电子电荷，是自由空间介电常数，d是平行极板间的距离，V0是下极板电压，是阻尼常数

13、，n0是y0处的电子密度。我们通过模型导航栏建立两个几何结构，几何结构1为二维结构，独立变量y，vy和z，偏微分方程通用形式的因变量为w。然后再增加一个几何结构2，使用独立变量x，y和z，偏微分方程通用形式的因变量为V。在几何结构1中，从y0到y1和vy2到vy8绘制一个矩形。定义以下常数：表1 常数表名 称表达式m9.1e-31k1.38e-23T300e1.6e-19e08.84e-12d0.001V00.3gamma1/1e-7n03e13vTsqrt(k*T/m)vCd*gammaalphae*V0/(m*vT*gamma*d)betae*d*d*n0/(V0*e0)在几何结构1中定义

14、以下标量表达式：表2 几何结构1中的标量表达式名 称表达式a2smooth0.1vy00.2w0(1/(2*vy0)*(flc2hs(vy+vy0-a,smooth)-flc2hs(vy-vy0-a,smooth)表3 几何结构1的子域设置0 -wvyF-(vT/vC)*vy*wy+w+vy*wvy-alpha*phix_extr*wvyInitw0在子域设置和边界设置中，在vy=-2，vy=8处设定w0，在y=0处设定w=w0。设定y1处为Dirichlet边界条件，其中r和g都为零。在几何结构2中从x0到x1绘制一条直线，并输入：表4 几何结构2的子域设置-VxF-beta*neInitV

15、0*x几何结构1的边界条件设定如下：表5 几何结构1边界条件位置类型系数y=0Dirichletg=0, r=w0-wy=1Dirichletg=0, r=0vy=-1Dirichletg=0, r=-wvy=10Dirichletg=0, r=-w对于几何结构2，设定x0处V0，x1处VV0。下面我们需要耦合两个几何结构。在几何结构1中设定以下映射耦合变量：表6 几何结构1的映射耦合变量名称表达式neWfluxvy*w打开Source Transformation菜单，选择General Transformation选项并设定x:y和y:vy。在destination选项卡中，分别选择Geo

16、metry：Geom 2，Level：Subdomain，Subdomain Selection：1，Destination Transformation x：x。在几何结构2中，如表7中定义导出耦合变量。选中General Transformation，并设定Source Transformation x：x。在destination选项卡中，选择几何结构1中的子域1，并设定Destination transformation x：y。在几何结构1的y0处和几何结构2的x0处需要加密网格。选择非线性求解器，设定Solution Form：weak and click Solve。表7 几何结

17、构2的导出耦合变量名称表达式phi_extrVx5附录弱解形式该模型必须使用弱形式求解。这是因为导出耦合变量对雅克比矩阵有非常大的影响。COMSOL为了考虑耦合影响必须使用弱形式。下面通过绘图来对比这两种方法（见第一章）。图9给出了使用通用形式求解时的雅克比稀疏矩阵。非零元素为219459个。第二张弱形式的图说明非零元素为1092160个。矩阵中有两带元素完全对应于导出耦合变量。当两个几何结构的交叉耦合影响占主导地位时，为了求解该模型，雅克比矩阵中的这些额外项就变得十分重要，否则牛顿法找不到合理的解。图9 使用通用解形式的雅克比稀疏矩阵图10使用弱解形式的雅克比稀疏矩阵6关于作者Ali Sha

18、jii博士现在是MKS Instruments公司高技术中心主任，MKS Instruments公司是半导体工业传输材料和反应活性气体（等离子体）设备的主要供应商。他1994年在MIT获得工程物理和应用数学博士学位。主要研究领域为反应和非反应流动力学、电动力学和复杂高非线性系统控制理论。目前已获得18项美国专利，在上述领域期刊中发表文章30余篇。Daniel Smith是MKS Instruments公司高技术中心研究员。2002年在圣安德鲁斯大学获得应用数学硕士学位，在曼彻斯特大学获得数值计算硕士学位。他的主要研究领域为粘性可压缩流动、多相流、计算电磁学和多物理场问题。参考文献1 H. Ri

19、sken, The Fokker-Planck Equation (Springer, 1996).2 J. R. Roth, Industrial Plasma Engineering (IOP Publishing, Ltd., 1995).3 M. A. Lieberman and A. J. Lichtenberg, Principles of Plasma Discharges and Materials Processing (Wiley, 1994).4 S. Harris, An Introductions to the theory of the Boltzmann Equa

20、tion (Dover, 2004).5 S. Ichimaru, Basic Principles of Plasma PhysicsA Statistical Approach (Addison-Wesley, 1980).6 J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, 2003).7 V. I. Kolobov, Fokker-Planck modeling of electron kinetics in plasmas and semiconductors, Computational Materials Science 28 (2003) 302-320.