CN104809286A - 一种等离子体中扩展坐标的完全匹配吸收边界的实现方法 - Google Patents

Abstract

本发明公开了一种等离子体中扩展坐标的完全匹配吸收边界的实现方法，包括：输入模型文件；初始化参数以及设置PML系数和吸收边界参数；更新计算整个计算区域磁场分量系数添加场源到电场分量系数更新计算整个计算区域的电场分量系数更新计算整个计算区域的电场分量系数更新计算整个计算区域的电子平均速度更新计算整个计算区域的电磁场分量系数的辅助变量；计算观测点处电磁场分量；将q+1赋值给q，并判断q是否达到预设值，若q未达到预设值，返回步骤3；若q达到预设值时，则结束。本发明的一种等离子体中扩展坐标的完全匹配吸收边界的实现方法，吸收边界能够很方便的与CFS因子中的参数相结合，对于低频和凋落波具有较好的吸收效果。

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方法被用于求解精细结构模型下的电磁场问题。

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

发明内容

本发明的目的是提供一种等离子体中扩展坐标的完全匹配吸收边界的实现方法，扩展坐标的完全匹配吸收边界能够很方便的与CFS因子中的参数相结合，对于低频和凋落波具有较好的吸收效果。

本发明所采用的技术方案是：一种等离子体中扩展坐标的完全匹配吸收边界的实现方法，包括以下步骤：

步骤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₄＝0，C₅＝2e/(ms)，C₆＝2；式中，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]，μ₀是空气中的磁导率，e，m分别是电子的电量和质量；

设置的参数具体包括：

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

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

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

α_η＝α_ηmax；

式中，η＝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]，其中m取值为4时边界的吸收效果最好，Δη取值范围为λ为源的波长；

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

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

C₄＝2en_e/(ε₀s)，C₅＝e/(0.5ms+mυ)，C₆＝s/(0.5s+υ)；

步骤3，更新计算整个计算区域磁场分量系数

步骤4，添加场源到电场分量系数更新计算整个计算区域的电场分量系数

场源的表达式为：

J_y(t)＝(t-t₀)/τ×exp(-(t-t₀)²/τ²)；

具体更新公式为：

步骤5，更新计算整个计算区域的电场分量系数Exq；具体按照以下公式进行更新计算：

步骤6，更新计算整个计算区域的电子平均速度ueqx,ueqy；具体按照以下公式进行更新计算：

$u_{\mathrm{ex}}^q{|}_{i+1/2,j}=-C_5{|}_{i+1/2,j}{E}_x^q{|}_{i+1/2,j}-C_6{|}_{i+1/2,j}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i+1/2,j};$

$u_{\mathrm{ex}}^q{|}_{i,j1/2}=-C_5{|}_{i,j1/2}{E}_x^q{|}_{i,j1/2}-C_6{|}_{i,j1/2}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i,j1/2};$

步骤7，更新计算整个计算区域的电磁场分量系数的辅助变量，即计算和具体按照以下公式进行更新计算：

式中，F_ζ＝E_x,E_y,H_z；η＝y,x；

步骤8，计算观测点处的电磁场分量，即计算E_x、E_y和H_z；具体按照以下公式更新计算：

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

步骤9，将q+1赋值给q，并判断q是否达到预设值，若q未达到预设值，则返回步骤3；若q达到预设值时，则结束。

本发明的特点还在于，

步骤3具体为：

首先，给出某些网格的磁场分量系数：

式中，C₃＝ε₀/μ₀，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

然后，根据步骤3.1中的公式，将整个计算区域的磁场分量系数写成矩阵方程形式：式中A表示矩阵系数，表示整个计算区域的第q阶的磁场分量系数，β^q-1(r)表示整个计算区域在q-1阶的变量，q-1阶的变量包括电磁场分量系数的和，整个计算区域的电子平均速度的和，辅助变量；

最后，通过矩阵求逆的方式求解矩阵方程得到整个计算区域的磁场分量系数

本发明的有益效果是：

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

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

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

附图说明

图1是本发明一种等离子体中扩展坐标的完全匹配吸收边界的实现方法的流程示意图；

图2是本发明实验中等离子体中电磁波传播的计算模型的示意图；

