時(shí)域有限差分法分析混響室中場的均勻性

黃 華 牛中奇 白 冰

(西安電子科技大學(xué)電子工程學(xué)院,陜西 西安710071)

1.引 言

從眾多研究者的工作可知,混響室可采用多種數(shù)值方法進(jìn)行仿真分析：矩量法(MOM)[1];時(shí)域有限差分法(FDTD)[2-6];有限元法(FEM)[7-8];FDTD與MOM混合法[9];FEM與FDTD混合法[10];平面波積分表示法[11];傳輸線矩陣法(TLM)[12]等。在經(jīng)典的仿真分析中,建立的混響室模型一般由混響室的內(nèi)部空間、攪拌器及天線三部分組成[13]。圖1所示即為混響室測試示意圖。

圖1 混響室測試示意圖

通常,仿真計(jì)算中耗費(fèi)的時(shí)間與混響室的體積、六個(gè)壁面所用導(dǎo)體的電導(dǎo)率、電磁波的頻率范圍以及受試設(shè)備自身的特性有關(guān),然而這會(huì)使得計(jì)算量很大,因而十分耗時(shí),有時(shí)要用微機(jī)實(shí)現(xiàn)仿真幾乎是不可能的。為克服這一困難,Franco Moglie把重疊平面波理論[14]用到了該問題的處理之中,其基本思路是,用重疊平面電磁波替代經(jīng)典仿真分析中當(dāng)攪拌器轉(zhuǎn)動(dòng)時(shí)產(chǎn)生的向不同方向傳播的電磁波,用在計(jì)算程序中設(shè)置的對(duì)源于不同方位的重疊波的重復(fù)計(jì)算次數(shù)等效經(jīng)典仿真分析方法中的攪拌器的離散步進(jìn)次數(shù),而在相應(yīng)的模型中只保留受試對(duì)象,如圖2所示。這樣便可模擬混響室內(nèi)真實(shí)的電磁環(huán)境。然而,在該文獻(xiàn)中并沒有分析重疊入射波數(shù)目與混響室內(nèi)場分布之間的相互關(guān)系,基于此,本文將對(duì)該問題進(jìn)行研究,尋求其中的規(guī)律性以便為工程實(shí)踐提供理論指導(dǎo)。

圖2 隨機(jī)重疊入射平面波

2.方法簡介

2.1 隨機(jī)平面電磁波的生成

要把重疊平面波理論用于混響室內(nèi)場的仿真分析,解決的關(guān)鍵問題是重疊平面波的生成。注意到攪拌器轉(zhuǎn)動(dòng)時(shí)所產(chǎn)生的電磁波的傳播方向是隨機(jī)變化的這一特征,本文用所處方位是隨機(jī)的一組等效源產(chǎn)生的所謂隨機(jī)平面電磁波對(duì)其進(jìn)行模擬。又由于本文將采用FDTD法對(duì)混響室內(nèi)的場進(jìn)行仿真計(jì)算,而在FDTD法中,入射波源是設(shè)置在總場區(qū)與散射場區(qū)的分界面即所謂的連接邊界面上的,因而,上述的一組等效源設(shè)置的具體位置應(yīng)在連接邊界上。同時(shí),構(gòu)成這一組源的各個(gè)子等效源的位置參數(shù)及各個(gè)子等效源輻射的電磁參數(shù)在FDTD程序中應(yīng)被設(shè)定,在每一暫態(tài)的迭代過程中只需調(diào)用即可。這些參數(shù)包括,每個(gè)子等效源在球坐標(biāo)系中的坐標(biāo)r、θ和φ以及由它們輻射出的電磁波的極化角α,即平面波的電場矢量與x軸的夾角,如圖3所示。圖中P點(diǎn)為一個(gè)子等效源所在點(diǎn),E、H和n分別為由該等效源所產(chǎn)生的電磁波的電場、磁場和傳播方向的單位矢。為了模擬攪拌器轉(zhuǎn)動(dòng)時(shí)所產(chǎn)生的電磁波的隨機(jī)性,設(shè)置在球坐標(biāo)系中的一組等效源的分布方位亦應(yīng)是隨機(jī)的,也就是說,一系列子波源設(shè)置點(diǎn)的方位角θ和φ是隨機(jī)的。但值得注意的是,θ和φ的隨機(jī)性要受到球坐標(biāo)中坐標(biāo)變量取值范圍的約束,同時(shí)在對(duì)θ和φ的隨機(jī)取值中還要防止等效源在某一區(qū)域內(nèi)的分布過于集中,以致帶來對(duì)混響室內(nèi)場隨機(jī)分布模擬的失實(shí),因而,也要設(shè)置約束條件來予以防止,鑒于以上考慮,設(shè)置以下約束條件：

