丁永平，王添文，王洪亮

（1.中國(guó)兵器工業(yè)新技術(shù)推廣研究所，北京 100089；2.中國(guó)船舶重工集團(tuán)公司第703研究所，哈爾濱 150036）

隨著電子信息科學(xué)的飛速發(fā)展，電子設(shè)備的結(jié)構(gòu)越來(lái)越復(fù)雜，所用器件的集成度、功率等級(jí)越來(lái)越高。這些變化給系統(tǒng)內(nèi)部以及系統(tǒng)間的電磁兼容性（Electromagnetic Compatible，EMC）帶來(lái)了許多問(wèn)題[1]。為了解決這一系列問(wèn)題，需要在系統(tǒng)設(shè)計(jì)的初期考慮電磁兼容方面的需求。

電磁兼容領(lǐng)域的數(shù)值分析方法可以幫助設(shè)計(jì)者建立起系統(tǒng)、設(shè)備或器件的數(shù)學(xué)模型，并通過(guò)仿真的方法，在設(shè)計(jì)初期就對(duì)電磁兼容性進(jìn)行規(guī)劃和設(shè)計(jì)。其實(shí)質(zhì)是針對(duì)現(xiàn)實(shí)中的電磁兼容問(wèn)題，運(yùn)用某種或幾種數(shù)學(xué)方法，在麥克斯韋方程的基礎(chǔ)上，建立該實(shí)際問(wèn)題的數(shù)學(xué)模型，并通過(guò)對(duì)模型解析解或數(shù)值解的求解來(lái)解決實(shí)際問(wèn)題的一門(mén)交叉學(xué)科。

1 電磁場(chǎng)數(shù)值方法可以解決的問(wèn)題

電磁兼容領(lǐng)域內(nèi)的數(shù)值分析方法可以解決的問(wèn)題包括以下幾類(lèi)。

1）基礎(chǔ)問(wèn)題的研究，主要包括輻射場(chǎng)的分布、傳播特性、材料的屏蔽效能以及場(chǎng)線(xiàn)耦合問(wèn)題等。

2）PCB 板級(jí)問(wèn)題的研究，主要包括信號(hào)完整性及電磁輻射問(wèn)題。

3）系統(tǒng)級(jí)電磁兼容問(wèn)題，主要是針對(duì)如汽車(chē)、飛機(jī)等大型系統(tǒng)進(jìn)行整體建模，并對(duì)其綜合電磁兼容特性進(jìn)行分析、改進(jìn)。

4）標(biāo)準(zhǔn)測(cè)試的建模及仿真，是利用數(shù)值分析的方法對(duì)各種電磁兼容標(biāo)準(zhǔn)測(cè)試過(guò)程進(jìn)行建模，在計(jì)算機(jī)內(nèi)對(duì)系統(tǒng)或設(shè)備進(jìn)行預(yù)測(cè)試仿真，其優(yōu)點(diǎn)是節(jié)約設(shè)計(jì)成本，縮短設(shè)計(jì)時(shí)間。

5）電磁生態(tài)效應(yīng)的研究，主要研究生物在電磁環(huán)境內(nèi)所受的影響，目前主要研究的熱點(diǎn)是手提電話(huà)和電磁爐對(duì)人體的影響。

2 常用的數(shù)值分析方法

2.1 矩量法（MoM法）

矩量法是由R.F.Harrington 于1968 年在“Field Computation by Moment Methods”中提出的[2]。它所面對(duì)的問(wèn)題主要有：天線(xiàn)設(shè)計(jì)、微波網(wǎng)絡(luò)、生物電磁學(xué)、PCB 板上的輻射效應(yīng)及微帶線(xiàn)研究等問(wèn)題。其核心思想是由自由空間格林函數(shù)的積分形式得到任意源激勵(lì)，產(chǎn)生表面電流分布，進(jìn)而得到電磁場(chǎng)的分布特性。

對(duì)于某個(gè)系統(tǒng)可以列出：

式中：L 為線(xiàn)性算子；g 為源函數(shù)；I 為未知的電流分布表達(dá)式。

再將整個(gè)系統(tǒng)根據(jù)一定的規(guī)則劃分為N個(gè)小區(qū)域，則其表面電流分布的表達(dá)式就可以表示為：

式中：αn為待定參數(shù)。

