王淑敏，陳曉寧，張海濤

(中國人民解放軍陸軍工程大學(xué) 國防工程學(xué)院，江蘇 南京 210007)

0 引言

雷電是一種常見的自然物理現(xiàn)象，伴隨著聲、光同時發(fā)生，全世界每天發(fā)生800萬次雷暴，也就是說每秒有100次左右的可能性發(fā)生雷擊飛機(jī)的事故[1]，雷擊飛機(jī)一旦發(fā)生，大電流、高電壓、強(qiáng)瞬態(tài)電磁場會對飛機(jī)造成嚴(yán)重的破壞，輕則設(shè)備儀器損毀，重則人員傷亡，因此，雷電對飛機(jī)的飛行安全存在著潛在的巨大威脅。

雷擊飛機(jī)時，雷電流從機(jī)體一個點流入從另一個點流出，這些附著位置稱為雷電附著點。對飛機(jī)雷電附著點的研究經(jīng)歷了如下幾個階段：20世紀(jì)末，對雷電附著的研究主要來自于飛機(jī)穿越雷暴區(qū)所得到的數(shù)據(jù)[2-3]；后來人們在實驗室進(jìn)行高壓放電試驗，對雷擊飛機(jī)過程進(jìn)行模擬[4-5]，雖然實驗室雷電附著點試驗?zāi)軌蜉^為真實地、直觀地模擬雷擊飛機(jī)的過程，但試驗對沖擊電壓發(fā)生器、高性能攝像機(jī)等設(shè)備的要求較高且危險性較大；隨著計算機(jī)技術(shù)的發(fā)展，各類電磁仿真軟件涌入了商業(yè)市場，得到廣泛地應(yīng)用，以仿真代替試驗[6-7]，可以不再受制于時間、空間、試驗條件的限制，但對計算機(jī)硬件配置的要求較高，且仿真參數(shù)的設(shè)置需要調(diào)試，不合適的參數(shù)會直接降低仿真結(jié)果準(zhǔn)確性；目前也有諸多學(xué)者運用各自編程的先導(dǎo)發(fā)展模型(Leader Progression Model,LPM)來模擬雷擊直升機(jī)過程[8-9]，但模型與實際雷電仍然存在一定的區(qū)別?？紤]到上述問題，本文提出了一種基于分形理論的三維先導(dǎo)發(fā)展模型，研究了如何能夠更加準(zhǔn)確地對雷擊主通道以及分支通道進(jìn)行描述，并將模型運用到某型直升機(jī)上計算了關(guān)鍵位置雷電附著頻率。

1 雷電先導(dǎo)發(fā)展模型的建立

LPM模型的建立基于長期以來人們對自然雷電的觀測和長間隙放電的發(fā)展，是對雷電發(fā)展過程的物理描述。本文利用基于分形理論的LPM模擬雷擊過程，主要包括兩個方面：負(fù)極性下行先導(dǎo)發(fā)展模型的建立、正極性上行先導(dǎo)發(fā)展模型的建立。

1.1 負(fù)極性下行先導(dǎo)發(fā)展模型

為了使負(fù)極性下行先導(dǎo)發(fā)展模型能夠更加貼近自然雷電的發(fā)展過程，需要對下行先導(dǎo)中的電荷分布以及發(fā)展角度進(jìn)行定義以使模型具備雷電的物理意義，并且運用了發(fā)展點概率選擇方法以使模型具備雷電的隨機(jī)性特征。整個模型的發(fā)展基于步進(jìn)的思想，即不斷地重復(fù)從已發(fā)展點向待發(fā)展點的步進(jìn)過程。

1.1.1下行先導(dǎo)發(fā)展模型的電荷分布

Cooray模型認(rèn)為在先導(dǎo)的發(fā)展過程中，電荷量主要集中于通道頭部位置，并按照近似線性的分布向云內(nèi)始發(fā)位置遞減[10]，因此本文利用線性分布的線電荷模型來描述先導(dǎo)通道,利用點電荷模型來描述先導(dǎo)頭部。

雷電通道可以大體分為主通道和分支通道。在負(fù)極性下行先導(dǎo)主通道內(nèi)，定義海拔高度為h位置處的電荷分布如式(1)所示。由于當(dāng)先導(dǎo)起始時，在云層中心位置處的先導(dǎo)通道電荷密度可近似為0，故b取0。另外，通過不同回?fù)綦娏鞣逯迪码姾擅芏扰c海拔高度的關(guān)系[8]可知，系數(shù)K與回?fù)綦娏鞣逯档年P(guān)系如式(2)所示。

