樊亞東 于建立, 詹清華 王建國 周 蜜 蔡 力 戚 為

（1.武漢大學(xué)電氣工程學(xué)院 武漢 430072 2.東北電力大學(xué)輸變電技術(shù)學(xué)院 吉林 132012 3.廣東電網(wǎng)公司 佛山 528000 4.國網(wǎng)寧波鄞州供電局 寧波 315100 ）

1 引言

雷電是造成架空電力線路故障和事故最主要的外部因素[1]，長期以來線路防雷的工作重點集中在直擊雷防護(hù)，對雷電感應(yīng)過電壓的研究相對較少[2]。隨著配電線路雷擊故障率在電網(wǎng)故障率中的比重越來越大，后者危害逐漸凸顯。配電線路走廊密度大且絕緣水平低，對雷電感應(yīng)過電壓敏感程度高，其故障率占雷擊總故障率的90%以上[2,3]。

傳統(tǒng)的雷電感應(yīng)過電壓常用解析法計算：國內(nèi)武漢高壓研究院在上世紀(jì) 80年代分析擬合出峰值的近似計算式[4]；日本和西方一些國家的學(xué)者研究并提出相應(yīng)計算公式，部分計算結(jié)果比觀測值小很多[4]；Chowdhuri[5]于 1966年提出較為詳細(xì)的解析計算式；Darveniza[5]于1971年提出考慮避雷線作用的近似計算公式。上述方法采用較多近似計算并不同程度忽略了部分影響因素（如大地電參數(shù)、回?fù)綦娏鞑ㄐ蔚龋?，計算精度無法保證。另外解析法計算對象多為線路感應(yīng)電壓峰值，無法得到波形，不能滿足越來越復(fù)雜的防雷工作。數(shù)值方法既能更多的還原物理原型又能充分考慮各影響因素，是計算雷電感應(yīng)過電壓的有力工具。而雷電電磁脈沖（LEMP）的計算是數(shù)值方法中不可或缺的環(huán)節(jié)，關(guān)于 LEMP的理論計算除回?fù)裟Ｐ蜕杏胁蛔阃?，傳統(tǒng)的近似計算方法也存在一些凾待解決的問題，如：未考慮大地電導(dǎo)率[6]（將大地作為良導(dǎo)體）；局限于大地電導(dǎo)率較高情況[7]；局限于距回?fù)敉ǖ垒^遠(yuǎn)（如大于100m）的場[8,9]；需要對近地面水平電場進(jìn)行修正[10]等。Yee[11]提出的時域有限差分（FDTD）法在計算LEMP方面不僅能考慮大地電參數(shù)的影響，而且在粗糙、不均勻反射面的計算環(huán)境中有著獨特的優(yōu)勢。

針對上述問題，本文采用數(shù)值方法模擬架空線雷電感應(yīng)過電壓，完善其從產(chǎn)生到發(fā)展的理論研究。重新推導(dǎo)了文獻(xiàn)[12]計算架空線雷電感應(yīng)過電壓的Agrawal耦合模型，提出一種基于多階 FDTD的數(shù)值計算方法。該方法采用一階FDTD計算LEMP，采用二階 FDTD[13]求解耦合電路方程，通過合理設(shè)置線路差分點實現(xiàn)兩過程有效結(jié)合。所編軟件適用于各種工況下架空線雷電感應(yīng)過電壓的計算。通過與火箭引雷試驗實測數(shù)據(jù)比較，驗證了本文方法的準(zhǔn)確性。

2 物理模型

計算雷電回?fù)敉ǖ乐車腖EMP及其在架空線上激勵產(chǎn)生的過電壓，需要建立回?fù)敉ǖ?、大地和架空線路之間的整體化模型，如圖1。在該模型中，將雷電回?fù)敉ǖ雷鳛閆坐標(biāo)軸，X-Y平面表示為地面，地面落雷點為坐標(biāo)系原點。對模型做如下設(shè)置：

圖1 模型結(jié)構(gòu)Fig.1 Scheme of model structure model

3 架空線雷電感應(yīng)過電壓計算方法

3.1 雷電回?fù)敉ǖ离娏?/h3>