同樣地，將源函數(shù)g劃分為N個(gè)小區(qū)域（微元）。

式中：βm為待定參數(shù)；Tm為小區(qū)域上的核函數(shù)。對(duì)于式（2）中的每一個(gè)小區(qū)域電流分布表達(dá)式，又可以看成是由被劃分后的源函數(shù)的每一個(gè)微元作用疊加的效果。即：

式中：lmn為待定系數(shù)。將式（2），（3），（4）整理得到：

這樣就將求解電磁場(chǎng)激勵(lì)下的電流分布問(wèn)題轉(zhuǎn)換為求解矩陣逆的問(wèn)題，在得到電流分布后就可以很容易地算出電磁場(chǎng)分布。

矩量法在解決無(wú)邊界輻射問(wèn)題、“細(xì)線(xiàn)”問(wèn)題和均勻電介質(zhì)問(wèn)題時(shí)表現(xiàn)得尤為突出。近年來(lái)矩量法還被應(yīng)用于解決復(fù)雜的電大尺寸問(wèn)題中。MLFMA（the multilevel fast multipole algorithm），ACA（the adaptive cross approximation algorithm），SVD（singular value decomposition technique）等算法也應(yīng)運(yùn)而生，降低了對(duì)CPU運(yùn)算能力和內(nèi)存的需求。矩量法還可以與解析方法(analytical methods)、傳輸線(xiàn)矩陣方法相結(jié)合解決屏蔽效能問(wèn)題和線(xiàn)纜輻射與受擾問(wèn)題。另外矩量法還可以與有限元法相結(jié)合解決非均勻材料結(jié)構(gòu)的感應(yīng)電流問(wèn)題。

2.2 時(shí)域有限差分法（FDTD法）

時(shí)域有限差分法是由Yee KS 于1966 年在“Numerical solution of initial boundary value problem involving Maxwell equation in isotropic media”中提出的[3]。它所面對(duì)的問(wèn)題主要有：傳輸線(xiàn)波導(dǎo)的傳輸、天線(xiàn)的接收檢測(cè)和輻射、耦合屏蔽和透入效應(yīng)、散射和逆散射、開(kāi)關(guān)過(guò)渡過(guò)程等。這種方法最大的特點(diǎn)是適合分析瞬態(tài)響應(yīng)問(wèn)題，特別是具有復(fù)雜幾何形狀和復(fù)雜環(huán)境的情形，但缺點(diǎn)是不適合分析低頻問(wèn)題[4]。其核心思想是建立在微分方程基礎(chǔ)上的。如圖1 所示，首先將空間按立方體分割，電磁場(chǎng)的6個(gè)分量在空間的取樣點(diǎn)分別放在立方體的邊沿和表面中心點(diǎn)上，電場(chǎng)與磁場(chǎng)分量在任何方向始終相差半個(gè)網(wǎng)格步長(zhǎng)。在時(shí)間上，Yee 也把電場(chǎng)分量與磁場(chǎng)分量差半個(gè)步長(zhǎng)取樣。使得利用一階導(dǎo)數(shù)的二階中心差分近似從Maxwell 方程獲得的FDTD 公式。

圖1 電場(chǎng)和磁場(chǎng)在空間和時(shí)間上的離散Fig. 1 Interleaving of E and H fields in space and time in the FDTD formulation

Maxwell旋度方程為：

式中：H，E分別為磁場(chǎng)和電場(chǎng)分量；D為電通密度；B為磁通密度；J為體電流密度。直角坐標(biāo)系中，式（7）可以寫(xiě)為：

式中：Hx，Hy，Hz，Ex，Ey，Ez分別為對(duì)應(yīng)方向上的磁場(chǎng)和電場(chǎng)分量；σ，μ分別為介質(zhì)電導(dǎo)率和磁導(dǎo)率；ε為介電常數(shù)。

令f（x，y，z，t）代表E和H在直角坐標(biāo)中的某一分量，i，j，k，n分別表示x軸、y軸、z軸以及時(shí)間軸的單位方向向量。在時(shí)間和空間中的離散取以下符號(hào)表示：

對(duì)于f（x，y，z，t）關(guān)于時(shí)間、空間的一階偏導(dǎo)數(shù)取中心差分近似，即

