CN104408256A - 一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法 - Google Patents

Abstract

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

Description

一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法

技术领域

本发明涉及数值仿真技术领域，特别涉及一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法。

背景技术

时域有限差分方法(FDTD)是求解麦克斯韦微分方程的直接时域方法，是电磁场数值计算方法中应用最广泛的计算方法之一。

然而，随着科学研究的深入和各种越来越广泛应用的需求，其算法本身受Courant Friedrichs Lewy(CFL)数值稳定性条件的限制的缺陷越来越凸显出来，在实际的应用中，所选取的时间步长是很短的，意味着在数值求解电磁场问题的时候，计算时间猛增，而且会使误差积累。为了克服这个缺点，一些计算数学中的方法被引入到计算电磁学中，无条件稳定的交替方向隐式(Alternating-Direction Implicit，ADI)FDTD方法、局部一维方法(Local One Dimensional method，LOD)FDTD方法和克兰克·尼克尔森(Crank-Nicolson，CN)FDTD方法相继被提出。这些方法都很好的消除了CFL数值稳定性条件的限制，从而时间步长的选择可以成倍地增加，计算时间也成倍的下降。

对于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域。利用双线性变换方法具有离散简单的优点。

由于人们对电磁波在生物体和水基物理领域的兴趣不断增加，Debye介质受到很大关注，因此Debye介质的CN-FDTD和CN-PML仿真亟需深入研究。

发明内容

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

为实现上述目的，一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法的具体方式如下：

将频域中麦克斯韦旋度方程修正为带有拉伸坐标算子的麦克斯韦方程。在截断Debye介质条件下，在SC-PML介质中的归一化麦克斯韦方程描述为：

$\mathrm{j\ωD}=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方向的偏导数。

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

截断一维Debye介质，SC-PML中的z方向极化，x方向传播的TEM波，归一化的麦克斯韦方程在频域中表示为

$\mathrm{j\ω}{D}_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方向的电导率。

将拉伸坐标变量由频域变换到z域，利用双线性变换关系

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

将变换到z域得

