本发明涉及数值仿真
技术领域：
，特别涉及一种截断二维Debye介质与Lorentz介质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稳定性条件限制的缺陷，提高CN-PML算法的计算效率和吸收效率而提出的一种截断二维Debye介质与Lorentz介质Crank-Nicolson完全匹配层实现算法。该算法结合Crank-Nicolson-Approximate-Decoupling(CNAD)求解思想，可以极大地改善CN-FDTD-PML算法的求解效率。一种截断二维Debye介质与Lorentz介质Crank-Nicolson完全匹配层实现算法，包括下列步骤：步骤1：将频域中二维TE模的麦克斯韦方程修正为带有拉伸坐标算子的麦克斯韦方程，并在直角坐标系中表示；针对色散介质项，利用辅助微分方程方法设置辅助变量；步骤2：根据频域和时域的映射变换关系，将直角坐标系中的二维修正的麦克斯韦方程变换到时域表示，同时基于辅助微分方程方法设置辅助变量；步骤3：基于Crank-Nicolson时域有限差分算法的时域展开形式，将时域形式的直角坐标系中二维麦克斯韦方程展开成时域有限差分的形式，同时也将辅助微分方程变换为时域有限差分的形式；步骤4：将时域有限差分形式整理成求解的形式，得到两个相互耦合的电场(Ex，Ey)迭代方程；步骤5：使用CNAD方法，将步骤4所得到的Ex迭代方程近似为可以高效求解的、系数为三对角矩阵的迭代方程；步骤6：利用步骤5所得到的迭代方程可以依次求解出两个电场分量值，然后将求解出的两个电场分量值代入到磁场的迭代方程中，求解出磁场分量，将求解出的电场值和磁场值代入到辅助变量的迭代方程中，求解出辅助变量的值；重复步骤6，从而在时间上迭代求解。采用CNAD算法和辅助微分方程方法可以有效地降低计算复杂度，优化计算过程，从而提高电磁场计算速度。附图说明：图1是本发明流程框图；图2是Debye色散介质下相对反射误差图；图3是Lorentz色散介质下相对反射误差图。具体实施方式：本发明的主旨是提出一种截断二维Debye介质与Lorentz介质Crank-Nicolson完全匹配层实现算法，利用CNAD算法和辅助微分方程方法减小计算量，从而提高电磁场计算速度。下面结合附图对本发明实施方式作进一步地详细描述。图1为本发明流程图，具体实现步骤如下：步骤1：将频域中二维TE模的麦克斯韦方程修正为带有拉伸坐标算子的麦克斯韦方程，并将频域中修正后的麦克斯韦方程在直角坐标系中表示；针对色散介质项，利用辅助微分方程方法设置辅助变量；二维TE波在线性色散介质中传播可以描述为jωϵr(ω)Ex=c0Sy-1∂yHz---(1)]]>jωϵr(ω)Ey=-c0Sx-1∂xHz---(2)]]>jωHz=-c0Sx-1∂xEy+c0Sy-1∂yEx---(3)]]>式中，c0是自由空间电磁波传播速度，Sη(η＝x，y)为PML复数拉伸坐标变量，在PML情况下，Sη可以表示为Sη=κη+σηjωϵ0---(4)]]>在CFS-PML情况下，Sη可以表示为Sη=κη+σηαη+jωϵ0---(5)]]>在色散介质中，εr(ω)可以表示为ϵr(ω)=ϵ∞+σjωϵ0+χ(ω)---(6)]]>式中，ε∞为无限频率介电常数，σ为电导率，χ(ω)为介质的电极化率。那么，式(1)与式(2)可以表示为jωϵ∞Ex+σϵ0Ex+jωPx=c0Sy-1∂yHz---(7)]]>jωϵ∞Ey+σϵ0Ey+jωPy=-c0Sx-1∂xHz---(8)]]>式中，Px和Py可由式(9)和式(10)得到Px＝χ(ω)Ex(9)Py＝χ(ω)Ey(10)步骤2：根据频域和时域的映射变换关系，将直角坐标系中的二维修正的麦克斯韦方程变换到时域表示，同时基于辅助微分方程方法设置辅助变量，即ϵ∞∂tEx+σϵ0Ex+δtPx=c0κy-1∂yHz-fxy---(11)]]>ϵ∞∂tEy+σϵ0Ey+∂tPy=-c0κx-1∂xHz+fyx---(12)]]>∂tHz=-c0κx-1∂xEy+c0κy-1∂yEx+gzx-gzy---(13)]]>式中，fxy、fyx、gzx和gzy为辅助变量。步骤3：基于Crank-Nicolson时域有限差分算法的时域展开形式，将时域形式的直角坐标系中二维麦克斯韦方程展开成时域有限差分形式，同时也将时域辅助微分方程变换为时域有限差分的形式，利用CN项将式(11)-(13)离散化，可得离散方程为Exi+1/2,jn+1=CeExi+1/2,jn+ChyjΓ1y(Hzi+1/2,j+1/2n+1)-CfΓ2(fxyi+1/2,jn+1)-CpΓ3(Pxi+1/2,jn+1)---(14)]]>Eyi,j+1,2n+1=CeEyi,j+1/2n-ChxiΓ1x(Hzi+1/2,j+1/2n+1)+CfΓ2(fyxi,j+1/2n+1)-CpΓ3(Pxi,j+1/2n+1)---(15)]]>Hzi+1/2,j+1/2n+1=Hzi+1/2,j+1/2n-Cexi+1/2Γ1x(Eyi,j+1/2n+1)-ΔthΓ2(gzyi+1/2,j+1/2n+1)+Ceyj+1/2Γ1y(Exi+1/2,jn+1)+ΔthΓ2(gzxi+1/2,j+1/2n+1)---(16)]]>式中，u＝2ε0ε∞，v＝σΔt，a＝c0Δth，Ce=u-vu+v,Chηk=ϵ0Δtu+v·κηk-1·c0Δη,Cf=ϵ0Δtu+v,Cp=2ϵ0u+v,]]>Ceηk=bηκηk-1,r1ηk=αηkκηk+σηkκηkϵ0Δth,r2ηk=σηkκηk2ϵ0,L1ηk=1-r1ηk1+r1ηk]]>Δη(η＝x，y)是空间单元尺寸，k(k＝i，j)是计算单元间的插入数值，为了清楚，Γη[*]为CN方法中的简写形式，如Γ1y(Hzi+1/2,j+1/2n+1)=Hzi+1/2,j+1/2n+1-Hzi+1/2,j-1/2n+1+Hzi+1/2,j+1/2n-Hzi+1/2,j-1/2n---(17)]]>对于简单的线性Debye和Lorentz色散介质，χ(ω)可以写为χ(ω)=d0e0+e1jω+e2(jω)2---(18)]]>式中，d0、e0、e1与e2为有理多项式的系数，式(9)和式(10)可以写为Pηn+1=Q0Pηn-Q1Pηn-1+Qe·[Eηn+1+Eηn+m(Eηn+Eηn-1)]---(19)]]>式中，η＝x，y，而m可以作为色散介质的开关变量，对应于Debye和Lorentz色散介质，m分别取0和1；步骤4：将时域有限差分形式整理成求解电场分量的形式，得到两个相互耦合的电场(Ex，Ey)迭代方程-Cy1jExi+1/2,j-1n+1+Cy2jExi+1/2,jn+1-Cy3jExi+1/2,j+1n+1=Cy1jExi+1/2,j-1n+Cy4jExi+1/2,jn+Cy3jExi+1/2,j+1n-Cy5i,jDxDy(Eyi,j+1/2n+1+Eyi,j+1/2n)+ρ1i,jn---(20)]]>-Cx1iEyi-1,j+1/2n+1+Cx2iEyi,j+1/2n+1-Cx3iEyi+1,j+1/2n+1=Cx1iEyi-1,j+1/2n+Cx4iEyi,j+1/2n+Cx3iEyi+1,j+1/2n-Cx5i,jDxDy(Exi+1/2,jn+1+Exi+1/2,jn)+ρ2i,jn---(21)]]>式中，和是n时刻已知的部分场量和辅助变量的简写形式，定义为ρ1i,jn=Cy6jDyHzi+1/2,j+1/2n-Cy7jfxyi+1/2,jn+Cy8i,jDygzxi+1/2,j+1/2n-Cy9jgzyi+1/2,j+1/2n+Cy10jgzyi+1/2,j-1/2n-Cp(Q0-1)Pxi+1/2,jn-CpQ1Pxi+1/2,jn-1---(22)]]>ρ2i,jn=-Cx6iDxHzi+1/2,j+1/2n+Cx7ifyxi,j+1/2n+Cx8i,jDxgzyi+1/2,j+1/2n-Cx9igzxi+1/2,j+1/2n+Cx10igzxi-1/2,j+1/2n-Cp(Q0-1)Pyi,j+1/2n-CpQ1Pyi,j+1/2n-1---(23)]]>步骤5：式(20)与式(21)仍然是一个比较复杂的矩阵，仍然需要很大的计算量，使用CNAD方法，将式(20)中的近似为通过整理得-Cy1jExi+1/2,j-1n+1+Cy2jExi+1/2,jn+1-Cy3jExi+1/2,j+1n+1=Cy1jExi+1/2,j-1n+Cy4jExi+1/2,jn+Cy3jExi+1/2,j+1n-2Cy5i,jDxDyEyi,j+1/2n+ρ1i,jn---(24)]]>步骤6：利用步骤5所得到的迭代方程可以依次求解出两个电场分量值，将求解出的两个电场分量值代入到磁场的迭代方程中，求解出磁场分量，将求解出的电场值和磁场值代入到辅助变量的迭代方程中，求解出辅助变量的值。循环步骤6，从而在时间上迭代求解。图2是本发明为Debye色散介质下相对反射误差图，图3为本发明Lorentz色散介质下相对反射误差图，为了验证所提方法，对本发明算法进行编程，通过计算机仿真得到图2和图3所示结果，其中，式中是传统FDTD算法能够保证数值稳定性的最大时间离散间隔。由图中可以看出，CN-PML的吸收性能在CFLN变化区域内可以接受，说明该算法具有无条件稳定性，仿真过程所需时间较传统算法仿真时间较短。以上所述仅为本发明的较佳实施例，并不限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。当前第1页1 2 3