二維MPSTD算法的研究

張建萍

（華北電力大學(xué) 電氣工程學(xué)院，河北 保定071003）

FPSTD算法于1997年最先被Q.H.Liu引入到電磁學(xué)中［1］，它的基本思想是：對(duì) Maxwell方程中時(shí)間微商利用中心差分或高階差分近似，對(duì)空間微分采用Fourier變換代替，并采用快速傅立葉變換（FFT）技術(shù)對(duì)離散傅立葉變換進(jìn)行計(jì)算。與FDTD方法相比，F(xiàn)PSTD算法由于采用Fourier變換及反變換來實(shí)現(xiàn)Maxwell方程中場量空間導(dǎo)數(shù)的計(jì)算，理論上精度可以無限高，所以不存在由于空間網(wǎng)格劃分帶來的誤差。在滿足Nyquist采樣定理的情況下，在每一個(gè)最小波長上只需要設(shè)置2個(gè)采樣點(diǎn)，即可實(shí)現(xiàn)對(duì)電磁問題的準(zhǔn)確分析。目前，該算法已經(jīng)被成功地應(yīng)用于求解空間簡單目標(biāo)的散射特性、建筑物內(nèi)電波場強(qiáng)分布和模擬探地雷達(dá)分析地下簡單目標(biāo)體的散射特性等。然而，由于FFT具有固有的周期性，F(xiàn)PSTD算法只適用于分析周期問題，如均勻介質(zhì)空間或者變化不大的非均勻介質(zhì)空間。若將其用于復(fù)雜媒質(zhì)電磁問題的分析時(shí)，如金屬物體和擁有不連續(xù)媒質(zhì)分布的問題時(shí)，由于Gibbs現(xiàn)象的影響，它的計(jì)算精度是較低的；其次，若將該算法用于復(fù)雜幾何形狀目標(biāo)問題的分析時(shí)，目標(biāo)外形的階梯近似使得該算法不能對(duì)目標(biāo)外形進(jìn)行準(zhǔn)確描述，造成了計(jì)算精度的下降。上述兩點(diǎn)使得該方法的應(yīng)用受到極大的約束。MPSTD算法于1996年被A.V.Kabakian引入到計(jì)算電磁學(xué)［2］。它的基本思想是：將整個(gè)計(jì)算域剖分成若干個(gè)與所分析電磁問題共形的曲邊四邊形子域（二維電磁問題）或曲面六面體子域（三維電磁問題），切比雪夫選配方法被單獨(dú)用在每個(gè)子域里進(jìn)行計(jì)算，子域之間通過子域分界面上的匹配邊界條件來實(shí)現(xiàn)信息的傳遞。已有的文獻(xiàn)表明：MPSTD算法每最小波長只需要設(shè)置π個(gè)采樣點(diǎn)即可實(shí)現(xiàn)對(duì)電磁問題的準(zhǔn)確分析。由于共形子域劃分的應(yīng)用和多域策略的引入，MPSTD算法比FPSTD算法具有更大的靈活性和適用性［3，4］。

由于PSTD算法采用傅立葉選配方法或切比雪夫選配方法來求解Maxwell方程中場量的空間導(dǎo)數(shù)，所以均具有譜域計(jì)算精度。同時(shí)，由于PSTD方法每最小波長只需要設(shè)置2個(gè)或π個(gè)采樣點(diǎn)，因而它對(duì)計(jì)算機(jī)的內(nèi)存需求大大降低，計(jì)算量也同時(shí)大大減少，已有大量文獻(xiàn)表明PSTD算法在許多應(yīng)用方面都顯得比FDTD更加有效。該算法作為一種新穎的電磁場數(shù)值方法，在計(jì)算電磁學(xué)中已經(jīng)得到越來越多的重視。

從以上的簡要描述中可以看出多區(qū)域偽譜時(shí)域算法是一種高階的、靈活的、有效的數(shù)值算法，具有很強(qiáng)的實(shí)用性。本文就二維多區(qū)域偽譜時(shí)域算法展開理論的研究［5－9］。

1 曲線坐標(biāo)下的Maxwell方程

由Maxwell方程

考慮一個(gè)各向同性，導(dǎo)電的不均勻媒質(zhì)，具有介電常數(shù)ε、磁導(dǎo)率μ和電導(dǎo)率σ。對(duì)于二維TMz極化的Maxwell方程可表達(dá)為

矩陣形式為：

式中，q＝（Hx，Hy，Ez）T。