帶入式（8）和（9）中，經(jīng)推導(dǎo)整理得到

式中：CA（m），CB（m）均為相應(yīng)的常數(shù)項(xiàng)[4]；Hx，Hy，Hz，Ex，Ey，Ez分別為對(duì)應(yīng)方向上的磁場(chǎng)和電場(chǎng)分量。同樣地，可以推導(dǎo)出的表達(dá)式。如果知道了初始條件和邊界條件，就可以由遞推公式得到空間任意復(fù)雜結(jié)構(gòu)的電磁場(chǎng)分布。

時(shí)域有限差分法未來(lái)發(fā)展方向主要有：時(shí)域有限差分法的穩(wěn)定性、精確性、誤差、網(wǎng)格效應(yīng)；進(jìn)一步發(fā)展邊界吸收條件完全的匹配層（PML）；與熱仿真、電路仿真、半導(dǎo)體器件仿真結(jié)合解決更多問(wèn)題；與矩量法、有限元法結(jié)合的混合算法等。

2.3 有限元法（FEM法）

有限元法是1969 年由Silveter 在“Finite element solution of homogeneous waveguide problems”將其推廣應(yīng)用于時(shí)諧電磁場(chǎng)問(wèn)題的[5]。它所面對(duì)的問(wèn)題主要有：波導(dǎo)設(shè)計(jì)、微型芯片設(shè)計(jì)、半導(dǎo)體器件設(shè)計(jì)、屏蔽效能分析、生物體對(duì)電磁能量的吸收等。其核心思想是將一個(gè)復(fù)雜連續(xù)介質(zhì)的求解區(qū)域分解為有限個(gè)形狀簡(jiǎn)單的子區(qū)域，作為原區(qū)域的等效域，從而把求解連續(xù)體的場(chǎng)變量問(wèn)題化簡(jiǎn)為求解有限個(gè)單元節(jié)點(diǎn)上的場(chǎng)變量值問(wèn)題。

選擇空間四面體為單位結(jié)構(gòu)，如圖2所示。在此單元內(nèi)，φe表示任意未知函數(shù)，它能夠近似為：

式中：ae，be，ce，de均為常數(shù)，e為未知變量。

可以證明，差值函數(shù)具有如下性質(zhì)：

圖2 有限元法四面體單位結(jié)構(gòu)Fig.2 Tetrahedron unit structure used in FEM

在利用變分原理和離散化方法建立了有限元矩陣方程后，就面臨著以結(jié)點(diǎn)值為未知數(shù)的矩陣方程的求解。將方程寫(xiě)為：

式中：A，b 均為由已知向量組成的矩陣；x 為待求的分布函數(shù)矩陣。當(dāng)上式右端的已知激勵(lì)矩陣b不為零時(shí)，為了確定方程的解，一般利用各種等效方法對(duì)矩陣A 求逆，其中最適用于有限元方法矩陣的是分解法。

有限元法未來(lái)發(fā)展方向主要有：通過(guò)設(shè)置時(shí)域終止邊界或與矩量法相結(jié)合將有限元法應(yīng)用于對(duì)輻射問(wèn)題和散射問(wèn)題的研究當(dāng)中。

2.4 傳輸線(xiàn)矩陣法（MTL法）

傳輸線(xiàn)矩陣法是由P.B.Johns 和R.L.Beurle 于1968 年在“Numerical solutions of 2-dimensional scattering problems using a transmission-line matrix”中提出的[6]。它所面對(duì)的問(wèn)題主要有：波導(dǎo)不連續(xù)性、散射問(wèn)題、信號(hào)完整性設(shè)計(jì)、電纜特性問(wèn)題分析、天線(xiàn)輻射原理。其核心思想是建立在電磁波傳遞特性與傳輸線(xiàn)中電壓、電流傳播特性相似性的基礎(chǔ)上發(fā)展而來(lái)的，現(xiàn)在已經(jīng)被推廣到了空間問(wèn)題中。

當(dāng)僅考慮一維傳輸線(xiàn)時(shí)，集總參數(shù)模型的等效電路如圖3所示。根據(jù)基爾霍夫電壓、電流定律，可以得到：

式中：V，I分別為電路中的電壓和電流；z為位置坐標(biāo)；t為時(shí)間坐標(biāo)；l，c分別為分布的電感和電容參數(shù)。經(jīng)過(guò)整理、去耦，可以得到：