回?fù)敉ǖ纼?nèi)雷電流的計算包括兩個主要過程：1）、通道基電流i(0,t)的確定；2）、通道內(nèi)回?fù)裟Ｐ偷拇_定。目前業(yè)內(nèi)應(yīng)用最廣泛的基電流函數(shù)有雙指數(shù)函數(shù)[14]：

和 Heidler函數(shù)[15]：

國內(nèi)外學(xué)者對雷電回?fù)裟Ｐ瓦M(jìn)行了大量的研究和討論，本文采用了應(yīng)用最多的工程模型[16]：

其中，I(z’,t)為任意高度z’和任意時間t的通道電流；u(t)是單位函數(shù)，t≥z’/vf時取 1 否則取 0；I(0,t)是通道基電流函數(shù)；p(z’)是電流衰減因子；v’是電流波傳輸速度；vf是回?fù)羲俣取?/p>

3.2 多階FDTD法計算雷電感應(yīng)過電壓

目前，解決外部電磁場激勵下傳輸線感應(yīng)電壓的問題主要包括 Taylor法[17]、Agrawal法[18]和Rachidi法[19]等，三種模型所采用的激勵場源不同。本文針對模型需求選擇Agrawal法，提出基于多階FDTD法計算線路感應(yīng)電壓。

該方法基本步驟如下：

1）選擇回?fù)裟Ｐ筒⑤斎牖負(fù)綦娏鲄?shù)（基電流、通道高度和回?fù)羲俣鹊龋?、架空線路參數(shù)（導(dǎo)線數(shù)、線徑和空間距離等參數(shù)）以及時間步長和總時間長度等基礎(chǔ)數(shù)據(jù)。

2）在計算空間內(nèi)選取合適的單位長度劃分空間網(wǎng)格，采用一階FDTD法求解Maxwell方程計算回?fù)敉ǖ乐車臻g的LEMP。

3）采用二階FDTD法處理Agrawal耦合方程，得到圖2所示電路模型的離散式，將步驟2）求得的空間水平電場和垂直電場代入電路方程的離散式中，迭代求解線路差分點的電壓和電流。

圖2 Agrawal 耦合模型電路圖Fig.2 Circuit of Agrawal coupling model

4）輸出每個時間步線路上“觀測點”的電壓和電流值，得到該點在LEMP激勵下的時域響應(yīng)。

3.2.1 FDTD法計算空間LEMP

雷電回?fù)敉ǖ揽梢暈檩S對稱的天線模型，產(chǎn)生的LEMP也具有軸對稱性。因而三維場可簡化為二維場，即含對稱軸在內(nèi)的半個剖面，模型如圖1。軸對稱情況下的電磁場具有TE模和 TM模兩組獨立的解，分別包括 Eφ、Hr、Hz分量和 Hφ、Er、Ez分量。天線輻射模型采用TM模，對應(yīng)柱坐標(biāo)系下的 Maxwell方程為[20]：

采用柱坐標(biāo)Yee網(wǎng)格的方程差分格式為[20]：

式中，Er、Ez和Hφ分別表示LEMP中的水平電場、垂直電場和水平磁場。ε、σ和μ0分別為傳播介質(zhì)的電容率、電導(dǎo)率和磁導(dǎo)率。n、Δt、Δz和Δr分別為時間步、單位時間步長、垂直和水平方向空間網(wǎng)格步長。本文旨在實現(xiàn)有限區(qū)域內(nèi)模擬無限空間的LEMP分布，須在計算區(qū)域邊界設(shè)置吸收邊界條件。工程上應(yīng)用性能較好的吸收邊界主要有Berenger提出的完全匹配層（PML）[20]、Chen提出的修正的完全匹配層（MPML）[20]以及Mur提出的單行波吸收條件[20]。本文對軸線和另外三個空間吸收邊界分別采用安培環(huán)路定理和Mur吸收條件處理Maxwell方程，一階FDTD法較為成熟不做重點介紹，故邊界上具體離散格式不再贅述。

3.2.2 場線耦合方程的二階FDTD離散方法