圖3 等效源的方位及其產(chǎn)生的電磁波的參數(shù)

1)θ的隨機(jī)取值應(yīng)在[0,π]之間,以便滿足所用坐標(biāo)系的約定。在此前提下,具體的隨機(jī)取值過程是,在 Fortran中首先由 random_element()在[0,1]之間抽取隨機(jī)數(shù),然后再乘以π弧度即可得到一系列θ的隨機(jī)值。

2)φ的隨機(jī)取值應(yīng)在[0,2π]之間,以滿足所用坐標(biāo)系的約定。對(duì)φ的隨機(jī)取值中應(yīng)特別值得注意的問題是,有時(shí)隨機(jī)取值的結(jié)果有可能出現(xiàn)某一區(qū)域的取值樣本過于密集,這種情況最有可能發(fā)生在高緯度即θ很小或θ很大的區(qū)域,對(duì)于發(fā)生這種情況的原因可作這樣的解釋：當(dāng)θ很小或θ很大時(shí),同球面上與之對(duì)應(yīng)的緯度圈的周長2πrsinθ將很小,因而在取同樣數(shù)量的φ樣本個(gè)數(shù)的情況下,其樣本點(diǎn)的密度便遠(yuǎn)大于其它緯度圈上的密度,因而若不加約束地按此種情況隨機(jī)設(shè)置等效源來模擬混響室內(nèi)的場分布,其結(jié)果將嚴(yán)重失實(shí)。為此,給予φ的隨機(jī)取值附加了這樣的約束：令φ=ψ/(2×sinθ),對(duì) ψ隨機(jī)取值,當(dāng)ψ的某一隨機(jī)值使得與之相應(yīng)的φ值落在[0,2π]之內(nèi),則保留之,否則則棄之重取,直至所取φ的個(gè)數(shù)滿足預(yù)設(shè)的等效源的個(gè)數(shù)為止。

3)α是線極化波的極化方向的表述,定義是電場矢量與x軸的夾角,稱其為極化角,其值在[0,2π]范圍內(nèi)。設(shè)各等效源產(chǎn)生的電磁波的電場幅度均是1,即Em=1 V/m.

上述約束條件中,θ和ψ是兩個(gè)相互獨(dú)立的隨機(jī)取樣變量,其中約束條件2)是為了使所選取的一系列等效源分布點(diǎn)P不至于過分密地集中在球面的兩極附近,以便保證之后采用FDTD法時(shí),設(shè)置在六個(gè)連接邊界面上的波源分布是較為均勻的,這樣才可較好地模擬在攪拌器轉(zhuǎn)動(dòng)時(shí)所產(chǎn)生的電磁波的傳播方向是隨機(jī)的效果。

2.2 隨機(jī)平面電磁波的引入

在仿真分析中,一是要對(duì)電磁波的特性進(jìn)行設(shè)置以模擬實(shí)際空間中實(shí)際存在的電磁波的特征,二是要把電磁波引入到仿真分析區(qū)域內(nèi)的各點(diǎn),以模擬電磁波在實(shí)際空間內(nèi)的傳播。一般說來對(duì)電磁波進(jìn)行設(shè)置的方法有兩種：一是解析方法,二是一維FDTD方法[15]。由于由一維FDTD方法引入的入射波在總場區(qū)內(nèi)的分布較為均勻且在散射場區(qū)泄漏較小,故本文選用第二種方法。在下面的仿真分析中,由程序設(shè)置的是連續(xù)的正弦平面電磁波,為了使引入到各點(diǎn)的正弦波能夠快速地達(dá)到穩(wěn)定,在開始的前半個(gè)周期加入了一個(gè)升余弦開關(guān)函數(shù),而在后半個(gè)周期仍為正弦函數(shù),如式(1)所示,其中,Em為電場振幅,f為電磁波的頻率,T為波的周期。

將Maxwell方程組中兩個(gè)旋度方程差分離散,再用一維FDTD法可得由電場表述磁場的如下表達(dá)式。

式中：μ0為真空中的磁導(dǎo)率;δ為討論空間的網(wǎng)格步長;Δt為計(jì)算中的時(shí)間步長;n為計(jì)算迭代的時(shí)間步數(shù);k為空間網(wǎng)格的位置。

2.3 程序?qū)崿F(xiàn)

本文用Fortran語言進(jìn)行計(jì)算編程。下面對(duì)主要步驟予以簡介。

2.3.1 電磁波的引入

