一种二维高精度迭代的非磁化等离子体中的实现方法与流程

本发明属于计算电磁学技术领域，具体涉及一种二维高精度迭代的非磁化等离子体中的实现方法。

背景技术：

时域有限差分(finite-differencetime-domain，fdtd)方法因其计算简单、容易实现等优点，被广泛用于色散媒质的电磁波传播的仿真中。但是，它的时间步长受柯西稳定性条件的限制，不能选取的较大，在多尺寸复杂精细结构模型中，fdtd方法计算速度慢，计算效率低。为了消除柯西稳定性条件的限制，人们提出了无条件稳定时域有限差分方法，比如：交替方向隐式(alternating-direction-implicit，adi)的时域有限差分(adi-fdtd)方法和基于加权拉盖尔多项式的时域有限差分(weighted-laguerre-polynomialsfinite-differencetime-domain，wlp-fdtd)方法。在这些方法中，adi-fdtd方法在使用较大的时间步长时会产生很大的色散误差，而wlp-fdtd方法既能消除柯西稳定性条件的限制，又能解决adi-fdtd方法在使用较大的时间步长时会产生很大的色散误差这个难题，因此wlp-fdtd方法可以高效的求解等离子体中的电磁问题。然而，这种wlp-fdtd方法在求解电磁场过程中，会产生一个大型的稀疏矩阵方程，直接求解此方程会使得计算较复杂，内存消耗较大，于是提出了一种因式分裂的wlp-fdtd方法，该方法在计算速度和计算效率上得到了较大的提高，但是由于该方法是通过添加微扰项，进行因式分裂得来的，计算中会产生分裂误差，为了减小分裂误差、提高计算精度，同时保证较快的计算速度，提出了一种迭代的加权拉盖尔多项式时域有限差分方法。

而由于计算机容量的限制，电磁场的计算只能在有限区域进行。为了能模拟开域电磁波传播过程，必须在计算区域的截断边界处给出吸收边界条件。有人提出了完全匹配层(perfectlymatchedlayer，pml)吸收边界，后来pml被广泛应用于计算区域的截断，而且被证明是非常有效的，但是研究发现这种传统pml对低频以及凋落波的吸收效果并不理想；使用带有复频率偏移(complexfrequencyshift，cfs)因子的pml(cfs-pml)吸收边界可以有效地改善传统pml对低频，凋落波与掠射情况的吸收效果。最近，有人提出了一种使用辅助微分方程的近似完全匹配吸收边界的wlp-fdtd方法，来解色散媒质中的电磁场问题，这种近似完全匹配吸收边界的吸收效果非常差、计算时存在误差，而且此算法计算时间长、内存消耗较大。

技术实现要素：

本发明的目的是提供一种二维高精度迭代的非磁化等离子体中的实现方法，计算速度快、内存消耗小、精度高，且对于低频和凋落波具有很好的吸收效果。

本发明所采用的技术方案是，一种二维高精度迭代的非磁化等离子体中的实现方法，按照以下步骤实施：

步骤1：输入模型文件；

步骤2：初始化参数和设置参数；

步骤3：添加场源到y方向上的电场分量系数中，设置电场分量系数记为初始场值

步骤4：更新计算整个计算区域的y方向上电场分量系数

步骤5：更新计算整个计算区域的x方向上电场分量系数

步骤6：将k+1赋值给k，并判断迭代次数k是否达到预设值，若未达到预设值，则返回步骤4，若达到预设值，则执行步骤7；

步骤7：更新计算整个计算区域的磁场分量系数

步骤8：更新计算整个计算区域的极化电流密度系数

步骤9：更新计算整个计算区域的电磁场分量系数的辅助变量；

步骤10：更新计算观测点处的电磁场分量；

步骤11：将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。

本发明的特点还在于：

步骤1具体为：

计算区域大小nx×ny，其中nx为x方向的网格数，ny为y方向的网格数；空间步长δη，η＝x，y，x为横坐标，y为纵坐标；时间步长δt；真空中的电导率σ，磁导率μ0，介电常数ε0；等离子体碰撞频率υ与等离子体频率ωp；等离子体在计算区域中的位置；吸收边界层数npml与相关参数κηmax，αηmax，σηmax；κηmax取整数，κηmax取值范围为[1，60]；αηmax取值范围为[0，1)；σηmax/σopt取值范围为(0，12]，σopt＝(m+1)/150πδh，m取值范围为[1，20]，δη取值范围为λ为源的波长；仿真计算时长tf；迭代次数k，k≥0且为整数；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数。

步骤2初始化的参数具体为：

