CN104809343B - 一种等离子体中使用电流密度卷积完全匹配层的实现方法 - Google Patents

Abstract

本发明公开了一种等离子体中使用电流密度卷积完全匹配层的实现方法，包括：输入模型文件；初始化参数、设置PML系数和吸收边界参数；分别更新计算整个计算区域y方向上和x方向上电场分量系数电场分量系数添加场源到磁场分量系数，并更新计算整个计算区域的磁场分量系数；更新计算整个计算区域的极化电流密度更新计算整个计算区域的电磁场分量系数的辅助变量；更新计算观测点处电磁场分量；将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，返回步骤3；若达到预设值，则结束。本发明的一种等离子体中使用电流密度卷积完全匹配层的实现方法，计算速度快，内存消耗小，且对于低频和凋落波具有很好的吸收效果。

Description

一种等离子体中使用电流密度卷积完全匹配层的实现方法

技术领域

本发明属于计算电磁学技术领域，涉及一种等离子体中使用电流密度卷积完全匹配层的实现方法。

背景技术

时域有限差分(Finite-difference time-domain，FDTD)方法因其实现简单，广泛用于色散媒质中对电磁波传播的仿真。但是，它的时间步长受柯西稳定性条件的限制，这限制了FDTD方法在精细结构模型中的应用。为了消除柯西稳定性条件的限制，人们提出了无条件稳定时域有限差分方法，比如：交替方向隐式(Alternating-Direction-Implicit，ADI)的时域有限差分(ADI-FDTD)方法和基于加权拉盖尔多项式的时域有限差分(Weighted-Laguerre-polynomials Finite-difference time-domain，WLP-FDTD)方法。在这些方法中，ADI-FDTD方法在使用较大的时间步长时会产生很大的色散误差，而WLP-FDTD方法既能消除柯西稳定性条件的限制，又能解决ADI-FDTD方法在使用较大的时间步长时会产生很大的色散误差这个难题，因此WLP-FDTD方法被用于求解精细结构模型下的电磁场问题。然而，这种传统的WLP-FDTD方法在求解精细结构的电磁场问题时，会产生一个大型的稀疏矩阵方程，直接求解此方程会使得计算较复杂，计算时间和内存消耗较大。

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

发明内容

本发明的目的是提供一种等离子体中使用电流密度卷积完全匹配层的实现方法，计算速度快，内存消耗小，且对于低频和凋落波具有很好的吸收效果。

本发明所采用的技术方案是：一种等离子体中使用电流密度卷积完全匹配层的实现方法，包括以下步骤：

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δη，η＝x或者η＝y，x为横坐标，y为纵坐标；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；等离子体中的碰撞频率υ与等离子体中的电子密度n_e；等离子体在计算区域中的位置；吸收边界层数NPML与相关参数κ_ηmax，α_ηmax，σ_ηmax；κ_ηmax取整数，κ_ηmax取值范围为[1，60]；α_ηmax取值范围为[0，1)；σ_ηmax/σ_opt取值范围为(0，12]；仿真计算时长T_f；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数；

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

初始化的参数具体包括：将整个计算区域的电磁场分量系数整个计算区域的极化电流密度整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度的和整个计算区域的辅助变量和其中F_ζ表示E_x,E_y,H_z，η＝x或者η＝y和拉盖尔多项式其中全部初始化为零；

PML系数(C_1η，C_2η，C₃，C₄，C₅，C₆)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1，C₃＝ε₀/μ₀，C₄＝2/(ε₀s)，C₅＝0，C₆＝2；式中，ε₀是真空中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]，μ₀是真空中的磁导率，e，m分别是电子的电量和质量；

设置的参数具体包括：

设置CFS-PML吸收边界的参数σ_η,κ_η,α_η；具体为：

σ_η＝σ_ηmax|η-η₀|^m/d^m；

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

α_η＝α_ηmax(d-|η-η₀|)/d；

式中，η＝x或者η＝y，η₀为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系数C_1η，C_2η和与等离子参数相关的系数C₅，C₆，；具体为：

C_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s)，C_2η＝(2α_η/ε₀s+1)；

