一种等离子体中使用辅助微分方程完全匹配层的实现方法

摘要

本发明公开了一种等离子体中使用辅助微分方程完全匹配层的实现方法，具体为：步骤1、输入模型文件；步骤2、初始化和设置步骤1模型文件参数；步骤3、利用步骤2参数计算电场分量系数步骤4、利用步骤2参数计算电场分量系数步骤5、利用步骤3和步骤4所得电场分量系数，计算磁场分量系数步骤6、利用步骤3和步骤4所得电场分量系数，计算中间变量系数步骤7、更新计算整个计算区域的电磁场分量系数的辅助变量；步骤8、更新计算观测点处的电磁场分量；步骤9、将q+1赋值给q，并判断q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。本发明计算速度快，内存消耗小，且对于低频和凋落波具有很好的吸收效果。

说明书

技术领域

本发明属于计算电磁学技术领域，具体涉及一种等离子体中使用辅助微 分方程完全匹配层的实现方法。

背景技术

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

而由于计算机容量的限制，电磁场的计算只能在有限区域进行。为了能 模拟开域电磁波传播过程，必须在计算区域的截断边界处给出吸收边界条 件。有人提出了完全匹配层(Perfectlymatchedlayer，PML)吸收边界，后 来PML被广泛应用于计算区域的截断，而且被证明是非常有效的，但是研 究发现这种传统PML对低频以及凋落波的吸收效果并不理想；使用带有复 频率偏移(Complexfrequencyshift，CFS)因子的扩展坐标的PML吸收边界 可以有效地改善传统PML对低频，凋落波与掠射情况的吸收效果。

发明内容

本发明的目的是提供一种等离子体中使用辅助微分方程完全匹配层的 实现方法，解决了现有求解磁化等离子体中的电磁问题时存在的计算速度 慢、对低频以及凋落波吸收效果不好的问题。

本发明所采用的技术方案是，一种等离子体中使用辅助微分方程完全匹 配层的实现方法，具体按照以下步骤：

步骤1、输入模型文件；

步骤2、初始化和设置步骤1模型文件的参数；

步骤3、添加场源到y方向上的电场分量系数中，并利用步骤2的参数 更新计算整个计算区域y方向上电场分量系数

步骤4、利用步骤2的参数更新计算整个计算区域的x方向上电场分量 系数

步骤5、利用步骤3和步骤4所得电场分量系数，更新计算整个计算区 域的磁场分量系数

步骤6、利用步骤3和步骤4所得电场分量系数，分别更新计算整个计 算区域的中间变量的系数

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

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

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

本发明的特点还在于：

步骤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]，σ_opt＝(m+1)/150πΔη，m取值范围为[1，20]， Δη取值范围为λ为源的波长；仿真计算时长T_f；加权拉盖尔多项 式的阶数q，其中q≥0且为整数；时间尺度因子s，其中s取值范围为 [10⁹,10¹³]；观测点；场源参数T_c，T_d。

步骤2初始化和设置步骤1模型文件的参数，具体为：

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

PML系数(C_1η，C_2η，C₃)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1，C₃＝ε₀/μ₀；

设置的参数具体包括：

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

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

式中，η＝x,y，η₀为PML层与非PML截面位置，d是PML吸收边界的厚度；

设置PML系数C_1η，C_2η；具体为：

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

步骤3具体为：

所添加场源的表达式为：

$J_y\left(t\right)=(t-T_c)/T_d\exp (-{(t-T_c)}^2/T_d^2)$

式中，T_c，T_d为场源参数；

步骤3.1、利用步骤2的参数计算电场分量系数在计算区域的方程为：

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

步骤3.2、使用追赶法求解步骤3.1的方程，得到整个计算区域y方向上 的电场分量系数

步骤4具体为：

步骤4.1、利用步骤2的参数计算电场分量系数在计算区域的方程为：

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

步骤5具体为：

步骤6具体更新公式为：