将整个计算区域的电磁场分量系数整个计算区域的极化电流密度系数整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度系数的和整个计算区域的辅助变量和其中表示表示拉盖尔多项式其中全部初始化为零；非磁化等离子体参数c5,c6初始化为c5＝0，c6＝2；pml系数c1η，c2η，c3，c4初始化为c1η＝2/ε0s，c2η＝1，c3＝e0/μ0，c4＝2/(ε0s)，式中，η＝x，y，x为横坐标，y为纵坐标，ε0是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]，μ0是空气中的磁导率。

设置的参数具体为：

设置带有cfs因子的sc-pml吸收边界的参数ση,κη,αη；具体为：

ση＝σηmax|η-η0|^m/d^m；

κη＝1+(κηmax-1)|η-η0|^m/d^m；

αη＝αηmax；

其中，η＝x,y，η0为pml层与非pml截面位置，d是pml吸收边界的厚度，κηmax取整数，κηmax取值范围为[1，60]；αηmax取值范围为[0，1)；σηmax根据σopt来设置，σηmax/σopt取值范围为(0，12]，σopt＝(m+1)/150πδη，m取值范围为[1，20]，δη取值范围为λ为源的波长；

设置pml系数c1η，c2η和与等离子参数相关的系数c5，c6；具体为：

c1η＝1/(κηαη+ση+0.5κηε0s)，c2η＝(2αη/ε0s+1)；

步骤3中所添加的场源的表达式为：

jy(t)＝(t-t0)/t×exp(-(t-t0)²/t²)

其中，t表示时间，单位为秒；t0,t为场源参数。

步骤4具体为：

步骤4.1：电场分量系数在计算区域的方程为：

其中，k表示迭代次数，i表示x轴方向上第i个计算网格的位置，j表示y轴方向上第j个计算网格的位置，i+1/2表示x轴方向上第i个半网格的位置，j表示y轴方向上第j个半网格的位置；c1x|i表示系数c1x在x轴方向上第i个网格处的值；表示第k+1次迭代x轴方向上第i-1个网格y轴方向上第j个半网格处的电场分量系数的值；

步骤4.2：使用追赶法求解步骤4.1的方程，得到整个计算区域的电场分量系数

步骤5具体为：

步骤5.1：电场分量系数在计算区域的方程为：

步骤5.2：使用追赶法求解步骤5.1的方程，得到整个计算区域的电场分量系数

步骤7具体更新方程为：

步骤8更新公式具体为：

步骤9更新公式具体为：

其中表示表示电场分量系数的和对y的导数、电场分量系数的和对x的导数、磁场分量系数的和对x或y的求导，表示辅助变量表示辅助变量的和。

步骤10更新公式具体为：

其中，u表示电磁场分量ex,ey,hz，r表示电磁场分量的位置，u^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s>0的扩展时间，是q阶拉盖尔多项式。

本发明的有益效果是：

1).本发明一种二维高精度迭代的非磁化等离子体中的实现方法，在直角坐标系下，通过用加权拉盖尔多项式表示电磁场分量，来解时域麦克斯韦方程，使得在更新计算整个计算区域的电磁场分量系数时不涉及到时间步长，只是在最后计算观测点处的电磁场分量时用到时间步长，因此计算过程中时间步长可以取得比柯西稳定性条件限制的时间步长更大；

2).本发明一种二维高精度迭代的非磁化等离子体中的实现方法，在求解非磁化等离子体中的电磁场分量系数时，使用迭代的因式分裂的方案将大型稀疏矩阵方程分裂成两个迭代的三对角矩阵方程，使得它在计算时比wlp-fdtd方法更简单、计算速度更快、内存消耗更少而且比无迭代的因式分裂的wlp-fdtd方法精度更高；

3).本发明一种二维高精度迭代的非磁化等离子体中的实现方法，在设置pml系数时，由于采用了cfs因子，并且通过调整cfs因子中的参数，可以使得该吸收边界对低频与凋落波的吸收更加有效；

4).本发明一种二维高精度迭代的非磁化等离子体中的实现方法，由于采用了复扩展坐标系，使得pml在实现时避免了场的分裂且与媒质无关。

附图说明

图1是本发明非磁化等离子体中的实现方法的流程示意图；

图2是本发明非磁化等离子体中点源辐射的计算模型的示意图；

图3是本发明的方法与传统的fdtd方法、因式分裂的wlp-fdtd方法在观测点处y方向上电场分量时域波形对比图；

图4是因式分裂的wlp-fdtd方法和本发明的方法在观测点处y方向上电场分量误差对比图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明的一种二维高精度迭代的非磁化等离子体中的实现方法，原理为：首先导出二维非磁化等离子体中，电磁场所满足的复扩展坐标系下的麦克斯韦方程，然后使用二维非磁化等离子体中迭代的因式分裂的wlp-fdtd方法推导出整个计算区域的电磁场分量系数和电流密度分量系数的更新方程，接着使用追赶法求得观测点处的电磁场分量系数，最后采用公式(20)求解观测点处的电磁场分量。