C₅＝2e²n_e/(ms+2mυ)，C₆＝s/(0.5s+υ)；

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

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

步骤5，添加场源到磁场分量系数中，并更新计算整个计算区域的磁场分量系数；

其中，所添加场源的表达式为：

I_mz(t)＝sin(2πf(t-t₀))×exp(-(t-t₀)²/τ²)；

式中，t₀，τ为场源参数；

具体更新公式为：

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

式中的q表示加权拉盖尔多项式的阶数，分别表示x方向上极化电流密度分量系数和y方向上极化电流密度分量系数；

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

步骤8，更新计算观测点处的电磁场分量，具体按照以下公式更新计算：

上式中U表示电磁场分量E_x,E_y,H_z，q表示加权拉盖尔多项式的阶数，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s＞0的扩展时间，是q阶拉盖尔多项式；

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

本发明的特征还在于，

步骤3具体为：

步骤3.1，给出电场分量系数在计算区域的方程，如下所示：

式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

步骤3.2，使用追赶法对整个计算区域的电场分量系数进行求解。

步骤4具体为：

步骤4.1，给出电场分量系数在计算区域的方程，如下所示：

步骤4.2，使用追赶法求解系数为三对角的整个计算区域的电场分量系数

本发明的有益效果是：

1).在直角坐标系下，通过用加权拉盖尔多项式表示电磁场分量，来解时域麦克斯韦方程，使得在更新计算整个计算区域的电磁场分量系数时不涉及到时间步长，只是在最后计算观测点处的电磁场分量时用到时间步长，因此计算过程中时间步长可以取得比柯西稳定性条件限制的时间步长更大；

2).在求解电磁场分量系数时，将大型稀疏矩阵方程分裂成两个三对角矩阵方程，使得它在计算时比传统的WLP-FDTD方法更简单、计算速度更快、内存消耗更少而且可以对大区域的电磁场问题进行求解；

3).在设置PML系数时，由于采用了CFS因子，并且通过调整CFS因子中的参数，可以使得该吸收边界对低频与凋落波的吸收更加有效；

4).由于采用了复扩展坐标系，使得PML在实现时避免了场的分裂且与媒质无关。

附图说明

图1是本发明一种等离子体中使用电流密度卷积完全匹配层的实现方法的流程示意图；

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

图3是本发明的方法与传统的FDTD方法在观测点处时域波形对比图；

图4是本发明实验中观测点的不同吸收边界相对反射误差。

具体实施方式

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

本发明的一种等离子体中使用电流密度卷积完全匹配层的实现方法，原理为：首先导出等离子体中，电磁场所满足的复扩展坐标系下的麦克斯韦方程，然后使用新的高速有效的电流密度卷积的加权拉盖尔多项式的时域有限差分方法(JEC-WLP-FDTD)推导出整个计算区域的电磁场分量系数和电流密度的更新方程，最后采用公式(15)的第一式求解观测点处的电磁场分量。

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

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

J＝-en_eu_e (4)

式中，H是磁场强度；E是电场强度；J是极化电流密度；ε₀，μ₀分别为真空中的介电常数和磁导率；n_e是电子密度；u_e是电子平均速度；υ是等离子体碰撞频率；e，m分别是电子的电量和质量。为修正后的微分算子，可以写成

s_x,s_y和s_z是坐标扩展变量，可以表示成

s_η＝k_η+σ_η/jωε₀ (6)

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

s_η＝k_η+σ_η/(α_η+jωε₀) (7)

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

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

然后，使用新的高速有效的JEC-WLP-FDTD方法推导出整个计算区域的电磁场分量系数和电流密度的更新方程；具体为：

为了计算方便，引入下面四个辅助变量；

将(7)代入(13)，然后利用jω→t的变换，可以得到

由于电磁场分量和其对时间的一阶偏导可以展开成一系列的电磁场分量系数与加权拉盖尔多项式的函数之和，公式如下：

式中，U表示电磁场分量E_x,E_y,H_z，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s＞0的扩展时间，是q阶拉盖尔多项式。将(15)式代入(14)式，然后应用加勒金的测试过程，即在等号两边同时乘以然后对时间进行积分，可以得到：