${\psi }_x^q{|}_{i+1/2,j}=\frac{-{\omega }_p^2{E}_x^q{|}_{i+1/2,j}-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_x^k{|}_{i+1/2,j}-{\upsilon }_{i+1/2,j}s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_x^k{|}_{i+1/2,j}}{0.25s^2+0.5{\upsilon }_{i+1/2,j}s}$

${\psi }_y^q{|}_{i,j+1/2}=\frac{-{\omega }_p^2{E}_y^q{|}_{i,j+1/2}-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_y^k{|}_{i,j+1/2}-{\upsilon }_{i,j+1/2}s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_y^k{|}_{i,j+1/2}}{0.25s^2+0.5{\upsilon }_{i,j+1/2}s}$

步骤7具体更新公式为：

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

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

本发明的有益效果是：

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

2)在求解电磁场分量系数时，将大型稀疏矩阵方程分裂成两个三对角 矩阵方程，使得它在计算时简单、计算速度更快、内存消耗更少；

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

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

附图说明

图1是本发明的流程示意图；

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

图3是本发明实施例中等离子体附近的非均匀网格示意图；

图4是本发明实施例中反射系数与透射系数振幅随频率的变化图；

图5是本发明实施例中反射系数与透射系数相位随频率的变化图。

具体实施方式

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

本发明一种等离子体中使用辅助微分方程完全匹配层的实现方法，磁化 等离子体中，电磁场所满足的复扩展坐标系下的麦克斯韦方程推导过程如 下：

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

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

$\widetilde{▿}=\hat{x}\frac{\partial }{s_x\partial x}+\hat{y}\frac{\partial }{s_y\partial y}+\hat{z}\frac{\partial }{s_z\partial 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\omega ϵE}}_x=\frac{1}{s_y}\frac{\partial H_z}{\partial y}---\left(5\right)$

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

${\mathrm{j\omega \mu }}_0{H}_z=\frac{1}{s_y}\frac{\partial E_x}{\partial y}-\frac{1}{s_x}\frac{\partial E_y}{\partial x}---\left(7\right)$

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

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

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

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

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

$\begin{array}{c}({\kappa }_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta }){\tilde{F}}_{\zeta \eta }+{\kappa }_{\eta }{ϵ}_0\frac{\partial {\tilde{F}}_{\zeta \eta }}{\partial t}={\alpha }_{\eta }\frac{\partial F_{\zeta }}{\partial \eta }+{ϵ}_0\frac{\partial }{\partial t}\left(\frac{\partial F_{\zeta }}{\partial \eta }\right)\\ (F_{\zeta }=E_x,E_y,H_z),(\eta =y,x)\end{array}---\left(9\right)$

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

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

式中，

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

频率形式的电位移矢量为

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

${ϵ}_r\left(\omega \right)=1-\frac{{\omega }_p^2}{{\omega }^2-j\omega \upsilon }---\left(15\right)$

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

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

即

令

则

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

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

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

${\psi }^q=\frac{-{\omega }_q^2{E}^q-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }^k-\upsilon s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }^k}{0.25s^2+0.5\upsilon s}---\left(23\right)$

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

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

${\psi }_x^q{|}_{i+1/2,j}=\frac{-{\omega }_p^2{E}_x^q{|}_{i+1/2,j}-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_x^k{|}_{i+1/2,j}-{\upsilon }_{i+1/2,j}s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_x^k{|}_{i+1/2,j}}{0.25s^2+0.5{\upsilon }_{i+1/2,j}s}---\left(24\right)$

${\psi }_y^q{|}_{i,j+1/2}=\frac{-{\omega }_p^2{E}_y^q{|}_{i,j+1/2}-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_y^k{|}_{i,j+1/2}-{\upsilon }_{i,j+1/2}s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_y^k{|}_{i,j+1/2}}{0.25s^2+0.5{\upsilon }_{i,j+1/2}s}---\left(25\right)$

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

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

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

将式(10)代入(26)、(27)和(7)中，应用加勒金测试过程分别得：