在求解二维非磁化等离子体中电磁波传播所满足的更新方程时，首先需要推导出复扩展坐标系下，pml中电磁场满足的麦克斯韦方程；

碰撞冷等离子体色散媒质中，扩展坐标下，麦克斯韦方程组和相关的联立方程为

式中，h是磁场强度；e是电场强度；j是极化电流密度；ε0、μ0分别为真空中的介电常数和磁导率；ωp是等离子体频率；υ是等离子体碰撞频率。为修正后的微分算子，可以写成

sx、sy和sz是坐标扩展变量，可以表示成

sη＝1+ση/jωε0(5)

加入cfs因子后，可以表示成

sη＝kη+ση/(αη+jωε0)(6)

其中η＝(x,y,z)，kη，ση和αη为pml的有关的参数。

应用扩展坐标的cfs-pml，仅考虑二维tez的情况，麦克斯韦方程组和相关的联立方程可化为：

下面我们引入几个辅助变量如下：

将(6)式分别代入式(12)-(15)，然后做如下变化于是可得下面几式：

我们可以知道

式中u代表ex、ey、hz，是加权拉盖尔多项式，p≥0；t≥0是p阶拉盖尔多项式。

将(20)、(21)代入(16)-(19)，然后应用加勒金的测试过程，可得

式中

s>0是时间尺度因子，q是加权拉盖尔多项式的阶数。

将(20)、(21)代入(7)—(11)，再应用加勒金的测试过程得到：

将(34)代入(31)，(35)代入(32)后得到

式中

将(36)、(37)和(33)式写成矩阵形式如下

式中

如果让那么(43)式可以写为

添加一个微扰项到上式，于是得到

为了强调迭代的特点，将(48)式重写为

上式可以分裂为下面两式

式中是一个非物理中间量，为了解(50)式，我们选择

式中

将(51)式代入(50)式，化简后得到

将上式扩展得到

将(54)式的第四式代入第二式和第五式，第一式和第四式代入第三式得到

对上式进行中心差分得到

上面五式中，i表示x轴方向上第i个计算网格的位置，j表示y轴方向上第j个计算网格的位置，i+1/2表示x轴方向上第i个半网格的位置，j表示y轴方向上第j个半网格的位置，k表示第k次迭代，c1x|i表示系数c1x在x轴方向上第i个网格处的值；表示第k+1次迭代x轴方向上第i-1个网格y轴方向上第j个半网格处的电场分量系数的值；在整个计算区域上，(56)式和(57)式可以写成三对角矩阵差分方程，与wlp-fdtd方法相比，这种迭代的因式分解的wlp-fdtd方法将大型稀疏矩阵方程的求解转变成两个三对角矩阵方程的求解，于是可以使用追赶法，非常简单的解得整个计算区域电磁场分量系数，最后通过公式(20)解得观测点的电磁场分量，而且这种方法添加了迭代的方案，使得其计算精度比因式分裂的wlp-fdtd要高。

本发明一种二维高精度迭代的非磁化等离子体中的实现方法，流程示意图如图1所示，按照以下步骤实施：

步骤1：输入模型文件，具体为：

计算区域大小nx×ny，其中nx为x方向的网格数，ny为y方向的网格数；空间步长δη，η＝x，y，x为横坐标，y为纵坐标；时间步长δt；真空中的电导率σ，磁导率μ0，介电常数ε0；等离子体碰撞频率υ与等离子体频率ωp；等离子体在计算区域中的位置；吸收边界层数npml与相关参数κηmax，αηmax，σηmax；κηmax取整数，κηmax取值范围为[1，60]；αηmax取值范围为[0，1)；σηmax/σopt取值范围为(0，12]，σopt＝(m+1)/150πδη，m取值范围为[1，20]，δη取值范围为λ为源的波长；仿真计算时长tf；迭代次数k，k≥0且为整数；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数。

步骤2：初始化参数和设置参数，具体为：

将整个计算区域的电磁场分量系数整个计算区域的极化电流密度系数整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度系数的和整个计算区域的辅助变量和其中表示表示拉盖尔多项式其中全部初始化为零；非磁化等离子体参数c5,c6初始化为c5＝0，c6＝2；pml系数c1η，c2η，c3，c4初始化为c1η＝2/ε0s，c2η＝1，c3＝ε0/μ0，c4＝2/(ε0s)，式中，η＝x，y，x为横坐标，y为纵坐标，ε0是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]，μ0是空气中的磁导率。

设置的参数具体为：

设置带有cfs因子的sc-pml吸收边界的参数ση,κη,αη；具体为：

