矩量法在電磁散射中的應用介紹

1、矩量法在電磁散射中的應用一矩量法在電磁散射問題中的應用電磁散射問題是電磁學中的一個重要研究領域，研究電磁波的散射機理以及計算其散射場強的大小與分布，具有十分重要的實際意義。矩量法作為一種有效的數(shù)值計算方法在其中有著廣泛的應用。但作為一種計算方法它也有著自己的缺陷，為了解決這些問題，人們提出了各種方 案，矩量法在這個過程中也獲得了很大的發(fā)展。MoM(Method of Moments)原本是一種近似求解線性算子方程的方 法，通過它可以將算子方程轉(zhuǎn)化為一矩陣方程，進而通過求解此矩陣方 程得到最終的近似解。MoM最早是由兩位數(shù)學家L. V. Kantorovich和V. I. Krylov提出的，后

2、來由 K. K. Mei引入計算電磁學，最終被 R. F. Harryington在其著作計算電磁場中的矩量法中加以系統(tǒng)描述。利用 矩量法求解電磁問題的主要優(yōu)點是：它嚴格地計算了各個子系統(tǒng)間的互耦，而算法本身又從根本上保證了誤差系統(tǒng)總體最小而不產(chǎn)生數(shù)值色 散。如今MoM被廣泛應用于計算電磁學中，雖然它不能處理電大尺寸 目標的電磁問題，但基于 MoM的各種加速方法仍受到極大重視，如多 層快速多極子方法 MLMFA等。電磁散射問題是電磁學中的一個重要研究領域，研究電磁波的散射機理以及計算其散射場強的大小與分布，具有十分重要的實際意義。在實際生活中，遇到的散射目標往往不僅具有復雜的幾何形狀，而且構(gòu)成

3、的材料也各不相同。因此對復雜目標的電磁散射特性進行快速、高效的 分析，具有重要的理論意義和實用價值。電磁散射問題只有在相對簡單的情況下才可以用嚴格的解析法來 求解，比如對極少數(shù)形狀規(guī)則的物體。對于電大物體，可以用高頻近似 方法，例如幾何光學法（GOb物理光學法（PO）、幾何繞射理論（GTD物 理繞射理論（PTD） 一致性幾何繞射理論（UTD）、復射線法（CT簿來求解散 射場。反之，對于電小物體，可以用準靜態(tài)場來進行分析。介乎這兩者 之間的物體，一般采用數(shù)值方法。目前分析電磁散射問題的數(shù)值方法主 要有微分方程法和積分方程法。微分方程法有有限差分法（FDM）、時域有限差分法（FDTD）頻域有BM差

4、分法（FDFD）時域平面波法（PWTD時 域多分辨分析法（MRTD有限元法（FEM傳輸線矩陣法（TLM）等，積分 方程法有表面積分方程法（SIEM）矩量法（MOM）、邊界元法（BEM）、體積 分方程法（VIEM）、快速多極子法（FMM）、時域積分方程法（IETD殍。這些 方法各有優(yōu)缺點，有的是為了避免矩陣求逆，有的是為了加快收斂，有 的是為了提高精度，還有的是為了減少貯存等。它們被廣泛應用于求解 復雜的工程電磁場問題中。應用微分方程法求解電磁散射問題時，由于散射體的外空間為無限大，需要人為設置截斷邊界使求解區(qū)域有限， 這種截斷邊界的引入會導致非物理的反射現(xiàn)象。矩量法是一種將連續(xù)方程離散化成代數(shù)

5、方程組的方法，經(jīng)常被看作數(shù)值精確解。它既適用于電磁場積分方程又適用于微分方程，由于已經(jīng)有有效的數(shù)值計算方法求解微分方程，所以目前矩量法大都用來求解積分方程。由于此方法能解決邊界比較復雜的一些問題，因而得到了廣泛的應用。需要注意的是，雖然矩量法中求解阻抗矩陣的表達式較為簡單，但其計算工作量很大， 對于以積分方程為基礎的離散方程，其系數(shù)矩陣通常為滿矩陣，所有元 素都需大量的數(shù)值計算。尤其隨著目標電尺寸的增大， 矩量法得到的系 數(shù)矩陣將迅速增大，這將給計算機內(nèi)存和 CPU帶來沉重的負擔。為了克 服這些困難，人們對傳統(tǒng)矩量法進行了一些改進，提出了一些快速、有 效的方法，如快速多極子方法（FMM）和多層

