雙時(shí)間步時(shí)域有限體積方法計(jì)算時(shí)變電磁場

雙時(shí)間步時(shí)域有限體積方法計(jì)算時(shí)變電磁場

許勇黃勇周鑄

(中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所,四川 綿陽 621000)

摘要雙時(shí)間步被引入到時(shí)域有限體積解算器這一直接求解麥克斯韋方程組的高精度數(shù)值方法中．雙時(shí)間步方法作為非定常時(shí)間推進(jìn)技巧,時(shí)間精度由物理時(shí)間步長決定,穩(wěn)定性要求由定常子迭代時(shí)間步長所滿足,從而放松了通常顯式時(shí)間格式和網(wǎng)格對(duì)物理時(shí)間步長的限制,達(dá)到節(jié)約計(jì)算量的目的．典型目標(biāo)電磁散射計(jì)算表明：通過對(duì)物理時(shí)間步長、最大子迭代步數(shù)、子迭代收斂判據(jù)的合理選取,雙時(shí)間步方法在保證計(jì)算精度的同時(shí),能提高計(jì)算效率．

關(guān)鍵詞雙時(shí)間步;時(shí)域有限體積法;電磁散射;雷達(dá)截面

中圖分類號(hào)O441.4

文獻(xiàn)標(biāo)志碼A

文章編號(hào)1005-0388(2015)04-0647-06

AbstractDual time stepping technique was introduced into the finite volume time domain (FVTD) method which is a high precision direct numerical solvers for time dependent electromagnetic fields controlled by Maxwell equations. As an unsteady time marching technique, the time precision is determined by physical time step, stability criterion is guaranteed by steady sub-iteration time step in this method, as a result, the physical time step was relaxed which ordinary limited by explicit time scheme and grids, the computational work is also decreased. The computation of canonical electromagnetic scattering problem demonstrated that the computational precision is conserved and efficiency is also increased; by the rational choose of physical time step, the maximum step and convergence criterion of sub-iteration in this technique.

收稿日期：2014-06-26

作者簡介

Dual time stepping FVTD method for computation of

time-dependent electromagnetic fields

XU YongHUANG YongZHOU Zhu

(ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchand

DevelopmentCenter,MianyangSichuan621000,China)

Key words dual time stepping; finite volume time domain method; electromagnetic scattering; radar cross section

資助項(xiàng)目: 裝備預(yù)先研究項(xiàng)目(51313010303)

聯(lián)系人： 周 鑄 E-mail: stephen000@sina.com

引言

時(shí)變電磁場滿足時(shí)域麥克斯韋方程組,隨著計(jì)算機(jī)技術(shù)的發(fā)展,直接求解該方程組成為可能．與歐拉方程相同的雙曲型數(shù)學(xué)特征促進(jìn)了計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)技術(shù)在電磁場計(jì)算中的應(yīng)用,其中時(shí)域有限差分(Finite-Difference Time-Domain,FDTD)法[1]和時(shí)域有限體積(Finite-Volume Time-Domain,FVTD)法[2-5]最為著名．

這些方法的時(shí)間推進(jìn)或采用二階中心差分或時(shí)空耦合Lax-Wendroff格式以及多步Runge-Kutta法,其共同點(diǎn)是時(shí)間計(jì)算的顯式格式．以Runge-Kutta方法為代表的顯式方法,既有編程簡便,又有易實(shí)現(xiàn)時(shí)間高精度的優(yōu)點(diǎn),對(duì)時(shí)域電磁場計(jì)算是可靠的時(shí)間離散方法．但是時(shí)間顯式方法有一個(gè)最大的缺陷,其時(shí)間步長受穩(wěn)定性的限制,整個(gè)計(jì)算空間必須采用統(tǒng)一的最小全局計(jì)算步長,為模擬幾何外形劇烈變化生成的貼體加密網(wǎng)格,會(huì)帶來很小的全局時(shí)間步長,從而使得到穩(wěn)定的時(shí)變電磁場需要較長的計(jì)算時(shí)間,特別是在求解多維非定常問題時(shí),計(jì)算時(shí)間的限制帶來計(jì)算量的顯著增加．另一方面,隱式計(jì)算方法雖能夠放寬計(jì)算步長的穩(wěn)定性限制,但伴隨而來的是時(shí)間精度下降和矩陣求逆運(yùn)算．