公式（19）即為傳輸線(xiàn)法中常用的電報(bào)方程。

圖3 傳輸線(xiàn)集總參數(shù)模型Fig.3 Lumped parameters model of transmission line

如圖4所示，將這一思想推廣到三維空間中，空間的每個(gè)節(jié)點(diǎn)由6個(gè)方向的小傳輸線(xiàn)段組成，每個(gè)方向上的傳輸線(xiàn)用相互垂直的2個(gè)電壓分量來(lái)描述[7]。在每個(gè)計(jì)算周期內(nèi)有12個(gè)電壓變化，產(chǎn)生12個(gè)反射電壓。而這一組反射電壓又稱(chēng)為臨近網(wǎng)格的激勵(lì)電壓，周而復(fù)始地重復(fù)這一過(guò)程，就可以模擬電磁波在空間中傳播的過(guò)程。在建模時(shí)需要注意，每個(gè)網(wǎng)格劃分最大尺寸是所關(guān)心最高頻率波長(zhǎng)的1/10。其過(guò)程可以用公式描述：

式中Vi是激勵(lì)電壓矩陣；Vr是反射電壓矩陣；S是散射矩陣；C是連接矩陣。

傳輸線(xiàn)矩陣法發(fā)展方向主要集中在模型的改進(jìn)方面，許多人提出了包含信息更多的SCN 模型。為降低對(duì)內(nèi)存的需求提出的ATLM，ARTLM 是2 種最新的方法。

圖4 對(duì)稱(chēng)濃縮節(jié)點(diǎn)Fig.4 Symmetrical condensed node（SCN）

2.5 局部元等效法（PEEC法）

局部元等效法是由A.Ruehli于1974年在“Equivalent circuit models for three-dimensional multiconductor systems”中提出的[8]。它所面對(duì)的問(wèn)題主要有：互聯(lián)機(jī)構(gòu)和封裝器件的電磁兼容分析、PCB 板級(jí)的各種結(jié)構(gòu)仿真問(wèn)題。其核心思想是將求解麥克斯韋問(wèn)題轉(zhuǎn)換成求解電路的支路電流和節(jié)點(diǎn)電位問(wèn)題。經(jīng)過(guò)等效以后，可以運(yùn)用電路的知識(shí)來(lái)分析電磁場(chǎng)問(wèn)題，一旦矩陣參數(shù)被確定，便可以通過(guò)求解等效的電流分布和電位分布，提供麥克斯韋方程的全波解決方案。

根據(jù)電磁學(xué)理論，電荷密度為ρ，電流密度為J的源在空間任意一點(diǎn)產(chǎn)生的電場(chǎng)強(qiáng)度可以表示為：

式中：σ為導(dǎo)體電導(dǎo)率；r為位置矢量；A和Φ分別為矢量磁位和標(biāo)量電位。

K個(gè)導(dǎo)體系統(tǒng)中：

式中：K（r，r′） 為積分核函數(shù)

將式（22），（23）代入式（21）中推導(dǎo)、整理。得到：

通過(guò)觀察可以看出式（24）中第1 項(xiàng)相當(dāng)于電阻項(xiàng)，第2 項(xiàng)相當(dāng)于電感項(xiàng)，第3 項(xiàng)相當(dāng)于常數(shù)項(xiàng)。則可以寫(xiě)成如下形式：

這樣就將電磁場(chǎng)的問(wèn)題等效成為電路的問(wèn)題。

3 驗(yàn)證數(shù)值分析模型有效性的方法

隨著各種數(shù)值分析方法在電磁兼容領(lǐng)域內(nèi)的應(yīng)用不斷增加，如何驗(yàn)證數(shù)值分析模型的有效性就成為一個(gè)值得關(guān)心的問(wèn)題[9]。針對(duì)這一問(wèn)題，提出幾種驗(yàn)證數(shù)值分析模型有效性的方法。

3.1 通過(guò)實(shí)驗(yàn)方法驗(yàn)證

