二維MPSTD算法的研究

張建萍

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

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

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

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

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

由Maxwell方程

考慮一個各向同性，導(dǎo)電的不均勻媒質(zhì)，具有介電常數(shù)ε、磁導(dǎo)率μ和電導(dǎo)率σ。對于二維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 特征變量子域匹配條件

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

其中

則特征矢量可得

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

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

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

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

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

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

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

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

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

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

4 完全匹配層

在時域方法中一個關(guān)鍵的問題就是采用吸收邊界條件去截斷無限大的計算區(qū)域，使得整個計算在有限區(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）是一個標(biāo)量系數(shù)；wη（η＝x，y）是衰減系數(shù)。為了簡單起見，令ax＝ay＝1。

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

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

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

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

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

5 時間積分方法

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

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

6 結(jié) 論

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

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

［5］姜永金，柴舜連，毛鈞杰.MPSTD算法子域分界面匹配條件的比較研究［J］.電波科學(xué)學(xué)報，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.