CN102722651B - 二维柱坐标完全匹配吸收边界的实现方法 - Google Patents

Abstract

本发明公开了一种二维柱坐标完全匹配吸收边界的实现方法，具体步骤如下：输入模型文件；对电场、磁场及相关辅助变量初始化为零；计算PML中的一维系数数组；更新计算整个计算区域的电位移矢量；更新计算整个计算区域的电场；更新计算电场场源；更新计算PML区域中电场辅助变量；更新计算整个计算区域的磁感应强度；更新计算整个计算区域的磁场；更新计算PML区域中磁场辅助变量；根据设定的仿真计算时长，判断是否继续更新；实际计算时间小于该设定值时，判断继续计算，返回至步骤3；否则输出计算得到的电场和磁场，结束。本发明能够很方便的与CFS参数相结合，能够更加有效的吸收向外传播的电磁波。

Description

二维柱坐标完全匹配吸收边界的实现方法

技术领域

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

背景技术

时域有限差分方法（FDTD,Finite-difference time-domain）作为一种电磁场时域数值计算方法已在电波传播、天线设计、电磁散射、电磁兼容等领域得到了广泛的应用。吸收边界的引入，使得有限区域的FDTD计算可以求解无限大媒质中电波的传播与散射问题。吸收边界的性能是影响FDTD方法性能和精度的主要因素，其中理想匹配层（PML,Perfectly matched layer）吸收边界是目前最常用的吸收边界。对于具有旋转对称结构的电磁问题，采用柱坐标系可以将3维问题简化成2维或者2.5维问题，可大幅度降低程序实现难度，同时提高计算效率。

经典的直角坐标系下的PML吸收边界有分裂场的PML(Split PML)和单轴各向异性PML（UPML:Uniaxial Anisotropic PML），但它们对低频以及凋落模的吸收效果并不理想；带有复频率偏移（CFS,Complex frequency shift）参数的PML（CFS-PML）吸收边界能够很好的改善PML对低频，凋落波与掠射情况的吸收效果；2000年J.Alan Rodeny与Stephen D.Gendney利用卷积关系，给出了直角坐标系实现CFS的卷积PML(convolution PML)吸收边界；2005年Xiaoting Dong Wen-Yan Yin与Yeow-Beng Gan利用双线性变换，给出了CFS在扩展直角坐标系与UPML中的应用。

现有的柱坐标系下PML吸收边界主要有：分裂场的PML和单轴各向异性UPML，近年柱坐标FDTD应用中的PML边界多采用复扩展坐标系下的UPML。但上述柱坐标PML吸收边界对低频以及凋落模的吸收效果很不理想。目前尚未见把CFS引入柱坐标PML吸收边界以改善PML对低频，凋落波与掠射波的吸收效果的研究报道和专利。

发明内容

本发明的目的是提供一种二维柱坐标完全匹配吸收边界的实现方法，能够很方便的与复频率偏移参数相结合，能够更加有效的吸收向外传播的电磁波。

本发明所采用的技术方案是，一种二维柱坐标完全匹配吸收边界的实现方法，其特征在于，具体步骤如下：

步骤1、输入模型文件：

所述模型文件内容包括：计算区域大小Nr×Nz，其中，Nr为柱坐标ρ方向网格数，Nz为柱坐标z方向网格数；空间步长Δζ，ζ=ρ,z；时间步长Δt；真空状态下，整个计算区域网格的电导率σ=0，磁导率μ₀，介电常数ε₀；PML吸收边界层数与相关参数κ_ζmax和α_ζζ，κ_ζmax,ζ=ρ,z为大于1的任一实数，α_ζζ，ζ=ρ,z，ζ=i,k,i+1/2,k+1/2为大于等于0小于1的任一实数，i，k表示空间网格坐标，i＝0,1,2…Nr，k=0,1,2…Nz；仿真计算时长；场源参数；

步骤2、根据步骤1输入的模型文件，对电场、磁场及相关辅助变量初始化为零；

其中，所述电场、磁场及相关辅助变量包括：柱坐标ρ方向二维电场数组E_ρ[Nr][Nz+1]，电位移矢量数组D_ρ[Nr][Nz+1]，柱坐标z方向二维电场数组E_z[Nr+1][Nz]，电位移矢量数组D_z[Nr+1][Nz]，柱坐标方向二维磁场数组柱坐标方向二维磁感应强度数组PML中电场辅助变量二维数组ψ_φρ[Nr+1][Nz]，ψ_φγ[Nr+1][Nz]，ψ_φz[Nr][Nz+1]，PML中磁场辅助变量二维数组Φ_ρz[Nr][Nz]和Φ_zρ[Nr][Nz]；