6、快速多極子方法（MLFMM）, 阻抗矩陣局部化（IML）方法和壓縮或稀疏化阻抗矩陣的小波分解法，（基于快速Fourier變換的CG-FFT法、稀疏矩陣規(guī)則網(wǎng)格（SMCG法和自適 應積分法（AIM）,來降低計算機內(nèi)存和計算量的需求。在這些快速分析 方法中，需要的計算量和內(nèi)存分別降為 D（NlogW）和0（N）, N為未知變量 數(shù)。但是，這些改進的方法仍然受傳統(tǒng)矩量法離散尺寸的限制。采用整 域基函數(shù)代替分域基函數(shù)可以降低矩量法系數(shù)矩陣的維數(shù)。然而，在絕大多數(shù)情況下，難以找到合適的整域基函數(shù)。為此，近年來，人們又相 繼開展了一些基于部分域（子域或塊）概念來降低矩陣維數(shù)的研究，如多 層矩量法（MMM）

7、、子域多層法（SMA卜合成基函數(shù)（SBF法等。這些方法 通過對問題的部分域進行分析來構(gòu)造宏基函數(shù)，宏基函數(shù)的域比傳統(tǒng)矩量法的大，因此可以降低未知變量數(shù)。這幾種方法是通過迭代的方式遞歸地修正互耦項來改進解的收斂性。 特征基函數(shù)法（CBFM尾近兩年提出 來的一種新的求解電磁散射問題的有效方法。二矩量法的基本原理矩量法是將算子方程轉(zhuǎn)化為矩陣方程，然后求解該矩陣方程的方法，由于在求解方程過程中，需要計算廣義矩量，所以我們稱這種方法稱為矩量法，其實質(zhì)是內(nèi)域基加權余量法。(一)離散化過程已知算子方程：Lf = gL程g已知，f唯一確定。設f = &fn (然是系數(shù)，fn是基函數(shù))。其 n中fj不一定是完備

8、的，n越大，不一定越好，只有fn完備，n越大越接近真值。方程離散為：- n Lf n = gn此過程的主要目的在于將算子方程轉(zhuǎn)為代數(shù)方程，具體步驟如下：(1)在算子L的定義域內(nèi)適當?shù)倪x擇一組基函數(shù)，fi,f2,f3.fn,他們應該是線性無關的。(2)將未知函數(shù)f(x)表示為該組基函數(shù)的線性組合，并取有限項近似，N即：f(x)= :nfn Tn 又= :nfn(2 1)n z1n 1(3)將公式(2-1)利用算子的線性，可將算子方程化為代數(shù)方程，即 ：NnL(fn)=g(2-2)n 1于是，待求f(x)的問題化解為求解fn的系數(shù)4的問題。(二)取樣檢驗過程為了使f(x)的近似函數(shù)fN(X)之間的

9、誤差極小，必須進行抽樣檢驗,在抽樣點上使加權平均誤差為零，從而確定未知系數(shù)可，這一過程的基本步驟如下：(1)在算子L的值域內(nèi)取適當?shù)倪x擇一組加權函數(shù)WI,W2, W3.Wn ,他們也是線性無關的。(2)將Wm與式(2-2)取內(nèi)積進行抽樣檢驗，因為要確定N個未知數(shù)，需要進行N次抽樣檢驗，則，ncL(fn),wm=(g,wm)(m =1,2,3. N)(2-3)(3)利用算子的線性和內(nèi)積的性質(zhì)，將(2-3)式轉(zhuǎn)化為矩陣方程，即：夕 n 比(fn),Wm)=Wm)(m=1,2,3N)(選配過程)(24)將它寫成矩陣形式：1-lmn 11 -n 1 = 0 1(2-5)于是，求解代數(shù)方程問題轉(zhuǎn)化為求解