Jameson在CFD領(lǐng)域提出了雙時(shí)間步(Dual Time Stepping)隱式迭代方法[6]．該方法通過在控制方程中引入定常的虛擬時(shí)間導(dǎo)數(shù)項(xiàng),在每個(gè)物理時(shí)間段內(nèi)對(duì)物理時(shí)間導(dǎo)數(shù)作線性化處理,從而使得物理時(shí)間推進(jìn)步長可以根據(jù)物理問題進(jìn)行選取而不受穩(wěn)定性的限制,計(jì)算的穩(wěn)定性要求由虛擬時(shí)間子迭代滿足,子迭代的定常計(jì)算可以采用局部時(shí)間步長加速收斂．從理論上來講,引入虛擬時(shí)間迭代過程的雙時(shí)間步法,只有在子迭代殘值趨于零時(shí)才接近時(shí)間高精度,但實(shí)際計(jì)算往往難以滿足這一點(diǎn),因此子迭代過程對(duì)非定常計(jì)算的影響是一個(gè)值得關(guān)注的問題．

雙時(shí)間步法在計(jì)算效率方面具有顯著的優(yōu)勢(shì),大量的CFD計(jì)算充分證實(shí)了這一點(diǎn)．但是雙時(shí)間步法求解時(shí)變電磁場時(shí),需對(duì)一些因素進(jìn)行必要的研究,包括物理時(shí)間迭代精度、物理時(shí)間步長、子迭代收斂標(biāo)準(zhǔn)以及最大子迭代步數(shù)的選擇,它們協(xié)同作用才能保證時(shí)間推進(jìn)的效率和計(jì)算精度．本文即在已有顯式多步Runge-Kutta電磁場FVTD解算器的基礎(chǔ)上,拓展應(yīng)用子迭代為多步Runge-Kutta的雙時(shí)間步法計(jì)算電磁散射場,極大縮短穩(wěn)定電磁場所需計(jì)算時(shí)間,并計(jì)算和分析以上因素對(duì)計(jì)算結(jié)果的影響規(guī)律．

1數(shù)值方法

1.1控制方程

時(shí)域麥克斯韋方程組的兩個(gè)旋度方程(法拉第電磁感應(yīng)定律和安培環(huán)路定律),在無源條件下的直角坐標(biāo)系守恒形式為

(1)

式中: Q=[BxByBzDxDyDz]T=[BD]T,

B是磁感應(yīng)強(qiáng)度矢量,D是電位移矢量; Fx=[0

-EzEy0Hz-Hy]T; Fy=[Ez0-Ex-Hz

0Hx]T; Fz=[-EyEx0Hy-Hx0]T.

對(duì)于復(fù)雜外形物體,在計(jì)算空間生成貼體多塊結(jié)構(gòu)網(wǎng)格,其坐標(biāo)變換為

k=k(x,y,z),k=ξ,η,ζ.

(2)

并得到曲線坐標(biāo)系下的麥克斯韋方程組守恒型

(3)

1.2空間積分方法

在計(jì)算中,守恒電磁場量在自由空間中使用散射場,可保持守恒方程形式,同式(1),從而避免在網(wǎng)格空間傳播入射場帶來的數(shù)值誤差,也有利于截?cái)嗫臻g無反射邊界條件的應(yīng)用,入射電磁場可由入射波解析給出,在完全導(dǎo)電壁處由邊界條件引入,這也相應(yīng)避免了傳統(tǒng)時(shí)域有限差分法中劃分總場區(qū)和散射場區(qū)的繁瑣做法．

有限體積法的空間精度體現(xiàn)在能否精確模擬守恒變量Q在網(wǎng)格單元分界面處的狀態(tài)變量,以得到相應(yīng)精確的分界面通量．這里采用Steger-Warming分裂進(jìn)行網(wǎng)格單元分界面通量計(jì)算

(4)

單元邊界左右流通量可統(tǒng)一寫為:

(5)

(6)

式中: S,S-為相似矩陣; Λ+,Λ-分別為正負(fù)特征值構(gòu)成的對(duì)角矩陣; QL,QR分別代表分界面處左右狀態(tài)變量,它們可采用MUSCL格式而達(dá)到最高三階精度:

