基于时域谱元法的短间隙气体放电数值仿真方法与流程
本发明属于短间隙气体放电数值模拟技术领域，特别是一种基于时域谱元法的短间隙气体放电数值仿真方法。
背景技术：
随着国民经济的不断发展，对能源需求的不断提高，超高压输电已经成为远距离输电的主要方式，电力系统的故障很大一部分是由气体绝缘介质击穿造成的，深入了解气体放电的过程对于保证电力系统安全有重要意义。此外气体放电会产生等离子体，等离子体与我们的生活、生产密不可分，在工业、农业、医疗、能源等方面都要涉及到等离子体处理技术。比如等离子体清洗、等离子体合成、等离子体冶炼、材料表面处理、聚合物薄膜制备等。近几年来随着计算机技术和等离子体模拟技术的发展。数值方法分析气体放电问题越来越引起各国科研人员的重视，成为当今研究的热点之一。流注放电模型是多种气体放电放电的初始阶段，是研究气体放电的一个重要的切入点。
通常认为流注放电是一个多尺度多场耦合的非线性动力学过程。描述流注放电的流体的控制方程包括粒子输运方程和与之耦合在一起的泊松方程。粒子输运方程求解的困难在于流注头部的带电粒子梯度很大，一般的数值算法很难处理，即使把网格划分很细，也会出现数值发散或者振荡的现象。通量传输校正技术(fct)是目前解决此类问题的主要方法，传统的通量传输校正只能用在显式迭代上，仿真时间受限于时间步长。
技术实现要素：
本发明的目的在于提供一种仿真时间短、精度高的基于时域谱元法的短间隙气体放电数值仿真方法。
实现本发明目的的技术解决方案为：一种基于时域谱元法的短间隙气体放电数值仿真方法，该方法结合隐式通量传输校正setd-fct技术，具体步骤如下：
第一步，建立短间隙气体的结构模型，使用曲六面体对模型进行离散，得到模型的结构信息；以描述短间隙气体放电的粒子连续性方程和泊松方程为控制方程；对电子密度、离子密度等未知量使用基于谱元法的gll基函数离散，进行伽辽金测试，从而得到高阶方案的表达式，该表达式的解称为高阶解；
第二步，对于高阶方案表达式中的刚度矩阵，通过消除对角线以外的负值元素，得到低阶方案的表达式；高阶方案的表达式减去低阶方案的表达式得到原始的反扩散通量；然后对反扩散通量进行预限制；计算校正因子；最终得到限定后的反扩散通量，在低阶方案的表达式上加入限定后的反扩散通量，从而得到粒子连续性方程最终的半离散形式；
第三步，针对粒子连续性方程最终的半离散形式，通过多步后向差分格式进行时间离散，在每个时间步上使用牛顿迭代进行求解，得到气体的电子密度、离子密度；
第四步，计算电场并更新输运参数，判断是否到达仿真时间：如达到则结束仿真；如未达到则转入第三步。
进一步地，第一步中所述对对流占优的粒子输运方程的电子密度、离子密度等未知量使用基于谱元法的gll基函数离散，进行伽辽金测试，从而得到高阶方案的表达式，具体为：
电子连续性方程如下：
其中，t代表时间，▽是拉普拉斯算子，ne为电子密度，νe为电子迁移速度，de为电子扩散速度，α为电离系数；
对电子连续性方程进行基函数展开，伽辽金测试得到：
其中，[tel]ij、[seu]ij、[sed]ij、[teu]ij分别为形成的矩阵，i、j表示矩阵的行和列，ni、nj分别为测试基函数和展开基函数；v是剖分单元的体积、δt为时间步长；
离子连续性方程的处理方式同上。
进一步地，步骤2所述对于高阶方案表达式中的刚度矩阵，通过消除对角线以外的负值元素，得到低阶方案的表达式，具体为：对于对流项▽·(νene)与扩散项▽2(dene)空间离散得到的矩阵，消除对角线以外的负值元素，从而构建出低阶方案；
以电子连续性方程为例，低阶方案表达式的构建方法如下：
其中，
离子连续性方程的处理方式同上。
进一步地，第三步所述，针对粒子连续性方程最终的半离散形式，通过多步后向差分格式进行时间离散，在每个时间步上使用牛顿迭代进行求解，得到气体的电子密度、离子密度，具体如下：
用4步的后向差分格式，进行时间离散，公式如下：
其中，n代表第n个时间步、δt为时间步长、limfn+4为限制后的反扩散通量；
式(7)中反扩散项为非线性项，在每个时间步采用牛顿迭代求解，公式如下：
其中，n代表第n个时间步，m代表每个时间步内的第m次牛顿迭代，limf(m)是第m次迭代时的反扩散通量。
本发明与现有技术相比，其显著优点为：(1)将一种可以用在隐式差分的通量传输校正方案(setd-fct)用于气体放电领域，从而使时间差分可以用隐格式展开，增大了时间步，减小了仿真时间；(2)同时时域谱元法作为一种高阶有限元方法，具有精度高、质量矩阵易求逆、曲六面体剖分单元可高度拟合待求物体模型等优点，剖分尺寸可以放大，允许使用较大的离散网格，从而减少未知量，提高仿真效率；(3)可以用在隐式差分上，使时间步长不在受限于剖分尺寸，从而可以放大时间步长，减少仿真时间。
下面结合附图对本发明作进一步详细描述。
附图说明
图1是布满电荷的圆盘模型示意图。
图2是气体放电电极配置图。
图3是不同时刻短间隙放电轴线上净电荷密度和电场强度的分布图，其中(a)为0.5ns流注通道内净电荷分布图，(b)为0.5ns流注通道内电场强度分布图，(c)为1.5ns流注通道内净电荷分布图，(d)为1.5ns流注通道内电场强度分布图，(e)为2.5ns流注通道内净电荷分布图，(f)为2.5ns流注通道内电场强度分布图。
具体实施方式
本发明使用可用于隐式时间差分格式的通量传输校正(fct)技术，对描短间隙气体放电物理过程的流体力学模型进行数值求解。其中，流体力学模型包含粒子输运方程和与之耦合在一起的泊松方程。求解粒子输运方程求解的难点，在于气体放电过程中，流注头部的带电粒子空间分布的梯度很大，一般的数值算法很难处理，即使在建模时以极细的网格对流注区域进行离散，在计算时也会出现数值发散或者数值振荡的现象。通量传输校正技术(fct)是解决这一难题的常用方法，但传统的fct技术使用显式时间差分格式，在时间上是有条件稳定的，本专利首次将一种可用于隐式时间差分格式的fct技术用于气体放电领域，并且结合时域谱元法求解该流体力学模型，从而实现对气体放电过程的精确仿真。由于该方法采用隐式时间差分格式，因此在时间上可以做到无条件稳定，从而极大地减少仿真时间。同时，作为一种高阶的有限元方法，时域谱元法还具有计算精度高的特点，并且允许使用较大的离散网格，从而减少未知量，提高仿真效率。
本发明首次将一种可以用在隐式差分的fct方案用于气体放电领域，并和时域谱元法结合，在时间迭代时采用多步后项差分格式展开时间偏导项，来求解描述流注放电的方程组，从而实现了对气体放电过程的仿真。
下面结合附图对本发明作进一步详细描述，结合图1图2说明隐式迭代的短间隙气体放电模型的数值分析方法，步骤如下：
第一步，建立短间隙气体的结构模型，使用曲六面体对模型进行离散，得到模型的结构信息；以描述短间隙气体放电的粒子连续性方程和泊松方程为控制方程；对电子密度、离子密度等未知量使用基于谱元法的gll基函数离散，进行伽辽金测试，从而得到高阶方案的表达式，该表达式的解称为高阶解；
(1)建立气体击穿的模型得到其结构的信息(包括曲六面体节点和单元信息)，设置仿真时需要的气压，背景电场，初始电荷条件，仿真时间等参数。
(2)对粒子输运方程进行伽辽金测试，对未知量基函数展开、得到高阶方案。
本发明求解电子连续性方程时，使用了可以用于隐式迭代的通量传输校正方案，其具体操作步骤如下：
以电子连续性方程为例，离子连续性方程相同处理：
其中，t代表时间，▽是拉普拉斯算子，ne为电子密度，νe为电子迁移速度，de为电子扩散速度，α为电离系数；
对其进行伽辽金测试可以得到
其中v为积分区域，i＝1,2···n。
未知量n进行基函数展开
将公式(3)、(4)、(5)代入(2)可以得到
其中，nj为展开基函数、ne为电子密度、ni为测试基函数。
将(6)进行化简可以得到
其中：
其中，nj为展开基函数、ne为电子密度、ni为测试基函数，[te]ij、[seu]ij、[sed]ij、[teu]ij分别为形成的矩阵，i、j表示矩阵的行和列，v是剖分单元的体积、δt为时间步长；
公式(7)就是使用setd-fct求解电子连续性方程时所构建出的高阶方案表达式，离子连续性方程的处理方式同上。
第二步，对于高阶方案表达式中的刚度矩阵，通过消除对角线以外的负值元素，得到低阶方案的表达式；高阶方案的表达式减去低阶方案的表达式得到原始的反扩散通量；然后对反扩散通量进行预限制；计算校正因子；最终得到限定后的反扩散通量，在低阶方案的表达式上加入限定后的反扩散通量，从而得到粒子连续性方程最终的半离散形式；
(1)构建低阶方案
对流项▽·(νene)与扩散项▽2(dene)空间离散得到的矩阵seu、sed，消除输运算子上负的偏离对角线的元素，构建低阶方案。具体为：
其中，
离子连续性方程的处理方式同上。
(2)原始的反扩散通量由公式(12)得到：
f＝-deune-dedne,
fij:＝0,iffij(nei-nej)＜0(13)
其中，f为原始的反扩散通量，反扩散通量为矩阵形式，其中每个元素按公式(12)和公式(13)填充，
(3)计算校正因子
使用zalesak通量限制器，计算校正因子αij。ni表示节点i上的解，si表示节点i以及其相邻节点的集合，用
对于局部最大或最小值，所有的反扩散通量必须完全取消：
节点i解的增量
其中,mi为矩阵te主对角线元素值，nei表示结点i上的电子密度值，
节点i可接受的通量为
其中，
节点i可接受的通量相比解增量的比例为：
其中，ri为节点i上可接受的通量和解增量的比值。
校正因子αij由如下给出：
第三步，针对粒子连续性方程最终的半离散形式，通过多步后向差分格式进行时间离散，在每个时间步上使用牛顿迭代进行求解，得到气体的电子密度、离子密度；
展开时间偏导项得到时间迭代的表达式，在低阶方案的右侧加入限定后的反扩散通量αijfij，隐式的时候反扩散通量是非线性项，需要用牛顿迭代求解下面的方程：
其中，n代表第n个时间步、∑αijfij为限制后的反扩散通量；
下面以4步的后向差分格式(吉尔格式)为例，进行时间离散，展开公式(19)
其中，其中，n代表第n个时间步、limfn+4为限制后的反扩散通量；
式(20)中反扩散项为非线性项在每个时间步采用牛顿迭代求解。牛顿迭代后每个时间步的表达式如(21)所示，
其中，n代表第n个时间步，m代表每个时间步内的第m次牛顿迭代，limf(m)是第m次迭代时的反扩散通量。
以上便是用无条件稳定的格式求解电子连续性方程方程的具体实施方式，对于正负离子连续性方程也可以采用类似的方法求解，这里不再赘述。
第四步，计算电场并更新输运参数，判断是否到达仿真时间：如达到则结束仿真；如未达到则转入第三步。
粒子连续性方程中输运参数是电场强度和气压的函数对于电场强度。电场强度可以通过求解泊松方程得到，特别的对于均匀场中的气体放电，电场强度的值可以用以下方法计算。
对于空间中半径为rd厚度为dz的圆盘上面布满电荷如图1所示，电荷密度为σ那么空间中的电场强度可以用下式来表示：
其中，e(z)为z方向上的电场强度，ε为介电常数，z为观察点距离阴极的距离，z'是该点到带电圆盘的距离。
考虑镜像电荷并对z'进行积分就能得到泊松场
加上外加电压产生的拉普拉斯电场el就能知道空间中的电场。
实施例1
为了验证本发明的正确性与有效性，下面分析短间隙气体放电中的双头流注的时空演变情况。
算例的电极配置情况如图2所示，其几何模型为：0.0001cm×0.0001cm×1cm的长方体，剖分尺寸0.0001cm，仿真的背景电场为50kv/cm，气压为760torr，仿真所用的输运参数如下：电离系数α＝5.7pexp(-259p/e)(cm-1)其中气压的单位为torr，电场的单位为v/cm；电子扩散系数为dz＝1800(cm2/s)；电子迁移率νe＝-2.9e5/p·e(cm2/s)。我们认为流注起始于极板中心位置的等离子团，带电粒子密度呈高斯分布。经过一个短暂的初始阶段后进入了一个稳态(steady-state)发展的传播过程，本算例通过观察轴线上电场和电子密度的分布描述了这一现象。
本专利中模拟的是氮气产生的双头流注，氮气是一种中性气体，在仿真过程中没有负离子出现因此不在考负离子连续性方程，并且由于仿真时间短，离子质量大，导致离子迁移和扩散速度远小于电子的迁移和扩散速度，因此可以忽略正离子的位移。在仿真初期，电离率远大于附着率，因此忽略了中的附着系数，同时用背景电离来代替光电离项。
分析：图3别描述了0.5ns、1.5ns、2.5ns氮气中的双头流注在轴线上的电子密度和电场强度的空间分布，其中(a)为0.5ns流注通道内净电荷分布图，(b)为0.5ns流注通道内电场强度分布图，(c)为1.5ns流注通道内净电荷分布图，(d)为1.5ns流注通道内电场强度分布图，(e)为2.5ns流注通道内净电荷分布图，(f)为2.5ns流注通道内电场强度分布图。从图3可以看出流注头部(净电荷空间梯度最大值)随着时间的发展而向前传播，传播速度在的数量级都在107～108cm/s之间，在流注头部存在0.1～0.2mm的静电荷层。流注内部的净电荷持续增加，这是因为流注通道内的电场没有被完全屏蔽(流注通道内部的电场值约为背景电场20％到70％)，导致流注内部气体继续被电离。流注中静电荷的存在导致空间电场发生畸变，使流注通道头部电场大于背景电场。此外从图3可以看出流注内部电压降小而外部电压降大，这是因为流注内部已经被电离形成了等离子体通道，而等离子体的电导率很大造成的。这些结论都和现有流注放电文献吻合。
图3还对比了显格式的龙格库塔差分和隐格式的吉尔差分的结果，并且给出了仿真耗时的对比，从图3可以看出二者电子密度和电场强度的分布曲线基本吻合，龙格库塔格式的时间步长要满足下面条件：
δx是剖分尺寸，ν是对流速度。从图3发现龙格库塔法在cfl≈0.1的情况下1.5ns、2.5ns时净电荷分布出现了振荡，时间步长不能继续放大，而吉尔差分法在cfl≈2.0时仍能保持结果的稳定，随着时间的推移出现差异的原因是时间步长放大时会影响反扩散项，但差异在可容忍的范围内。本算例中龙格库塔法的cfl≈0.1，时间步长为1×10-13s，用时225.7min。而吉尔差分方法的cfl≈2，时间步长为20×10-13s用时28.1min。对比仿真消耗的时间可以发现，使用隐式差分时，虽然破坏了质量矩阵的块对角特性，并且由于引入了非线性项每一个时间步都要进行牛顿迭代求解，但在保证精度和稳定性的情况下仿真时间缩短了约87％。这都说明了本专利提出的算法的准确性和有效性。
- 还没有人留言评论。精彩留言会获得点赞！