10、矩陣方程的問題。(三)矩陣求逆一旦求得了矩陣方程，通過常規(guī)的矩陣求逆或求解線性方程組，就可以得到矩陣方程的解。如果 Ln】是非奇異的，則它的逆矩陣hmnL是存在的，則矩陣方程的解為：除二lmn I 6m 】(26)將求得的展開系數(shù) 4帶入到(2-1)式中，就可以得到原來算子方程的 近似解為：NfN X Cnfn(x)(2-7)n 1以上所述是矩量法求解算子方程的基本過程，在矩量法的所有應用 中，通常都要遵循這個統(tǒng)一的過程。由于 f(x)和fN(x)之間存在近似， 故算子方程的左端近似值與右端精確值之間存在如下關系： N鞏x)= ncL(fn(x) -g(x)(2-8)n 1三基函數(shù)與權函數(shù)的選

11、擇矩量法的求解過程是簡單的， 求解步驟是統(tǒng)一的，應用起來比較方 便。求解精度取決于很多因素，例如離散化的成都、基函數(shù)和權函數(shù)的 選取和矩陣的求解過程等，其中基函數(shù)與權函數(shù)的選擇尤其重要，可以直接影響到求解的精度?；瘮?shù)可以分為全域基和分域基，權函數(shù)可以分為權函數(shù)、分域權和點匹配，它們之間的不同組合便形成不同的方法。下面簡單的對此進 行介紹。1、 全域基函數(shù)：全域基函數(shù)是指算子 L定義域內(nèi)的全域上存在一組基函數(shù)。他們應該滿足邊界條件。在矩量法求解的離散 化過程中，選擇全域基函數(shù)作為展開函數(shù)，實質(zhì)上將未知函 數(shù)表示為全域存在的若干個離散化基函數(shù)的線性組合。2、 分域基函數(shù)：分域基函數(shù)不是在算子 L

12、定義域的全域上存在的，而僅僅存在于算子定義域的各個分域上的函數(shù)。選擇分 域基函數(shù)作為未知函數(shù)的展開函數(shù)，在矩量法求解的離散化 過程中是一種區(qū)域離散，即未知函數(shù)表示為各個分域上存在 的函數(shù)線性組合。3、 全域權：在算子L的值域內(nèi)選擇一組權函數(shù) W1,W2，W3.Wm,它是一組在L值域內(nèi)的全域上存在的權函數(shù)。在此方法中，如 果選擇Wn=fn,即權函數(shù)等于基函數(shù)，則稱此方法為伽略金 法。4、 點匹配:全域基、全域權的伽略金發(fā)是一種常用的求解方法。但是，如果算子本身是復雜的積分算子，或者選擇了比較復 雜的基函數(shù)，由于內(nèi)積運算本身又是積分運算，從而使矩陣 元素lmn和gm的形成十分困難。這時為了簡化，我