(7)

(8)

式(7)～(8)中:φ是限制器,本文中取φ=1;κ=1/3是三階精度格式的控制參數(shù);和Δ分別是后差和前差算符．

1.3雙時(shí)間步法

對(duì)于時(shí)空分別積分的半離散線方法,顯式時(shí)間計(jì)算方面采用常微分方程求解中的Runge-Kutta法:

Qn+k/m=Qn-λαkR(Qn+(k-1)/m),k=1,m,

(9)

λ=Δt/V,

(10)

αk=1/m-k+1．

(11)

顯式方法穩(wěn)定性要求全網(wǎng)格空間采用相同物理時(shí)間步長,會(huì)加大計(jì)算量;雙時(shí)間步方法則采用時(shí)間雙重循環(huán),外循環(huán)為物理時(shí)間推進(jìn),內(nèi)循環(huán)為虛擬時(shí)間迭代,相應(yīng)雙時(shí)間步法控制方程組修改為

(12)

(13)

式中:

(14)

Q*是Qn+1(n是物理時(shí)間步)的近似,物理時(shí)間導(dǎo)數(shù)采用后向差分為二階時(shí)間精度．式(13)類似式(9)作顯式子迭代:

k=1,m,

(15)

上標(biāo)l、l+1代表虛擬定常子迭代過程中的兩個(gè)相鄰時(shí)間層．

在數(shù)值模擬過程中,需要對(duì)內(nèi)迭代收斂過程進(jìn)行判斷和控制,使得內(nèi)迭代過程在保證一定時(shí)間精度的前提下盡快收斂．Arnone指出[8],式(15)在物理時(shí)間步小于或與虛擬時(shí)間步同量級(jí)時(shí)變得不穩(wěn)定．Melson[9]進(jìn)一步明確不穩(wěn)定性來源于殘差中3Q*/(2Δt)的顯式計(jì)算,物理時(shí)間步Δt越小不穩(wěn)定越明顯,解決之道就是隱式化該項(xiàng)并作線性化處理,最后得到每層子迭代的Runge-Kutta關(guān)系為

(16)

除了式(11)對(duì)應(yīng)最大時(shí)間精度系數(shù)外,可以針對(duì)空間格式調(diào)節(jié)系數(shù)來增加最大時(shí)間步和改善穩(wěn)定性．表1即為對(duì)應(yīng)一階空間格式,步數(shù)分別為3、4、5的Runge-Kutta最大CFL數(shù)和相應(yīng)系數(shù)．

表1　多步Runge-Kutta系數(shù)和最大CFL數(shù)

由此得到簡諧波入射FVTD計(jì)算目標(biāo)雷達(dá)截面(Radar Cross-Section,RCS)步驟:首先經(jīng)時(shí)間推進(jìn)在整個(gè)網(wǎng)格區(qū)域建立穩(wěn)定的電磁振蕩,然后在空間Stratton-Chu積分面處逐周期進(jìn)行傅里葉變換,獲得收斂的頻域電磁散射場,進(jìn)而由積分方程得到目標(biāo)雙站RCS分布[7]．

2計(jì)算結(jié)果與分析

雙時(shí)間步方法時(shí)間精度由物理時(shí)間步長確定,為二階時(shí)間精度;采用越小的物理時(shí)間步長時(shí)間精度就越有保證,采用較大的物理時(shí)間間隔有利于提高計(jì)算效率,但物理時(shí)間步長也不能夠無節(jié)制增加．圖1比較了不同物理時(shí)間步長對(duì)圓柱RCS結(jié)果的影響,k為波數(shù),a為圓柱半徑,計(jì)算網(wǎng)格為51×46,其中物面網(wǎng)格選取每波長20個(gè)網(wǎng)格點(diǎn),遠(yuǎn)場邊界在3波長外,輻向網(wǎng)格壁面加密,平均約每波長15個(gè)網(wǎng)格點(diǎn)．

(a) TM波

(b) TE波 圖1　物理時(shí)間步長對(duì)圓柱(ka=2)RCS的影響

