基于時(shí)域多分辨分析方法的地閃雷電電磁脈沖計(jì)算

江志東 周璧華 劉亞文 楊 波

（解放軍理工大學(xué) 電磁環(huán)境效應(yīng)與電光工程國家級重點(diǎn)實(shí)驗(yàn)室，江蘇 南京210007）

引 言

雷電是發(fā)生在大氣中的一種瞬態(tài)大電流、高電壓和強(qiáng)電磁輻射的天氣現(xiàn)象.雷電電磁脈沖（Lightning Electromagnetic Pulse，LEMP）的計(jì)算對于架空電纜耦合評估、電氣和電子設(shè)備防護(hù)以及雷電間接效應(yīng)的研究至關(guān)重要.目前計(jì)算地閃放電回?fù)敉ǖ垒椛涞腖EMP有三類方法［1］，即精確解析法、近似法和數(shù)值計(jì)算法.其中時(shí)域有限差分（Finite Difference-Time Domain，F(xiàn)DTD）法具有大地電參數(shù)引入方便、一次計(jì)算可得到所有場點(diǎn)信息、得出的時(shí)域波形可直接用于耦合計(jì)算及便于編程實(shí)現(xiàn)等優(yōu)點(diǎn).自文獻(xiàn)［2］首次將FDTD法引入到地閃放電回?fù)敉ǖ缼资追秶鷥?nèi)LEMP的計(jì)算后，F(xiàn)DTD法廣泛應(yīng)用于地閃LEMP計(jì)算［3］、LEMP線纜耦合電壓計(jì)算［4］、其他算法精度評估［5］等方面.

盡管FDTD法具有上述諸多優(yōu)點(diǎn)，但由于數(shù)值色散和數(shù)值穩(wěn)定性的影響，實(shí)際應(yīng)用時(shí)易受計(jì)算機(jī)內(nèi)存和速度的限制.如由于色散誤差的限制，空間網(wǎng)格必須足夠小，通常每個波長至少要取10個網(wǎng)格單元.由于穩(wěn)定條件的限制，時(shí)間步長取值不能太大，計(jì)算時(shí)間勢必增大.盡管當(dāng)前計(jì)算平臺發(fā)展迅速，計(jì)算效率仍是需要考慮的問題.

Krumpholz等［6］在1996年首次提出時(shí)域多分辨分析法（Multi-resolution time-domain，MRTD），該方法將矩量法（Moment of method，MoM）和多分辨分析（Multi-resolution analysis，MRA）結(jié)合來離散麥克斯韋方程組，在時(shí)域?qū)㈦妶雠c磁場分量分別展開成一系列空間上的尺度函數(shù)和時(shí)間上的脈沖函數(shù)之和.與FDTD法相比，MRTD法不僅具有更好的線性色散特性［7］，其最突出特點(diǎn)是對空間進(jìn)行劃分網(wǎng)格后每個波長采樣兩個點(diǎn)即可達(dá)到FDTD每個波長采10個點(diǎn)的精度，即Nyquist定理的極限，從而可以大大減少計(jì)算網(wǎng)格數(shù)，節(jié)省內(nèi)存和計(jì)算時(shí)間.

為提高LEMP時(shí)域計(jì)算效率、減少內(nèi)存消耗及克服經(jīng)典FDTD法的數(shù)值色散問題，本文嘗試將MRTD法引入到地閃LEMP的計(jì)算中，并與傳統(tǒng)的FDTD法計(jì)算結(jié)果對比.理論和模擬實(shí)驗(yàn)結(jié)果表明，MRTD法在保持計(jì)算精度的前提下能夠提高計(jì)算效率、節(jié)省計(jì)算資源.

1 基于MRTD的LEMP計(jì)算模型

