本发明涉及数值仿真
技术领域：
，特别涉及一种基于Douglas-Gunn(DG)算法的二维色散介质Crank-Nicolson完全匹配层实现算法。
背景技术：
：时域有限差分方法(FDTD)作为一种计算电磁方法被广泛地应用于各种时域的电磁仿真计算中，如天线、射频电路、光学器件和半导体等。FDTD具有广泛的适用性、适合并行计算、计算程序通用性等特点。然而，随着科学研究的深入和各种越来越广泛应用的需求，其算法本身受CourantFriedrichsLewy(CFL)数值稳定性条件的限制的缺陷越来越明显。算法本身所受数字稳定性条件限制：在计算过程中时间步长和空间步长必须满足CFL约束条件，即Δt≤(c(Δx)-2+(Δy)-2+(Δz)-2)-1]]>式中，Δt为计算时间步长，c为FDTD计算域中介质的光速，Δx、Δy和Δz为三维空间步长。在实际计算中，空间离散步长和时间步长相对波长和周期都非常小，所以必然会在计算电大尺寸目标时出现资源不足的情况，导致FDTD的计算效率很低。因此为了消除CFL条件的限制，无条件稳定的交替方向隐式(AlternatingDirectionImplicit，ADI)FDTD方法、局部一维(LocalOneDimension，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提出的完全匹配层(PerfectlyMatchedLayer，PML)是目前应用较广的吸收边界条件，PML可以理解为：通过在FDTD区域截断边界处设置一种特殊介质层，该层介质的波阻抗与相邻介质波阻抗完全匹配，从而使入射波无反射地穿过分界面而进入PML层，PML层是有耗介质，最后将电磁波吸收。目前常用的PML吸收边界主要有拉伸坐标变换完全匹配层(SC-PML)和单轴各项异性完全匹配层(UPML)。由于人们对电磁波与色散介质相互作用问题的兴趣不断增加，因此涉及色散介质的CN-FDTD和CN-PML仿真亟需深入研究。技术实现要素：本发明的目的是针对FDTD算法受到CFL稳定性条件限制的缺陷，提高色散介质的PML算法的计算效率和吸收效率而提出的基于DG算法与辅助微分方程方法及CN-FDTD的SC-PML算法。该算法可以极大地改善CN-FDTD-PML算法的求解效率。基于DG算法的二维色散Crank-Nicolson完全匹配层实现算法，包括下列步骤：步骤1：将频域中麦克斯韦方程修正为带有拉伸坐标算子的麦克斯韦方程，并在直角坐标系中表示；针对色散介质项，利用辅助微分方程方法设置辅助变量；步骤2：根据频域和时域的映射变换关系，将直角坐标系中的二维修正的麦克斯韦方程变换到时域表示，同时基于辅助微分方程方法设置辅助变量；步骤3：基于Crank-Nicolson时域有限差分算法的时域展开形式，将时域形式的直角坐标系中二维麦克斯韦方程展开成时域有限差分的形式，同时也将辅助微分方程变换为时域有限差分的形式；步骤4：将时域有限差分形式整理成求解的形式，结果产生一组电场和磁场的耦合隐式方程，经整理后得到系数矩阵为块三对角矩阵形式的隐式迭代方程；步骤5：使用Crank-NicolsonDouglas-Gunn(CNDG)方法，将步骤4所得到的电场迭代方程近似分解为可以高效求解的、系数为三对角矩阵的两个迭代方程；步骤6：利用步骤5所得到的迭代方程求解出电场值，然后将求解出的电场值代入到磁场的迭代方程中，求解出磁场分量，将求解出的电场值和磁场值代入到辅助变量的迭代方程中，求解出辅助变量的值；重复步骤6，从而在时间上迭代求解。采用CNDG方法可以有效地将系数矩阵为块三对角矩阵形式的电场迭代方程近似分解为可以高效求解的、系数为三对角矩阵的两个迭代方程，降低计算复杂度，提高计算效率，对FDTD算法有指导意义。附图说明：图1是本发明流程框图；图2是本发明在不同CFLN条件下截断Debye色散介质的时域相对反射误差特性曲线；图3是本发明在不同CFLN条件下截断Debye色散介质的频域相对反射误差特性曲线；图4是本发明在不同CFLN条件下截断Lorentz色散介质的时域相对反射误差特性曲线；图5是本发明在不同CFLN条件下截断Lorentz色散介质的频域相对反射误差特性曲线。具体实施方式：本发明的主旨是提出一种基于DG算法的二维色散介质Crank-Nicolson完全匹配层实现算法，利用Douglas-Gunn求解思想极大地提高了电磁场计算速度。下面结合附图对本发明实施方式作进一步地详细描述。图1为本发明流程图，具体实现步骤如下：步骤1：将频域中麦克斯韦方程修正为带有拉伸坐标算子的麦克斯韦方程，并将频域中修正后的麦克斯韦方程在直角坐标系中表示；针对色散介质项，利用辅助微分方程方法设置辅助变量；TM波在线性色散介质中传播可以描述为-jωHx=c0Sy-1∂yEz---(1)]]>jωHy=c0Sx-1∂xEz---(2)]]>jωϵr(ω)Ez=c0Sx-1∂xHy-c0Sy-1∂yHx---(3)]]>式中，c0是自由空间电磁波传播速度，Sη(η＝x，y)为PML复数拉伸坐标变量。在PML情况下，Sη可以表示为Sη=κη+σηjωϵ0---(4)]]>在CFS-PML情况下，Sη可以表示为Sη=κη+σηαη+jωϵ0---(5)]]>色散介质中，εr(ω)可以表示为ϵr(ω)=ϵ∞+σjωϵ0+χ(ω)---(6)]]>式中，ε∞为无限频率介电常数，σ为电导率，χ(ω)为介质的电极化率。式(3)可以表示为jωϵ∞Ez+σϵ0Ez+jωPz=c0Sx-1∂xHy-c0Sy-1∂yHx---(7)]]>式中，Pz可由下式得到Pz＝χ(ω)Ez(8)步骤2：根据频域和时域的映射变换关系，将直角坐标系中的二维修正的麦克斯韦方程变换到时域表示，同时基于辅助微分方程方法设置辅助变量，即∂tHx=-c0κy-1∂yEz+gy---(9)]]>∂tHy=c0κx-1∂xEz-gx---(10)]]>ϵ∞∂tEz+σϵ0·Ez+∂tPz=c0κx-1∂xHy-c0κy-1∂yHx-fx+fy---(11)]]>式中，fx、fy、gx和gy为辅助变量。步骤3：基于Crank-Nicolson时域有限差分算法的时域展开形式，将时域形式的直角坐标系中二维麦克斯韦方程展开成时域有限差分形式，同时也将时域辅助微分方程变换为时域有限差分的形式，利用CN项将式(9)-(11)离散化，可得离散方程为Hxi,j+1/2n+1=Hxi,j+1/2n-Ceyj+1/2Γ1y(Ezi,jn+1)+ΔthΓ2(gyi,j+1/2n+1)---(12)]]>Hyi+1/2,jn+1=Hyi+1/2,jn+Cexi+1/2Γ1x(Ezi,jn+1)-ΔthΓ2(gxi+1/2,jn+1)---(13)]]>Ezi,jn+1=CeEzi,jn+ChxiΓ1x(Hyi+1/2,jn+1)-ChyjΓ1y(Hxi,j+1/2n+1)-(Δth/ϵ∞)Γ2(fxi,jn+1)+(Δth/ϵ∞)Γ2(fyi,jn+1)-ϵ∞-1Γ3(Pzi,jn+1)---(14)]]>式中，Δth=Δt2,u=2ϵ0ϵ∞,v=σΔt,Ce=u-vu+v,Ceηk+1/2=c0Δthκηk+1/2Δη,]]>Chηk=Δthϵ∞κηkΔη,r1ηk=αηkκηk+σηkκηkϵ0Δth,r2ηk=c0Δthσηkϵ0κηk2Δx,L1ηk=1-r1ηk1+r1ηk,]]>Δη(η＝x，y)是空间单元尺寸，k(k＝i，j)是计算单元间的插入数值，为了清楚，Гη[*]为CN方法中的简写形式，如Γ1x(Ezi,jn+1)=Ezi+1,jn+1-Ezi,jn+1+Ezi+1,jn-Ezi,jn---(15)]]>对于简单的线性Debye和Lorentz色散介质，χ(ω)可以写为χ(ω)=d0e0+e1jω+e2(jω)2---(16)]]>式中，d0、e0、e1与e2为有理多项式的系数，式(8)可以离散为n+1时间步长公式Pzi,jn+1=Q0Pzi,jn-Q1Pzi,jn-1+Qe·[Ezi,jn+1+Ezi,jn+m(Ezi,jn+Ezi,jn-1)]---(17)]]>式中的m可以作为色散介质的开关变量；对应于Debye和Lorentz色散介质，m分别取0和1；步骤4：将时域有限差分的形式整理成求解的形式，结果产生一组电场和磁场的耦合方程，这是一组隐式方程，将这组方程去耦，整理后获得左边为块三对角矩阵形式的系数电场隐式迭代方程(1-D2x-D2y)Ezi,jn+1=(ϵ∞Ce-(1+m)Qeϵ∞+Qe+D2x+D2y)Ezi,jn+ρ---(18)]]>式中，D2x定义为D2x(Ezi,jn+1)=Cxi+Ezi+1,jn+1-(Cxi++Cxi-)Ezi,jn+1+Cxi-Ezi-1,jn+1---(19)]]>D2y具有类似的定义；ρ是n时刻已知的场量和辅助变量的简写形式，定义为ρ=ξxi·(Hyi+1/2,jn-Hyi-1/2,jn)-ζyj·(Hxi,j+1/2n-Hxt,j-1/2n)-Sxi+gxi+1/2,jn+Sxi-gxt-1/2,jn-Syj+gyi,j+1/2n+Syj-gyi,j-1/2n-θxifxi,jn-Q0-1ϵ∞+QePzi,jn+Q1ϵ∞+QePzi,jn-1-mQeϵ∞+QeEzi,jn-1+θyjfyt,jn---(20)]]>步骤5：式(18)仍然是一个比较复杂的矩阵，仍然需要很大的计算量，使用CNDG方法，在式(18)的左边和右边分别添加与整理得(1-D2x)Ezi,j*=(ϵ∞Ce-(1+m)Qeϵ∞+Qe+D2x+2D2y)Ezi,jn+ρ---(21)]]>(1-D2y)Ezi,jn+1=Ezi,j*-D2yEzi,jn---(22)]]>步骤6：利用步骤5所得到的迭代方程求解出电场值，然后将求解出的电场值代入到磁场的迭代方程中，求解出磁场分量，将求解出的电场值和磁场值代入到辅助变量的迭代方程中，求解出辅助变量的值。重复步骤6，从而在时间上迭代求解。图2是本发明在不同CFLN条件下截断Debye色散介质的时域相对反射误差特性曲线，图4是本发明在不同CFLN条件下截断Lorentz色散介质的时域相对反射误差特性曲线，为了验证所提方法，对本发明算法进行编程，通过计算机仿真得到附2和图4所示结果，其中，式中是传统FDTD算法能够保证数值稳定性的最大时间离散间隔。由图中可以看出，CN-PML的吸收性能并不随CFLN的增大而变化；图3是本发明在不同CFLN条件下截断Debye色散介质的频域相对反射误差特性曲线，图5是本发明在不同CFLN条件下截断Lorentz色散介质的频域相对反射误差特性曲线，可以看出PML吸收性能均不随CFLN的增大而变化，说明该算法具有无条件稳定性，仿真过程所需时间较传统算法仿真时间较短。以上所述仅为本发明的较佳实施例，并不限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。当前第1页1 2 3