如2.2節(jié)所述,電磁波是采用一維FDTD法引入的,且對(duì)電磁波的種類及其表達(dá)式作了說明?，F(xiàn)在介紹如何用Fortran程序予以實(shí)現(xiàn)。首先,在連續(xù)的正弦波上等間隔地采樣10000個(gè),然后把這10000個(gè)樣本數(shù)據(jù)由一維FDTD法逐次推進(jìn)到討論空間的每個(gè)網(wǎng)格內(nèi),每一時(shí)刻各網(wǎng)格內(nèi)均占有10個(gè)采樣數(shù)據(jù)。這樣設(shè)置好后,再將由各個(gè)角度入射的波的傳播方向投影到直角坐標(biāo)系的x、y、z方向上。當(dāng)把如式(1)所示的、其參數(shù)Em、f和T為固定取值的、來至P1點(diǎn)等效源所產(chǎn)生的電磁波在混響室內(nèi)的場分布計(jì)算完以后,再對(duì)入射波參數(shù)進(jìn)行更新,即對(duì)P2點(diǎn)的等效源產(chǎn)生的電磁波在混響室內(nèi)的場分布進(jìn)行計(jì)算,直至把Pn點(diǎn)的等效源產(chǎn)生的場計(jì)算完畢。之后對(duì)場進(jìn)行矢量疊加可得總場。

2.3.2 歸一化設(shè)置

由于同一電磁波的電場與磁場的振幅的數(shù)值之比為波阻抗,而對(duì)于空氣而言波阻抗η=377Ω,這表明,電場與磁場的量值不在同一數(shù)量級(jí),為了在計(jì)算中獲得同樣的準(zhǔn)確度,需對(duì)這兩個(gè)量進(jìn)行歸一化處理[16],而當(dāng)計(jì)算程序結(jié)束進(jìn)行數(shù)據(jù)存儲(chǔ)時(shí)再將其還原。本文是對(duì)電場進(jìn)行歸一化處理的,即程序處理時(shí),令?E=E/η而H不變。另外,計(jì)算時(shí),取討論空間中沿各坐標(biāo)方向的步長相等,即Δx=Δy=Δz=δ,并令討論空間內(nèi)的電導(dǎo)率σ=0,這樣Maxwell方程離散式前的系數(shù)就只剩下了 Δt/(ε0δ)和Δt/(μ0δ)兩項(xiàng)。由于

故有

將式(4)代入式(2),有

若將上式中的因子[en(k+1)-en(k)]/η記為(即歸一化為),則式(5)可寫為

將其它的離散式進(jìn)行同樣的處理。再將歸一化電場?E=E/η代入式(1),有

上面是對(duì)電場歸一化處理的具體步驟,存儲(chǔ)數(shù)據(jù)時(shí)再用E=η?E對(duì)電場進(jìn)行還原。當(dāng)然,在方法上也可將磁場歸一化為?H=ηH而E不變,處理后可得與上述類似的表達(dá)式,只是在數(shù)據(jù)存儲(chǔ)時(shí)需對(duì)磁場進(jìn)行還原。

3.仿真計(jì)算及結(jié)果分析

作為有實(shí)際背景的事例,計(jì)算中將討論空間的網(wǎng)格數(shù)取為90×90×90,其中沿各坐標(biāo)方向,連接邊界與吸收邊界的距離為10個(gè)網(wǎng)格。取電磁波頻率為915 MHz;空間步長為0.025 m;時(shí)間步長為Δt=4.5×10-11s,Δt這樣的取值是為了滿足數(shù)值穩(wěn)定性的要求。對(duì)計(jì)算數(shù)據(jù)的最終處理是將由Fortran程序計(jì)算出的數(shù)據(jù)導(dǎo)入Matlab程序中完成的。

3.1 對(duì)平面電磁波正弦特征的檢驗(yàn)

2.2節(jié)中已經(jīng)說過,在本文中由程序設(shè)置連續(xù)平面正弦電磁波以對(duì)在正弦電磁波隨機(jī)入射時(shí)混響室內(nèi)的場分布進(jìn)行仿真,但為了使推進(jìn)到各網(wǎng)格的正弦波快速地達(dá)到穩(wěn)定狀態(tài)以滿足仿真要求的條件,因而作了式(1)那樣的處理。那么這樣處理之后電磁波還會(huì)保持正弦電磁波那樣的特點(diǎn)嗎?這就需要檢驗(yàn)。圖4和圖5分別為當(dāng)將正弦電場和磁場離散為200個(gè)取樣點(diǎn)后,再用式(1)處理的結(jié)果。由圖可以看出,其波仍保持著好的正弦特征,且電場幅度和磁場幅度之比正好等于波阻抗。