MRTD算法從麥克斯韋方程組出發(fā)，將相應(yīng)的場分量按選定的基函數(shù)展開，并將時(shí)間微分近似為中心差分，空間微分近似為空間點(diǎn)場值加權(quán)的平均，從而建立差分迭代方程.目前各種具有有限支撐域的小波尺度函數(shù)如Haar小波［8］、以及（Cohen Daubechies Feauveau，CDF）雙正交小波［9］，Daubechies正交小波［10］等的MRTD算法已陸續(xù)被提出和應(yīng)用于諧振腔、微帶線和散射體的計(jì)算.

鑒于Daubechies正交小波尺度函數(shù)的緊支撐性和近似的移位內(nèi)插特性，能夠很好地平衡計(jì)算精度和計(jì)算復(fù)雜度之間的矛盾，本文選擇Daubechies小波的尺度函數(shù)作為展開基.

根據(jù)矩量法和多分辨分析特性，將麥克斯韋旋度方程組中電場和磁場分量利用空間和時(shí)間基函數(shù)展開，以磁場分量Hz為例有如下展式

將磁場分量Hz的展開式即式（1）代入麥克斯韋旋度方程，并采用Galerkin采樣方法［9］，得到場展開系數(shù)H1的MRTD迭代公式

a（v）的值參考文獻(xiàn)［11］，根據(jù)a（v）的對稱關(guān)系a（-v-1）＝-a（v），只需要考慮0≤v≤2時(shí)的情形.

1.1 二維柱坐標(biāo)下LEMP計(jì)算模型

近似認(rèn)為地閃放電通道垂直于地面且四周圍無窮空間，則地閃放電電流形成的場具有對稱性，可在二維柱坐標(biāo)系下求解，地閃放電通道附近的LEMP網(wǎng)格劃分如圖1所示.

圖1 二維柱坐標(biāo)下MRTD計(jì)算域網(wǎng)格劃分

二維柱坐標(biāo)系下MRTD差分方程如下，其格式和FDTD類似，不同點(diǎn)在于任意時(shí)刻每個場點(diǎn)值由相應(yīng)的展開系數(shù)和相鄰場點(diǎn)值相乘后加權(quán)求得，即式（4）、（5）和（6）中的而展開系數(shù)的個數(shù)l與小波尺度函數(shù)的緊支撐性有關(guān)，本文采用的是db2，故LS為3.由此，得到基于Daubechies小波尺度函數(shù)的MRTD迭代公式：

式（4）、（5）和（6）中σ和ε分別表示自由空間或者大地中的電導(dǎo)率和介電常數(shù)，εrg為大地的相對介電常數(shù).空氣中σ＝0，ε＝ε0，大地中σ＝σg，εg＝ε0εrg.

1.2 激勵源的引入及吸收邊界處理

由于地閃回?fù)綦娏魑挥诙S柱坐標(biāo)系中的z軸上，在圖1的LEMP計(jì)算模型中，z軸上只有Ez一個分量，因此可以通過Ez在MRTD差分網(wǎng)格中引入激勵源.

考慮到所計(jì)算問題的軸對稱性，為減少計(jì)算量，只計(jì)算包含回?fù)敉ǖ涝趦?nèi)的半個剖面內(nèi)的場，為此需要對軸線上的Ez作特殊處理.根據(jù)安培環(huán)路定理，在雷電流放電通道高度范圍外不存在雷電流，在雷電放電通道高度范圍內(nèi)，采用如下差分格式

式（7）中，Ⅰ0，j＋1／2為距離地面高（j＋1／2）Δz處的電流元.

由于MRTD場展開方式的特殊性，需對特殊點(diǎn)進(jìn)行處理.在計(jì)算空間某一點(diǎn)的電磁場分量系數(shù)時(shí)，關(guān)聯(lián)系數(shù)a（v）決定了要考慮該分量周圍電磁場系數(shù)的個數(shù).由于基函數(shù)相互有重疊，導(dǎo)致邊界條件復(fù)雜，設(shè)計(jì)滿足特定應(yīng)用的匹配層吸收邊界是MRTD應(yīng)用的一大難點(diǎn).本文采用PML完全匹配層吸收邊界，吸收層外側(cè)采用PEC進(jìn)行截?cái)?，并利用鏡像原理對其進(jìn)行處理.二維圓柱坐標(biāo)情況下自由空間及有耗介質(zhì)中PML內(nèi)各場量的處理方式參見文獻(xiàn)［12］.