步骤3、计算PML中的一维系数数组：C_hφρi0,C_hφρi1,C_hφρi2，C_hφγi0,C_hφγi1,C_hφγi2，C_hφzk0,C_hφzk1,C_hφzk2，C_eρzk0,C_eρzk1,C_eρzk2，C_ezρi0,C_ezρi1,C_ezρi2；

具体过程为： $C_{\mathrm{h\φ\ρi}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}}},$

$C_{\mathrm{h\φ\ρi}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\ρi}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}-2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\γi}0}=-\frac{\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})-2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}}}{\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}}},$

$C_{\mathrm{h\φ\γi}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\γi}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}-2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φzk}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}}},$

$C_{\mathrm{h\φzk}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}}+2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}})},$

$C_{\mathrm{h\φzk}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}}-2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}})},$

$C_{\mathrm{e\ρzk}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_z+{\mathrm{\σ}}_{\mathrm{zk}+1/2})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_z+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{e\ρzk}1}=\frac{{\mathrm{\Δt\α}}_{\mathrm{zk}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_{\mathrm{zk}+1/2}+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{e\ρzk}2}=\frac{{\mathrm{\Δt\α}}_{\mathrm{zk}+1/2}-2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_{\mathrm{zk}+1/2}+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{ez\ρi}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

$C_{\mathrm{ez\ρi}1}=\frac{{\mathrm{\Δt\α}}_{\mathrm{\ρi}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

$C_{\mathrm{ez\ρi}2}=\frac{{\mathrm{\Δt\α}}_{\mathrm{\ρi}+1/2}-2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}};$

其中，

$A_{\mathrm{\ρi}}=1+({\mathrm{\κ}}_{\mathrm{\ρ}\max }-1)\frac{{({\mathrm{\ρ}}_i-{\mathrm{\ρ}}_0)}^{m+1}}{(m+1)d^m}\frac{1}{{\mathrm{\ρ}}_i},$

$B_{\mathrm{\ρi}}=\frac{{\mathrm{\σ}}_{\max }{({\mathrm{\ρ}}_i-{\mathrm{\ρ}}_0)}^{m+1}}{(m+1)d^m}\frac{1}{{\mathrm{\ρ}}_i};$

ρ₀表示柱坐标ρ方向PML层靠近FDTD区的界面的位置坐标，m为多项式的阶数，d为PML层厚度，ρ_i表示柱坐标ρ方向第i个网格的坐标；

κ_ζζ，ζ=ρ,z，ζ=i,k,i+1/2,k+1/2表示ζ网格点处的κ_ζ的取值，κ_ζ的计算公式为：

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

σ_ζζ，ζ=ρ,z，ζ=i,k,i+1/2,k+1/2表示ζ网格点处的σ_ζ的取值，σ_ζ的计算公式为：

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

σ_ζmax＝Cσ_opt＝C(m+1)/150πΔζ，

C为大于0任一实数；

步骤4、更新计算整个计算区域的电位移矢量；

其中，电位移矢量在z方向的分量的更新公式为：

$D_z{|}_{i,k+1/2}^{n+1}=D_z{|}_{i,k+1/2}^n+\mathrm{\Δt}({\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n-1/2}+C_{\mathrm{h\φ\ρi}1}\frac{{\∂H}_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}+{\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n-1/2}+C_{\mathrm{h\φ\γi}1}\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}),$

电位移矢量在ρ方向的分量的更新公式为：

$D_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}=D_{\mathrm{\ρ}}{|}_{i+1/2,k}^n-\mathrm{\Δt}({\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n-1/2}+C_{\mathrm{h\φzk}1}\frac{{\∂H}_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,k}^{n+1/2}),$

n表示离散时间步，

$\frac{{\∂H}_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{\mathrm{\Δ\ρ}},$

$\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}+H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{2i\·\mathrm{\Δ\ρ}},$

$\frac{{\∂H}_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,k}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}}{\mathrm{\Δz}};$

步骤5、更新计算整个计算区域的电场；

步骤6、更新计算电场场源；

步骤7、更新计算PML区域中电场辅助变量；

其中，电场辅助变量的更新计算公式为：

${\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n+1/2}=C_{\mathrm{h\φ\ρi}0}{\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n-1/2}+(C_{\mathrm{h\φ\ρi}0}{C}_{\mathrm{h\φ\ρi}1}+C_{\mathrm{h\φ\ρi}2})\frac{{\∂H}_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2},$

