CN104794289A - 一种扩展直角坐标系下完全匹配吸收边界的实现方法 - Google Patents

Abstract

本发明公开了一种扩展直角坐标系下完全匹配吸收边界的实现方法，包括以下步骤：输入模型文件；初始化参数以及设置参数；添加场源到电场分量系数中，并更新计算整个计算区域的y方向上电场分量系数；更新计算整个计算区域的x方向上电场分量系数；更新计算整个计算区域的磁场分量系数；更新计算整个计算区域的电磁场分量系数的辅助变量；更新计算观测点处的电磁场分量；将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回；若达到预设值，则结束。本发明的一种扩展直角坐标系下完全匹配吸收边界的实现方法，计算速度快，且对低频与凋落波的吸收更加有效。

Description

一种扩展直角坐标系下完全匹配吸收边界的实现方法

技术领域

本发明属于计算电磁学技术领域，涉及一种扩展直角坐标系下完全匹配吸收边界的实现方法。

背景技术

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

另外，由于计算机容量的限制，电磁场的计算只能在有限区域进行。为了能模拟开域电磁波传播过程，在计算区域的截断边界处必须给出吸收边界条件。现有的直角坐标系下吸收边界主要有：Mur吸收边界，分裂场的PML(Split-field PML)，单轴各向异性PML(Uniaxial Anisotropic PML，UPML)。以上三种吸收边界均可以应用于WLP-FDTD方法的电磁场计算中，但上述吸收边界，对低频以及凋落波的吸收效果都不理想。

发明内容

本发明的目的是提供一种扩展直角坐标系下完全匹配吸收边界的实现方法，计算速度快，且对于低频和凋落波具有很好的吸收效果。

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

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δζ，ζ＝x，y；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；吸收边界层数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_h＝E_x,E_y,H_z，ζ＝y，x)、拉盖尔多项式(其中)初始化为零；

初始化PML系数(C_1ζ,C_2ζ,C_3ζ)，具体为：

C_1ζ＝1/(1+0.5ε₀s)

C_2ζ＝1

C_3ζ＝ε₀/μ₀

式中，ζ＝x，y，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

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

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

C_3ζ＝ε₀/μ₀+2α_ζ/(μ₀s)；

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

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

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

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

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

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

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

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

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

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

本发明的特点还在于，

步骤3具体为：首先给出电场分量系数在计算区域的方程，如下所示：

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

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

步骤4具体为：首先给出电场分量系数在计算区域的方程，如下所示：

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

本发明的有益效果是：

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

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

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

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

附图说明

图1是一种扩展直角坐标系下完全匹配吸收边界的实现方法的流程示意图；

图2是本发明实验中的计算模型的结构示意图；

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

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

具体实施方式

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

本发明一种扩展直角坐标系下完全匹配吸收边界的实现方法，所依据的原理为：首先导出复扩展直角坐标系下，PML中电磁场所满足的麦克斯韦方程；然后利用一种新的高速有效的WLP-FDTD方法推导出电磁场分量系数的更新方程；最后采用公式(10)中的第一式求解观测点处的电磁场分量。

在求电磁场分量系数更新方程时，首先需要推导出复扩展直角坐标系下，PML中电磁场满足的麦克斯韦方程；

复扩展直角坐标系下，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是坐标扩展变量，可以表示成

s_ζ＝k_ζ+σ_ζ/jωε₀    (3)

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

s_ζ＝k_ζ+σ_ζ/(α_ζ+jωε₀)    (4)

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

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

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

$\mathrm{j\ω}{\mathrm{\ϵ}}_0{E}_y=\frac{1}{s_x}\frac{\∂H_z}{\∂x}-J_y---\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方向的磁场。

然后，求出电磁场分量系数的更新方程；

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