采用坐標(biāo)變換將坐標(biāo)系（x，y）中具有曲邊界的子區(qū)域映射成坐標(biāo)系（ξ，η）中的單位正方形。引入坐標(biāo)變換關(guān)系為：ξ＝ξ（x，y），η＝η（x，y），式（6）變?yōu)椋?/p>

式中，q＝（Hx，Hy，Ez）T

且ξx＝?ξ／?x，ξy＝?ξ／?y，ηx＝?η／?x，ηy＝?η／?y。

2 切比雪夫選配方法

由于切比雪夫選配方法具有十分優(yōu)越的近似特性，因此該方法被廣泛地應(yīng)用于偏微分方程的求解中。K階切比雪夫多項(xiàng)式定義為：

TN（x）＝cos（Ncos－1（x）） 其中｜x｜≤1

在坐標(biāo)系（ξ，η）中，定義單位正方形上的Chebyshev－Gauss－Lobatto配置點(diǎn)為：

采用張量積Chebyshev－Lagrange多項(xiàng)式，q（ξ，η）為

式中，gi（α）為Lagrange差值多項(xiàng)式

式中，α＝（ξ，η）。

于是場量q（ξ，η）定義在網(wǎng)格點(diǎn)上的空間偏導(dǎo)數(shù)為：

3 特征變量子域匹配條件

為了獲得整個(gè)區(qū)域的解，就需要采用匹配條件將每個(gè)子區(qū)域中的解聯(lián)系起來。特征變量法是一種有效地施加邊界條件的方式。將式（7）中的矩陣A對(duì)角化可得

其中

則特征矢量可得

式中，特征變量R－代表沿－ξ傳播的特征波的幅度；特征變量R＋代表沿＋ξ傳播的特征波的幅度；特征變量R0代表不傳播的特征波的幅度。對(duì)角矩陣Λ中對(duì)角線上的元素對(duì)應(yīng)著特征波傳播的速度。矩陣B的特征變量類似于矩陣A的特征變量。

一般情況下在兩個(gè)子區(qū)域之間有兩類交界面：一類是兩個(gè)子區(qū)域具有不同的媒質(zhì)；另一類是兩個(gè)子區(qū)域具有相同的媒質(zhì)。對(duì)這兩類情況的處理方法是不同的。首先考慮兩個(gè)相鄰的子區(qū)域具有相同媒質(zhì)的情況。假定ξ軸是從子區(qū)域1指向子區(qū)域2。在這種境況下，在交界面上匹配特征變量可得

式中，上標(biāo)（1）和（2）分別定義為子區(qū)域1和子區(qū)域2。式（17）～式（19）分別代表的含義如下：

（1）入射的特征變量波是由相鄰子區(qū)域的出射的特征變量波所決定；

（2）不傳播的特征變量波保持連續(xù)；

（3）出射的特征變量波不改變地進(jìn)入下一子區(qū)域。

對(duì)于兩個(gè)相鄰的子區(qū)域具有不同媒質(zhì)的情況，在交界面上施加物理邊界條件。對(duì)于介質(zhì)交界面，保持切向電磁場連續(xù)：

式中，下標(biāo)t表示交界面處的切向場分量。

對(duì)于理想導(dǎo)體，切向電場和法向磁場為零：

式中，下標(biāo)n表示交界面處的法向場分量。

4 完全匹配層

在時(shí)域方法中一個(gè)關(guān)鍵的問題就是采用吸收邊界條件去截?cái)酂o限大的計(jì)算區(qū)域，使得整個(gè)計(jì)算在有限區(qū)域中進(jìn)行且在吸收邊界處不引入反射。這里采用一種系統(tǒng)的方法去推導(dǎo)良態(tài)提出的（well－posed）有耗媒質(zhì)理想匹配層（PML）公式。

根據(jù)復(fù)坐標(biāo)構(gòu)架技術(shù)引入復(fù)坐標(biāo)構(gòu)架變量

式中，aη（η＝x，y）是一個(gè)標(biāo)量系數(shù)；wη（η＝x，y）是衰減系數(shù)。為了簡單起見，令ax＝ay＝1。

下面推導(dǎo)有耗媒質(zhì)PML公式：

將式（3）～（5）變成頻域形式

經(jīng)過適當(dāng)運(yùn)算，變成時(shí)域形式得：

則式（27）～（29）變?yōu)椋?/p>