物理時(shí)間步長決定了傅里葉變換中取樣周期的點(diǎn)數(shù)目,在保證精度前提下的更大取值有利于節(jié)約計(jì)算量,其選取與目標(biāo)電尺度無關(guān),跟電磁場空間散射特性相關(guān),其選取由典型二維、三維計(jì)算驗(yàn)證選取．由圖1知二維問題用入射電磁波周期無量綱化的物理時(shí)間步長Δt=0.01計(jì)算精度已可以,建議三維問題選擇Δt≈0.001量級(jí),隨著物理時(shí)間間隔的增大,計(jì)算結(jié)果與參考值差別逐漸增大,在過大的物理時(shí)間步長條件下,雙時(shí)間步方法的計(jì)算精度變低．

雙時(shí)間方法非定常計(jì)算的時(shí)間精度還受到每一物理時(shí)間步上子迭代步數(shù)制約．迭代步數(shù)多,收斂控制嚴(yán)格,時(shí)間精度才有保證．真實(shí)時(shí)間步小,則子迭代可以取得少一些,數(shù)值計(jì)算中為防止場變化劇烈處、以及子迭代CFL數(shù)取得過大等情況時(shí),殘差很有可能無法下降到收斂判據(jù)值,陷入死循環(huán)而給出最大子迭代步數(shù),另一方面對(duì)子迭代收斂進(jìn)行判斷和控制,在保證一定時(shí)間精度的前提下盡快結(jié)束迭代．本文子迭代收斂判據(jù)使用電場和磁場在計(jì)算空間最大幅度差的絕對(duì)值．

圖2是TM波入射圓柱(ka=10,網(wǎng)格為201×46)最大子迭代步數(shù)的影響,圖3是TE波入射圓柱(ka=10)子迭代收斂判據(jù)的影響,可見子迭代殘差下降2個(gè)量級(jí)已能滿足要求,這與CFD中殘差下降2～3個(gè)量級(jí)一致．

圖2　最大子迭代步數(shù)對(duì)圓柱(ka=10)RCS的影響

圖3　子迭代收斂判據(jù)對(duì)圓柱(ka=10)RCS的影響

以Runge-Kutta方法為代表的顯式方法容易實(shí)現(xiàn)時(shí)間求解的高精度,對(duì)非定常計(jì)算來說是理論上最為可靠的時(shí)間離散方法,以下針對(duì)雙時(shí)間步方法與顯式Runge-Kutta方法,就時(shí)間精度、結(jié)果好壞及計(jì)算效率作一些比較．

采用完全導(dǎo)電金屬球電磁散射為例,其電尺寸為ka=2,a為球半徑．計(jì)算網(wǎng)格為61×31×46,時(shí)間迭代都是四步Runge-Kutta推進(jìn),其他計(jì)算條件相同．圖4同一空間點(diǎn)分別用顯式Runge-Kutta法和雙時(shí)間步方法計(jì)算,其散射電磁場量的時(shí)間振蕩歷程、空間兩者幾無區(qū)別,說明物理時(shí)間后向差分、二階精度的雙時(shí)間方法精度滿足傅里葉變換需要．

圖5比較了顯式和雙時(shí)間步計(jì)算的RCS雙站分布,可見后者在前向(θ=180°)與解析解吻合更好．圖6則是這兩種方法逐周期進(jìn)行傅里葉變換所得頻域值收斂歷程比較,在達(dá)到同樣的收斂標(biāo)準(zhǔn)(0.001),雙時(shí)間步收斂略快,這對(duì)大網(wǎng)格量和存在很小體積網(wǎng)格單元的情況,能大大提高計(jì)算效率．因?yàn)榇藭r(shí)顯式Δt≈10-4量級(jí)甚至更小,相反雙時(shí)間步方法則根據(jù)時(shí)間精度選擇Δt,并可用局部時(shí)間步等技巧加速虛擬時(shí)間子迭代的收斂．

圖4　時(shí)變散射場振蕩歷程比較(金屬球,ka=2)

圖5　雙時(shí)間步計(jì)算金屬球雙站(RCS) 分布比較(ka=2)

圖6　電磁場頻域收斂歷程比較(金屬球,ka=2)

3結(jié)論