單導(dǎo)線 Agrawal模型等值電路耦合電路方程為[18]：

Vs為導(dǎo)線散射電壓；ξg為大地阻抗的瞬態(tài)值，表示大地?fù)p耗的影響；符號?表示卷積積分；Ex(x,h,t)為差分點處沿導(dǎo)線方向水平電場值；L和C分別為導(dǎo)線單位長電感電容。對上式進(jìn)行向量擴(kuò)展，可得多導(dǎo)線模型的耦合矩陣方程[18]：

可見，耦合方程從單導(dǎo)線模型擴(kuò)展到多導(dǎo)線模型將各變量擴(kuò)展為變量矩陣。導(dǎo)線上總感應(yīng)電壓Vi(x,t)理解為散射電壓和入射電壓之和[18]：，其中入射電壓，h為第 i根導(dǎo)線對地i高度，Ez(x,z,t)為差分點處垂直電場。為便于理解將模型簡化處理，在兩端設(shè)置邊界條件如下：

R0和RL分別為線路始端和末端接地電阻。文獻(xiàn)[12,13]采用一階FDTD對式（10-11）離散處理，本文從提高數(shù)值穩(wěn)定性和計算精度的角度對式（10-11）采用二階FDTD法處理。分別在空間域x和時間域 t內(nèi)作變形處理，式（10）和式（11）可改寫為：

o(Δt3)為展開式余項。為方便使用FDTD處理方程，將待求變量和激勵源變量寫成空間步x和時間步t上更加直觀的離散格式：

其中，Δx指線路差分點間單位長度，k表示沿線空間步數(shù)。將式（18-19）中散射電壓、電流和水平電場對空間步長x和時間步長t的偏微分分別做一階、二階差分替換，即得式（18-19）的二階FDTD離散形式:

圖3 一階FDTD離散差分圖Fig.3 Scheme of first order FDTD discrete difference

對于一階 FDTD差分和二階 FDTD差分的區(qū)別，可通過圖3和圖4得到清晰直觀的理解：

圖4 二階FDTD離散差分圖Fig.4 Scheme of second order FDTD discrete difference

由圖3可見一階FDTD法電壓計算與電流計算非同步進(jìn)行，差分格式中散射電壓節(jié)點和電流節(jié)點交替分布（理解為取節(jié)點電壓和兩節(jié)點之間電流為求解變量），若離散電壓節(jié)點共有kmax個，則離散電流節(jié)點共有kmax-1個；圖4的二階FDTD差分格式中離散電壓節(jié)點和電流節(jié)點重合為同一點，該離散格式的數(shù)值穩(wěn)定性更高且能很好的簡化線路不連續(xù)點（如避雷線接地點）的處理。按照上述過程對式（12-13）進(jìn)行相同方法的推導(dǎo)和差分代換，得到考慮大地?fù)p耗的多導(dǎo)線模型二階FDTD離散方程：

3.2.3 不連續(xù)點處理

線路中每級桿塔的接地位置和有裝設(shè)避雷器的位置在計算中作為不連續(xù)點對待，該節(jié)點電壓可采用圖5模型單獨計算：節(jié)點k的相鄰兩節(jié)點（k-2、k-1、k+1、k+2）在第 n+1時間步的散射電壓和電流以及阻抗函數(shù)Γ和該處垂直電場 Ez為已知量。節(jié)點k的散射電壓表達(dá)式為：其中，為接地點處分路阻抗中的電流，計算式為：；Γ為的函數(shù)表示該阻抗上的壓降，分路阻抗若為線性電阻，則Γ =RgIg，可得：。對多導(dǎo)線情況可參照上文對單導(dǎo)線的擴(kuò)展方法：

圖5 線路不連續(xù)點模型圖Fig.5 Scheme of discontinuity point model

若導(dǎo)線數(shù)為N，則矩陣[Γ]為：

4 計算方法驗證

4.1 LEMP的計算與實測結(jié)果對比

LEMP的準(zhǔn)確計算是研究線路雷電感應(yīng)過電壓的基礎(chǔ)。為檢驗本文FDTD法計算LEMP，對以下火箭引雷試驗進(jìn)行計算對比。