$S_x\left(z\right)=b_x\·\frac{1-z^{-1}}{1-a_x{z}^{-1}}---\left(8\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为计算时间步长。

将式(4)和式(5)变换到z域，采取分子消元法后，再基于Crank-Nicolson时域有限差分算法的时域展开形式，以及根据z域和时域的映射关系，将z域形式的直角坐标系中一维麦克斯韦旋度方程展开成时域有限差分的形式，可得

$\begin{array}{c}{D}_z^{n+1}\left(i\right)=a_x\left(i\right)D_z^n\left(i\right)+c_x\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(9\right)$

$\begin{array}{c}{H}_y^{n+1}(i+\frac{1}{2})=a_x(i+\frac{1}{2})H_y^n(i+\frac{1}{2})\\ +c_x(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(10\right)$

式中，c_x(i)＝c₀Δtb_x(i)/(2Δx)。

根据Debye介质的色散关系，电位移矢量和电场分量有如下关系

$D_z=\widehat{\mathrm{\ϵ}}\left(\mathrm{\ω}\right)E_z---\left(11\right)$

式中，为

$\widehat{\mathrm{\ϵ}}\left(\mathrm{\ω}\right)={\mathrm{\ϵ}}_{\∞}+\frac{{\mathrm{\ϵ}}_s-{\mathrm{\ϵ}}_{\∞}}{1+\mathrm{j\ω}{t}_0}---\left(12\right)$

式中，ε_∞是在频率为无穷大时的相对介电常数，ε_s是在频率为零时的相对介电常数，t₀为弛豫时间，推导电位移矢量和电场分量的关系并利用双线性关系将表达式变换到z域表示；

$D_z^{n+1}\left(i\right)=p_1{D}_z^n\left(i\right)+p_2{E}_z^{n+1}\left(i\right)p_3{E}_z^n\left(i\right)---\left(13\right)$

式中，p₁＝(1-τ)/(1+τ)，p₂＝(ε_∞+ε_sτ)/(1+τ)，p₃＝(ε_∞-ε_sτ)/(1+τ)，τ＝Δt/t₀。将色散关系的迭代方程(13)代入到电场的迭代方程中(9)

$\begin{array}{c}{p}_2{E}_z^{n+1}\left(i\right)=(a_x\left(i\right)-p_1)D_z^n\left(i\right)+p_3{E}_z^n\left(i\right)\\ +c_{1x}\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(14\right)$

将这组电场和磁场的隐式方程进行去耦，即将磁场分量的方程代入到电场分量的方程中，得到等式左边为三对角矩阵形式的系数的电场显式迭代方程。

$\begin{array}{c}-{\mathrm{\ξ}}_{1x}^{-}\left(i\right)E_z^{n+1}(i-1)+{\mathrm{\ξ}}_{2x}\left(i\right)E_z^{n+1}\left(i\right)-{\mathrm{\ξ}}_{1x}^{+}\left(i\right)E_z^{n+1}(i+1)={\mathrm{\φ}}_x\left(i\right)D_z^n\left(i\right)\\ +{\mathrm{\ξ}}_{1x}^{-}\left(i\right)E_z^n(i-1)+{\mathrm{\ξ}}_{3x}\left(i\right)E_z^n\left(i\right)+{\mathrm{\ξ}}_{1x}^{+}\left(i\right)E_z^n(i+1)\\ +{\mathrm{\β}}_x^{+}\left(i\right)H_y^n(i+\frac{1}{2})-{\mathrm{\β}}_x^{-}\left(i\right)H_y^n(i-\frac{1}{2})\end{array}---\left(15\right)$

式中，

${\mathrm{\ξ}}_{1x}^{\±}\left(i\right)=c_x\left(i\right)\·c_x(i\±\frac{1}{2});{\mathrm{\ξ}}_{2x}\left(i\right)=p_2+{\mathrm{\ξ}}_{1x}^{+}\left(i\right)+{\mathrm{\ξ}}_{1x}^{-}\left(i\right);{\mathrm{\φ}}_x\left(i\right)=(a_x\left(i\right)-p_1);$

${\mathrm{\ξ}}_{3x}\left(i\right)=p_3-{\mathrm{\ξ}}_{1x}^{+}\left(i\right)-{\mathrm{\ξ}}_{1x}^{-}\left(i\right);{\mathrm{\β}}_x^{\±}\left(i\right)=c_x\left(i\right)\·(1+a_x(i\±\frac{1}{2}))$

利用追赶法求解系数为三对角矩阵的电场迭代方程，得到电场分量。再将求解出的电场值代入到磁场的迭代方程中，求解出磁场分量的值。返回到式(15)进行时间步的迭代循环。

附图说明：

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

图1本发明流程框图

图2空间模型图

图3相对反射误差图

具体实施方式：

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

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

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

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

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

式中，代表在观察点处的时域离散电场分量E_z(t)的仿真解，代表参考解，代表参考解在仿真时间内的最大值，FT{·}代表傅里叶变换。参考解是通过将FDTD计算域向外延伸到10000Δx，即计算域足够大，在仿真时间内在观察点处不存在外边界反射波，得到的反射误差如图3所示，从仿真结果来看，在频域中反射系数在CNLN变化的情况下保持在有效的范围内，该算法具有无条件稳定性。

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

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

Claims (5)

1.一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法，其特征在于具体设置步骤：

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

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

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

步骤4：根据双线性变换方法的z域和频域的变换关系，将拉伸坐标变量的z域表示式代入到原方程中，采用消元法进行推导变换；

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

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

步骤7：根据Debye介质的色散关系，推导电位移矢量和电场分量的关系并利用双线性关系将表达式变换到z域表示，将关系表达式代入电场的隐式方程中，整理为关于电场的隐式方程；

步骤8：将这组隐式方程进行去耦，即将磁场分量的方程代入到电场分量的方程中；

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

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

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

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

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

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

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

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

利用映射关系 $\mathrm{j\ω}\⇒\frac{2}{\mathrm{\Δt}}\frac{1-z^{-1}}{1+z^{-1}},$ 得到

$S_x\left(z\right)=\mathrm{ZT}\left[S_x^{-1}\right]=-b_x\·\left[\frac{1-z^{-1}}{1-a_x{z}^{-1}}\right]---\left(3\right)$

式中， $b_x=\frac{1}{(1+\frac{\mathrm{\Δt}{\mathrm{\σ}}_x}{2{\mathrm{\ϵ}}_0})},a_x=\frac{(1-\frac{\mathrm{\Δt}{\mathrm{\σ}}_x}{2{\mathrm{\ϵ}}_0})}{(1+\frac{\mathrm{\Δt}{\mathrm{\σ}}_x}{2{\mathrm{\ϵ}}_0})},$ Δt为计算时间步长，ZT[·]表示双线性变换。

3.由权利1所述的截断一维Debye介质Crank-Nicolson完全匹配层实现算法，其特征在于：步骤5，基于Crank-Nicolson时域有限差分算法的时域展开形式：

$\begin{array}{c}{S}_z^{n+1}\left(i\right)=a_x\left(i\right)D_z^n\left(i\right)+c_x\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(4\right)$

$\begin{array}{c}{H}_y^{n+1}(i+\frac{1}{2})=a_x(i+\frac{1}{2})H_y^n(i+\frac{1}{2})\\ +c_x(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(5\right)$ 式中，c_x(i)＝c₀Δtb_x(i)/(2Δx)。

4.由权利1所述的一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法，其特征在于：步骤7，电位移矢量和电场分量有如下关系

$D_z=\widehat{\mathrm{\ϵ}}\left(\mathrm{\ω}\right)E_z---\left(6\right)$

式中，为

$\widehat{\mathrm{\ϵ}}\left(\mathrm{\ω}\right)={\mathrm{\ϵ}}_{\∞}+\frac{{\mathrm{\ϵ}}_s-{\mathrm{\ϵ}}_{\∞}}{1+\mathrm{j\ω}{t}_0}---\left(7\right)$

式中，ε_∞是在频率为无穷大时的相对介电常数，ε_s是在频率为零时的相对介电常数，t₀为弛豫时间，推导电位移矢量和电场分量的关系并利用双线性关系将表达式变换到z域表示

$D_z^{n+1}\left(i\right)=p_1{D}_z^n\left(i\right)+p_2{E}_z^{n+1}\left(i\right)-p_3{E}_z^n\left(i\right)---\left(8\right)$

式中，p₁＝(1-τ)/(1+τ)，p₂＝(ε_∞+ε_sτ)/(1+τ)，p₃＝(ε_∞-ε_sτ)/(1+τ)，τ＝Δt/t₀。

5.由权利1所述的一种截断一维Debye介质Crank-Nicolson完全匹配层实现算法，其特征在于：步骤9，将代入磁场后的场量的方程进行整理，整理后获得等式左边为三对角矩阵形式的系数的电场显式迭代方程为

$\begin{array}{c}-{\mathrm{\ξ}}_{1x}^{-}\left(i\right)E_z^{n+1}(i-1)+{\mathrm{\ξ}}_{2x}\left(i\right)E_z^{n+1}\left(i\right)-{\mathrm{\ξ}}_{1x}^{+}\left(i\right)E_z^{n+1}(i+1)={\mathrm{\φ}}_x\left(i\right)D_z^n\left(i\right)\\ +{\mathrm{\ξ}}_{1x}^{-}\left(i\right)E_z^n(i-1)+{\mathrm{\ξ}}_{3x}\left(i\right)E_z^n\left(i\right)+{\mathrm{\ξ}}_{1z}^{+}\left(i\right)E_z^n(i+1)\\ +{\mathrm{\β}}_x^{+}\left(i\right)H_y^n(i+\frac{1}{2})-{\mathrm{\β}}_x^{-}\left(i\right)H_y^n(i-\frac{1}{2})\end{array}---\left(9\right)$

式中

${\mathrm{\ξ}}_{1x}^{\±}\left(i\right)=c_x\left(i\right)\·c_x(i\±\frac{1}{2});{\mathrm{\ξ}}_{2x}\left(i\right)=p_2+{\mathrm{\ξ}}_{1x}^{+}\left(i\right)+{\mathrm{\ξ}}_{1x}^{-}\left(i\right);{\mathrm{\φ}}_x\left(i\right)=(a_x\left(i\right)-p_1);$

${\mathrm{\ξ}}_{3x}\left(i\right)=p_3-{\mathrm{\ξ}}_{1x}^{+}\left(i\right)-{\mathrm{\ξ}}_{1x}^{-}\left(i\right);{\mathrm{\β}}_x^{\±}\left(i\right)=c_x\left(i\right)\·(1+a_x(i\±\frac{1}{2})).$

再将求解出的电场值代入到权利3中的公式(5)中求解出磁场的值。