雙時(shí)間步時(shí)域有限體積方法計(jì)算時(shí)變電磁場,徹底放松了顯式格式穩(wěn)定性對(duì)物理時(shí)間步的限制,使得時(shí)間精度由物理時(shí)間間隔決定,穩(wěn)定性由定常虛擬時(shí)間子迭代(可以應(yīng)用加速收斂技巧)保證,該方法可以由添加很少幾段程序由一個(gè)顯式計(jì)算程序修改而來,這樣既能保證計(jì)算精度又能大大提高計(jì)算效率的方法,值得在時(shí)域電磁場數(shù)值計(jì)算中推廣應(yīng)用．

致謝:衷心感謝中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所江雄、劉剛等人在雙時(shí)間方法應(yīng)用上提供的幫助,使得該算法能在時(shí)變電磁場的計(jì)算中得到推廣．

參考文獻(xiàn)

[1]YEE K S. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media [J]. IEEE Trans on AP, 1966, 14(5):302-307.

[2]SHANKAR V, MOHAMMADIAN A H, HALL W F. A time-domain, finite-volume treatment for the Maxwell equations [J]. Electromagnetics, 1990, 10:127-145.

[3]SHANG J S, GAITONDE D. Scattered electromagnetic field of a reentry vehicle [C]//AIAA 94-0231, 1994.

[4] 路占波, 張安, 候新宇，等. 采用FDTD/FVTD混合算法分析蝶形微帶天線[J]. 電波科學(xué)學(xué)報(bào), 2007, 22(3): 527-532.

LU Zhanbo, ZHANG An, HOU Xinyu, et al. Analysis of bow-tie microstrip antenna by hybrid FDTD/FVTD algorithm[J]. Journal of Radio Science, 2007, 22(3): 527-532. (in Chinese)

[5]鄧聰, 尹文祿, 柴舜連, 等. 兩種適用時(shí)域有限體積法的截?cái)噙吔鏪J]. 電波科學(xué)學(xué)報(bào), 2009, 24(3): 537-540.

DENG Cong, YIN Wenlu, CHAI Shunlian, et al. Two boundary conditions for the truncation of cell-centered FVTD algorithm on unstructured lattices [J]. Chinese Journal of Radio Science, 2009, 24(3):537-540. (in Chinese)

[6]JAMESON A. Time dependent calculations using multigrid with applications to unsteady flows past airfoils and wings [C]//AIAA 91-1596, 1991.

[7]許勇, 樂嘉陵. 基于CFD的電磁散射數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2004, 22(2): 185-189.

XU Yong, LE Jialing. CFD based numerical simulation of electromagnetic scattering [J]. Acta Aerodynamica Sinica, 2004, 22(2):185-189. (in Chinese)

[8]ARNONE A, LIOU M S, POVINELLI L A. Multigrid time accurate integration of the navier-stokes equations [C]//AIAA 93-3361, 1993.

[9]MELSON N D, SANETRIK M D, ATKINS H L. Time-accurate Navier-Stokes calculations with multigrid acceleration [C]//Proc 6th Copper Mountain Conf. on Multigrid Methods, 1993: 423-439.

許勇(1971-)，男，四川人，中國空氣動(dòng)力研究與發(fā)展中心副研究員，研究領(lǐng)域?yàn)橛?jì)算電磁學(xué)，研究方向包括再入體電磁特性、電磁場數(shù)值計(jì)算方法.

黃勇(1970-)，男，四川人，中國空氣動(dòng)力研究與發(fā)展中心研究員．研究方向?yàn)橛?jì)算空氣動(dòng)力學(xué)、飛行器氣動(dòng)多目標(biāo)優(yōu)化設(shè)計(jì)．

周鑄 (1973-),男,重慶人,中國空氣動(dòng)力研究與發(fā)展中心研究員.研究方向?yàn)橛?jì)算空氣動(dòng)力學(xué)、復(fù)雜流動(dòng)分析.

胡榮旭,吳振森,張金鵬. 中國近海蒸發(fā)波導(dǎo)反演中最佳雷達(dá)參數(shù)分析[J]. 電波科學(xué)學(xué)報(bào)，2015,30(4):653-660. doi: 10.13443/j.cjors. 2014080501

HU Rongxu, WU Zhensen, ZHANG Jinpeng. Analysis for best radar parameter in inversion of evaporation duct above China sea areas[J]. Chinese Journal of Radio Science,2015,30(4):653-660. (in Chinese). doi: 10.13443/j.cjors. 2014080501