[发明专利]利用混合的三棱柱—四面体网格实现DGTD中PML的方法

权利要求书

1.利用混合的三棱柱—四面体网格实现DGTD中PML的方法，其特征在于，包括以下步骤：

S1、根据实际情况选取计算模型，根据具体模型设计PML吸收层，并分别建立计算域中PML区和非PML区的麦克斯韦方程；

S2、进行网格划分，并保证周期边界主面和从面上的网格匹配；具体实现方法为：对Ω_TOT区域采用四面体网格划分；

对Ω_PML区域采用三棱柱网格划分，具体操作方法为：在进行Ω_TOT区域的网格划分时在Ω_TOT和Ω_PML的交界面上形成多个三角形，在进行Ω_PML区域的网格划分时将会以该交界面上的三角形作为三棱柱的下底面，然后根据所需要的网格密度设定棱柱的高，如果该区域的三棱柱网格是直三棱柱那么该三棱柱的上底面和下底面完全一致，如果是共形三棱柱则该三棱柱的上底面的形状和下底面形状一致；当划分完一层之后所形成的新交界面必须保证和Ω_TOT和Ω_PML的交界面形状一致；然后重复上述步骤直到将PML区域划分完毕；

S3、选择基函数，将电场强度矢量和磁场强度矢量在不同区域的网格采用不同的基函数进行展开；具体实现方法为：

根据步骤S1建立的麦克斯韦方程得到以下公式：

公式中的ke、kh、ve、vh分别取1/2、1/2、1/(2*ym)、1/(2*zm)，ym为导纳、zm为阻抗；Φ^m为单元的基函数，E^m和H^m为本单元的电场和磁场，E^m+和H^m+为相邻单元的电场和磁场，dv为体积分单元，ds为面积分单元；τ_m表示的是四面体网格；n为边界面上的单位外法向矢量，Γ为相交面；将PML区域和非PML区域的相对磁导率和相对介电质常数统一表示为和

对公式(3)中的电场和磁场采用棱边基函数进行离散：

其中和为t时刻的待求解的系数，P为每个网格单元的矢量棱边基函数的个数；表示第m个网格单元的第i条矢量棱边基函数，x,y,z表示是坐标参数，t表示的是在t时刻；

S4、运用DGTD算法先进行空间离散再进行时间离散得到其时间离散方程；具体实现方法为：

非PML区空间离散和时间离散如下：

PML区空间离散和时间离散如下：

M表示质量矩阵，S表示刚度矩阵，表示PML区中的质量矩阵，F_ke,F_ve,F_kh,F_vh分别表示通量矩阵，M,S,F_ke,F_ve,F_kh,F_vh分别通过公式(9)计算；Δt为时间步长；Μ^m、J^m、和表示的是第m个网格单元的中间参量，无具体物理意义；dt为时间积分微元；表示的是第m个网格单元的电场系数在第n*dt时刻的值，表示第m个网格单元的电场系数在第(n+1)*dt时刻的值，表示第m+1个网格单元的电场系数在第n*dt时刻的值；表示第m个网格单元的磁场系数在第(n+1/2)*dt时刻的值，表示第m个网格单元的磁场系数在第(n-1/2)*dt时刻的值，表示第m+1个网格单元的磁场系数在第(n-1/2)*dt时刻的值，表示第m+1个网格单元的磁场系数在第(n+1/2)*dt时刻的值；

在公式(9)中，v＝{ve,vh}，k＝{ke,kh}，qq'表示矩阵的第q行第q'列；由公式(11)求解；δ_m,δ_e分别表示非PML区中的相对磁导率和相对电导率；

其中：

其中，(x,y,z)为坐标参量，(x₀,y₀,z₀)为交界面出的坐标值，d为PML层厚度，ε_r为相对电介质常数；

S5、集成单元矩阵，根据条件得到所求的单元矩阵；具体实现方法为：在非PML区时处理M和S矩阵将四面体的每个棱边基函数按照公式(9)处理，得到所求的单元矩阵；对于F矩阵按照公式(9)进行处理，得到所求的单元矩阵；对于PML区和非PML区交界面的时候需要对公式(5)中的F⁺矩阵进行特殊处理，由于分界面两端的基函数表现形式不一样，对面积部分的积分采用高斯积分来进行处理；

根据公式(1)和公式(2)在交界面有以下形式的公式部分需要进行特殊处理：

设在PML区和非PML区的交界面为Γ，根据公式(9)和公式(13)看出需要对公式(8)中的进行特殊处理：设属于PML区的时候其棱边基函数为Φ_iΦ_jΦ_k，属于非PML区的时候其棱边基函数为将与Φ_iΦ_jΦ_k形成对应关系，设Φ_iΦ_jΦ_k分别对应于将公式(9)修改为：

上述公式中Φ_q代表的是Φ_iΦ_jΦ_k而代表的是等号左端部分代表的是已知公式的变形；当这些矩阵求解完毕之后同时需要对相交面Γ上的E^m+和H^m+进行调整有：

公式(15)说明的是将相邻单元的第a,b,c条棱边上的电场磁场系数，在与公式(14)的矩阵相乘时调整到位置i,j,k上，同时变换矩阵维数；

S6、进行时间迭代得到每个时刻电场磁场值。