CN104375975A - 基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法 - Google Patents

Abstract

本发明提出了一种基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法，属于数值仿真技术领域，该方法的目的是将计算机有限的内存空间仿真成无限空间。本发明的技术特征是：通过利用双线性变换方法将复数拉伸坐标变量由频域变换到z域，然后利用Crank-Nicolson时域有限差分方法将麦克斯韦方程在时域进行离散，推导出电场的显式迭代方程，最后求解出电磁场分量的值，本发明具有无条件稳定性，提高电磁场计算速度和节约内存的优点。

Description

基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法

技术领域

本发明涉及数值仿真技术领域，特别涉及一种基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法。

背景技术

近年来，时域有限差分方法(FDTD)作为一种计算电磁方法被广泛地应用于各种时域的电磁仿真计算中，如天线，射频电路，光学器件和半导体等等。FDTD具有广泛的适用性、适合并行计算、计算程序通用性等特点。

然而，随着科学研究的深入和各种越来越广泛应用的需求，其算法本身受Courant FriedrichsLewy(CFL)数值稳定性条件的限制的缺陷越来越凸显出来。算法本身所受数字稳定性条件限制：在计算过程中时间步长和空间步长必须满足CFL约束条件，即其中，Δt为计算时间步长，c₀为计算介质中的光速，Δx、Δy和Δz为三维空间步长。在实际计算中，空间离散步长和时间步长相对波长和周期都非常小，所以必然会在计算电大尺寸目标时出现资源不足的情况，导致FDTD的计算效率很低。因此为了消除CFL条件的限制，无条件稳定的交替方向隐式(Alternating-Direction Impolicit，ADI)FDTD方法、局部一维(LocalOne Dimension，LOD)FDTD方法和克兰克·尼克尔森(Crank-Nicolson，CN)FDTD方法相继被提出。

对于ADI-FDTD算法和LOD-FDTD算法虽然在一定程度上克服了稳定性条件限制，但算法的计算精度过低，性能并不理想，其原因是由于当时间步长增大后，导致的数值色散增大，进而导致算法的误差较大。2004年，G.Sun等人采用Crank-Nicolson差分格式对麦克斯韦方程进行离散化处理，即CN-FDTD，算法在时间步长取值远大于稳定性条件(如20倍)仍能保持良好的稳定精度，展现出更好的适用性，并且CN-FDTD算法是一种更加简便的无条件稳定的方法，将前面两种算法中所需的2个运算过程简化到1个运算过程，从而大大降低了运算资源，因此学者们一致认为CN-FDTD具有更广阔的发展前景。

由于计算机内存空间的限制，数值计算只能在有限的区域内进行，为了能模拟开放或者半开放区域的电磁辐射和散射等问题，在计算区域的截断边界处必须设置吸收边界条件，以便用有限的网格空间模拟开放的无限空间，来解决任意介质内的电磁波传播以及各种电磁问题。由Berenger提出的完全匹配层(Perfectly Matched Layer，PML)是目前应用较广的吸收边界条件，PML可以理解为：通过在FDTD区域截断边界处设置一种特殊介质层，该层介质的波阻抗与相邻介质波阻抗完全匹配，从而使入射波无反射地穿过分界面而进入PML层，PML层是有耗介质，最后将电磁波吸收。目前常用的PML吸收边界主要有拉伸坐标变换完全匹配层(SC-PML)和单轴各项异性完全匹配层(UPML)。2002年，Ramadan利用双线性变换(Bilinear-Transform)方法将复数拉伸坐标变量变换到z域。利用双线性变换方法具有离散简单的优点。

发明内容

本发明的目的是针对FDTD算法受到CFL稳定性条件限制的缺陷，提高CN-PML算法的计算效率和吸收效率而提出的基于双线性变换方法和CN-FDTD的SC-PML算法。

为实现上述目的，基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法的具体步骤方法如下：

步骤1：将频域中麦克斯韦旋度方程修正为带有拉伸坐标算子的麦克斯韦方程。

真空条件下，在SC-PML介质中的归一化麦克斯韦方程描述为：

$\mathrm{j\ωE}=c_0{\▿}_s\×H---\left(1\right)$

$\mathrm{j\ωH}=-c_0{\▿}_s\×E---\left(2\right)$

式中，c₀是真空中的光速，为拉伸坐标算子，定义为