将(28)、(29)和(30)写成矩阵的形式如下：

$\left(\begin{array}{c}{W}_E^q\\ W_H^q\end{array}\right)=\left(\begin{array}{c}{V}_E^{q-1}\\ V_H^{q-1}\end{array}\right)+\left(\begin{array}{cc}0& D_H\\ D_E& 0\end{array}\right)\left(\begin{array}{c}{W}_E^q\\ W_H^q\end{array}\right)---\left(31\right)$

式中

$W_E^q=\left(\begin{array}{c}{E}_x^q\\ E_y^q\end{array}\right),W_H^q=\left[H_z^q\right]---\left(32\right)$

$D_H=\left(\begin{array}{c}{C}_4{C}_{1y}{C}_{2y}{D}_y\\ -C_4{C}_{1x}{C}_{2x}{D}_x\end{array}\right),D_E=\left[\begin{array}{cc}{C}_3{C}_{2y}{C}_{1y}{D}_y& -C_3{C}_{2x}{C}_{1x}{D}_x\end{array}\right]---\left(33\right)$

式中

$P\_{\mathrm{temp}}_x=\frac{-s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_x^k+0.5s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_x^k}{0.25s+0.5\upsilon }---\left(35\right)$

$P\_{\mathrm{temp}}_y=\frac{-s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_y^k+0.5s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_y^k}{0.25s+0.5\upsilon }---\left(36\right)$

如果让那么式(31)写为：

$(I-A-B)W^q=V_{EH}^{q-1}---\left(37\right)$

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

$(I-A)(I-B)W^q=V_{EH}^{q-1}+{\mathrm{ABV}}_{EH}^{q-1}---\left(38\right)$

上式可以分裂为下面两式

$\left(\begin{array}{c}(I-A)W^{*q}=(I+B)V_{EH}^{q-1}\\ (I-B)W^q=W^{*q}-{\mathrm{BV}}_{EH}^{q-1}\end{array}\right)---\left(39\right)$

式中是非物理中间量，为解式(39)，令

$A=\left(\begin{array}{cc}0& D_{Ha}\\ D_{Ea}& 0\end{array}\right),B=\left(\begin{array}{cc}0& D_{Hb}\\ D_{Eb}& 0\end{array}\right)---\left(40\right)$

$\begin{array}{c}{D}_{Ha}={\left(\begin{array}{cc}0& -C_4{C}_{1x}{C}_{2x}{D}_x\end{array}\right)}^T,D_{Hb}={\left(\begin{array}{cc}{C}_4{C}_{1y}{C}_{2y}{D}_y& 0\end{array}\right)}^T\\ D_{Ea}=\left[\begin{array}{cc}0& -C_3{C}_{2x}{C}_{1x}{D}_x\end{array}\right],D_{Eb}=\left[\begin{array}{cc}{C}_3{C}_{2y}{C}_{1y}{D}_y& 0\end{array}\right]\end{array}---\left(41\right)$

将式(40)代入式(39)得到

$\left(\begin{array}{cc}I& -D_{Ha}\\ -D_{Ea}& I\end{array}\right)\left(\begin{array}{c}{W}_E^{*q}\\ W_H^{*q}\end{array}\right)=\left(\begin{array}{cc}I& D_{Hb}\\ D_{Eb}& I\end{array}\right)\left(\begin{array}{c}{V}_E^{q-1}\\ V_H^{q-1}\end{array}\right)$

$\left(\begin{array}{cc}I& -D_{Hb}\\ -D_{Eb}& I\end{array}\right)\left(\begin{array}{c}{W}_E^q\\ W_H^q\end{array}\right)=\left(\begin{array}{c}{W}_E^{*q}\\ W_H^{*q}\end{array}\right)-\left(\begin{array}{cc}0& D_{Hb}\\ D_{Eb}& 0\end{array}\right)\left(\begin{array}{c}{V}_E^{q-1}\\ V_H^{q-1}\end{array}\right)$