實(shí)驗(yàn)是最常用的一種驗(yàn)證方法，采用這種方法時(shí)，首先要注意實(shí)驗(yàn)的各種參數(shù)設(shè)置應(yīng)與仿真保持一致。其次，需要在仿真過(guò)程中考慮實(shí)際測(cè)量時(shí)的各種限制，比如仿真時(shí)觀測(cè)點(diǎn)探頭的輸入阻抗應(yīng)與測(cè)試設(shè)備相同。當(dāng)用到天線(xiàn)時(shí)，仿真用的是一個(gè)近乎完美的天線(xiàn)，而測(cè)量用到的實(shí)際天線(xiàn)的各種參數(shù)會(huì)隨著頻率和方位的變化而變化。最后，應(yīng)注意實(shí)驗(yàn)本身的準(zhǔn)確性，在普通環(huán)境下做的測(cè)試，其準(zhǔn)確性通常不夠。即使是在一些商業(yè)的電磁兼容專(zhuān)業(yè)實(shí)驗(yàn)室中測(cè)試得到的結(jié)果通常也會(huì)有±6 dB的誤差。產(chǎn)生誤差的因素很多，主要是測(cè)量設(shè)備自身的誤差、測(cè)試電纜敷設(shè)的方式、天線(xiàn)因數(shù)的變化、實(shí)驗(yàn)室墻面反射系數(shù)的變化等。

在某些情況下，必須考慮到測(cè)試環(huán)境的精度問(wèn)題。當(dāng)測(cè)試數(shù)據(jù)與仿真結(jié)論不符時(shí)，應(yīng)該考慮到是否為測(cè)量結(jié)果不準(zhǔn)確，不應(yīng)該完全相信測(cè)量結(jié)果而判定仿真結(jié)果錯(cuò)誤。

3.2 通過(guò)多種仿真技術(shù)驗(yàn)證

另一種流行的驗(yàn)證方式就是通過(guò)2 種或2 種以上不同的數(shù)值分析方法來(lái)對(duì)同一個(gè)電磁兼容問(wèn)題進(jìn)行建模[10]。同樣需要注意，對(duì)于每種建模方法要保證各種參數(shù)設(shè)置是相同的，這樣的比較才是有效的。前面提到的數(shù)值分析的方法在采用的時(shí)候，要盡量選取機(jī)理差異較大的兩種進(jìn)行對(duì)比，通常FDTD，F(xiàn)EM，TLM 方法被認(rèn)為是與體積相關(guān)的分析方法，而MoM，PEEC 法被認(rèn)為是與表面相關(guān)的分析方法。在驗(yàn)證時(shí)，最好從每一類(lèi)中都選取一個(gè)來(lái)進(jìn)行對(duì)比。這樣的選取方法要比使用多種基于相同機(jī)理的仿真軟件來(lái)驗(yàn)證的方法得到的結(jié)論更為可信。

3.3 通過(guò)中間結(jié)論驗(yàn)證

數(shù)值分析方法提供了比較直觀的中間結(jié)論，比如空間電場(chǎng)、磁場(chǎng)的分布情況，導(dǎo)體表面電流的分布情況等等。這些結(jié)論不能通過(guò)實(shí)驗(yàn)的方法觀察到，但可以很好地幫助驗(yàn)證模型的有效性。在MoM 法和PEEC 法中的穩(wěn)態(tài)電流分布就是一個(gè)非常重要的中間結(jié)論。其中一個(gè)結(jié)論是，在兩個(gè)連續(xù)微元內(nèi)電流是否有突變，在導(dǎo)線(xiàn)的末端電流是否為零，電流的流向是否正確等等，這些都可以幫助判斷模型的正確性。另一個(gè)重要的中間結(jié)論是，在用FDTD 法，PEEC 法，TLM 法時(shí)電磁場(chǎng)的特性，是否在通過(guò)屏蔽層時(shí)明顯減弱，是否與理論分析的分布情況一致，也都是判定模型有效性的有力依據(jù)。

3.4 通過(guò)標(biāo)準(zhǔn)問(wèn)題和定量問(wèn)題驗(yàn)證

所謂標(biāo)準(zhǔn)問(wèn)題和定量問(wèn)題是指那些已經(jīng)被實(shí)驗(yàn)和仿真論證過(guò)，具有定量結(jié)論的標(biāo)準(zhǔn)化問(wèn)題，比如一些半波偶極子的輻射問(wèn)題等。在采用一種新的數(shù)值分析方法時(shí)，可以先針對(duì)這些已經(jīng)有結(jié)論的標(biāo)準(zhǔn)問(wèn)題、定量問(wèn)題進(jìn)行建模和仿真，用仿真得到的數(shù)值或波形與已有的結(jié)論進(jìn)行對(duì)比，從而也可以驗(yàn)證這種新型數(shù)值分析方法的有效性。