ρ(h)=K(Hc-h)+b

(1)

K×Hc=-1.5×10-5×Ip

(2)

式中，ρ(h)為通道內(nèi)電荷密度，C/m；Hc為雷雨云高度，m；h為先導(dǎo)頭部海拔高度值，m；K、b為系數(shù)；Ip為回?fù)綦娏鞣逯?，A。

利用點電荷來模擬先導(dǎo)頭部，并將單位積分長度內(nèi)先導(dǎo)通道看成帶電直線，則在負(fù)極性下行先導(dǎo)頭部流注區(qū)域內(nèi)，運用高斯定理積分可求解得下行先導(dǎo)頭部電荷表達(dá)式為：

(3)

式中，ε為電介質(zhì)的介電常數(shù)，F(xiàn)/m；Q為積分區(qū)域內(nèi)電荷總量；h為先導(dǎo)頭部海拔高度值，m；Es為流注區(qū)平均場強(qiáng)。

對于先導(dǎo)的分支通道，考慮到其與主通道具有相似的電荷密度分布規(guī)律，因此可以采用在主通道電荷分布的基礎(chǔ)上加入修正系數(shù)C來描述，如式(4)、(5)所示，其中C取0.312 5[8]。

(4)

(5)

1.1.2基于角度約束的下行先導(dǎo)待發(fā)展點選擇方法

為了確保下行先導(dǎo)的發(fā)展方向始終朝向地面，需要進(jìn)一步明確角度約束條件，分析文獻(xiàn)[8]中提到的兩次雷電觀測中下行先導(dǎo)分支通道角度提取情況：發(fā)生于2008年7月22日，編號CG20080722；發(fā)生于2011年7月31日，編號20110731，由于雷電的主通道始終約與地面垂直，暫時不予考慮。對兩次雷擊中各分支通道方向與鉛錘方向的夾角度數(shù)進(jìn)行統(tǒng)計，如圖1所示，可知小于60°的通道數(shù)約占總通道數(shù)的4/5，因此本文定義角度約束條件為：先導(dǎo)已發(fā)展點和先導(dǎo)待發(fā)展點的連線方向與鉛錘方向的夾角度數(shù)范圍為0°～120°，并利用上述角度約束條件篩選出待發(fā)展點的坐標(biāo)。

圖1 雷電先導(dǎo)通道夾角統(tǒng)計圖

1.1.3基于概率選擇的下行先導(dǎo)發(fā)展點確定方法

(6)

(7)

鑒于流注向先導(dǎo)轉(zhuǎn)化過程與大氣環(huán)境有關(guān)，在不同的壓強(qiáng)、溫度、濕度等大氣參數(shù)下Eth應(yīng)該取不同的值，通過將大氣參數(shù)轉(zhuǎn)化為包含海拔高度的物理量，得到臨界擊穿電場Eth與先導(dǎo)頭部海拔高度h的關(guān)系：

(8)

式中，Eth0為負(fù)極性下行先導(dǎo)在標(biāo)準(zhǔn)大氣環(huán)境下的電場擊穿閾值，取750 kV/m[13]。

另外，為了確定最終的發(fā)展點，運用了遺傳算法(Genetic Algorithm,GA)中的輪盤賭法進(jìn)行非確定性搜索，保證了發(fā)展概率大的待發(fā)展點被選擇的幾率大，在滿足模型物理意義的基礎(chǔ)上使先導(dǎo)模型的發(fā)展具有了隨機(jī)性，即每一次計算得到的雷擊路徑都不相同。

1.1.4負(fù)極性下行先導(dǎo)發(fā)展模型的大氣電學(xué)驗證

分形維數(shù)(fractal dimension)的概念來源于幾何學(xué)中的分形理論，反映了物體或分形體占有空間的有效及自相似性程度。大量觀測結(jié)果表明，自然雷電的分形維數(shù)D大體上介于1.1～1.4之間[14]，為了評估本文建立的負(fù)極性下行先導(dǎo)隨機(jī)發(fā)展模型的準(zhǔn)確性，本文利用計盒維數(shù)法對先導(dǎo)的分形維數(shù)進(jìn)行了檢驗，具體方法為：在三維xoy坐標(biāo)系下，將模擬的雷電下行先導(dǎo)分別投影到兩個相互垂直的平面上，運用計盒維數(shù)法對兩個二維圖形進(jìn)行分形維數(shù)計算，后取平均值得到三維先導(dǎo)模型的分形維數(shù)D。