${\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n+1/2}=C_{\mathrm{h\φ\γi}0}{\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n-1/2}+(C_{\mathrm{h\φ\γi}0}{C}_{\mathrm{h\φ\γi}1}+C_{\mathrm{h\φ\γi}2})\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2},$

${\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n+1/2}=C_{\mathrm{h\φzk}0}{\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n-1/2}+(C_{\mathrm{h\φzk}0}{C}_{\mathrm{h\φzk}1}+C_{\mathrm{h\φzk}2})\frac{{\∂H}_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,,k}^{n+1/2},$

步骤8、更新计算整个计算区域的磁感应强度；

其中，磁感应强度的更新计算公式为：

其中， $\frac{{\∂E}_{\mathrm{\ρ}}}{\∂z}{|}_{i+1/2,k+1/2}^n=\frac{E_{\mathrm{\ρ}}{|}_{i+1/2,k+1}^n-E_{\mathrm{\ρ}}{|}_{i+1/2,k}^n}{\mathrm{\Δz}},\frac{{\∂E}_z}{\∂\mathrm{\ρ}}{|}_{i+1/2,k+1/2}^n=\frac{E_z{|}_{i+1,k+1/2}^n-E_z{|}_{i-1,k+1/2}^n}{\mathrm{\Δ\ρ}};$

步骤9、更新计算整个计算区域的磁场；

步骤10、更新计算PML区域中磁场辅助变量；

其中，磁场辅助变量的更新计算公式为：

${\mathrm{\Φ}}_{\mathrm{\ρz}}{|}_{i+1/2,k+1/2}^n=C_{\mathrm{e\ρzk}0}{\mathrm{\Φ}}_{\mathrm{\ρz}}{|}_{i+1/2,k+1/2}^{n-1}+(C_{\mathrm{e\ρzk}0}{C}_{\mathrm{e\ρzk}1}+C_{\mathrm{e\ρzk}2})\frac{{\∂E}_{\mathrm{\ρ}}}{\∂z}{|}_{i+1/2,k+1/2}^n,$

${\mathrm{\Φ}}_{\mathrm{z\ρ}}{|}_{i+1/2,k+1/2}^n=C_{\mathrm{ez\ρi}0}{\mathrm{\Φ}}_{\mathrm{z\ρ}}{|}_{i+1/2,k+1/2}^{n-1}+(C_{\mathrm{ez\ρi}0}{C}_{\mathrm{ez\ρi}1}+C_{\mathrm{ez\ρi}2})\frac{{\∂E}_z}{{\∂}_{\mathrm{\ρ}}}{|}_{i+1/2,k+1/2}^n,$

步骤11、根据步骤1设定的仿真计算时长，判断是否继续更新；实际计算时间小于该设定值时，继续计算，返回至步骤3；否则输出步骤5计算得到的电场和步骤9计算得到的磁场，结束。

步骤5中，真空中，电场在z方向分量的更新公式为：

$E_z{|}_{i,k+1/2}^{n+1}=\frac{D_z{|}_{i,k+1/2}^{n+1}}{{\mathrm{\ϵ}}_0},$

电场在ρ方向分量的更新公式为：

$E_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}=\frac{D_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}}{{\mathrm{\ϵ}}_0}.$

步骤9中，真空中，磁场的更新计算公式为：

$H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}=\frac{B_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}}{{\mathrm{\μ}}_0}.$

本发明二维柱坐标完全匹配吸收边界的实现方法的有益效果是：步骤3中计算PML中的一维系数数组时，系数数组的计算方法中包含了CFS参数α_ζζ，因此，本发明方法能够更加有效的吸收向外传播的电磁波。

附图说明

图1是本发明实施例1的计算模型示意图；

图2是本发明实施例1计算得到的观测点处的z方向电场强度；

图3是采用现有解析公式计算得到的观测点处的z方向电场强度。

具体实施方式

本发明二维柱坐标完全匹配吸收边界的实现方法，具体步骤如下：

步骤1、输入模型文件。

模型文件内容包括：计算区域大小Nr×Nz，其中，Nr为柱坐标ρ方向网格数，Nz为柱坐标z方向网格数；空间步长Δζ，ζ=ρ,z；时间步长Δt；真空状态下，整个计算区域网格的电导率σ=0，磁导率μ₀，介电常数ε₀；PML吸收边界层数与相关参数κ_ζmax和α_ζζ，κ_ζmax,ζ＝ρ,z为大于1的任一实数，α_ζζ，ζ=ρ,z，ζ=i,k,i+1/2,k+1/2为大于等于0小于1的任一实数，i，k表示空间网格坐标，i＝0,1,2…Nr，k=0,1,2…Nz；仿真计算时长；场源参数。