${\▿}_s=\hat{x}{S}_x^{-1}{\∂}_x+\hat{y}{S}_y^{-1}{\∂}_y+\hat{z}{S}_z^{-1}{\∂}_z---\left(3\right)$

式中，是关于x，y，z方向的偏导数。

步骤2：将修正后的频域中一维麦克斯韦旋度方程在直角坐标系中表示。

一维真空中，SC-PML中的z方向极化、x方向传播的TEM波，归一化的麦克斯韦方程在频域中表示为

$\mathrm{j\ω}{E}_z=c_0{S}_x^{-1}\frac{\∂}{\∂x}{H}_y---\left(4\right)$

$\mathrm{j\ω}{H}_y=c_0{S}_x^{-1}\frac{\∂}{\∂x}{E}_z---\left(5\right)$

式中，拉伸坐标变量S_x表示为

$S_x^{-1}=\frac{1}{1+\frac{{\mathrm{\σ}}_x}{\mathrm{j\ω}{\mathrm{\ϵ}}_0}}---\left(6\right)$

式中，σ_x是PML层中沿x方向的电导率。

步骤3：根据频域和z域的映射关系，将直角坐标系中的一维麦克斯韦方程变换到z域表示。

$\frac{1-z^{-1}}{\mathrm{\Δt}}{E}_z=c_0\·S_x\left(z\right)\frac{\∂}{\∂x}{H}_y---\left(7\right)$

$\frac{1-z^{-1}}{\mathrm{\Δt}}{H}_y=c_0\·S_x\left(z\right)\frac{\∂}{\∂x}{E}_z---\left(8\right)$

式中，利用了映射关系，S_x(z)为的z域表示。

步骤4：根据双线性变换方法，将拉伸坐标变量由频域变换到z域，并代入到一维麦克斯韦方程的z域表达式，同时设置辅助变量。

根据双线性变换关系

$\mathrm{j\ω}\⇒\frac{2}{\mathrm{\Δt}}\frac{1-z^{-1}}{1+z^{-1}}---\left(9\right)$

将变换到z域得

$S_x\left(z\right)=b_x\·\frac{1-z^{-1}}{1-a_x{z}^{-1}}---\left(10\right)$

式中， $b_x=\frac{1}{1+\mathrm{\Δt}{\mathrm{\σ}}_x/\left(2{\mathrm{\ϵ}}_0\right)},a_x=b_x(1-\frac{\mathrm{\Δt}{\mathrm{\σ}}_x}{2{\mathrm{\ϵ}}_0}),$ Δt为计算时间步长。

将式(10)代入到式(7)和式(8)中得

$\frac{1-z^{-1}}{\mathrm{\Δt}}{E}_z=c_0\·b_x\·\frac{1-z^{-1}}{1-a_x{z}^{-1}}\frac{\∂}{\∂x}{H}_y---\left(11\right)$

$\frac{1-z^{-1}}{\mathrm{\Δt}}{H}_y=c_0\·b_x\·\frac{1-z^{-1}}{1-a_x{z}^{-1}}\frac{\∂}{\∂x}{E}_z---\left(12\right)$

定义两个辅助变量

$F_{\mathrm{zx}}=\frac{1}{1-a_x{z}^{-1}}\frac{\∂}{\∂x}{H}_y=a_x{z}^{-1}{F}_{\mathrm{zx}}+\frac{\∂}{\∂x}{H}_y---\left(13\right)$

$G_{\mathrm{yx}}=\frac{1}{1-a_x{z}^{-1}}\frac{\∂}{\∂x}{E}_z=a_x{z}^{-1}{G}_{\mathrm{yx}}+\frac{\∂}{\∂x}{E}_z---\left(14\right)$

则式(11)和式(12)变为

$\frac{1-z^{-1}}{\mathrm{\Δt}}{E}_z=c_0\·b_x\·(F_{\mathrm{zx}}-z^{-1}{F}_{\mathrm{zx}})---\left(15\right)$

$\frac{1-z^{-1}}{\mathrm{\Δt}}{H}_y=c_0\·b_x\·(G_{\mathrm{yx}}-z^{-1}{G}_{\mathrm{yx}})---\left(16\right)$

将式(13)和式(14)代入到式(15)和式(16)中得