图3是实验中利用本发明的方法计算等离子体中电波传播的反射和透射系数幅度的曲线图；

图4是实验中利用本发明的方法计算等离子体中电波传播的反射和透射系数相位的曲线图。

具体实施方式

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

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

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

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

$\widetilde{\▿}\×H={\mathrm{\ϵ}}_0\frac{\∂E}{\∂t}+J---\left(1\right)$

$\widetilde{\▿}\×E=-{\mathrm{\μ}}_0\frac{\∂H}{\∂t}---\left(2\right)$

$\frac{\∂u_e}{\∂t}=-\frac{e}{m}E-\mathrm{\υ}{u}_e---\left(3\right)$

J＝-en_eu_e       (4)

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

$\widetilde{\▿}=\hat{x}\frac{\∂}{s_x\∂x}+\hat{y}\frac{\∂}{s_y\∂y}+\hat{z}\frac{\∂}{s_z\∂z}---\left(5\right)$

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

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

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

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

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

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

$\frac{1}{s_y}\frac{\∂H_z}{\∂y}={\mathrm{\ϵ}}_0\frac{\∂E_x}{\∂t}-e{n}_e{u}_{\mathrm{ex}}---\left(8\right)$

$-\frac{1}{s_x}\frac{\∂H_z}{\∂x}={\mathrm{\ϵ}}_0\frac{\∂E_y}{\∂t}-{\mathrm{en}}_e{u}_{\mathrm{ey}}---\left(9\right)$

${\mathrm{\μ}}_0=\frac{\∂H_z}{\∂t}=\frac{1}{s_y}\frac{\∂E_x}{\∂y}-\frac{1}{s_x}\frac{\∂E_y}{\∂x}---\left(10\right)$

$\frac{\∂u_{\mathrm{ex}}}{\∂t}=-\frac{e}{m}{E}_x-\mathrm{\υ}{u}_{\mathrm{ex}}---\left(11\right)$

$\frac{\∂u_{\mathrm{ey}}}{\∂t}=-\frac{e}{m}{E}_y-\mathrm{\υ}{u}_{\mathrm{ey}}---\left(12\right)$

然后，使用JEC-WLP-FDTD方法推导出整个计算区域的电磁场分量系数和电子平均速度的更新方程；

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

${\tilde{H}}_{\mathrm{zy}}=\frac{1}{s_y}\frac{\∂H_z}{\∂y},{\tilde{H}}_{\mathrm{zx}}=\frac{1}{s_x}\frac{\∂H_z}{\∂x},{\tilde{E}}_{\mathrm{xy}}=\frac{1}{s_y}\frac{\∂E_x}{\∂y},{\tilde{E}}_{\mathrm{yx}}=\frac{1}{s_x}\frac{\∂E_y}{\∂x}---\left(13\right)$

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

$\begin{array}{c}({\mathrm{\κ}}_{\mathrm{\η}}{\mathrm{\α}}_{\mathrm{\η}}+{\mathrm{\σ}}_{\mathrm{\η}}){\tilde{F}}_{\mathrm{\ζ\η}}+{\mathrm{\κ}}_{\mathrm{\η}}{\mathrm{\ϵ}}_0\frac{\∂{\tilde{F}}_{\mathrm{\ζ\η}}}{\∂t}={\mathrm{\α}}_{\mathrm{\η}}\frac{\∂F_{\mathrm{\ζ}}}{\∂\mathrm{\η}}+{\mathrm{\ϵ}}_0\frac{\∂}{\∂t}\left(\frac{\∂F_{\mathrm{\ζ}}}{\∂\mathrm{\η}}\right)\\ (F_{\mathrm{\ζ}}=E_x,E_y,H_z),(\mathrm{\η}=y,x)\end{array}---\left(14\right)$

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

上式中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)式，再应用加勒金的测试过程得到：

$u_{\mathrm{ex}}^q=-C_5{E}_x^q-C_6\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k---\left(22\right)$

$u_{\mathrm{ey}}^q=-C_5{E}_y^q-C_6\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ey}}^q---\left(23\right)$

上式中

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

C₃＝ε₀/μ₀，C₄＝2en_e/(ε₀s)      (25)

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

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