如图1所示，本实施例用于真空中电偶极子的辐射的计算，场源为z方向电流源，采用微分高斯脉冲激励。其更新公式为：

$E_z^{n+1}\left(\mathrm{0,25}\right)=E_z^{n+1}\left(\mathrm{0,25}\right)-2\×{10}^{-10}\×\left(\frac{4\mathrm{\Δt}}{\mathrm{\Δ\ρ}\·\mathrm{\Δ\ρ}\·\mathrm{\Δz}\·{\mathrm{\ϵ}}_0}\right)\×(t-3\mathrm{\τ})/{\mathrm{\τ}}^2\×\exp (-{(t-3\mathrm{\τ})}^2/{\mathrm{\τ}}^2),$ τ=2ns。

计算区域大小Nr×Nz（50×50），横坐标为柱坐标ρ方向，纵坐标为柱坐标z方向，PML吸收边界层数为10个网格，A点为电场场源所在网格点，B点为电场观测点。空间步长Δζ，ζ=ρ,z，Δρ=Δz=5cm。时间步长Δt=83.333ps。仿真计算时长为100ns。α_ζζ＝0.8×10^-3。

步骤2、根据步骤1输入的模型文件，对电场、磁场及相关辅助变量初始化为零。

电场、磁场及相关辅助变量包括：柱坐标ρ方向二维电场数组E_ρ[Nr][Nz+1]，电位移矢量数组D_ρ[Nr][Nz+1]，柱坐标z方向二维电场数组E_z[Nr+1][Nz],电位移矢量数组D_z[Nr+1][Nz]，柱坐标方向二维磁场数组柱坐标方向二维磁感应强度数组PML中电场辅助变量二维数组ψ_φρ[Nr+1][Nz]，ψ_φγ[Nr+1][Nz]，ψ_φz[Nr][Nz+1]，PML中磁场辅助变量二维数组Φ_ρz[Nr][Nz]和Φ_zρ[Nr][Nz]。

步骤3、计算PML中的一维系数数组：C_hφρi0,C_hφρi1,C_hφρi2，C_hφγi0,C_hφγi1,C_hφγi2，C_hφzk0,C_hφzk1,C_hφzk2，C_eρzk0,C_eρzk1,C_eρzk2，C_ezρi0,C_ezρi1,C_ezρi2。

具体过程为： $C_{\mathrm{h\φ\ρi}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}}},$

$C_{\mathrm{h\φ\ρi}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\ρi}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}-2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\γi}0}=-\frac{\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})-2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}}}{\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}}},$

$C_{\mathrm{h\φ\γi}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\γi}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}-2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φzk}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}}},$

$C_{\mathrm{h\φzk}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}}+2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}})},$

$C_{\mathrm{h\φzk}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}}-2{\mathrm{\ϵ}}_0)}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}})},$