注意到上面的PML方程是不分裂場形式的。除了一些低階項(xiàng)以外，它和原始的Maxwell方程一樣都是對(duì)稱的雙曲系統(tǒng)，因此這個(gè)PML是 Well－posed。當(dāng)σ＝0時(shí)，有耗媒質(zhì)的ABC就變?yōu)闊o耗的ABC；當(dāng)wx＝wy＝0時(shí)，PML方程就退化為最初的 Maxwell方程。

5 時(shí)間積分方法

采用五步、四階、低存儲(chǔ)的Runge－Kutta方法進(jìn)行時(shí)間積分。將式（7）表示為

定義qn＝q（tn），tn＝nΔt，其中Δt是時(shí)間步間隔，則對(duì)于2 N存儲(chǔ)，M步K 階Runge－Kutta方法可表示為

6 結(jié) 論

本文分析了便于復(fù)雜媒質(zhì)求解的MPSTD技術(shù)，探討了MPSTD所使用的切比雪夫選配方法，子域匹配，完全匹配層（PML）。MPSTD算法的實(shí)現(xiàn)歸結(jié)為以下五個(gè)步驟：（1）計(jì)算模型的建立。首先，根據(jù)待分析問題的特征，選擇合理的坐標(biāo)系（本論文選擇直角坐標(biāo)系），建立準(zhǔn)確的計(jì)算模；其次，根據(jù)計(jì)算模型的特征，進(jìn)行合理的子域剖分，并根據(jù)每個(gè)子域的特點(diǎn)，選擇每個(gè)子域剖分的網(wǎng)格數(shù)目；最后，選擇合理的子域到單位正方形之間的坐換關(guān)系式；（2）時(shí)間步推進(jìn)。選擇合理的時(shí)間步方法處理Maxwell微分方程組中關(guān)于時(shí)間的偏微分計(jì)算（本論文選擇Runge－Kutta法）；（3）空間導(dǎo)數(shù)的計(jì)算。在各子域內(nèi)部采用切比雪夫譜選配方法計(jì)算Maxwell微分方程組中的空間導(dǎo)數(shù)，并采用匹配邊界條件處理子域間分界面上的場（本論文選擇基于特征變量CV的子域分界面匹配條件）；（4）吸收邊界處理。MPSTD算法和其它數(shù)值計(jì)算方法一樣，在計(jì)算“開放”電磁系統(tǒng)時(shí)，必須人為設(shè)定一個(gè)吸收邊界進(jìn)行截?cái)啵ū菊撐倪x擇的是PML吸收邊界條件）；（5）計(jì)算結(jié)果的提取和輸出。根據(jù)仿真目的的不同，選擇相應(yīng)的參數(shù)提取方式，實(shí)現(xiàn)算法的后處理。

［1］Liu Q H.The PSTD Algorithm：a Time－domain Method Requiring Only Two Cells Per Wavelength［J］.Microwave Opt Technol Lett，1997，15（3）：158－165.

［2］Kabakian A V.A three－dimensional spectral collocation time－domain solver for electromagnetic wave scattering［C］.New York：Aircraft Industries Association of America，1997.

［3］楊 虎.多區(qū)域時(shí)域偽譜算法在電磁分析中的理論和應(yīng)用研究［D］.長沙：國防科學(xué)技術(shù)大學(xué)，2006.

［4］史 琰.高階時(shí)域電磁計(jì)算方法［D］.西安：西安電子科技大學(xué)，2005.

［5］姜永金，柴舜連，毛鈞杰.MPSTD算法子域分界面匹配條件的比較研究［J］.電波科學(xué)學(xué)報(bào)，2010，（1）：18－20.

［6］Yang B，Hesthaven J S.Multidomain pseudospectral computation of Maxwell＇s equations in 3Dgeneral curvilinear coordinates［J］.Appl.Numerical Math，2000，33：281－289.

［7］Fan G X，Liu Q H，Hesthaven J S.Multidomain pesudospectral time－domain simulations of scattering by objects buried in lossy media［J］.IEEE Transactions on Geoscience and Remote Sensing，2002，40（6）：1366－1373.

［8］Fan G X，Liu Q H.A Well－Posed PML Absorbing Boundary Condition for Lossy Media［C］.IEEE Antennas and Propagation Society International Symposium，2001.

［9］Fan G X，Liu Q H.A Strongly Well－Posed PML in Lossy Media［C］.IEEE Antennas and wireless propagation letters，2003，（2）：97－100.