$\frac{1-z^{-1}}{\mathrm{\Δt}}{E}_z=c_0\·b_x\·((a_x-1)z^{-1}{F}_{\mathrm{zx}}+\frac{\∂}{\∂x}{H}_y)---\left(17\right)$

$\frac{1-z^{-1}}{\mathrm{\Δt}}{H}_y=c_0\·b_x\·((a_x-1)z^{-1}{G}_{\mathrm{yx}}+\frac{\∂}{\∂x}{E}_z)---\left(18\right)$

步骤5：基于Crank-Nicolson时域有限差分算法的时域展开形式，以及根据z域和时域的映射关系，将z域形式的直角坐标系中一维麦克斯韦旋度方程展开成时域有限差分的形式，同时也将z域形式的辅助变量变换为时域有限差分的形式。

由式(17)和式(18)，可得

$\begin{array}{c}\frac{E_z^{n+1}\left(i\right)-E_z^n\left(i\right)}{\mathrm{\Δt}}=c_0\·b_x\left(i\right)\·(a_x\left(i\right)-1)F_{\mathrm{zx}}^n\left(i\right)\\ +\frac{c_0\·b_x\left(i\right)}{2\mathrm{\Δx}}[H_y^{n+1}(i+\frac{1}{2})-H_y^{n+1}(i-\frac{1}{2})+H_y^n(i+\frac{1}{2})-H_y^n(i-\frac{1}{2})]\end{array}---\left(19\right)$

$\begin{array}{c}\frac{H_y^{n+1}(i+\frac{1}{2})-H_y^n(i+\frac{1}{2})}{\mathrm{\Δt}}=c_0\·b_x(i+\frac{1}{2})\·(a_x(i+\frac{1}{2})-1)G_{\mathrm{yx}}^n\\ +c_0\·b_x(i+\frac{1}{2})\frac{1}{2\mathrm{\Δx}}[E_z^{n+1}(i+1)-E_z^{n+1}\left(i\right)+E_z^n(i+1)-E_z^n\left(i\right)]\end{array}---\left(20\right)$

辅助变量的离散方程为

$F_{\mathrm{zx}}^{n+1}\left(i\right)=a_x\left(i\right)F_{\mathrm{zx}}^n\left(i\right)$

$+\frac{1}{2\mathrm{\Δx}}[H_y^{n+1}(i+\frac{1}{2})-H_y^{n+1}(i-\frac{1}{2})+H_y^n(i+\frac{1}{2})-H_y^n(i-\frac{1}{2})]---\left(21\right)$