佛羅里達(dá)州國際雷電研究測試中心（ICLRT）在2000年7月11日人工引雷方法測得Flash S0022,stroke 1[21]。實測15m處地面垂直電場如圖6a。按照實測雷電流波形采用一個 Heidler函數(shù)與一個雙指數(shù)函數(shù)疊加擬合基電流，F(xiàn)DTD法計算電場如圖6b。

圖6 實測電場與計算電場對比Fig.6 The measured electric field and calculation

由以上算例可見，本文FDTD法計算電場結(jié)果與實測電場基本吻合，證明LEMP算法正確性。

4.2 雷電感應(yīng)過電壓的計算與實測結(jié)果對比

1 數(shù)值計算與火箭引雷試驗對比

2008年8月12日在廣東火箭引雷試驗中對220V架空民用線感應(yīng)電壓進(jìn)行了觀測。試驗布置如下：兩根水平架空導(dǎo)線長 1 500m，高 3m，引流桿距導(dǎo)線垂直距離約 40m且距導(dǎo)線一端（A端）約為1 250m，觀測點距 A端約 1 300m。本文采用兩個Heidler函數(shù)疊加擬合首次回?fù)艋娏?，如圖7。采用本文方法對觀測點感應(yīng)電壓進(jìn)行仿真計算,圖8為本文仿真計算感應(yīng)電壓與實測結(jié)果相對比。

圖7 回?fù)敉ǖ赖撞侩娏鱂ig.7 Base current waveform of the return channel

圖8 觀測點感應(yīng)過電壓Fig.8 Induced over-voltage of observation point

2 數(shù)值計算與自然閃電觀測結(jié)果對比

2008年7月31日在上述觀測系統(tǒng)中測得了由自然閃電引起的感應(yīng)電壓。在記錄的1秒鐘內(nèi)共有7次回?fù)暨^程，本文對第三次回?fù)舻挠^測點電壓波形進(jìn)行數(shù)值模擬，圖9為實測電壓波形和數(shù)值計算電壓波形對比結(jié)果。

圖9 觀測點感應(yīng)過電壓Fig.9 Induced over-voltage of observation point

由以上兩算例可見，本文多階FDTD法計算觀測點雷電感應(yīng)過電壓結(jié)果與實測數(shù)據(jù)吻合較好，證明該計算方法的正確性。

5 計算方法應(yīng)用

本文多階FDTD法可計算不同導(dǎo)線數(shù)在多工況下的雷電感應(yīng)過電壓。為檢驗該方法的計算性能，對雷電流部分參數(shù)和兩種導(dǎo)線排布方式的影響進(jìn)行計算分析。

5.1 雷電流幅值Im的影響

以圖1所示簡單結(jié)構(gòu)為計算基礎(chǔ)。模型參數(shù)設(shè)置為：單架空線長600m高3m，雷擊點距線路垂直距離 50m且與線路兩端距離相等，觀測點距一端400m。采用 Heidler函數(shù)擬合基電流，采用 MTLL回?fù)裟Ｐ?，大地電?dǎo)率σ取5×10-3S/m。

基電流波形及回?fù)魠?shù)相同，基電流幅值Im取值分別為11kA、16kA和21kA，所對應(yīng)的觀測點電壓Um1、Um2和Um3的波形如圖10所示。

圖10 基電流幅值的影響Fig.10 The influence of base current amplitude

由圖10可知，雷電流幅值 Im對線路感應(yīng)電壓的影響主要體現(xiàn)在幅值變化，Im從 11kA增大到16kA和21kA，提高比例分別為45%和91%，對應(yīng)觀測點感應(yīng)電壓幅值從20kV提高到32kV和45kV，提高比例分別為60%和125%?？梢奍m的提高明顯引起了感應(yīng)電壓幅值提高，該過程對電壓波形影響很小。

5.2 雷電流波頭陡度的影響

基電流幅值及回?fù)魠?shù)相同，其他輸入設(shè)置同5.1，波頭時間分別為 0.1μs、0.4μs和 1μs，所對應(yīng)的觀測點感應(yīng)電壓如圖11。

