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

Abstract

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

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对低频，凋落波与掠射情况的吸收效果。有人将CFS-PML吸收边界应用到普通媒质的WLP-FDTD方法中，分析了CFS-PML和无CFS因子的PML的吸收效果，但是这种方法应用于等离子体中，则会比较复杂，计算量比较大。

发明内容

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

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

步骤1，输入模型文件；

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

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

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

PML系数(C_1η,C_2η)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1；式中，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

设置的参数具体包括：

设置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_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s)；

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

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

场源的表达式为：

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

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

式中， $C_4=\frac{0.25s^2+0.5\mathrm{\υs}}{0.25s^2+0.5\mathrm{\υs}+{\mathrm{\ω}}_p^2};$

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

${\mathrm{\ψ}}_x^p{|}_{i+1/2,j}=\frac{-{\mathrm{\ω}}_p^2{E}_x^p{|}_{i+1/2,j}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}_x^k{|}_{i+1/2,j}-{\mathrm{\υ}}_{i+1/2,j}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}_x^k{|}_{i+1/2,j}}{0.25s^2+0.5{\mathrm{\υ}}_{i+1/2,j}s};$

${\mathrm{\ψ}}_y^p{|}_{i,j+1/2}=\frac{-{\mathrm{\ω}}_p^2{E}_y^p{|}_{i,j+1/2}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}_y^k{|}_{i,j+1/2}-{\mathrm{\υ}}_{i,j+1/2}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}_y^k{|}_{i,j+1/2}}{0.25s^2+0.5{\mathrm{\υ}}_{i,j+1/2}s};$

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

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

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

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

本发明的特点还在于，

步骤3具体为：

步骤3.1，按照以下公式计算某些网格的磁场分量系数：

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

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

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

本发明的有益效果是：

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

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

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

附图说明

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

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

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

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

图5是本发明实验二中观测点的时域反射误差的曲线图。

具体实施方式：

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

本发明的一种二维等离子体中扩展坐标的完全匹配吸收边界的实现方法，原理为：首先导出等离子体中，电磁场所满足的复扩展坐标系下的麦克斯韦方程，然后使用WLP-FDTD方法推导出整个计算区域的电磁场分量系数更新方程，最后采用公式(10)中的第一式求解整个计算区域的电磁场分量。

等离子体中，电磁场所满足的复扩展坐标系下的麦克斯韦方程推导过程如下：

在复扩展直角坐标系下，PML中电磁场满足的麦克斯韦方程如下：

其中，表示电场矢量，表示磁场矢量，j为虚数单位，ω为角频率，μ为介质的磁导率，ε为介质的介电常数，为修正后的微分算子，可以写成

$\widetilde{\▿}=\hat{x}\frac{\∂}{S_x\∂x}+\hat{y}\frac{\∂}{S_y\∂y}+\hat{z}\frac{\∂}{S_z\∂z}---\left(2\right)$

上式中的S_x,S_y和S_z是坐标扩展变量，可以用下面的公式3统一表示：

S_η＝k_η+σ_η/jωε₀         (3)

在上面的公式(3)中加入CFS因子后，表达式为：

S_η＝k_η+σ_η/(α_η+jωε₀)         (4)

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

本发明中，仅考虑无耗色散媒质中二维横电波的情况，因此，复扩展直角坐标下的麦克斯韦方程可以写成：

$\mathrm{j\ω\ϵ}{E}_x=\frac{1}{s_y}\frac{\∂H_z}{\∂y};---\left(5\right)$

$\mathrm{j\ω\ϵ}{E}_y=-\frac{1}{s_x}\frac{\∂H_z}{\∂x};---\left(6\right)$

$\mathrm{j\ω}{\mathrm{\μ}}_0{H}_z=\frac{1}{s_y}\frac{\∂E_x}{\∂y}-\frac{1}{s_x}\frac{\∂E_y}{\∂x};---\left(7\right)$

其中E_x,E_y分别表示x,y方向的电场，H_z分别表示z方向的磁场；

电磁场分量系数的更新方程推导过程如下：

首先引入四个辅助变量来求解公式(5)、(6)和(7)中的电磁场分量系数更新方程，分别为：