運用上文的建模方法模擬負(fù)極性下行先導(dǎo)的對地放電過程，如圖2所示。經(jīng)過大量的統(tǒng)計計算，得到本文建立的負(fù)極性下行先導(dǎo)發(fā)展模型分形維數(shù)分布在1.2～1.37之間，如圖3所示，由此說明了模型能夠較好地模擬自然雷電的發(fā)展過程。

圖2 負(fù)極性下行先導(dǎo)對地放電

圖3 先導(dǎo)模型分形維數(shù)

1.2 正極性上行先導(dǎo)發(fā)展模型

Becerra、Cooray模型是通過不斷迭代計算先導(dǎo)前方幾米范圍內(nèi)的電勢分布來計算某位置處是否能成為雷擊附著點，被廣泛應(yīng)用于模擬雷電上行先導(dǎo)的起始和發(fā)展過程[15]。將該模型應(yīng)用到本雷擊直升機(jī)附著點的判斷中，具體實現(xiàn)過程如下：

初始電暈起始后，認(rèn)為電暈區(qū)域幾米范圍內(nèi)電場強(qiáng)度恒定，直升機(jī)雷擊關(guān)鍵位置處附近的空間電勢分布U發(fā)生畸變，可以用直線來近似表示為:

U(0)=kUz+U0

(9)

式中，上標(biāo)“(i)”表示迭代到第i步，i=0,1,2,…；kU為直線的斜率，kV/m；U0為直線的截距，kV；z表示距離，m。

則此時初始電暈區(qū)域的電荷量ΔQ以及電暈區(qū)的頭部長度ls可以表示為：

(10)

(11)

式中，Es為正極性流注區(qū)平均場強(qiáng)，取4.5×105kV/m；KQ為系數(shù)，表征電暈區(qū)電荷總量與電勢分布畸變關(guān)系，取4×10-11C/(V·m)。

當(dāng)初始電暈區(qū)域的電荷量積累到了臨界值1 μC時，判定為不穩(wěn)定上行先導(dǎo)開始發(fā)展，此時先導(dǎo)頭部的電勢表示為：

(12)

式中，lL(i)為發(fā)展到第i步時上行先導(dǎo)長度；E∞為穩(wěn)態(tài)時先導(dǎo)通道的平均場強(qiáng)；x0為初始狀態(tài)時的先導(dǎo)長度。

不穩(wěn)定上行先導(dǎo)頭部流注放電，計算前方流注區(qū)的頭部位置和流注區(qū)電荷量：

(13)

(14)

流注向先導(dǎo)轉(zhuǎn)化，不穩(wěn)定上行先導(dǎo)繼續(xù)發(fā)展，計算不穩(wěn)定上行先導(dǎo)先導(dǎo)發(fā)展長度如式(15)所示，若達(dá)到臨界值2 m時，則判定不穩(wěn)定上行先導(dǎo)轉(zhuǎn)化為了穩(wěn)定上行先導(dǎo)，該先導(dǎo)起始位置即為雷電附著點。

(15)

式中，qL為流注向先導(dǎo)轉(zhuǎn)化時每單位長度的電荷量。

2 關(guān)鍵位置附著情況計算方法

雷擊飛機(jī)的過程可以分為兩種方式：自然雷電附著的情形和飛機(jī)觸發(fā)雷電的情形，其中，相比于飛機(jī)觸發(fā)雷電來說，自然雷電附著時的雷電流幅值高，對飛機(jī)造成的危害更大，因此，本文對自然雷電附著的物理過程進(jìn)行仿真，并根據(jù)標(biāo)準(zhǔn)[16]要求，按照棒-板電極的試驗布局方法對先導(dǎo)起始位置和仿真空間進(jìn)行設(shè)定。具體方法為：在以直升機(jī)為球心的球面上設(shè)置放電電極位置，并以30°為增量將球面劃分為72個位置，但由于直升機(jī)結(jié)構(gòu)具有左右對稱性，故可將其簡化為37個位置點，如圖4所示，其中，本文的電極位置即下行先導(dǎo)起始位置。另外，為了滿足電極位置與直升機(jī)模型中心之間的距離為最大機(jī)身尺寸的3倍，球半徑取45 m。

圖4 放電起始位置方位圖

本文試驗對象某型直升機(jī)，是一種中型通用直升機(jī)。按照真機(jī)尺寸進(jìn)行三維等比例建模，整個模型主要包括了前設(shè)備艙、玻璃、駕駛艙、座艙、過渡段、動力艙、尾梁艙、尾旋翼和起落架等。

