[發明專利]利用混合的三棱柱—四面體網格實現DGTD中PML的方法

有效

• 申請號: 201711471544.6
• 申請日: 2017-12-29
• 公開（公告）號: CN108229000B
• 公開（公告）日: 2020-06-09
• 發明（設計）人: 徐立;唐鵬飛;楊中海;李斌
• 申請（專利權）人: 電子科技大學
• 主分類號: G06F30/23
• 分類號: G06F30/23
• 代理公司: 成都虹盛匯泉專利代理有限公司 51268
• 代理人: 王偉

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、進行時間迭代得到每個時刻電場磁場值。