2 數(shù)值仿真及分析

2.1 數(shù)值計(jì)算結(jié)果對比

數(shù)值計(jì)算采用的地閃回?fù)敉ǖ滥Ｐ蜑椋╩odified transmission line model with linear current decay with height，MTLL）模型，回?fù)敉ǖ栏叨仍O(shè)為7 200m，回?fù)羲俣葹?.1×108m／s.通道基電流表達(dá)采用雙指數(shù)函數(shù)，分別選用1.2／50μs和0.25／100 μs波形作為首次回?fù)艉秃罄m(xù)回?fù)衾纂娏鞯牡湫筒ㄐ?大地電參數(shù)σg＝0.001S／m，εrg＝10.計(jì)算空間3 600m×2 400m.

如圖2所示，對于首次回?fù)艉秃罄m(xù)回?fù)鬖EMP的計(jì)算，F(xiàn)DTD法分別采用1m×1m和0.5m×0.5m的空間網(wǎng)格剖分，其計(jì)算的場分量作為標(biāo)準(zhǔn)參考值.MRTD法分別采用3m×3m和2m×2m的空間網(wǎng)格剖分，計(jì)算結(jié)果表明了MRTD法的有效性.同時(shí)，不同網(wǎng)格剖分對于垂直電場的影響不大，水平電場隨著網(wǎng)格的變大，在波形的初始階段和峰值處出現(xiàn)一定程度的振蕩.由于水平電場易受大地電參數(shù)的影響，計(jì)算空間網(wǎng)格的增大使得大地電參數(shù)對計(jì)算結(jié)果的影響更加明顯.

圖2 距回?fù)敉ǖ?20m處垂直電場和水平電場

如圖3所示，與圖2結(jié)果類似，不同的空間網(wǎng)格剖分對垂直電場分量的計(jì)算結(jié)果沒有影響.水平電場分量的計(jì)算結(jié)果整體保持一致，隨著網(wǎng)格的變大出現(xiàn)一定程度的振蕩.

圖3 距回?fù)敉ǖ?80m處垂直電場和水平電場

2.2 討論與分析

為分析MRTD算法中差分近似后的數(shù)值色散特性，采用類似于FDTD的分析方法，將平面波方程帶入MRTD迭代方程，得到差分近似后相速和光速之比的各項(xiàng)異性的關(guān)系，如圖4所示.

圖4 差分近似后相速和光速之比的各項(xiàng)異性

圖4中λ0／δ為單個波長內(nèi)的采樣點(diǎn)數(shù)，當(dāng)λ0／δ為參變量時(shí)，從圖中相速和光速之比vφ／c與平面波傳播方向角φ之間的關(guān)系可以看出，MRTD算法單個波長內(nèi)采樣點(diǎn)取3時(shí)的各項(xiàng)異性程度優(yōu)于FDTDF法采樣點(diǎn)為15時(shí)的情形，故MRTD法在取較大網(wǎng)格時(shí)仍然能夠保持較高的精度和較低的數(shù)值色散誤差.

3 結(jié) 論

本文將MRTD算法引入地閃LEMP的計(jì)算，選用Daubechies小波尺度函數(shù)作為基函數(shù)，采用MTLL地閃回?fù)裟Ｐ驮诙S柱坐標(biāo)系下計(jì)算了地閃放電通道不同距離處LEMP場分量.數(shù)值仿真結(jié)果證實(shí)了MRTD算法的有效性，與FDTD算法相比，在保證計(jì)算精度的前提下，采用較大計(jì)算空間網(wǎng)格，減少計(jì)算機(jī)內(nèi)存開銷.目前，采用MRTD法的主要難點(diǎn)在于，無論如何選擇展開基底，其時(shí)間穩(wěn)定性條件比FDTD法都要苛刻.此外，雖然MRTD算法計(jì)算大目標(biāo)時(shí)能夠以幾何級數(shù)的量級節(jié)省內(nèi)存，但由于基函數(shù)互相有重疊，導(dǎo)致邊界條件復(fù)雜化，設(shè)計(jì)滿足特定應(yīng)用的匹配吸收邊界是另一難點(diǎn).下一步可考慮優(yōu)化吸收邊界匹配層的設(shè)計(jì)，以便將MRTD法更好的應(yīng)用于地閃LEMP數(shù)值分析、線纜耦合研究等相關(guān)領(lǐng)域.