$C_{\mathrm{e\ρzk}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_z+{\mathrm{\σ}}_{\mathrm{zk}+1/2})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_z+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{e\ρzk}1}=\frac{{\mathrm{\Δt\α}}_{\mathrm{zk}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_{\mathrm{zk}+1/2}+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{e\ρzk}2}=\frac{{\mathrm{\Δt\α}}_{\mathrm{zk}+1/2}-2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_{\mathrm{zk}+1/2}+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{ez\ρi}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

$C_{\mathrm{ez\ρi}1}=\frac{{\mathrm{\Δt\α}}_{\mathrm{\ρi}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

$C_{\mathrm{ez\ρi}2}=\frac{{\mathrm{\Δt\α}}_{\mathrm{\ρi}+1/2}-2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

其中，

$A_{\mathrm{\ρi}}=1+({\mathrm{\κ}}_{\mathrm{\ρ}\max }-1)\frac{{({\mathrm{\ρ}}_i-{\mathrm{\ρ}}_0)}^{m+1}}{(m+1)d^m}\frac{1}{{\mathrm{\ρ}}_i},$

$B_{\mathrm{\ρi}}=\frac{{\mathrm{\σ}}_{\max }{({\mathrm{\ρ}}_i-{\mathrm{\ρ}}_0)}^{m+1}}{(m+1)d^m}\frac{1}{{\mathrm{\ρ}}_i};$

ρ₀表示柱坐标ρ方向PML层靠近FDTD区的界面的位置坐标，m为多项式的阶数，取m=4。d为PML层厚度，d＝10Δζ。ρ_i表示柱坐标ρ方向第i个网格的坐标；k_ζmax=10。

κ_ζζ，ζ=ρ,z，ζ=i,k,i+1/2,k+1/2表示ζ网格点处的κ_ζ的取值，κ_ζ的计算公式为：

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

σ_ζζ，ζ=ρ,z，ζ=i,k,i+1/2,k+1/2表示ζ网格点处的σ_ζ的取值，σ_ζ的计算公式为：

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

σ_ζmax=Cσ_opt=C(m+1)/150πΔζ，

C为大于0任一实数。本实施例中C=1。

步骤4、更新计算整个计算区域的电位移矢量。

电位移矢量在z方向的分量的更新公式为：

$D_z{|}_{i,k+1/2}^{n+1}=D_z{|}_{i,k+1/2}^n+\mathrm{\Δt}({\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n-1/2}+C_{\mathrm{h\φ\ρi}1}\frac{{\∂H}_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}+{\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n-1/2}+C_{\mathrm{h\φ\γi}1}\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}),$

电位移矢量在ρ方向的分量的更新公式为：

$D_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}=D_{\mathrm{\ρ}}{|}_{i+1/2,k}^n-\mathrm{\Δt}({\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n-1/2}+C_{\mathrm{h\φzk}1}\frac{{\∂H}_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,k}^{n+1/2}),$

n表示离散时间步，

$\frac{{\∂H}_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{\mathrm{\Δ\ρ}},$

$\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}+H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{2i\·\mathrm{\Δ\ρ}},$

$\frac{{\∂H}_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,k}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}}{\mathrm{\Δz}}.$

步骤5、更新计算整个计算区域的电场。

步骤5中，真空中，电场在z方向分量的更新公式为：

$E_z{|}_{i,k+1/2}^{n+1}=\frac{D_z{|}_{i,k+1/2}^{n+1}}{{\mathrm{\ϵ}}_0},$

电场在ρ方向分量的更新公式为：

$E_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}=\frac{D_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}}{{\mathrm{\ϵ}}_0}.$

步骤6、更新计算电场场源。

在采用微分高斯脉冲激励时，电场场源为z方向的电流源，且电场场源位于网格A点（0，25）处，此时，电场场源更新公式为：

$E_z^{n+1}\left(\mathrm{0,25}\right)=E_z^{n+1}\left(\mathrm{0,25}\right)-2\×{10}^{-10}\×\left(\frac{4\mathrm{\Δt}}{\mathrm{\Δ\ρ}\·\mathrm{\Δ\ρ}\·\mathrm{\Δz}\·{\mathrm{\ϵ}}_0}\right)\×(t-3\mathrm{\τ})/{\mathrm{\τ}}^2\×\exp (-{(t-3\mathrm{\τ})}^2/{\mathrm{\τ}}^2)\left(1\right),$

其中，τ＝2ns。

步骤7、更新计算PML区域中电场辅助变量。

电场辅助变量的更新计算公式为：

${\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n+1/2}=C_{\mathrm{h\φ\ρi}0}{\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n-1/2}+(C_{\mathrm{h\φ\ρi}0}{C}_{\mathrm{h\φ\ρi}1}+C_{\mathrm{h\φ\ρi}2})\frac{{\∂H}_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2},$

${\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n+1/2}=C_{\mathrm{h\φ\γi}0}{\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n-1/2}+(C_{\mathrm{h\φ\γi}0}{C}_{\mathrm{h\φ\γi}1}+C_{\mathrm{h\φ\γi}2})\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2},$

${\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n+1/2}=C_{\mathrm{h\φzk}0}{\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n-1/2}+(C_{\mathrm{h\φzk}0}{C}_{\mathrm{h\φzk}1}+C_{\mathrm{h\φzk}2})\frac{{\∂H}_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,,k}^{n+1/2},$

步骤8、更新计算整个计算区域的磁感应强度。

磁感应强度的更新计算公式为：

其中，

步骤9、更新计算整个计算区域的磁场。

真空中，磁场的更新计算公式为：

$H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}=\frac{B_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}}{{\mathrm{\μ}}_0}.$

步骤10、更新计算PML区域中磁场辅助变量。

磁场辅助变量的更新计算公式为：