${\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的变换，可以得到四组方程，这里仅给出第一个方程

$({\mathrm{\κ}}_y{\mathrm{\α}}_y+{\mathrm{\σ}}_y){\tilde{H}}_{\mathrm{zy}}+{\mathrm{\κ}}_y{\mathrm{\ϵ}}_0\frac{\∂{\tilde{H}}_{\mathrm{zy}}}{\∂t}={\mathrm{\α}}_y\frac{\∂H_z}{\∂y}+{\mathrm{\ϵ}}_0\frac{\∂}{\∂t}\left(\frac{\∂H_z}{\∂y}\right)---\left(9\right)$

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

上式中U表示电磁场分量E_x,E_y,H_z，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s>0的扩展时间，是q阶拉盖尔多项式。将(10)代入(5)、(6)和(7)中，然后让方程两边同乘以可以得到

上面三式中，b＝2ε₀/μ₀，q是加权拉盖尔多项式的阶数，D_x和D_y分别是沿x和y方向上的微分算子，和是q阶电磁场分量系数，C_iζ,i＝1,2,3；ζ＝x,y是与坐标网格有关的PML系数，计算式为：

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

C_2ζ＝1+2α_ζ/(ε₀s)    (14)

C_3ζ＝ε₀/μ₀+2α_ζ/(μ₀s)

和是拉盖尔域中电磁场分量和辅助变量的低阶和，公式如下：

(F_h＝E_x,E_y,H_z),(ζ＝y,x)

上式中的辅助变量的计算式为

将(11)、(12)和(13)写成一个矩阵形式的方程，如下：

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

式中

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

D_H＝[C_2yC_1yD_y -C_2xC_1xD_x]^T    (19)

D_E＝[C_3yC_1yD_y -C_3xC_1xD_x]

令 $W^q={\left[\begin{array}{cc}{W}_E^q& W_H^q\end{array}\right]}^T,V_{\mathrm{EH}}^{q-1}={\left[\begin{array}{cc}{V}_E^{q-1}& V_H^{q-1}\end{array}\right]}^T,A+B=\left[\begin{array}{cc}0& D_H\\ D_E& 0\end{array}\right],$ 则(17)式变为

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

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

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

上式可以分裂为下面两式

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

式中 $W^{*q}={\left[\begin{array}{cc}{W}_E^{*q}& W_H^{*q}\end{array}\right]}^T={\left[\begin{array}{ccc}{E}_x^{*q}& E_y^{*q}& H_z^{*q}\end{array}\right]}^T$ 是一个非物理中间量，为了解(23)式，使上式中I为单位对角矩阵，A和B的等式如下

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

式中

D_Ha＝[0 -C_2xC_1xD_x]^T,D_Hb＝[C_2yC_1yD_y 0]^T    (25)

D_Ea＝[0 -C_3xC_1xD_x],D_Eb＝[C_3yC_1yD_y 0]

将(24)式代入(23)式化简后得到

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

将上式展开后得到

$(1-C_{2y}{C}_{1y}{D}_y{C}_{3y}{C}_{1y}{D}_y)E_x^q=E_x^{*q}-C_{2y}{C}_{1y}{D}_y{C}_{3x}{C}_{1x}{D}_x{E}_y^{*q}---\left(29\right)$

$E_y^q=E_y^{*q}---\left(30\right)$

将(30)式代入(28)式和(31)式，(27)代入(29)得到

对上面三式进行中心差分，离散化后，得到

上面三式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；在整个计算区域上，(35)式和(36)式可以写成三对角矩阵差分方程，与传统WLP-FDTD方法相比，这种新的高速有效的WLP-FDTD方法将大型稀疏矩阵方程的求解转变成两个三对角矩阵方程的求解，于是可以使用追赶法，非常简单的解得整个计算区域电磁场分量系数，最后通过(10)式解得观测点的电磁场分量。

本发明一种扩展直角坐标系下完全匹配吸收边界的实现方法，包括以下步骤：

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δζ，ζ＝x，y；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；吸收边界层数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_h＝E_x,E_y,H_z，ζ＝y，x)、拉盖尔多项式(其中)初始化为零；

初始化PML系数(C_1ζ,C_2ζ,C_3ζ)；具体为：

C_1ζ＝1/(1+0.5ε₀s)

C_2ζ＝1

C_3ζ＝ε₀/μ₀

式中，ζ＝x，y，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

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

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

C_3ζ＝ε₀/μ₀+2α_ζ/(μ₀s)；

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

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

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

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

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

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

式中，i表示横坐标上的第i个计算网格，j表示纵坐标上的第j个计算网格；由于方程的左边有三个电场分量系数因此可以使用追赶法对整个计算区域的电场分量系数进行求解；

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

步骤4具体为：

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

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

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

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

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

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

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

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

实验：点源辐射的计算

按照本发明的方法步骤进行实施，实验中整个计算区域为50×50网格，网格大小为1cm×1cm，即Δx＝Δy＝1cm。四个边界采用8层网格的PML吸收边界，计算中所加的源位于网格(25，25)，所加场源的表达式如下：

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

其中，T_c＝5ns,T_d＝1ns。观测点位于(40，40)网格处。时间步长Δt＝23.5ps,加权拉盖尔多项式的阶数q＝120，时间扩展因子s＝2.5×10¹⁰，整个仿真时间为T_f＝17ns，PML吸收边界参数κ_ζmax＝27，σ_ζmax＝0.6×σ_opt，α_ζmax＝0.027。采用本发明方法计算的观测点处的电场分量E_y与采用传统FDFD方法计算的结果参见图3。从图3中可见，传统FDTD方法与本发明方法计算结果一致，验证了本发明方法的正确性。图4为观测点的不同吸收边界相对反射误差，其计算公式可以表示为：

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

其中，E_pml为当存在SC-PML吸收边界时，观测点的时域波形，E_ref(t)为参考波形，max|E_ref(t)|为参考波形绝对值的最大值。由图4可知，带有CFS因子的SC-PML吸收边界最大的反射误差为-50dB，它比分裂场的PML吸收边界的吸收效果更好。因此，本发明的方法吸收边界的吸收效果比分裂场的更好。

Claims (3)

1.一种扩展直角坐标系下完全匹配吸收边界的实现方法，其特征在于，包括以下步骤：

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_x×N_y，其中N_x为x方向的网格数，N_y为y方向的网格数；空间步长Δζ，ζ＝x，y；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；吸收边界层数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，ζ＝y，x)、拉盖尔多项式其中初始化为零；

初始化PML系数(C_1ζ,C_2ζ,C_3ζ)；具体为：

C_1ζ＝1/(1+0.5ε₀s)

C_2ζ＝1

C_3ζ＝ε₀/μ₀；

式中，ζ＝x，y，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

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

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

C_3ζ＝ε₀/μ₀+2α_ζ/(μ₀s)；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.根据权利要求1或2所述的一种扩展直角坐标系下完全匹配吸收边界的实现方法，其特征在于，所述步骤4具体为：

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

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