［1］RAKOV V A，RACHIDI F.Overview of recent progress in lightning research and lightning protection［J］.IEEE Trans Electromagn Compat，2009，51（3）：428-442.

［2］YANG Chunshan，ZHOU Bihua.Calculation methods of electromagnetic fields very close to lightning［J］.IEEE Trans Electromagn Compat，2004，46（1）：133-141.

［3］楊 波，周璧華，孟 鑫.地閃雷電電磁脈沖在大地中的分布研究［J］.物理學(xué)報(bào)，2010，59（12）：8978-8985.YANG Bo，ZHOU Bihua，MENG Xin.Distribution of cloud-to-ground lightning electromagnetic pulse fields under the ground［J］.Acta Physica Sinica，2010，59（12）：8978-8985.（in Chinese）

［4］YANG Bo，ZHOU Bihua，CHEN Bin，et al.Numerical study of lightning-induced currents on buried cables and shield wire protection method［J］.IEEE Trans Electromagn Compat，2012，54（2）：323-331.

［5］ZHANG Qilin，LI Dongshuai，Zhang Yuanyuan，et al.On the accuracy of wait＇s formula along a mixed propagation path within 1km from the lightning channel［J］.IEEE Trans Electromagn Compat，2012，54（5）：1042-1047.

［6］KRUMPHOLZ M，KATEHI L P B.MRTD：new time-domain schemes based on multiresolution analysis［J］.IEEE Transactions on Microwave Theory and Techniques，1996，44（4）：555-571.

［7］代少玉，吳振森，徐仰彬.用基于Daubechies尺度函數(shù)的時(shí)域多分辨分析計(jì)算電磁散射［J］.物理學(xué)報(bào)，2007，56（2）：786-790.DAI Shaoyu，WU Zhensen，XU Yangbin.Using the MRTD based on Daubechies scaling function to solve the problem of electromagnetic scattering［J］.Acta Physica Sinica，2007，56（2）：786-790.（in Chinese）

［8］FUJII M，HOEFER W J R.Multiresolution analysis similar to the FDTD method-derivation and application［J］.IEEE Transactions on Microwave Theory and Techniques，1998，46（12）：2463-2475.

［9］CHEONG Y W，LEE Y M，RA K H，et al.Wavelet-Galerkin scheme of time-dependent inhomogeneous electromagnetic problems［J］.IEEE Microwave Guid Wave Lett，1999，9（8）：297-299.

［10］KOVVALI N，LIN W，CARIN D L.Order of accuracy analysis for multiresolution time-domain using Daubechies bases［J］.Microw Opt Technol Lett，2005，45（4）：290-293.

［11］高強(qiáng)業(yè)，周建江，曹群生.MRTD方法的色散特性分析和電磁散射應(yīng)用［J］.南京航空航天大學(xué)學(xué)報(bào)，

2010，42（2）：191-197. GAO Qiangye，ZHOU Jianjiang，CAO Qunsheng.Dispersion property analysis and electromagnetic scattering application for MRTD Method［J］.Journal of Nanjing University of Aeronautics ＆Astronautics，2010，42（2）：191-197.（in Chinese）［12］LIU Yawen，CHEN Yiwang，CHEN Bin，et al.A cylindrical MRTD algorithm with PML and Quasi-PML［J］.IEEE Transactions on Microwave Theory and Techniques，2013，61（3）：1006-1017.