${\mathrm{\Φ}}_{\mathrm{\ρz}}{|}_{i+1/2,k+1/2}^n=C_{\mathrm{e\ρzk}0}{\mathrm{\Φ}}_{\mathrm{\ρz}}{|}_{i+1/2,k+1/2}^{n-1}+(C_{\mathrm{e\ρzk}0}{C}_{\mathrm{e\ρzk}1}+C_{\mathrm{e\ρzk}2})\frac{{\∂E}_{\mathrm{\ρ}}}{\∂z}{|}_{i+1/2,k+1/2}^n,$

${\mathrm{\Φ}}_{\mathrm{z\ρ}}{|}_{i+1/2,k+1/2}^n=C_{\mathrm{ez\ρi}0}{\mathrm{\Φ}}_{\mathrm{z\ρ}}{|}_{i+1/2,k+1/2}^{n-1}+(C_{\mathrm{ez\ρi}0}{C}_{\mathrm{ez\ρi}1}+C_{\mathrm{ez\ρi}2})\frac{{\∂E}_z}{{\∂}_{\mathrm{\ρ}}}{|}_{i+1/2,k+1/2}^n,$

步骤11、根据步骤1设定的仿真计算时长，判断是否继续更新；实际计算时间小于该设定值时，继续计算，返回至步骤3；否则输出步骤5计算得到的电场和步骤9计算得到的磁场，结束。

如图2所示，是本发明实施例1计算的观测点B处（25，10）的z方向电场强度。如图3所示，采用现有解析方法计算得到的观测点处的z方向电场强度。从图2与图3对比可见，本发明方法计算与解析计算结果一致，验证了本发明方法有效的实现了二维柱坐标完全匹配吸收边界。

Claims (3)

1.一种二维柱坐标完全匹配吸收边界的实现方法，其特征在于，具体步骤如下：

步骤1、输入模型文件：

所述模型文件内容包括：计算区域大小Nr×Nz，其中，Nr为柱坐标ρ方向网格数，Nz为柱坐标z方向网格数；空间步长Δζ，ζ＝ρ,z；时间步长Δt；真空状态下，整个计算区域网格的电导率σ＝0，磁导率μ₀，介电常数ε₀；PML吸收边界层数与相关参数κ_ζmax和κ_ζmax为大于1的任一实数，ζ＝ρ,z，为大于等于0小于1的任一实数，ζ＝ρ,z，i，k表示空间网格坐标，i＝0,1,2…Nr，k＝0,1,2…Nz；仿真计算时长；场源参数；

步骤2、根据步骤1输入的模型文件，对电场、磁场及相关辅助变量初始化为零；

其中，所述电场、磁场及相关辅助变量包括：柱坐标ρ方向二维电场数组E_ρ[Nr][Nz+1]，电位移矢量数组D_ρ[Nr][Nz+1]，柱坐标z方向二维电场数组E_z[Nr+1][Nz],电位移矢量数组D_z[Nr+1][Nz]，柱坐标方向二维磁场数组柱坐标方向二维磁感应强度数组PML中电场辅助变量二维数组ψ_φρ[Nr+1][Nz]，ψ_φγ[Nr+1][Nz]，ψ_φz[Nr][Nz+1]，PML中磁场辅助变量二维数组Φ_ρz[Nr][Nz]和Φ_zρ[Nr][Nz]；

步骤3、计算PML中的电场，电位移矢量，磁场，磁感应强度以及辅助变量更新所需要用到的一维系数数组：C_hφρi0,C_hφρi1,C_hφρi2，C_hφγi0,C_hφγi1,C_hφγi2，C_hφzk0,C_hφzk1,C_hφzk2，C_eρzk0,C_eρzk1,C_eρzk2，C_ezρi0,C_ezρi1,C_ezρi2，分别用于电场辅助变量和磁场辅助变量的计算公式更新，其中各系数中符号含义为：下标ρ表示柱坐标中的ρ方向，下标z表示柱坐标中的z方向，下标i表示ρ方向第i个网格，下标k表示z方向第k个网格，下标0、1、2用于区别一个辅助变量中更新的三个系数；

具体过程为： $C_{\mathrm{h\φ\ρi}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}}},$

$C_{\mathrm{h\φ\ρi}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0\mathrm{})}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\ρi}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0\mathrm{})}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+{\mathrm{\σ}}_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\γi}0}=-\frac{\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})-2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}}}{\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}}},$