将公式(21)～(23)、(27)和(28)进行中心差分，离散化后分别得到

$u_{\mathrm{ex}}^q{|}_{i+1/2,j}=-C_5{|}_{i+1/2,j}{E}_x^q{|}_{i+1/2,j}-C_6{|}_{i+1/2,j}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i+1/2,j}---\left(30\right)$

$u_{\mathrm{ex}}^q{|}_{i,j1/2}=-C_5{|}_{i,j1/2}{E}_x^q{|}_{i,j1/2}-C_6{|}_{i,j1/2}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i,j1/2}---\left(31\right)$

将(32)式和(33)式代入(29)式得到

将上式写成矩阵方程形式，具体为：

$\left[A\right]\left[H_z^q\left(r\right)\right]=\left[{\mathrm{\β}}^{q-1}\left(r\right)\right];$

式中，A表示矩阵系数，表示整个计算区域的第q阶的磁场分量系数，β^q-1(r)表示整个计算区域在q-1阶的变量，q-1阶的变量包括电磁场分量系数的和，整个计算区域的电子平均速度的和，辅助变量；通过矩阵求逆的方式求解矩阵方程得到整个计算区域的磁场分量系数然后往回迭代就可以得到整个计算区域的电场分量系数接着就可以得到整个计算区域的电子平均速度最后应用(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/σ_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₄＝0，C₅＝2e/(ms)，C₆＝2；式中，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]，μ₀是空气中的磁导率，e，m分别是电子的电量和质量；

设置的参数具体包括：

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

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

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

α_η＝α_ηmax；

式中，η＝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]，其中m取值为4时边界的吸收效果最好，Δη取值范围为λ为源的波长；

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

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

C₄＝2en_e/(ε₀s)，C₅＝e/(0.5ms+mυ)，C₆＝s*(0.5s+υ)；

步骤3，更新计算整个计算区域磁场分量系数

首先给出某些网格的磁场分量系数，如下：

式中，C₃＝ε₀/μ₀，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

然后将整个计算区域的磁场分量系数写成矩阵方程形式：式中A表示矩阵系数，表示整个计算区域的第q阶的磁场分量系数，β^q-1(r)表示整个计算区域在q-1阶的变量，q-1阶的变量包括电磁场分量系数的和，整个计算区域的电子平均速度的和，辅助变量；通过矩阵求逆的方式求解矩阵方程得到整个计算区域的磁场分量系数

步骤4，添加场源到电场分量系数更新计算整个计算区域的电场分量系数

场源的表达式为：

J_y(t)＝(t-t₀)/τ×exp(-(t-t₀)²/τ²)；

具体更新公式为：

步骤5，更新计算整个计算区域的电场分量系数具体按照以下公式进行更新计算：

步骤6，更新计算整个计算区域的电子平均速度具体按照以下公式进行更新计算：

$u_{\mathrm{ex}}^q{|}_{i+1/2,j}=-C_5{|}_{i+1/2,j}{E}_x^q{|}_{i+1/2,j}-C_6{|}_{i+1/2,j}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i+1/2,j};$

$u_{\mathrm{ex}}^q{|}_{i,j1/2}=-C_5{|}_{i,j1/2}{E}_x^q{|}_{i,j1/2}-C_6{|}_{i,j1/2}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i,j1/2};$

步骤7，更新计算整个计算区域的电磁场分量系数的辅助变量，即计算和具体按照以下公式进行更新计算：

式中，F_ζ＝E_x,E_y,H_z,η＝y,x；

步骤8，计算观测点处的电磁场分量，即计算E_x、E_y和H_z；具体按照以下公式更新计算：

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

步骤9，将q+1赋值给q，并判断q是否达到预设值，若q未达到预设值，则返回步骤3；若q达到预设值时，则结束。

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

实验：等离子体中电磁波传播的计算

按照本发明的方法步骤进行实验过程，实验中加权拉盖尔多项式的阶数设置为300，然后计算电磁场的反射和透射系数幅度以及相位。如图2所示，整个计算区域为100×50网格，网格大小为150μm×150μm，即Δx＝Δy＝150μm；等离子体放置在x方向上第30至69的网格上，其参数n_e＝1.0231×10¹⁹m^-3，υ＝2×10¹⁰rad/s；四个边界均采用10层网格的PML吸收边界，PML吸收边界的参数计算公式为：