3.2 隨機(jī)入射平面電磁波的等效源分布

本文中,分別將產(chǎn)生隨機(jī)入射正弦平面電磁波的等效源的個(gè)數(shù)設(shè)置為 1、10、20、100、200 、300、400和500個(gè),以便從結(jié)果中看出場分布與源的個(gè)數(shù)間的關(guān)系。在此,僅將源為500個(gè)時(shí)它們?cè)谇蛎嫔系姆植际居趫D6中。當(dāng)然,對(duì)它們具體位置的設(shè)置是按照2.1節(jié)中的約束條件選取的。其中,圖6(a)為觀察者位于球坐標(biāo)θ=28°,φ=60°處的俯視圖,圖6(b)為觀察者位于θ=32°,φ=23°處的俯視圖。

從圖中可以看出,波源在球面上的分布是隨機(jī)均勻的,是滿足統(tǒng)計(jì)規(guī)律要求的。

3.3 混響室內(nèi)場分布仿真及數(shù)據(jù)分析

在此,以對(duì)尺寸為2.1 m×2.3 m×3.5 m的混響室的工作區(qū)域按90×90×90個(gè)網(wǎng)格劃分后進(jìn)行場分布仿真計(jì)算為例。對(duì)該區(qū)域內(nèi)場的均勻性與哪些因素有關(guān)進(jìn)行討論。具體做法是,分別選擇隨機(jī)入射電磁波的個(gè)數(shù)為 1、10、20 、100 、200 、300 、400 和500個(gè),計(jì)算區(qū)域內(nèi)不同剖面層中的場分布,其結(jié)果以圖示形式給出。限于文章篇幅,文中僅給出了yoz平面內(nèi)的第45層,xoz平面內(nèi)的第56層和xoy平面內(nèi)的第68層的電場分布于圖7的(a1)→(h1),(a2)→(h2)和(a3)→(h3)中。圖中N表示重疊入射波的個(gè)數(shù),網(wǎng)格從1→10、80→91代表吸收邊界到連接邊界的網(wǎng)格數(shù)。在yoz、xoz和xoy三個(gè)平面中的第45層、第56層和第68層的這三個(gè)層面是任意選取的。

其中層數(shù)的劃分是對(duì)90×90×90的計(jì)算區(qū)域進(jìn)行的,每個(gè)坐標(biāo)方向上均可化為90層,如圖8所示,虛線正方體為計(jì)算區(qū)域。

圖8 計(jì)算區(qū)域?qū)訑?shù)劃分示意圖

基于上述仿真結(jié)果并按照IEC61000-4-21標(biāo)準(zhǔn)中規(guī)定的計(jì)算方法就可以計(jì)算出混響室工作區(qū)域中的場均勻度,如表1所示。當(dāng)入射波的個(gè)數(shù)為1時(shí),電場強(qiáng)度的標(biāo)準(zhǔn)偏差為3.40 dB,它不滿足標(biāo)準(zhǔn)中所提出的3 dB要求,并且可從圖7的(a1)、(a2)和(a3)直觀看出,此時(shí)場均勻性較差,它們的電場強(qiáng)度起伏較大,且包絡(luò)為正弦規(guī)律,這就說明單個(gè)波入射的情況下混響室內(nèi)場的均勻性較差;但隨著入射電磁波數(shù)目的增加,比如,當(dāng)入射波的個(gè)數(shù)大于或等于10時(shí),標(biāo)準(zhǔn)偏差在小于3 dB的前提下逐漸減小,這表明場的均勻性越來越好,但與此同時(shí),計(jì)算時(shí)間亦會(huì)逐漸增加,對(duì)計(jì)算機(jī)硬件的要求也同步提高。因此,對(duì)那些復(fù)雜的、要求精度高的問題在單個(gè)計(jì)算機(jī)不能處理或其消耗的時(shí)間相當(dāng)長的情況下可采用并行算法進(jìn)行求解。上述結(jié)果也提示人們,在對(duì)研究對(duì)象進(jìn)行數(shù)值分析之前必須要依據(jù)問題的復(fù)雜性及其求解的精度要求,適當(dāng)?shù)剡x取入射波的個(gè)數(shù),以達(dá)到既節(jié)約計(jì)算時(shí)間又滿足求解精度的目的。

表1 頻率為915 MHz時(shí)場的均勻性隨入射波個(gè)數(shù)的變化

4.結(jié) 論