$\left(\begin{array}{c}{W}_E^{*q}-D_{Ha}{W}_H^{*q}=V_E^{q-1}+D_{Hb}{V}_H^{q-1}\\ -D_{Ea}{W}_E^{*q}+W_H^{*q}=D_{Eb}{V}_E^{q-1}+V_H^{q-1}\\ W_E^q-D_{Hb}{W}_H^q=W_E^{*q}-D_{Hb}{V}_H^{q-1}\\ -D_{Eb}{W}_E^q+W_H^q=W_H^{*q}-D_{Eb}{V}_E^{q-1}\end{array}\right)---\left(42\right)$

$\left(\begin{array}{c}(I-D_{Ha}{D}_{Ea})W_E^{*q}=(I+D_{Ha}{D}_{Eb})V_E^{q-1}+(D_{Hb}+D_{Ha})V_H^{q-1}\\ (I-D_{Hb}{D}_{Eb})W_E^q=(I+D_{Hb}{D}_{Ea})W_E^{*q}\\ W_H^q=D_{Eb}{W}_E^q+D_{Ea}{W}_E^{*q}+V_H^{q-1}\end{array}\right)---\left(43\right)$

$\left(\begin{array}{cc}1-C_4{C}_{1y}{C}_{2y}{D}_y{C}_3{C}_{2y}{C}_{1y}{D}_y& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}{E}_x^q\\ E_y^q\end{array}\right)=\left(\begin{array}{cc}1& -C_4{C}_{1y}{C}_{2y}{D}_y{C}_3{C}_{2x}{C}_{1x}{D}_x\\ 0& 1\end{array}\right)\left(\begin{array}{c}{E}_x^{*q}\\ E_y^{*q}\end{array}\right)---\left(45\right)$

将上式扩展得到

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

对上式进行中心差分得到

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

本发明一种等离子体中使用辅助微分方程完全匹配层的实现方法，流程

如图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]，σ_opt＝(m+1)/150πΔη，m取值范围为[1，20]， Δη取值范围为λ为源的波长；仿真计算时长Tf；加权拉盖尔多项 式的阶数q，其中q≥0且为整数；时间尺度因子s，其中s取值范围为 [10⁹,10¹³]；观测点；场源参数T_c，T_d。

步骤2初始化和设置步骤1模型文件的参数，具体为：

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

PML系数(C_1η，C_2η，C₃)初始化为C_1η＝1/(1+0.5ε₀s)，C_2η＝1，C₃＝ε₀/μ₀；

设置的参数具体包括：

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

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

式中，η＝x,y，η₀为PML层与非PML截面位置，d是PML吸收边界的厚度；

设置PML系数C_1η，C_2η；具体为：

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

步骤3添加场源到y方向上的电场分量系数中，并利用步骤2的参数更 新计算整个计算区域y方向上电场分量系数具体为：

所添加场源的表达式为：

$J_y\left(t\right)=(t-T_c)/T_d\exp (-{(t-T_c)}^2/T_d^2)$

式中，T_c，T_d为场源参数；

步骤3.1、利用步骤2的参数计算电场分量系数在计算区域的方程为：

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

步骤3.2、使用追赶法求解步骤3.1的方程，得到整个计算区域y方向上 的电场分量系数

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

步骤4.1、利用步骤2的参数计算电场分量系数在计算区域的方程为：

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

步骤5利用步骤3和步骤4所得电场分量系数，更新计算整个计算区域 的磁场分量系数具体为：

步骤6利用步骤3和步骤4所得电场分量系数，分别更新计算整个计算 区域的中间变量的系数具体更新公式为：

${\psi }_x^q{|}_{i+1/2,j}=\frac{-{\omega }_p^2{E}_x^q{|}_{i+1=2,j}-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_x^k{|}_{i+1/2,j}-{\upsilon }_{i+1/2,j}s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_x^k{|}_{i+1/2,j}}{0.25s^2+0.5{\upsilon }_{i+1/2,j}s}$