式中，

C_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s) (18)

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

将(15)式代入(8)～(12)式，再应用加勒金的测试过程得到：

上式中

C_2η＝(2α_η/ε₀s+1) (24)

C₅＝2e²n_e/(ms+2mυ)，C₆＝s/(0.5s+υ) (26)

将(22)式代入(19)式，(23)式代入(20)式后分别得到

将(27)、(28)和(21)式写成矩阵形式如下

式中

如果让则(29)式可以写为

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

上式可以分裂为下面两式

式中是一个非物理中间量，为了解(35)式，令矩阵A和矩阵B为：

式中：

D_Ea＝[0 -C₃C_2xC_1xD_x],D_Eb＝[C₃C_2yC_1yD_y 0] (37)

将(36)式代入(35)式化简后得到

将上式扩展得到

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

对上式进行中心差分得到

上面五式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

由于在整个计算区域上，(43)式和(44)式可以写成三对角矩阵差分方程，于是可以使用追赶法，解得整个计算区域电磁场分量系数，最后通过(15)式解得观测点的电磁场分量。

本发明一种等离子体中使用电流密度卷积完全匹配层的实现方法，具体实施过程如图1所示，包括以下步骤：

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δη，η＝x或者η＝y，x为横坐标，y为纵坐标；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；等离子体中的碰撞频率υ与等离子体中的电子密度n_e；等离子体在计算区域中的位置；吸收边界层数NPML与相关参数κ_ηmax，α_ηmax，σ_ηmax；κ_ηmax取整数，κ_ηmax取值范围为[1，60]；α_ηmax取值范围为[0，1)；σ_ηmax/σo_pt取值范围为(0，12]；仿真计算时长T_f；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数；

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

初始化的参数具体包括：将整个计算区域的电磁场分量系数整个计算区域的极化电流密度整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度的和整个计算区域的辅助变量和其中F_ζ表示E_x,E_y,H_z，η＝x或者η＝y和拉盖尔多项式其中全部初始化为零；

PML系数(C_1η，C_2η，C₃，C₄，C₅，C₆)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1，C₃＝ε₀/μ₀，C₄＝2/(ε₀s)，C₅＝0，C₆＝2；式中，ε₀是真空中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]，μ₀是真空中的磁导率，e，m分别是电子的电量和质量；

设置的参数具体包括：

设置CFS-PML吸收边界的参数σ_η,κ_η,α_η；具体为：

σ_η＝σ_ηmax|η-η₀|^m/d^m；

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

α_η＝α_ηmax(d-|η-η₀|)/d；

式中，η＝x或者η＝y，η₀为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系数C_1η，C_2η和与等离子参数相关的系数C₅，C₆；具体为：

C_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s)，C_2η＝(2α_η/ε₀s+1)；

C₅＝2e²n_e/(ms+2mυ)，C₆＝s/(0.5s+υ)；

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

电场分量系数在计算区域的方程，如下所示：

式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格。

由于上式方程的左边有三个电场分量系数因此可以使用追赶法对整个计算区域的电场分量系数进行求解；

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

电场分量系数在计算区域的方程，如下所示：

根据上式，我们可以使用追赶法求解系数为三对角的整个计算区域的电场分量系数

步骤5，添加场源到磁场分量系数中，并更新计算整个计算区域的磁场分量系数；

其中，所添加场源的表达式为：

I_mz(t)＝sin(2πf(t-t₀))×exp(-(t-t₀)²/τ²)；

式中，t₀，τ为场源参数；

具体更新公式为：

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

式中的q表示加权拉盖尔多项式的阶数，分别表示x方向上极化电流密度分量系数和y方向上极化电流密度分量系数；

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

步骤8，更新计算观测点处的电磁场分量，具体按照以下公式更新计算：

上式中U表示电磁场分量E_x,E_y,H_z，q表示加权拉盖尔多项式的阶数，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s＞0的扩展时间，是q阶拉盖尔多项式；

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

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

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

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