${\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(8\right)$

将公式(4)代入公式(8)中，然后利用jω→t的变换，可以得到：

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

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

将公式(10)代入公式(9)中，然后使用加勒金的测试过程，即在等号两边同时乘以然后对时间进行积分，可以得到：

式中，

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

频率形式的电位移矢量为

其中，若忽略离子的影响，非磁化等离子体的相对介电常数为

${\mathrm{\ϵ}}_r\left(\mathrm{\ω}\right)=1-\frac{{\mathrm{\ω}}_p^2}{{\mathrm{\ω}}^2-\mathrm{j\ω\υ}}---\left(15\right)$

上式中ω_p是等离子体频率，υ是等离子体碰撞频率；

将(15)式代入(14)式得到

即

令

则

使用变化，则(19)式和(16)式变为

电磁场分量对时间的二阶偏导为：

将(10)和(22)代入(20)，然后使用加勒金的测试过程，化简后可以得到

${\mathrm{\ψ}}^p=\frac{-{\mathrm{\ω}}_p^2{E}^p-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}^k-\mathrm{\υs}\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}^k}{0.25s^2+0.5\mathrm{\υs}}---\left(23\right)$

上式中ψ可表示为ψ_x或ψ_y,E可表示为E_x或E_y。

将公式(23)进行中心差分得到中间变量更新方程：

${\mathrm{\ψ}}_x^p{|}_{i+1/2,j}=\frac{-{\mathrm{\ω}}_p^2{E}_x^p{|}_{i+1/2,j}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}_x^k{|}_{i+1/2,j}-{\mathrm{\υ}}_{i+1/2,j}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}_x^k{|}_{i+1/2,j}}{0.25s^2+0.5{\mathrm{\υ}}_{i+1/2,j}s}---\left(24\right)$

${\mathrm{\ψ}}_y^p{|}_{i,j+1/2}=\frac{-{\mathrm{\ω}}_p^2{E}_y^p{|}_{i,j+1/2}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}_y^k{|}_{i,j+1/2}-{\mathrm{\υ}}_{i,j+1/2}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}_y^k{|}_{i,j+1/2}}{0.25s^2+0.5{\mathrm{\υ}}_{i,j+1/2}s}---\left(25\right)$

接下来将(21)式代入(5)式和(6)式得到

${\mathrm{\ϵ}}_0\frac{\∂E_x}{\∂t}-{\mathrm{\ϵ}}_0\frac{\∂{\mathrm{\ψ}}_x}{\∂t}=\frac{1}{s_y}\frac{\∂H_z}{\∂y}---\left(26\right)$

${\mathrm{\ϵ}}_0\frac{\∂E_y}{\∂t}-{\mathrm{\ϵ}}_0\frac{\∂{\mathrm{\ψ}}_y}{\∂t}=-\frac{1}{s_x}\frac{\∂H_z}{\∂x}---\left(27\right)$

将(10)式和(23)式代入(26)式、(27)式以及将(10)式代入(7)式，然后应用加勒金的测试过程，最后进行中心差分得到

上面三式中

C_2η＝(2α_η/ε₀s+1)，C₃＝ε₀/μ₀       (31)

最后，将(24)式代入(28)式,(25)式代入(29)式得到电场分量系数更新方程：

式中

$C_4=\frac{0.25s^2+0.5\mathrm{\υs}}{0.25s^2+0.5\mathrm{\υs}+{\mathrm{\ω}}_p^2}---\left(34\right)$

将(32)式和(33)式代入(30)式可得到：

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

本发明一种二维等离子体中扩展坐标的完全匹配吸收边界的实现方法，具体实施过程如图1所示，包括以下步骤：

步骤1，输入模型文件；

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

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

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

PML系数(C_1η,C_2η)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1；式中，ε₀是空气中的介电常数，s为时间尺度因子，取值为[10⁹，10¹³]；

设置的参数具体包括：

设置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_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s)；

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

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

其中，场源的表达式为：

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

磁场分量系数按照以下公式进行更新计算：

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

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

在求解整个计算区域的磁场分量系数时，根据上式，将整个计算区域的电磁场分量系数写成矩阵方程形式：式中A表示矩阵系数，表示整个计算区域的第p阶的磁场分量系数，β^p-1(r)表示整个计算区域在p-1阶的变量，p-1阶的变量包括电磁场分量系数的和，中间变量系数的和辅助变量；最后通过矩阵求逆的方式解得矩阵方程中整个计算区域的磁场分量系数；

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

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

步骤5，更新计算整个计算区域的中间变量的系数具体更新公式为：

${\mathrm{\ψ}}_x^p{|}_{i+1/2,j}=\frac{-{\mathrm{\ω}}_p^2{E}_x^p{|}_{i+1/2,j}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}_x^k{|}_{i+1/2,j}-{\mathrm{\υ}}_{i+1/2,j}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}_x^k{|}_{i+1/2,j}}{0.25s^2+0.5{\mathrm{\υ}}_{i+1/2,j}s};$

${\mathrm{\ψ}}_y^p{|}_{i,j+1/2}=\frac{-{\mathrm{\ω}}_p^2{E}_y^p{|}_{i,j+1/2}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}(p-k){\mathrm{\ψ}}_y^k{|}_{i,j+1/2}-{\mathrm{\υ}}_{i,j+1/2}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{\mathrm{\ψ}}_y^k{|}_{i,j+1/2}}{0.25s^2+0.5{\mathrm{\υ}}_{i,j+1/2}s};$

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