σ_η＝ση_max|η-η₀|^m/d^m     (35)

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

α_η＝α_ηmax       (37)

σ_opt＝(m+1)/150πΔζ       (38)

计算中所加的源位于x＝20的线上，表达式如下：

J_y(t)＝(t-t₀)/τ×exp(-(t-t₀)²/τ²)     (39)

其中t₀＝0.05ns,τ＝0.01ns。时间步长Δt＝1ps，加权拉盖尔多项式的阶数q＝500，时间扩展因子s＝5.18×10¹¹，整个仿真时间为T_f＝1ns，PML吸收边界参数κ_ηmax＝1，σ_ηmax＝σ_opt，α_ηmax＝0。采用本发明方法计算的反射和透射系数振幅的计算结果参见图3，从图3中可见，解析解的计算结果与本发明方法计算结果和FDTD方法的计算结果一致，验证了本发明方法的正确性。采用本发明方法计算的反射和透射系数相位的计算结果参见图4，从图4中可见，解析解的计算结果与本发明方法计算结果一致，FDTD方法反射系数相位与解析解的结果吻合的不好，验证了本发明方法比FDTD方法更好。

Claims (2)

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，η＝x，y)和拉盖尔多项式(其中)全部初始化为零；

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

设置的参数具体包括：

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

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

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

α_η＝α_ηmax；

式中，η＝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]，其中m取值为4时边界的吸收效果最好，Δη取值范围为λ为源的波长；

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

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

C₄＝2en_e/(ε₀s)，C₅＝e/(0.5ms+mυ)，C₆＝s/(0.5s+υ)；

步骤3，更新计算整个计算区域磁场分量系数H_z ^q；

步骤4，添加场源到电场分量系数更新计算整个计算区域的电场分量系数

场源的表达式为：

J_y(t)＝(t-t₀)/τ×exp(-(t-t₀)²/τ²)；

具体更新公式为：

步骤5，更新计算整个计算区域的电场分量系数具体按照以下公式进行更新计算：

步骤6，更新计算整个计算区域的电子平均速度具体按照以下公式进行更新计算：

$u_{\mathrm{ex}}^q{|}_{i+1/2,j}=-C_5{|}_{i+1/2,j}{E}_x^q{|}_{i+1/2,j}-C_6{|}_{i+1/2,j}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ex}}^k{|}_{i+1/2,j};$

$u_{\mathrm{ey}}^q{|}_{i,j+1/2}=-C_5{|}_{i,j+1/2}{E}_y^q{|}_{i,j+1/2}-C_6{|}_{i,j+1/2}\underset{k=0,q>0}{\overset{q-1}{\mathrm{\Σ}}}{u}_{\mathrm{ey}}^k{|}_{i,j+1/2};$

步骤7，更新计算整个计算区域的电磁场分量系数的辅助变量，即计算和具体按照以下公式进行更新计算：

式中，F_ζ＝E_x,E_y,H_z；η＝y,x；

步骤8，计算观测点处的电磁场分量，即计算E_x、E_y和H_z；具体按照以下公式更新计算：

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

步骤9，将q+1赋值给q，并判断q是否达到预设值，若q未达到预设值，则返回步骤3；若q达到预设值时，则结束。

2.根据权利要求1所述的一种等离子体中扩展坐标的完全匹配吸收边界的实现方法，其特征在于，所述步骤3按照以下步骤实施：

首先，给出某些网格的磁场分量系数，如下：

式中，C₃＝ε₀/μ₀，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；

然后，将整个计算区域的磁场分量系数写成矩阵方程形式，如下：

$\left[A\right]\left[H_z^q\left(r\right)\right]=\left[{\mathrm{\β}}^{q-1}\left(r\right)\right];$

式中，A表示矩阵系数，表示整个计算区域的第q阶的磁场分量系数，β^q-1(r)表示整个计算区域在q-1阶的变量，q-1阶的变量包括电磁场分量系数的和，整个计算区域的电子平均速度的和，以及辅助变量；

最后，通过矩阵求逆的方式求解矩阵方程得到整个计算区域的磁场分量系数