13、們可以利用 5函數(shù)的選擇性將6函數(shù)作為權函數(shù)將使內(nèi)積計算得以簡 化。四算例問題：計算帶電導體板的電容。正方形板，邊長為 2a,位于z = 0。 空間中任意一點電位 (x,y,z)=/dxr dy,其中 *4二；0RR = (x-x)2 +(y-y)2 +(z-z)2。在板上 二V 為常數(shù)，有 V =dxdy-,其中電荷分布 仃(x,y)為待求，。4 二；0RR=J(x-x)2 +(y-y)2 ,電容為：C=9 = l/dx fdyQ(x,y計算當a =4時方板的歸一化電容V V 2a1 .選取分域基：方板各邊N等分，2b =紅，其中M =N2個子域。定義: N1 n小2fn = ,其中在嶼上為

14、1,否則為0。將。(x,y2 fn帶入0n 1a , a ，二 X , yV = dx dy -中，有 一 4二;0Rn =N2Vlmnfnn 1其中 lmn = ：dx ：dy _aa有而Wr/ lmn是g上單位振幅的均勻電荷密度在的中心處產(chǎn)生的電位。因此有1 n zN2CnS =、lmnSV n 4m,n物體的電容式其各小塊電容的總和加上每一對小塊間的互電容。, a , a ,二 X , yL f =L 二 x ,y = dx dy 2, 2 =V4 二；。；(x-x ) (y-y )2.取Wm =8(x-xm)5(y-ym)lmn N2 N2 L n N2 1 = V L2Ln h21

15、= lmn V1 .C = V LSnn = ：Sn 1 M |lJmn M M 1 L 13. lmn的計算：lnn：由本身面上的單位電荷密度在其中心處產(chǎn)生的電位(自作用單元)b b12b-lnn = dx dy-22 =ln 124 二 0 ；x y，0而互作用單元：dxdy6nl 二 lmn 一一. .2、2、2 .、2.S1 4 0 ；(x -xm) (y _ ym)4 _0,(xnxm)，Hn - Ym)經(jīng)過計算模擬可以發(fā)現(xiàn)當N取值越大時，結(jié)果將趨近于一個定值。當N為20時，計算得到歸一化電容為40.29PF非常接近理論計算值40PF,誤差為0.725%,在誤差允許范圍內(nèi)。程序：#p

16、ragma onceclass Martix(private:const int M;const int N;bool flag;double *array;public:Martix(double *_array, int _M, int _N) :M(_M), N(_N)(-array = _array;flag = false;Martix(int _M, int _N) :M(_M), N(_N)(flag = true;array = new double *M;for (int i = 0; i M; i+)arrayi = new doubleN;for (int i = 0;

17、i M; i+)(for (int j = 0; j N; j+)(arrayij = 0;M artix Gauss();void show();void swap(int i, int j);int Max(int i);M artix operator+(Martix &T);M artix operator-(Martix &T);M artix operator*(Martix &T);M artix operator*(double t);friend Martix operator*(double t, const Martix &T)Martix sum(T.M, T.N);f

18、or (int i = 0; i T.M; i+) (for (int j = 0; j array;) int GetM()(return this-M;)int GetN()(return this-N;)bool GetFlag()(return this-flag;)Martix(););#include Martix.h#include#includeusing namespace std;void Martix:show()(for (int i = 0; i M; i+)(cout fixed right;for (int j = 0; j N; j+) (cout setw(1

19、0) setprecision(6) arrayij )cout endl;Martix Martix:Gauss() (Martix sum(M, N);for (int i = 0; i N; i+)sum.arrayii = 1;for (int i = 0; i N; i+)if (abs(arrayii) = 0.00001&iM - 1)(if (abs(arrayMax(i)i) = 0.00001)continue;int count = Max(i);swap(i, count);sum.swap(i, count);)else if (abs(arrayii) 0.0000

20、1&i = N - 1)(cout 沒有逆矩陣 endl;exit(0);)/*cout 變換前show();*/double temp = arrayii;for (int j = 0; j N; j+)(arrayij = arrayij / temp;sum.array皿=sum.arrayij / temp;)for (int k = i + 1; k = 0.00001)(for (int j = 0; j N; j+)(arraykj = arraykj - temp*arrayij;sum.arraykj = sum.arraykj - temp*sum.arrayij;)dou

21、ble t = abs(array00);for (int i = 1; i arrayii);if (abs(t) = 0.0001) (cout 沒有逆矩陣= 0; i-)(for (int j = 0; j N; j+)(for (int k = i + 1; k N; k+)(sum.arrayij = sum.arrayij - arrayik * sum.arraykj; )return sum;)void Martix:swap(int i, int j)(double temp = 0;for (int k = 0; k M, this-N);for (int i = 0; i

22、 M; i+)(for (int j = 0; j arrayij;) return sum;)Martix Martix:operator+(Martix &T)(Martix sum(T.M, T.N);if (M != T.M | N != T.N) (cout 錯誤 endl;)for (int i = 0; i M; i+)(for (int j = 0; j arrayij + T.arrayij;)return sum;)Martix Martix:operator-(Martix &T)(M artix sum(T.M, T.N);if (M != T.M | N != T.N

23、)(cout 錯誤 endl;)for (int i = 0; i M; i+)(for (int j = 0; j arrayij - T.arrayij;)return sum;)Martix Martix:operator*(Martix &T)(M artix sum(T.M, T.N);if (N != T.M)(cout 錯誤 endl;)for (int i = 0; i M; i+)(for (int j = 0; j T.N; j+)(for (int k = 0; k arrayik * T.arraykj;return sum;)Martix Martix:operato

24、r*(double t)(Martix sum(this-M, this-N);for (int i = 0; i M; i+)(for (int j = 0; j N; j+)( sum.arrayij = this-arrayij * t;) return sum;)int Martix:Max(int i)(int count = i;double temp = abs(arrayii);for (int k = i + 1; k N; k+) (if (temp abs(arrayki)temp = arrayki;)while (temp != abs(arraycounti)(count+;) return count;)Martix:-Martix()()#include#include#include#includeMartix.h using namespace std;const double pi = 3.1415926;const double eps = 1 / (36 * pi)*pow(10, -9);double a = 2.0;int