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

Abstract

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

Description

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

技术领域

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

背景技术

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

而由于计算机容量的限制，电磁场的计算只能在有限区域进行。为了能模拟开域电磁波传播过程，必须在计算区域的截断边界处给出吸收边界条件。2006年，有人将二阶Mur吸收边界应用到直角坐标下传统WLP-FDTD方法的电磁场计算中。同年，有人将单轴各向异性完全匹配吸收边界(uniaxial anisotropic perfectly matched layer absorbingboundary condition,UPML)应用到直角坐标传统WLP-FDTD方法的电磁场计算中，分析了边界的吸收效果。后来，有人将分裂场的完全匹配层(split-field PML)吸收边界用到了柱坐标下的传统WLP-FDTD方法的计算中。有人将split-field PML应用到直角坐标下新的高速有效的WLP-FDTD方法中。上述吸收边界，对低频以及凋落波的吸收效果并不理想。

发明内容

本发明的目的是提供一种扩展柱坐标系下完全匹配吸收边界的实现方法，计算速度快，且对低频与凋落波的吸收更加有效。

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

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_ρ×N_z，其中N_ρ为ρ方向的网格数，N_z为z方向的网格数；空间步长Δζ，ζ＝ρ，z；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；吸收边界层数N与相关参数κ_ζmax，σ_ζmax，α_ζmax；κ_ζmax取整数，κ_ζmax取值范围为[1，60]；α_ζmax取值范围为[0，1)；σ_ζmax/σ_opt取值范围为(0，12]；仿真计算时长T_f；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数；

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

初始化的参数具体包括：将整个计算区域的电磁场分量系数整个计算区域的电磁场分量系数的和整个计算区域的辅助变量( 其中ζ＝ρ，z)、拉盖尔多项式(其中)初始化为零；

PML系数初始化为 和式中，ζ＝ρ，z，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

设置的参数包括：

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

σ_ζ＝σ_ζmax|ζ-ζ₀|^m/d^m

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

α_ζ＝α_ζmax|ζ-ζ₀|/d

式中ζ＝ρ,z，ζ₀为PML层与非PML截面位置，d是PML吸收边界的厚度，κ_ζmax取整数，κ_ζmax取值范围为[1，60]；α_ζmax取值范围为[0，1)；σ_ζmax根据σ_opt来设置，σ_ζmax/σ_opt取值范围为(0，12]；σ_opt＝(m+1)/15πNΔζ，N为PML层数，m取值范围为[1，20]，Δζ取值范围λ为源的波长；

设置PML系数，具体为按照以下公式设置；

式中

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

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

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

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

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

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

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

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

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

本发明的特点还在于，

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

首先给出电场分量系数在轴线上的方程和非轴线上的方程，如下所示：

轴线上电场分量系数计算公式为：

非轴线上电场分量系数计算公式为：

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

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

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

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

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

本发明的有益效果是：

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

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

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

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

附图说明

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

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

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

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

具体实施方式

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

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

导出复扩展柱坐标系下，PML中电磁场满足的麦克斯韦方程的过程如下：

在复扩展柱坐标系下，PML中电磁场满足的麦克斯韦方程为

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

其中和是复扩展坐标，可以表示为

其中和为扩展复坐标，ρ₀和z₀是PML层与非PML截面位置，s_ρ和s_z是坐

标扩展变量，可以表示成

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

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

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

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

经过复坐标扩展之后，在复坐标下和是连续的(保证了在非PML与PML连接处不会发生反射)，同时在PML还有衰减项(虚数部分)。将公式(3)和(4)带入公式(2)中，又可以表示为：

本发明仅考虑二维关于的横电波情况，于是复扩展柱坐标下的麦克斯韦方程可以写成

其中E_ρ,E_z分别表示ρ,z方向的电场，分别表示方向的磁场。

导出非轴线上电磁场分量系数的更新方程如下：

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

将(6)代入(11)，然后利用jω→t的变换，可以得到五组方程，这里给出第一个方程

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

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

上面三式中，q是加权拉盖尔多项式的阶数，D_ρ和D_z分别是沿ρ和z方向上的微分算子，和是q阶电磁场分量系数，也是辅助变量，计算公式如下

上式中的计算式为

式中和是与坐标网格有关的 PML系数。(11)式中的第三个辅助变量的计算式为

上式中的计算式为

和是PML层系数，式中

将(14)、(15)和(16)写成一个矩阵形式的方程，如下：

式中

令则(21)式变为

上式中I是单位对角矩阵，A和B的等式如下

式中

通过添加微扰项方程(25)可以分裂成两个方程，如下所示：

式中将A,B,W^*q,W^q代入上式得到

将上式展开后得到

将(33)代入(31)，(34)代入(32)和(35)，得到

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

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

导出轴线上电磁场分量系数的更新方程如下：

从(39)式可以看出：在轴向上，i＝0，计算中会出现奇点，但该奇点并非麦克斯韦方程本身固有的，而是由微分方程转变为差分方程后带来的，可以通过下式来解决

对上式进行积分，然后将(13)式代入上式中，让方程两边同乘以最后进行中心差分可以得到下式

将上式与(14)式和(16)式写成一个矩阵方程，然后使用上述求解非轴线上更新方程的方法可以得到(36)式在轴线上的差分方程为

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

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_ρ×N_z，其中N_ρ为ρ方向的网格数，N_z为z方向的网格数；空间步长Δζ，ζ＝ρ，z；时间步长Δt；真空中的电导率σ，磁导率μ₀，介电常数ε₀；吸收边界层数N与相关参数κ_ζmax，σ_ζmax，α_ζmax；κ_ζmax取整数，κ_ζmax取值范围为[1，60]；α_ζmax取值范围为[0，1)；σ_ζmax/σ_opt取值范围为(0，12]；仿真计算时长T_f；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数；

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

