二維MPSTD算法的研究

張建萍

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

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

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

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

1 曲線坐標下的Maxwell方程

由Maxwell方程

考慮一個各向同性，導電的不均勻媒質，具有介電常數ε、磁導率μ和電導率σ。對于二維TMz極化的Maxwell方程可表達為

矩陣形式為：

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

采用坐標變換將坐標系（x，y）中具有曲邊界的子區域映射成坐標系（ξ，η）中的單位正方形。引入坐標變換關系為：ξ＝ξ（x，y），η＝η（x，y），式（6）變為：

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

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

2 切比雪夫選配方法

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

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

在坐標系（ξ，η）中，定義單位正方形上的Chebyshev－Gauss－Lobatto配置點為：

采用張量積Chebyshev－Lagrange多項式，q（ξ，η）為

式中，gi（α）為Lagrange差值多項式

式中，α＝（ξ，η）。

于是場量q（ξ，η）定義在網格點上的空間偏導數為：

3 特征變量子域匹配條件

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

其中

則特征矢量可得

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

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

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

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

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

（3）出射的特征變量波不改變地進入下一子區域。

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

式中，下標t表示交界面處的切向場分量。

對于理想導體，切向電場和法向磁場為零：

式中，下標n表示交界面處的法向場分量。

4 完全匹配層

在時域方法中一個關鍵的問題就是采用吸收邊界條件去截斷無限大的計算區域，使得整個計算在有限區域中進行且在吸收邊界處不引入反射。這里采用一種系統的方法去推導良態提出的（well－posed）有耗媒質理想匹配層（PML）公式。

根據復坐標構架技術引入復坐標構架變量

式中，aη（η＝x，y）是一個標量系數；wη（η＝x，y）是衰減系數。為了簡單起見，令ax＝ay＝1。

下面推導有耗媒質PML公式：

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

經過適當運算，變成時域形式得：

則式（27）～（29）變為：

注意到上面的PML方程是不分裂場形式的。除了一些低階項以外，它和原始的Maxwell方程一樣都是對稱的雙曲系統，因此這個PML是 Well－posed。當σ＝0時，有耗媒質的ABC就變為無耗的ABC；當wx＝wy＝0時，PML方程就退化為最初的 Maxwell方程。

5 時間積分方法

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

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

6 結 論

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

［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］楊 虎.多區域時域偽譜算法在電磁分析中的理論和應用研究［D］.長沙：國防科學技術大學，2006.

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

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