圖11 基電流波頭陡度的影響Fig.11 The influence of base current wave steepness

由圖11可知，雷電流波頭陡度降低，波頭時間從0.1μs增大到0.4μs和1μs，對應(yīng)的感應(yīng)電壓幅值從 39kV降低到38kV和 35kV，降低比例為 3%和10%，波頭時間從 1.8μs增大到 2.0μs和 2.5μs，分別增大11%和39%?？梢娎纂娏鞑^陡度降低，引起了感應(yīng)電壓波頭時間較明顯增大，陡度降低，而幅值則降低很小。

5.3 導(dǎo)線排布方式影響

基電流和回?fù)魠?shù)相同，三根架空導(dǎo)線水平布置(與雷擊點距離分別為 50m、51m 和 52m，高度3m)和垂直布置（與雷擊點距離 50m，高度分別為3m、4m和5m）時各導(dǎo)線觀測點電壓波形如圖12。

由圖12可知，水平排布的導(dǎo)線觀測點感應(yīng)電壓幅值均為30kV左右，相差不足3%；而垂直排布的導(dǎo)線觀測點感應(yīng)電壓則相差明顯，由下而上分別約為30kV、45kV和60kV，中間和最高導(dǎo)線比最低導(dǎo)線分別提高了 50%和 100%?？梢妰煞N排布方式的“屏蔽效應(yīng)”差異大，垂直排布的導(dǎo)線感應(yīng)電壓幅值差明顯高于水平排布。兩種排布方式對感應(yīng)電壓波形影響不大。

圖12 導(dǎo)線排布方式的影響Fig.12 The influence of conductor configuration

6 結(jié)論

提出一種多階FDTD的架空線雷電感應(yīng)過電壓計算方法，該方法采用一階FDTD計算LEMP，采用二階 FDTD推導(dǎo)場線耦合方程，比單純的一階FDTD法數(shù)值穩(wěn)定性更高，得到結(jié)論如下：

1）與火箭引雷試驗和自然閃電觀測數(shù)據(jù)對比驗證了所提方法的正確性和實用性。該方法可實現(xiàn)不同導(dǎo)線數(shù)在各種工況下雷電感應(yīng)過電壓的準(zhǔn)確計算，且能更簡捷的處理線路上的不連續(xù)點。

2）通過計算得到：雷電流幅值對線路感應(yīng)電壓幅值影響明顯，對電壓波形影響小；雷電流波頭陡度對感應(yīng)電壓幅值影響較小，但對波形有較大影響。多導(dǎo)線之間有“屏蔽效應(yīng)”，該效應(yīng)受導(dǎo)線排布方式影響明顯，垂直排布的導(dǎo)線感應(yīng)電壓差較水平排布明顯要大。

[1] 李瑞芳,吳廣寧,曹曉斌,等.雷電流幅值概率計算公式[J].電工技術(shù)學(xué)報,2011,26(4): 161-167.Li Ruifang,Wu Guangning,Cao Xiaobin,et al.Formula for probability of lightning current amplitude[J].Transctions of China Electrotechnical Society,2011,26(4): 161-167.

[2] 熊小伏,方偉陽,程韌俐,等.基于實時雷擊信息的輸電線強(qiáng)送決策方法[J].電力系統(tǒng)保護(hù)與控制,2013,41(19): 7-11.Xiong Xiaofu,Fang Weiyang,Cheng Renli,et al.Forced power supply decision of transmission lines based on real-time lightning information[J].Power System Protection and Control,2013,41(19): 7-11.

[3] 李瑞芳,吳廣寧,曹曉斌,等.輸電線路雷電繞擊率的三維計算方法[J].電工技術(shù)學(xué)報,2009,24(10):134-138.Li Ruifang,Wu Guangning,Cao Xiaobin,et al.Three-Dimensional Calculation Method on Shielding Failure Rate of Transmission Lines[J].Transctions of China Electrotechnical Society,2009,24(10): 134-138.

[4] 谷定燮.500kV輸變電工程設(shè)計中雷電過電壓問題[J].高電壓技術(shù),2000,26(6): 60-62.Gu Dingxie.The Research on Lightning Overvoltages of 500kV Transmission Engineering Desgn[J].High Voltage Engineering,2000,26(6): 60-62.