初始化的参数具体包括：将整个计算区域的电磁场分量系数整个计算区域的电磁场分量系数的和整个计算区域的辅助变量( 其中)、拉盖尔多项式(其中)初始化为零；

PML系数初始化为 和式中，ζ＝ρ，z，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

设置的参数包括：

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

σ_ζ＝σ_ζmax|ζ-ζ₀|^m/d^m

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

α_ζ＝α_ζmax|ζ-ζ₀|/d

式中ζ＝ρ,z，ζ₀为PML层与非PML截面位置，d是PML吸收边界的厚度，κ_ζmax取整数，κ_ζmax取值范围为[1，60]；α_ζmax取值范围为[0，1)；σ_ζmax根据σ_opt来设置，σ_ζmax/σ_opt取值范围为(0，12]；σ_opt＝(m+1)/15πNΔζ，N为PML层数，m取值范围为[1，20]，Δζ取值范围λ为源的波长；

设置PML系数，具体为按照以下公式设置；

式中

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

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

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

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

首先给出电场分量系数在轴线上的方程和非轴线上的方程，如下所示：

其中，轴线上电场分量系数计算公式为：

非轴线上电场分量系数计算公式为：

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

通过观察上面两式，发现非轴线上电场分量系数方程的左边有三个电场分量系数轴线上电场分量系数方程的左边虽然只有两个电场分量系数但是可以认为另一个量为零，因此整个计算区域的电场分量系数可以写成一个矩阵系数为三对角的矩阵方程，然后使用追赶法对整个计算区域的电场分量系数进行求解；

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

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

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

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

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

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

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

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

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

实验：带有理想导体(PEC)薄板的点源辐射的计算

按照本发明的方法步骤进行实施，薄板放置在第24行0～29列网格处，添加PEC薄板的目的是为了产生凋落波，整个计算区域为50×50网格，网格大小为1cm×1cm，即Δρ＝Δz＝1cm。左边界采用轴对称边界，其余三个

边界采用15层网格的PML吸收边界，PML吸收边界的参数计算公式为：

σ_ζ＝σ_ζmax|ζ-ζ₀|^m/d^m

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

α_ζ＝α_ζmax|ζ-ζ₀|/d

σ_opt＝(m+1)/(15πNΔζ)

式中ζ＝ρ,z，ζ₀为PML层与非PML截面位置，N为PML层数，d是PML吸收边界的厚度，m是一个常数。计算中所加的源位于网格(0，25)，表达式如下：

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

其中E_pml为存在吸收边界时，观测点的时域波形，E_ref(t)为参考波形，max|E_ref(t)|为参考波形绝对值的最大值。从图4可以发现，带有CFS因子的SC-PML吸收边界最大的反射误差为-60dB。它比没有CFS因子的SC-PML吸收边界最大的反射误差要小20多个dB。因此，可以通过调整CFS因子来改善边界的吸收效果。

Claims (3)

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

步骤1，输入模型文件；

输入的模型文件具体为：计算区域大小N_ρ×N_z，其中N_ρ为ρ方向的网格数，N_z为z方向的网格数；空间步长Δζ，ζ＝ρ或ζ＝z；时间步长Δt；真空中的电导率σ，磁导率μ₀，ε₀表示真空中的介电常数；吸收边界层数N与相关参数κ_ζmax，σ_ζmax，α_ζmax；κ_ζmax取整数，κ_ζmax取值范围为[1，60]；α_ζmax取值范围为[0，1)；σ_ζmax/σ_opt取值范围为(0，12]；σ_opt＝(m+1)/150πΔζ，仿真计算时长T_f；加权拉盖尔多项式的阶数q，q≥0且为整数；时间尺度因子s，s取值范围为[10⁹，10¹³]；观测点；场源参数；

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

初始化的参数具体包括：将整个计算区域的电磁场分量系数整个计算区域的电磁场分量系数的和整个计算区域的辅助变量 拉盖尔多项式初始化为零，其中F_η＝E_ρ或F_η＝E_z或ζ＝ρ或ζ＝z，

PML系数初始化为 和式中，ζ＝ρ，z，ε₀是空气中的介电常数，s为时间尺度因子，取值范围为[10⁹，10¹³]；

设置的参数包括：

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

σ_ζ＝σ_ζmax|ζ-ζ₀|^m/d^m

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

α_ζ＝α_ζmax|ζ-ζ₀|/d

式中ζ＝ρ,z，ζ₀为PML层与非PML截面位置，d是PML吸收边界的厚度，κ_ζmax取整数，κ_ζmax取值范围为[1，60]；α_ζmax取值范围为[0，1)；σ_ζmax根据σ_opt来设置，σ_ζmax/σ_opt取值范围为(0，12]；σ_opt＝(m+1)/15πNΔζ，N为PML层数，m取值范围为[1，20]，Δζ取值范围λ为源的波长；

设置PML系数，具体为按照以下公式设置；

式中

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

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

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

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

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

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

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

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

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

2.根据权利要求1所述的一种扩展柱坐标系下完全匹配吸收边界的实现方法，其特征在于，所述步骤3中更新计算整个计算区域z方向上电场分量系数的过程具体为：

首先给出电场分量系数在轴线上的方程和非轴线上的方程，如下所示：

轴线上电场分量系数计算公式为：

非轴线上电场分量系数计算公式为：

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

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

3.根据权利要求1或2所述的一种扩展柱坐标系下完全匹配吸收边界的实现方法，其特征在于，所述步骤4中更新计算整个计算区域ρ方向上电场分量系数的过程具体为：

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

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