基于時域多分辨分析方法的地閃雷電電磁脈沖計算

江志東 周璧華 劉亞文 楊 波

（解放軍理工大學 電磁環境效應與電光工程國家級重點實驗室，江蘇 南京210007）

引 言

雷電是發生在大氣中的一種瞬態大電流、高電壓和強電磁輻射的天氣現象.雷電電磁脈沖（Lightning Electromagnetic Pulse，LEMP）的計算對于架空電纜耦合評估、電氣和電子設備防護以及雷電間接效應的研究至關重要.目前計算地閃放電回擊通道輻射的LEMP有三類方法［1］，即精確解析法、近似法和數值計算法.其中時域有限差分（Finite Difference-Time Domain，FDTD）法具有大地電參數引入方便、一次計算可得到所有場點信息、得出的時域波形可直接用于耦合計算及便于編程實現等優點.自文獻［2］首次將FDTD法引入到地閃放電回擊通道幾十米范圍內LEMP的計算后，FDTD法廣泛應用于地閃LEMP計算［3］、LEMP線纜耦合電壓計算［4］、其他算法精度評估［5］等方面.

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

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

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

1 基于MRTD的LEMP計算模型

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

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

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

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

a（v）的值參考文獻［11］，根據a（v）的對稱關系a（-v-1）＝-a（v），只需要考慮0≤v≤2時的情形.

1.1 二維柱坐標下LEMP計算模型

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

圖1 二維柱坐標下MRTD計算域網格劃分

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

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

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

由于地閃回擊電流位于二維柱坐標系中的z軸上，在圖1的LEMP計算模型中，z軸上只有Ez一個分量，因此可以通過Ez在MRTD差分網格中引入激勵源.

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

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

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

2 數值仿真及分析

2.1 數值計算結果對比

數值計算采用的地閃回擊通道模型為（modified transmission line model with linear current decay with height，MTLL）模型，回擊通道高度設為7 200m，回擊速度為1.1×108m／s.通道基電流表達采用雙指數函數，分別選用1.2／50μs和0.25／100 μs波形作為首次回擊和后續回擊雷電流的典型波形.大地電參數σg＝0.001S／m，εrg＝10.計算空間3 600m×2 400m.

如圖2所示，對于首次回擊和后續回擊LEMP的計算，FDTD法分別采用1m×1m和0.5m×0.5m的空間網格剖分，其計算的場分量作為標準參考值.MRTD法分別采用3m×3m和2m×2m的空間網格剖分，計算結果表明了MRTD法的有效性.同時，不同網格剖分對于垂直電場的影響不大，水平電場隨著網格的變大，在波形的初始階段和峰值處出現一定程度的振蕩.由于水平電場易受大地電參數的影響，計算空間網格的增大使得大地電參數對計算結果的影響更加明顯.

圖2 距回擊通道120m處垂直電場和水平電場

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

圖3 距回擊通道480m處垂直電場和水平電場

2.2 討論與分析

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

圖4 差分近似后相速和光速之比的各項異性

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

3 結 論

本文將MRTD算法引入地閃LEMP的計算，選用Daubechies小波尺度函數作為基函數，采用MTLL地閃回擊模型在二維柱坐標系下計算了地閃放電通道不同距離處LEMP場分量.數值仿真結果證實了MRTD算法的有效性，與FDTD算法相比，在保證計算精度的前提下，采用較大計算空間網格，減少計算機內存開銷.目前，采用MRTD法的主要難點在于，無論如何選擇展開基底，其時間穩定性條件比FDTD法都要苛刻.此外，雖然MRTD算法計算大目標時能夠以幾何級數的量級節省內存，但由于基函數互相有重疊，導致邊界條件復雜化，設計滿足特定應用的匹配吸收邊界是另一難點.下一步可考慮優化吸收邊界匹配層的設計，以便將MRTD法更好的應用于地閃LEMP數值分析、線纜耦合研究等相關領域.

［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］.物理學報，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尺度函數的時域多分辨分析計算電磁散射［J］.物理學報，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］高強業，周建江，曹群生.MRTD方法的色散特性分析和電磁散射應用［J］.南京航空航天大學學報，

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.