ση＝σηmax|η-η0|^m/d^m；

κη＝1+(κηmax-1)|η-η0|^m/d^m；

αη＝αηmax；

其中，η＝x,y，η0为pml层与非pml截面位置，d是pml吸收边界的厚度，κηmax取整数，κηmax取值范围为[1，60]；αηmax取值范围为[0，1)；σηmax根据σopt来设置，σηmax/σopt取值范围为(0，12]，σopt＝(m+1)/150pδη，m取值范围为[1，20]，δη取值范围为λ为源的波长；

设置pml系数c1η，c2η和与等离子参数相关的系数c5，c6；具体为：

c1η＝1/(κηαη+ση+0.5κηε0s)，c2η＝(2αη/ε0s+1)；

步骤3：添加场源到y方向上的电场分量系数中，设置电场分量系数记为初始场值

步骤3中所添加的场源的表达式为：

jy(t)＝(t-t0)/t×exp(-(t-t0)²/t²)

其中，t表示时间，单位为秒；t0,t为场源参数。

步骤4：更新计算整个计算区域的y方向上电场分量系数具体为：

步骤4.1：电场分量系数在计算区域的方程为：

其中，k表示迭代次数，i表示x轴方向上第i个计算网格的位置，j表示y轴方向上第j个计算网格的位置，i+1/2表示x轴方向上第i个半网格的位置，j表示y轴方向上第j个半网格的位置；c1x|i表示系数c1x在x轴方向上第i个网格处的值；表示第k+1次迭代x轴方向上第i-1个网格y轴方向上第j个半网格处的电场分量系数的值。

步骤4.2：使用追赶法求解步骤4.1的方程，得到整个计算区域的电场分量系数

步骤5：更新计算整个计算区域的x方向上电场分量系数具体为：

步骤5.1：电场分量系数在计算区域的方程为：

步骤5.2：使用追赶法求解步骤5.1的方程，得到整个计算区域的电场分量系数

步骤6：将k+1赋值给k，并判断迭代次数k是否达到预设值，若未达到预设值，则返回步骤4，若达到预设值，则执行步骤7；

步骤7：更新计算整个计算区域的磁场分量系数具体更新方程为：

步骤8：更新计算整个计算区域的极化电流密度系数更新公式具体为：

步骤9：更新计算整个计算区域的电磁场分量系数的辅助变量，更新公式具体为：

其中表示表示电场分量系数的和对y的导数、电场分量系数的和对x的导数、磁场分量系数的和对x或y的求导，表示辅助变量表示辅助变量的和。

步骤10：更新计算观测点处的电磁场分量，更新公式具体为：

其中，u表示电磁场分量ex,ey,hz，u^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s>0的扩展时间，是q阶拉盖尔多项式。

步骤11：将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。

下面通过实验对本发明的效果进行说明：

实验：非磁化等离子体中点源辐射的计算

以图2所示的等离子体中点源辐射的计算模型为例，按照本发明的方法步骤进行实施，实验中整个计算区域为50×50网格，网格大小为0.15mm×0.15mm，即δx＝δy＝0.15mm，等离子体充满整个计算区域，其参数ωp＝1.261×10¹¹rad/s,υ＝2×10¹⁰rad/s。四个边界采用10层网格的pml吸收边界，计算中所加的源位于网格(25，25)，所加场源的表达式如下：

jy(t)＝(t-t0)/t×exp(-(t-t0)²/t²)(61)

其中，t0＝0.05ns,t＝0.01ns。观测点位于(38，38)网格处。时间步长δt＝0.25ps，加权拉盖尔多项式的阶数q＝245，时间扩展因子s＝1.256×10¹²，整个仿真时间为tf＝0.75ns，迭代次数k＝2，pml吸收边界参数κηmax＝1，σηmax＝σopt，αηmax＝0。采用本发明方法计算的观测点处的电场分量ey与采用传统fdfd方法和因式分裂的wlp-fdtd方法计算的结果参见图3。采用因式分裂的wlp-fdtd方法和本发明的方法计算观测点处y方向上电场分量误差结果参见图4。从图3可以看出传统fdtd方法与本发明方法计算结果一致，验证了本发明方法的正确性。更进一步的我们可以从图4清晰地发现本发明方法的计算精度相比于因式分裂的wlp-fdtd方法的计算精度要高很多，尤其是在仿真时间的后期。在运行计算时间上，传统fdtd方法需要4.382秒，因式分裂的wlp-fdtd方法需要0.437秒，本发明专利中的计算方法需要0.904秒，由此可知本发明专利中的计算方法相对于fdtd方法计算速度得到较大的提高，虽然计算速度比因式分裂的wlp-fdtd方法要稍慢点但精度比因式分裂的wlp-fdtd方法要高。