[5] 周詩健,孫景群譯.雷電/上下卷[M].北京: 水利電力出版社,1982.R.H.Golde.Lightning[M].Zhou Shijian,Sun Jingqun translated.Water Resources and Electric Power Press,1982.

[6] 馮雷,周佩白.高層建筑遭受雷擊時感應(yīng)電場的計算[J].高電壓技術(shù),2001,27(1): 52-54.Feng Lei,Zhou Peibai.Computation of the Induced Electric Field in a High Building[J].High Voltage Engineering,2001,27(1): 52-54.

[7] C.A.Nucci,F.Rachidi,M.V.Ianoz,et al.Lightning induced voltages on overhead lines.IEEE Trans.Electromag.Compat,1993,35(1): 75-83.

[8] V.Cooray.Effects of propagation on the return stroke radiation fields.Radio Sci,1987,22: 757-768.

[9] V.Cooray.Some consideration on the “Cooray-Rubinstein” formulation used in deriving the horizontal electric field of lightning return strokes over finitely conducting ground.IEEE Trans.Electromg.Compat.2002,44(4): 560-565.

[10] 邊凱,陳維江,李成榕,等.架空配電線路雷電感應(yīng)過電壓計算研究[J].中國電機(jī)工程學(xué)報,2012,32(31): 191-199.Bian Kai,Chen Weijiang,Li Chengrong,et al.Calculation of lightning induced overvoltage on overhead distribution lines[J].Proceedings of the CSEE,2012,32(31): 191-199.

[11] K.S.Yee.Numerical solution of initial boundary value problems involving Maxwell equations in isotropic media.IEEE Trans.Antennas Propagat.1966,14(3): 302~307.

[12] He-Ming Ren,Bi-Hua Zhou,Vladimir A.Rakov et al.Analysis of Lightning-Induced Voltages on Overhead Lines Using a 2-D FDTD Method and Agrawal Coupling Model[J].IEEE Trans on Electromagnetic Compativility,2008,50(3): 651-659.

[13] M.Paolone.Lightning Electromagnetic Field Coupling to Overhead Lines: Theory,Numerical Simulations,and Experimental Validation [J].IEEE Transactions on Electromagnetic Compatibility,2009,51(3): 532-547.

[14] R.H.Golde.Lightning surges on overhead distribution lines caused by indirect and direct lightning strokes[J].Trans.Amer.Inst.Elec.Engrs.,1954,73,437-447.

[15] Heidler F,Cvetic J M,Stanic B V.Calculation of lightning current parameters[J].IEEE Transactions on Power Delivery,1999,14(2): 399-404.

[16] V.Cooray.On the concepts used in return stroke models applied in engineering practice.IEEE Transactions on Electromagnetic Compatibility,2003,45(1):101-108.

[17] C.D.Taylor,R.S.Satterwhite,C.W.Harrison,et al.The response of a terminated two-wire transmission line excited by a nonuniform elecromagnetic field[J].IEEE Trans.Electromagn.Compat,1965,13(6):987-989.

[18] A.K.Agrawal,H.J.Price,S.H.Gurbaxani.Transient response of multiconductor transmission lines excited by a nonuniform electromagnetic field[J].IEEE Transactions on Electromagnetic Compatibility,1980,22(2):119-129.

[19] F.Rachidi.Formulation of field-to-transmission line coupling equations in terms of magnetic excitation Field[J].IEEE Transactions on Electromagnetic Compatibility,1993,35(3): 404-407.

[20] 葛德彪,閆玉波.電磁波時域有限差分法(第三版)[M].西安: 西安電子科技大學(xué)出版社,2011:29-31.Ge Debiao,Yan Yubo.Finite-difference time-domain method for electromagnetic waves[M].Xian: xidian university press,2011: 29-31.

[21] M.Miki,V.A.Rakov.et al.Electric fields near triggered lightning channels measured with Pockels sensors.J.Geophys.Res.,2002,107: 10.1029/2001JD001087.