$\begin{array}{c}{G}_{\mathrm{yx}}^{n+1}(i+\frac{1}{2})=a_x(i+\frac{1}{2})G_{\mathrm{yx}}^n(i+\frac{1}{2})\\ +\frac{\∂}{2\mathrm{\Δx}}[E_z^{n+1}(i+1)-E_z^{n+1}\left(i\right)+E_z^n(i+1)-E_z^n\left(i\right)\end{array}---\left(22\right)$

步骤6：将时域有限差分的形式整理成求解的形式。结果产生一组电场和磁场耦合的方程，是一组隐式方程。

$\begin{array}{c}{E}_z^{n+1}\left(i\right)=E_z^n\left(i\right)+m_{1x}\left(i\right)F_{\mathrm{zx}}^n\left(i\right)\\ +m_{2x}\left(i\right)[H_y^{n+1}(i+\frac{1}{2})-H_y^{n+1}(i-\frac{1}{2})+H_y^n(i+\frac{1}{2})-H_y^n(i-\frac{1}{2})]\end{array}---\left(23\right)$

$\begin{array}{c}{H}_y^{n+1}(i+\frac{1}{2})=H_y^n(i+\frac{1}{2})+m_{1x}(i+\frac{1}{2})G_{\mathrm{yx}}^n(i+\frac{1}{2})\\ +m_{2x}(i+\frac{1}{2})[E_z^{n+1}(i+1)-E_z^{n+1}\left(i\right)+E_z^n(i+1)-E_z^n\left(i\right)]\end{array}---\left(24\right)$

式中，m_1x(i)＝c₀Δt·b_x(i)·(a_x(i)-1)，

步骤7：将这组隐式方程进行去耦，即将磁场分量的方程代入到电场分量的方程中，同时将辅助变量的显式方程组代入，得到

$E_z^{n+1}\left(i\right)=E_z^n\left(i\right)+m_{1x}\left(i\right)F_{\mathrm{zx}}^n\left(i\right)$

$+m_{2x}\left(i\right)\left[\begin{array}{c}{H}_y^n(i+\frac{1}{2})\\ +m_{1x}(i+\frac{1}{2})G_{\mathrm{yx}}^n(i+\frac{1}{2})\\ +m_{2x}(i+\frac{1}{2})[E_z^{n+1}(i+1)-E_z^{n+1}\left(i\right)+E_z^n(i+1)-E_z^n\left(i\right)]\\ -H_y^n(i-\frac{1}{2})\\ -m_{1x}(i-\frac{1}{2})G_{\mathrm{yx}}^n(i-\frac{1}{2})\\ -m_{2x}(i-\frac{1}{2})[E_z^{n+1}\left(i\right)-E_z^{n+1}(i-1)+E_z^n\left(i\right)-E_z^n(i-1)]\\ +H_y^n(i+\frac{1}{2})-H_y^n(i-\frac{1}{2})\end{array}\right]---\left(25\right)$

步骤8：将被代入磁场后的电场的方程进行整理，整理后获得等式左边为三对角矩阵形式的系数的电场显式迭代方程。

$\begin{array}{c}-m_{2x}\left(i\right)m_{2x}(i-\frac{1}{2})E_z^{n+1}(i-1)\\ +[1+m_{2x}\left(i\right)(m_{2x}(i-\frac{1}{2})+m_{2x}(i+\frac{1}{2}))]E_z^{n+1}\left(i\right)\\ -m_{2x}\left(i\right)m_{2x}(i+\frac{1}{2})E_z^{n+1}(i+1)\\ =m_{2x}\left(i\right)m_{2x}(i-\frac{1}{2})E_z^n(i-1)\\ +[1-m_{2x}\left(i\right)(m_{2x}(i-\frac{1}{2})+m_{2x}(i+\frac{1}{2}))]E_z^n\left(i\right)+m_{2x}\left(i\right)m_{2x}(i+\frac{1}{2})E_z^n(i+1)\\ +m_{1x}\left(i\right)F_{\mathrm{zx}}^n\left(i\right)+m_{2x}\left(i\right)m_{1x}(i+\frac{1}{2})G_{\mathrm{yx}}^n(i+\frac{1}{2})-m_{2x}\left(i\right)m_{1x}(i-\frac{1}{2})G_{\mathrm{yx}}^n(i-\frac{1}{2})\\ +2m_{2x}\left(i\right)(H_y^n(i+\frac{1}{2})-H_y^n(i-\frac{1}{2}))\end{array}---\left(26\right)$

步骤9：利用追赶法求解系数为三对角矩阵的电场迭代方程，得到电场分量。

步骤10：将求解出的电场值代入到磁场的迭代方程中，求解出磁场分量。

步骤11：将求解出的电场和磁场值代入到辅助变量的迭代方程中，求解出辅助变量的值，返回到步骤9，循环步骤9-11，从而在时间上迭代求解。

附图说明：

为了更清楚地说明发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造劳动性的前提下，还可以根据附图获得其他的附图。

图1本发明流程框图

图2空间模型图

图3相对反射误差图

具体实施方式：

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

步骤1：如图1，根据算法流程图用FORTRAN语言编写程序，并在程序中给出具体算例，设定计算空间模型的电磁参数，其中包括计算区域大小：L+2×N，其中L为真空长度(单位：元胞)，N为完全匹配层的厚度(单位：元胞)；空间步长Δx；时间步长为其中为传统FDTD的时间步长，CFLN为CN-FDTD的时间步长相对于传统FDTD的时间步长的倍数。

步骤2：如图2，根据空间模型图设置激励源位置和激励源的类型，本实施例以微分高斯脉冲作为激励源。

步骤3：在空间中设置观察点(距离PML层一个Δx)，将空间中电磁波的数据写入外部文件作为仿真解，并通过扩大计算域的方法(一维情况下，将计算域向两侧分别扩大10000Δx，使计算域足够大，使在既定的计算时间内没有反射波出现)，在所设置的观察点不变的情况下获得参考解，将仿真解和参考解代入到反射误差的公式(27)中计算相对反射误差，来确定完全匹配层对电磁波的吸收效果。

$R_{\mathrm{dB}}\left(t\right)=20{\log }_{10}\left(\right|\frac{E_z^T\left(t\right)-E_z^R\left(t\right)}{E_{z\_\max }^R\left(t\right)}\left|\right)---\left(\text{27}\right)$

式中，代表在观察点处的时域离散电场分量E_z(t)的仿真解，代表参考解，代表参考解在仿真时间内的最大值。参考解是通过将FDTD计算域向外延伸到10000Δx，即计算域足够大，在仿真时间内在观察点处不存在外边界反射波。所得到的反射误差如图3所示，结果显示相对反射误差的最大值随着CNLN的增加保持不变，说明该算法具有无条件稳定性。

本领域技术人员可以理解附图只是一个优先实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法，其特征在于具体设置步骤： 

步骤1：将频域中麦克斯韦旋度方程修正为带有拉伸坐标算子的麦克斯韦方程； 

步骤2：将频域中修正后的一维麦克斯韦旋度方程在直角坐标系中表示； 

步骤3：根据频域和z域的映射变换关系，将直角坐标系中的一维麦克斯韦方程变换到z域表示； 

步骤4：根据双线性变换方法，将拉伸坐标变量由频域变换到z域，并代入到一维麦克斯韦方程的z域表达式，同时设置辅助变量； 

步骤5：基于Crank-Nicolson时域有限差分算法的时域展开形式，以及根据z域和时域的映射关系，将z域形式的直角坐标系中一维麦克斯韦旋度方程展开成时域有限差分的形式，同时也将z域形式的辅助变量变换为时域有限差分的形式； 

步骤6：将时域有限差分形式的方程整理成求解的形式，结果产生一组电场和磁场耦合的方程，是一组隐式方程； 

步骤7：将这组隐式方程进行去耦，即将磁场分量的方程代入到电场分量的方程中，同时将辅助变量的显式方程组代入； 

步骤8：将代入磁场和辅助变量后的电场分量的迭代方程进行整理，整理后获得等式左边为三对角矩阵形式的系数的电场显式迭代方程； 

步骤9：利用追赶法求解系数为三对角矩阵的电场迭代方程，得到电场分量的值； 

步骤10：将求解出的电场值代入到磁场的迭代方程中，求解出磁场分量的值； 

步骤11：将求解出的电场值和磁场值代入到辅助变量的迭代方程中，求解出辅助变量的值；返回到步骤9，循环步骤9-11，从而在时间上迭代求解。 

2.由权利1所述的基于双线性变换的一维真空Crank-Nicolson完全匹 配层实现算法，其特征在于：步骤4，利用双线性变换方法将拉伸坐标变量由频域变换到z域过程： 

对于一维PML介质空间，z方向极化、x方向传播的麦克斯韦方程中的拉伸坐标变量的频域表达式为 

式中，σ_x是PML层中沿x方向的电导率；S_x的倒数为 

利用映射关系得到 

式中，Δt为计算时间步长，ZT[·]表示双线性变换。 

3.由权利1所述的基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法，其特征在于：步骤5，基于Crank-Nicolson时域有限差分算法的时域展开形式： 

式中，m_1x(i)＝c₀Δt·b_x(i)·(a_x(i)-1)， 和分别为辅助变量，其迭代方程为 

。

4.由权利1所述的基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法，其特征在于：步骤8，将代入磁场和辅助变量后的电场分量的迭代方程进行整理，整理后获得等式左边为三对角矩阵形式的系数的电场显式迭代方程为 

式中， 

a_x1(i)＝κ²b_x(i)·b_x(i-1/2)； 

a_x2(i)＝1+κ²b_x(i)·(b_x(i-1/2)+b_x(i+1/2))； 

a_x3(i)＝κ²b_x(i)b_x(i+1/2)； 

a_x4(i)＝1-κ²b_x(i)·(b_x(i-1/2)+b_x(i+1/2))； 

a_x5(i)＝c₀·Δt·κ·b_x(i)·b_x(i+1/2)·(a_x(i+1/2)-1)； 

a_x6(i)＝c₀·Δt·κ·b_x(i)·b_x(i-1/2)·(a_x(i-1/2)-1)； 

a_x7＝2κb_x(i)； 

a_x8＝c₀Δtb_x(i)·(a_x(i)-1)。 

再将求解的电场值代入到权利3所述的公式(5)中求得磁场值，将求得的电场值和磁场值代入到由权利3所述的辅助变量的迭代方程中，可求得辅助变量的值。