典型的直升機(jī)雷電附著點位置有主旋翼翼尖、尾旋翼翼尖、機(jī)頭[16]，這幾處是進(jìn)行直升機(jī)雷電防護(hù)設(shè)計的關(guān)鍵位置?；谏衔牡哪Ｐ徒⒎椒?，在MATLAB中模擬雷擊直升機(jī)過程，在不同負(fù)極性下行先導(dǎo)起始位置(電極位置)的情況下，計算這些關(guān)鍵位置成為雷電附著點的頻率，如圖5所示。

3 計算結(jié)果與分析

依據(jù)上文的方法模擬雷擊直升機(jī)過程，研究直升機(jī)主旋翼、尾旋翼和機(jī)頭三個關(guān)鍵位置處的雷電附著情況，在圖4中的每個電極位置下各仿真10次，共1 110次，附著次數(shù)結(jié)果列于表1，另分別繪制了同緯度、同經(jīng)度電極位置雷電附著次數(shù)統(tǒng)計圖，如圖6所示。

圖5 關(guān)鍵位置雷電附著判斷流程圖

分析表1和圖6可知：

(1)在37個電極位置的放電試驗中，主旋翼翼尖、尾旋翼翼尖和機(jī)頭分別被雷電附著182、158、129次，共469次，主旋翼、尾旋翼和機(jī)頭分別占49.2%、42.7%、 34.9%，表明三個關(guān)鍵位置被雷電附著的可能性最大的為主旋翼，其次是尾旋翼，可能性最小的是機(jī)頭。

(2)當(dāng)雷電起始于11號電極位置時，擊中三個關(guān)鍵位置的概率最大，達(dá)到了66.7%，其次是當(dāng)雷電起始于10號電極位置時，擊中概率為60%。對于主旋翼來說，當(dāng)雷電起始于11號電極位置時，其被附著的可能性最大；對于尾旋翼來說，當(dāng)雷電起始于9、10號電極位置時，其被附著的可能性最大；對于機(jī)頭來說，當(dāng)雷電起始于29號電極位置時，其被附著的可能性最大。

(3)對于圖6(a)，當(dāng)電極位置處于同一緯度上時，隨著經(jīng)度的增加，三個關(guān)鍵位置處的雷電附著次數(shù)規(guī)律為：主旋翼附著次數(shù)呈現(xiàn)中間大兩邊小的趨勢；尾旋翼附著次數(shù)基本呈現(xiàn)減小的趨勢，除了當(dāng)雷電起始于21號電極位置時，附著次數(shù)為極大值；機(jī)頭附著次數(shù)呈現(xiàn)增大的趨勢。分析可知，基本符合關(guān)鍵點位置與電極位置距離越近，附著次數(shù)越大的規(guī)律。

(4)對于圖6(b)，在相同經(jīng)度下，隨著電極位置緯度的增加，雷電在三個關(guān)鍵點處的附著次數(shù)呈現(xiàn)先增加后減少，在90°緯線的時候為極小值，后又先增加再減少，大體呈現(xiàn)對稱趨勢。其中，當(dāng)電極位置處于同一經(jīng)度上時，機(jī)頭位置處的雷電附著次數(shù)具有非常明顯的對稱性，且隨著經(jīng)度的增加而附著次數(shù)逐漸增多。在圖中，三條附著次數(shù)曲線相交較少，隨著經(jīng)度的增加，附著次數(shù)最大的依次為尾旋翼、尾旋翼、主旋翼、主旋翼、主旋翼、機(jī)頭、機(jī)頭。

4 結(jié)論

本文通過討論雷電先導(dǎo)發(fā)展模型主通道、分支通道的電荷分布以及發(fā)展規(guī)律，建立了具備自然雷電物理意義與隨機(jī)性特征先導(dǎo)的三維雷電先導(dǎo)發(fā)展模型。為了驗證模型的有效性，運用計盒維數(shù)法對分形維數(shù)進(jìn)行了計算，計算結(jié)果均分布在1.2～1.37之間，說明模型符合自然雷電的幾何特征。另外，依據(jù)相關(guān)雷擊附著點試驗標(biāo)準(zhǔn)，將模型運用到某型直升機(jī)上研究了關(guān)鍵位置雷電附著情況，結(jié)果表明直升機(jī)三個關(guān)鍵位置處的雷電附著概率由大到小依次為主旋翼(49.2%)、 尾旋翼(42.7%)、機(jī)頭(34.9%)，綜上所述，本文提出的雷電先導(dǎo)發(fā)展模型能夠為直升機(jī)雷電附著點的判斷提供一種快速的數(shù)值計算方法。

表1 直升機(jī)關(guān)鍵位置處雷電附著情況

圖6 關(guān)鍵位置處雷電附著次數(shù)統(tǒng)計圖