I_mz(t)＝sin(2πf(t-t₀))×exp(-(t-t₀)²/τ²) (48)

其中，t₀＝144ps,τ＝48ps,f＝10.5GHz。观测点位于(38，38)网格处。时间步长Δt＝1.587ps,加权拉盖尔多项式的阶数q＝250，时间扩展因子s＝1.15×10¹²，整个仿真时间为T_f＝0.8ns，PML吸收边界参数κ_ηmax＝6，σ_ηmax＝0.8×σ_opt，α_ηmax＝0.62。采用本发明方法计算的观测点处的磁场分量H_z与采用传统FDFD方法计算的结果参见图3。从图3中可见，传统FDTD方法与本发明方法计算结果一致，验证了本发明方法的正确性。图4为观测点的不同吸收边界相对反射误差，其计算公式可以表示为：

其中，H_pml为当存在SC-PML吸收边界时，观测点的时域波形，H_ref(t)为参考波形，max|H_ref(t)|为参考波形绝对值的最大值。由图4可知，带有CFS因子的SC-PML吸收边界最大的反射误差为-74dB，它比没有CFS因子的SC-PML吸收的吸收效果更好。因此，我们可以调节CFS因子的参数来使边界的吸收效果更好。

Claims (1)

1.一种等离子体中使用电流密度卷积完全匹配层的实现方法，其特征在于，包括以下步骤：

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δη，η＝x或者η＝y，x为横坐标，y为纵坐标；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；等离子体中的碰撞频率υ与等离子体中的电子密度n_e；等离子体在计算区域中的位置；吸收边界层数NPML与相关参数κ_ηmax，α_ηmax，σ_ηmax；κ_ηmax取整数，κ_ηmax取值范围为[1，60]；α_ηmax取值范围为[0，1)；σ_ηmax/σ_opt取值范围为(0，12]；仿真计算时长T_f；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数；

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

初始化的参数具体包括：将整个计算区域的电磁场分量系数整个计算区域的极化电流密度整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度的和整个计算区域的辅助变量和拉盖尔多项式全部初始化为零；其中F_ζ表示E_x,E_y,H_z，其中

PML系数(C_1η，C_2η，C₃，C₄，C₅，C₆)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1，C₃＝ε₀/μ₀，C₄＝2/(ε₀s)，C₅＝0，C₆＝2；式中，e，m分别是电子的电量和质量；

设置的参数具体包括：

设置CFS-PML吸收边界的参数σ_η,κ_η,α_η；具体为：

σ_η＝σ_ηmax|η-η₀|^m/d^m；

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

α_η＝α_ηmax(d-|η-η₀|)/d；

式中，η₀为PML层与非PML截面位置，d是PML吸收边界的厚度，σ_ηmax根据σ_opt来设置，σ_opt＝(m+1)/150πΔη，m取值范围为[1，20]，Δη取值范围为λ为源的波长；

设置PML系数C_1η，C_2η和与等离子参数相关的系数C₅，C₆，具体为：

C_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s)，C_2η＝(2α_η/ε₀s+1)；

C₅＝2e²n_e/(ms+2mυ)，C₆＝s/(0.5s+υ)；

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

具体为：

步骤3.1，给出电场分量系数在计算区域的方程，如下所示：

式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

步骤3.2，使用追赶法对整个计算区域的电场分量系数进行求解；

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

步骤4.1，给出电场分量系数在计算区域的方程，如下所示：

步骤4.2，使用追赶法求解系数为三对角的整个计算区域的电场分量系数

步骤5，添加场源到磁场分量系数中，并更新计算整个计算区域的磁场分量系数；

其中，所添加场源的表达式为：

I_mz(t)＝sin(2πf(t-t₀))×exp(-(t-t₀)²/τ²)；

式中，t₀，τ为场源参数；

具体更新公式为：

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

式中的分别表示x方向上极化电流密度分量系数和y方向上极化电流密度分量系数；

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

步骤8，更新计算观测点处的电磁场分量，具体按照以下公式更新计算：

上式中U表示电磁场分量E_x,E_y,H_z，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s>0的扩展时间；

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