把重疊波理論運(yùn)用到混響室的仿真分析中,使得計(jì)算時(shí)間只與測試體的特性相關(guān),這樣不僅能使仿真過程中混響室內(nèi)的電磁環(huán)境與實(shí)際混響室內(nèi)的電磁環(huán)境更接近,而且與傳統(tǒng)的混響室仿真分析相比還可以大大縮短計(jì)算時(shí)間;同時(shí),分析比較了不同重疊入射平面波情況下混響室內(nèi)場分布的均勻度,得出了欲使混響室內(nèi)場分布達(dá)到預(yù)設(shè)均勻度要求時(shí),所應(yīng)選取的最佳重疊入射平面波的數(shù)目,而本文工作中所表明的混響室內(nèi)這種場分布的預(yù)期精度與應(yīng)選取重疊入射平面波個(gè)數(shù)間的關(guān)系是文獻(xiàn)[14]沒有提及和體現(xiàn)的,從而為混響室的設(shè)計(jì)提供了可靠的理論指導(dǎo)。

[1] LEGNAB T,FREYER G,HATFIELD M,et al.Verification of fields applied to an EUT in a Reverberation chamber using numerical modeling[C]//IEEE Int.Symp.Electrom.Compatibility,1988,1：28-33.

[2] BONNET P,VERNET R,GIRARD S,et al.FDTD modeli-ng of reverberation chamber[J].Electronics Letters,2005,41(20)：1101-1102.

[3] BAI L,WANG L,WANG B,et al.Reverberation chamber modeling using FDTD[C]//IEEE Int.Symp.Electromagn.Compatibility,1999,1：7-11.

[4] CHUNG S Y,RHEE J G,RHEE H J,et al.Field unitormity characteristics of an asymmetric structure reverberation chamber by FDTD method[C]//IEEE Int.Symp.Electromagnetic Compatibility,2001,1：429-434.

[5] HARIM A K.FDTD analysis of electromagnetic fields in a reverberation chamber[J].IEICE Trans.Communication,1998,E81-B(10)：1946-1950.

[6] KOUVELIOTIS N K,TRAKADAS P T,CAPSALIS C N.Examination of field uniformity in vibration intrinsic reverberation chamber using the FDTD method[J].Electro.Letters,2002,38(3)：109-110.

[7] BUNTING C F.Statistical characterization and the simulation of a reverberation chamber using finite-element techniques[J].IEEE Trans.Electromagnetic Compatibility,2002,41(1)：214-221.

[8] BUNTING C F,MOELLER K J,REDDY C J,et al.A two-dimentional finite-element analysis of reverberation chambers[J].IEEE Trans.Electromagnetic Compatibility,1999,41(4)：280-289.

[9] HOEPPE F,GINESTE P,DEMOULIN B,et al.Numerical predictions applied to mode-stirred reverberation chambers[C]//Proc.Reverberation Chamber,Anechoic Chamber and OATS Users Meeting.Seattle,WA,2001.

[10] SURIANO C R,THIELE G A,SURIANO J R.Low frequency behavior of a reverberation chamber with monopole antenna[C]//IEEE Int.Symp.Electromagnetic Compatibility,2002,2：645-650.

[11] HILL D A.Plane wave integral representation for fields in reverberation chambers[J].IEEE Tran.E-lectromagnetic Compatibility,1998,40(3)：209-217.

[12] PETIRSCH M,SCHWAB A J.Investigation of the field uniformity of a mode-stirred chamber using diffusers based on acoustic theory[J].IEEE Transaction electromagnetic Compatibility,1999,41(4)：446-451.

[13] 沈遠(yuǎn)茂,石 丹,高攸綱,等.利用多天線源攪拌改善混響室場均勻性的分析[J].電波科學(xué)學(xué)報(bào),2009,24(4)：682-686.SHEN Yuanmao,SHI Dan,GAO Yongang,et al.Improvement in field uniformity introduced by multiple-antenna in source-stirred reverberation chamber[J].Chinese Journal of Radio Science,2009,24(4)：682-686.(in Chinese)

[14]MOGLIEF,PASTORE A P.FDTDanalysis of plane wave superposition to simulate susceptibility tests in reverberation chambers[J].IEEE Trans.Electromagn.Compatibility,2006,48(1)：195-202.

[15] 葛徳彪,石守元,朱之偉.一種新的FDTD入射場設(shè)置方法[J].微波學(xué)報(bào),1995,11(3)：187-190.GE Debiao,SHI Shouyuan,ZHU Zhiwei.A new FDTD scheme for introducing incident fields[J].Journal of Microwaves,1995,11(3)：187-190.(in Chinese)

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