$C_{\mathrm{h\φ\γi}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0\mathrm{})}{(\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φ\γi}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}}+2{\mathrm{\ϵ}}_0\mathrm{})}{(\mathrm{\Δt}(A_{\mathrm{\ρi}}{\mathrm{\α}}_{\mathrm{\ρi}}+B_{\mathrm{\ρi}})+2{\mathrm{\ϵ}}_0{A}_{\mathrm{\ρi}})},$

$C_{\mathrm{h\φzk}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}}},$

$C_{\mathrm{h\φzk}1}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}}+2{\mathrm{\ϵ}}_0\mathrm{})}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}})},$

$C_{\mathrm{h\φzk}2}=\frac{(\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}}+2{\mathrm{\ϵ}}_0\mathrm{})}{(\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}}{\mathrm{\α}}_{\mathrm{zk}}+{\mathrm{\σ}}_{\mathrm{zk}})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}})},$

$C_{\mathrm{e\ρzk}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_z+{\mathrm{\σ}}_{\mathrm{zk}+1/2})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_z+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{e\ρzk}1}=\frac{\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_{\mathrm{zk}+1/2}+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{e\ρzk}2}=\frac{\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{zk}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{zk}+1/2}{\mathrm{\α}}_{\mathrm{zk}+1/2}+{\mathrm{\σ}}_{\mathrm{zk}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{zk}+1/2}},$

$C_{\mathrm{ez\ρi}0}=-\frac{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})-2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

$C_{\mathrm{ez\ρi}1}=\frac{\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

$C_{\mathrm{ez\ρi}2}=\frac{\mathrm{\Δt}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+2{\mathrm{\ϵ}}_0}{\mathrm{\Δt}({\mathrm{\κ}}_{\mathrm{\ρi}+1/2}{\mathrm{\α}}_{\mathrm{\ρi}+1/2}+{\mathrm{\σ}}_{\mathrm{\ρi}+1/2})+2{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{\mathrm{\ρi}+1/2}},$

其中，A_ρi，B_ρi为计算所用的中间变量，其计算公式为：

$A_{\mathrm{\ρi}}=1+({\mathrm{\κ}}_{\mathrm{\ρ}\max }-1)\frac{{({\mathrm{\ρ}}_i-{\mathrm{\ρ}}_0)}^{m+1}}{(m+1)d^m}\frac{1}{{\mathrm{\ρ}}_i},$

$B_{\mathrm{\ρi}}=\frac{{\mathrm{\σ}}_{\max }{({\mathrm{\ρ}}_i-{\mathrm{\ρ}}_0)}^{m+1}}{(m+1)d^m}\frac{1}{{\mathrm{\ρ}}_i};$

ρ₀表示柱坐标ρ方向PML层靠近FDTD区的界面的位置坐标，m为多项式的阶数，d为PML层厚度，ρ_i表示柱坐标ρ方向第i个网格的坐标；

表示ζ方向网格点处的κ_ζ的取值，κ_ζ的计算公式为：

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

表示ζ方向网格点处的σ_ζ的取值，σ_ζ的计算公式为：

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

σ_ζmax＝Cσ_opt＝C(m+1)/150πΔζ，

C为大于0任一实数；

步骤4、更新计算整个计算区域的电位移矢量；

其中，电位移矢量在z方向的分量的更新公式为：

$D_z{|}_{i,k+1/2}^{n+1}=D_z{|}_{i,k+1/2}^n+\mathrm{\Δt}({\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n-1/2}+C_{\mathrm{h\φ\ρi}1}\frac{\∂H_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}+{\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n-1/2}+C_{\mathrm{h\φ\γi}1}\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}),$

电位移矢量在ρ方向的分量的更新公式为：

$D_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}=D_{\mathrm{\ρ}}{|}_{i+1/2,k}^n-\mathrm{\Δt}({\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n-1/2}+C_{\mathrm{h\φzk}1}\frac{\∂H_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,k}^{n+1/2}),$

n表示离散时间步，

$\frac{\∂H_{\mathrm{\φ}}}{\∂\mathrm{\φ}}{|}_{i,k+1/2}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{\mathrm{\Δ\ρ}},$

$\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{2i\·\mathrm{\Δ\ρ}},$

$\frac{\∂H_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,k}^{n+1/2}=\frac{H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\mathrm{\φ}}{|}_{i-1/2,k+1/2}^{n+1/2}}{\mathrm{\Δz}};$

步骤5、更新计算整个计算区域的电场；

步骤6、更新计算电场场源；

步骤7、更新计算PML区域中电场辅助变量；

其中，电场辅助变量的更新计算公式为：