步骤7，计算观测点处的电磁场分量，即计算E_x、E_y和H_z；具体计算公式为：

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

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

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

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

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

σ_η＝σ_ηmax|η-η₀|^m/d^m          (36)

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

α_η＝α_ηmax           (38)

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

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

I_mz(t)＝(t-t₀)/τ×exp(-(t-t₀)²/τ²)      (40)

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

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

以图4所示的等离子体中点源辐射的计算模型为例，按照本发明的方法步骤进行实验过程，实验二中设置加权拉盖尔多项式的阶数为300。如图4所示，整个计算区域为32×32网格，网格大小为150μm×150μm，即Δx＝Δy＝150μm，时间步长Δt＝0.25ps，等离子体充满整个计算区域，其参数ω_p＝1.803×10¹¹rad/s,υ＝2×10¹⁰rad/s，四个边界均采用5层网格的PML吸收边界。计算中所加的源位于(16，16)，其中t₀＝10ps,τ＝3ps。时间扩展因子s＝4×10¹¹，PML吸收边界参数κ_ηmax＝1，σ_ηmax＝1.15×σ_opt，α_ηmax＝0.3。观测点为(25，25)，图5为观测点的时域反射误差，其计算公式可以表示为：

$H_{\mathrm{err}}\left(t\right)=20\×{\log }_{10}\left(\frac{|H_{\mathrm{pml}}\left(t\right)-H_{\mathrm{ref}}\left(t\right)}{\max \left|H_{\mathrm{ref}}\left(t\right)\right|}\right)---\left(41\right)$

其中，H_pml为当存在CFS-PML吸收边界时，观测点的时域波形，H_ref(t)为参考波形，max|H_ref(t)|为参考波形绝对值的最大值。由图5可知，CFS-PML最大的反射误差为-68dB。因而，本发明的方法吸收效果较好。

Claims (2)

1.一种二维等离子体中扩展坐标的完全匹配吸收边界的实现方法，其特征在于，包括以下步骤：

步骤1，输入模型文件；

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

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

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

PML系数(C_1η,C_2η)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1；式中，ε₀是空气中的介电常数，s为时间尺度因子，取值为[10⁹，10¹³]；

设置的参数具体包括：

设置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_1η＝1/(κ_ηα_η+σ_η+0.5κ_ηε₀s)；

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

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

所述场源的表达式为：

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

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

式中， $C_4=\frac{0.25s^2+0.5\mathrm{\υs}}{0.25s^2+0.5\mathrm{\υs}+{\mathrm{\ω}}_p^2};$

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

${{\mathrm{\ψ}}_x^p|}_{i+1/2,j}=\frac{{-{\mathrm{\ω}}_p^2{E}_x^p|}_{i+1/2,j}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{(p-k){\mathrm{\ψ}}_x^k|}_{i+1/2,j}-{\mathrm{\υ}}_{i+1/2,j}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{{\mathrm{\ψ}}_x^k|}_{i+1/2,j}}{0.25s^2+0.5{\mathrm{\υ}}_{i+1/2,j}s};$

${{\mathrm{\ψ}}_y^p|}_{i,j+1/2}=\frac{{-{\mathrm{\ω}}_p^2{E}_y^p|}_{i,j+1/2}-s^2\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{(p-k){\mathrm{\ψ}}_y^k|}_{i,j+1/2}-{\mathrm{\υ}}_{i,j+1/2}s\underset{k=0,p>0}{\overset{p-1}{\mathrm{\Σ}}}{{\mathrm{\ψ}}_y^k|}_{i,j+1/2}}{0.25s^2+0.5{\mathrm{\υ}}_{i,j+1/2}s};$

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

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

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

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

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

步骤3.1，按照以下公式计算某些网格的磁场分量系数：

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

步骤3.2，根据步骤3.2中的公式，将整个计算区域的磁场分量系数写成矩阵方程形式： $\left[A\right]\left[H_z^p\left(r\right)\right]=\left[{\mathrm{\β}}^{p-1}\left(r\right)\right],$ 式中A表示矩阵系数，表示整个计算区域的第p阶的磁场分量系数，β^p-1(r)表示整个计算区域在p-1阶的变量，p-1阶的变量包括电磁场分量系数的和，中间变量系数的和，辅助变量；

步骤3.3，通过矩阵求逆的方式求解矩阵方程 $\left[A\right]\left[H_z^p\left(r\right)\right]=\left[{\mathrm{\β}}^{p-1}\left(r\right)\right],$ 得到整个计算区域的磁场分量系数