${\psi }_y^q{|}_{i,j+1/2}=\frac{-{\omega }_p^2{E}_y^q{|}_{i,j+1/2}-s^2\underset{k=0,q>0}{\overset{q-1}{\Sigma }}(q-k){\psi }_y^k{|}_{i,j+1/2}-{\upsilon }_{i,j+1/2}s\underset{k=0,q>0}{\overset{q-1}{\Sigma }}{\psi }_y^k{|}_{i,j+1/2}}{0.25s^2+0.5{\upsilon }_{i,j+1/2}s}$

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

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

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

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

实施例

选取厚度为6mm的非均匀等离子体平板，计算其对电磁波的反射与透 射系数。仿真模型如图2所示。计算区域为180×36个网格，沿着x方向从第 31到第150的网格为非均匀等离子体，其中等离子体1占据了第31到第50 的网格以及第131到150的网格，厚度皆为1.5毫米，其等离子体参数为 ω_p＝1.80327×10¹¹rad/s，υ＝2×10¹⁰rad/s，等离子体2占据了第51到第130 的网格，厚度为3毫米，其等离子体参数为ω_p＝5.6985×10¹⁰rad/s， υ＝2×10¹⁰rad/s。使用渐变的非均匀网格对求解区域进行网格划分，非均匀 等离子体附近的网格划分如图3所示，沿x方向等离子体2的网格尺寸为 37.5μm，等离子体1的网格尺寸为75μm，其他计算区域网格尺寸为150μm， 最小网格大小为37.5μm×37.5μm。源加到电场分量E_y上，其表达式为

$J_y\left(t\right)=(t-T_c)/T_d\exp (-{\left(t-T_c\right)}^2/T_d^2)$

式中t_c＝0.05ns,t_d＝0.01ns。WLP-FDTD方法中，我们将时间步长设为 Δt＝3.5355ps，仿真总时间T_f＝424.26ps，WLPs的阶数为100，时间尺度因子 为s＝5.18×10¹¹。对于传统FDTD方法，时间步长受柯西稳定性条件的限制， 我们选择Δt＝62.5fs，时间步数为6789。两个观测点放置在左边和右边的等 离子体平板处，采用本专利中的发明方法、传统FDTD方法和WLP-FDTD 方法计算出观测点处的时域电场值，然后将其转换到频域上计算出反射和透 射系数，图4和图5分别显示出不同方法的振幅和相位，从图可以看出它们 吻合的非常好，验证了本发明方法的的正确性。两种WLP-FDTD方法的计 算效率如表1所示：

表1

与WLP-FDTD方法相比，本专利中的方法的计算时间大大减小。在求 解等离子体中的电磁场传播过程中，两种WLP-FDTD方法所用到的变量的 数目是一致的，但是由于WLP-FDTD方法需要求解大型稀疏矩阵方程，内 存消耗需要更大，因此本专利中的方法比WLP-FDTD所需内存更少。

本发明一种等离子体中使用辅助微分方程完全匹配层的实现方法法，在 直角坐标系下，通过用加权拉盖尔多项式表示电磁场分量，来解时域麦克斯 韦方程，使得在更新计算整个计算区域的电磁场分量系数时不涉及到时间步 长，只是在最后计算观测点处的电磁场分量时用到时间步长，因此计算过程 中时间步长可以取得比柯西稳定性条件限制的时间步长更大；在求解电磁场 分量系数时，将大型稀疏矩阵方程分裂成两个三对角矩阵方程，使得它在计 算时比WLP-FDTD方法更简单、计算速度更快、内存消耗更少而且可以对 大区域的电磁场问题进行求解；在设置PML系数时，由于采用了CFS因子， 并且通过调整CFS因子中的参数，可以使得该吸收边界对低频与凋落波的吸 收更加有效；由于采用了复扩展坐标系，使得PML在实现时避免了场的分 裂且与媒质无关。