3.5 通過(guò)收斂性和參數(shù)變化特性驗(yàn)證

收斂性是指采用FDTD法、MoM法、FEM法時(shí)需要對(duì)分析域進(jìn)行網(wǎng)格劃分，網(wǎng)格的大小一般是所關(guān)心最高頻率波長(zhǎng)的1/10。但這種劃分方法往往不夠，需要進(jìn)一步縮小網(wǎng)格的尺寸，如果得到的結(jié)論出現(xiàn)錯(cuò)誤或不收斂，則說(shuō)明模型是無(wú)效的。如果不具備更強(qiáng)的計(jì)算能力，則應(yīng)考慮采取其它驗(yàn)證方法。

在建模過(guò)程中往往包含一些比較重要的參數(shù)，通過(guò)這些參數(shù)的改變，相應(yīng)的仿真結(jié)論會(huì)有所變化。這樣就可以改變這些參數(shù)，然后運(yùn)用理論分析和以往經(jīng)驗(yàn)對(duì)仿真結(jié)果的變化趨勢(shì)進(jìn)行判斷，如果結(jié)論與判斷不同則說(shuō)明仿真模型也是無(wú)效的。

4 結(jié)論

1）介紹了電磁兼容領(lǐng)域內(nèi)的數(shù)值分析方法可以解決的問(wèn)題有基礎(chǔ)理論的研究、PCB 板級(jí)問(wèn)題的研究、系統(tǒng)級(jí)電磁兼容的研究、標(biāo)準(zhǔn)測(cè)試的建模及仿真以及電磁生態(tài)效應(yīng)的研究。

2）介紹了包括矩量法、時(shí)域有限差分法、有限元法、傳輸線(xiàn)矩陣法、局部元等效法等幾種常用的數(shù)值分析方法，分析了各種方法的適用范圍和核心思想。

3）對(duì)如何驗(yàn)證數(shù)值分析方法有效性進(jìn)行了分析，提出了通過(guò)實(shí)驗(yàn)、多種仿真技術(shù)、中間結(jié)論、標(biāo)準(zhǔn)問(wèn)題和定量問(wèn)題以及收斂性和參數(shù)變化特性等5種靈活、有效的驗(yàn)證方法。

[1]李明，朱中文，蔡偉勇. 電磁兼容技術(shù)研究現(xiàn)狀與趨勢(shì)[J].認(rèn)證與標(biāo)準(zhǔn)電磁兼容，2007，8（7）：61—64.

[2]HARRINGTON R. Field Computation by Moment Methods[M]. New York：Macmillan publishing company，1968：20—38.

[3]YEE K. Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media[J]. IEEE Transaction on Antennas Propagation，1966，14（5）：302—307.

[4]葛德彪，閆玉波.電磁波時(shí)域有限差分法[M].西安：西安電子科技大學(xué)出版社，2005：13—20.

[5]SILVESTER P. Finite Element Solution of Homogeneous Waveguide Problems [J]. Alta Frequenza，1969，5（5）：6—16.

[6]JOHNS P，BEURLE R.Numerical Solutions of 2-dimensional Scattering Problems Using a Transmission-line Matrix[J].IEEE Procesure，1971，118（5）：1203—1208.

[7]HEINZ D，CHRISTIAN S，HERMANN S.Numerical Electromagnetic Field Analysis for EMC Problems[J].IEEE Transactions on Electromagnetic Compatibility，2007，49（2）：253-262.

[8]RUEHLI A.Equivalent Circuit Models for Three-dimensional Multiconductor Systems[J].IEEE Transaction on Microwave Theory Technology，1974，22（3）：216—221.

[9]IEEE Project 1597.1，Standard for Validation of Computational Electromagnetics（CEM）Computer Modeling and Simulation[S].

[10]DUFFY A，COLEBY D，WOOLFSON M.Progress in Quantifying Validation Data[J].IEEE Int Symp on EMC，2003，18（22）：323—328.