${\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n+1/2}=C_{\mathrm{h\φ\ρi}0}{\mathrm{\ψ}}_{\mathrm{\φ\ρ}}{|}_{i,k+1/2}^{n-1/2}+(C_{\mathrm{h\φ\ρi}0}{C}_{\mathrm{h\φ\ρi}1}+C_{\mathrm{h\φ\ρi}2})\frac{\∂H_{\mathrm{\φ}}}{\∂\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2},$

${\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n+1/2}=C_{\mathrm{h\φ\γi}0}{\mathrm{\ψ}}_{\mathrm{\φ\γ}}{|}_{i,k+1/2}^{n-1/2}+(C_{\mathrm{h\φ\γi}0}{C}_{\mathrm{h\φ\γi}1}+C_{\mathrm{h\φ\γi}2})\frac{H_{\mathrm{\φ}}}{\mathrm{\ρ}}{|}_{i,k+1/2}^{n+1/2},$

${\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n+1/2}=C_{\mathrm{h\φzk}0}{\mathrm{\ψ}}_{\mathrm{\φz}}{|}_{i+1/2,k}^{n-1/2}+(C_{\mathrm{h\φzk}0}{C}_{\mathrm{h\φzk}1}+C_{\mathrm{h\φzk}2})\frac{\∂H_{\mathrm{\φ}}}{\∂z}{|}_{i+1/2,,k}^{n+1/2},$

步骤8、更新计算整个计算区域的磁感应强度；

其中，磁感应强度的更新计算公式为：

其中， $\frac{{\∂E}_{\mathrm{\ρ}}}{\∂z}{|}_{i+1/2,k+1/2}^n=\frac{E_{\mathrm{\ρ}}{|}_{i+1/2,k+1}^n-E_{\mathrm{\ρ}}{|}_{i+1/2,k}^n}{\mathrm{\Δz}},\frac{\∂E_z}{\∂\mathrm{\ρ}}{|}_{i+1/2,k+1/2}^n=\frac{E_z{|}_{i+1,k+1/2}^n-E_z{|}_{i-1,k+1/2}^n}{\mathrm{\Δ\ρ}};$

步骤9、更新计算整个计算区域的磁场；

步骤10、更新计算PML区域中磁场辅助变量；

其中，磁场辅助变量的更新计算公式为：

${\mathrm{\Φ}}_{\mathrm{\ρz}}{|}_{i+1/2,k+1/2}^n=C_{\mathrm{e\ρzk}0}{\mathrm{\Φ}}_{\mathrm{\ρz}}{|}_{i+1/2,k+1/2}^{n-1}+(C_{\mathrm{e\ρzk}0}{C}_{\mathrm{e\ρzk}1}+C_{\mathrm{e\ρzk}2})\frac{\∂E_{\mathrm{\ρ}}}{\∂\mathrm{\ρ}}{|}_{i+1/2,k+1/2}^n,$

${\mathrm{\Φ}}_{\mathrm{z\ρ}}{|}_{i+1/2,k+1/2}^n=C_{\mathrm{ez\ρi}0}{\mathrm{\Φ}}_{\mathrm{z\ρ}}{|}_{i+1/2,k+1/2}^{n-1}+(C_{\mathrm{ez\ρi}0}{C}_{\mathrm{ez\ρi}1}+C_{\mathrm{ez\ρi}2})\frac{\∂E_z}{\∂\mathrm{\ρ}}{|}_{i+1/2,k+1/2}^n,$

步骤11、根据步骤1设定的仿真计算时长，判断是否继续更新；实际计算时间小于该设定值时，继续计算，返回至步骤3；否则输出步骤5计算得到的电场和步骤9计算得到的磁场，结束。

2.按照权利要求1所述的二维柱坐标完全匹配吸收边界的实现方法，其特征在于，步骤5中，真空中电场在z方向分量的更新公式为：

$E_z{|}_{i,k+1/2}^{n+1}=\frac{D_z{|}_{i,k+1/2}^{n+1}}{{\mathrm{\ϵ}}_0},$

电场在ρ方向分量的更新公式为：

$E_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}=\frac{D_{\mathrm{\ρ}}{|}_{i+1/2,k}^{n+1}}{{\mathrm{\ϵ}}_0}.$

3.按照权利要求1所述的二维柱坐标完全匹配吸收边界的实现方法，其特征在于，步骤9中，真空中磁场的更新计算公式为：

$H_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}=\frac{B_{\mathrm{\φ}}{|}_{i+1/2,k+1/2}^{n+1/2}}{{\mathrm{\μ}}_0}.$