空間分辨率對上行先導(dǎo)起始過程數(shù)值模擬的影響

裴曉芳，丁潔，郭秀峰,,4*，紀梓煜，陳佳祥，張超凡，范陽

(1.南京信息工程大學(xué)電子與信息工程學(xué)院，南京 210044；2.無錫學(xué)院大氣與遙感學(xué)院，無錫 214105；3.南京信息工程大學(xué)江蘇省大氣環(huán)境與裝備技術(shù)協(xié)同創(chuàng)新中心，南京 210044；4.中國氣象科學(xué)研究院災(zāi)害天氣國家重點實驗室，北京 100081)

隨著社會的發(fā)展，由雷電造成的人身安全問題及經(jīng)濟財產(chǎn)損失越來越多，對此合理的規(guī)避由雷電造成的損失顯得尤為重要[1-5]。在雷暴期間，建筑物上能否始發(fā)上行先導(dǎo)是判斷其是否能被閃電擊中的一個標準。因此，無論在雷電物理研究還是雷電防護領(lǐng)域，物體上如何發(fā)生上行先導(dǎo)是一個關(guān)鍵的問題。隨著計算機技術(shù)的發(fā)展，目前對于上行先導(dǎo)的研究，除了對閃電的直接觀測和實驗室實驗?zāi)M之外，數(shù)值模擬成為當(dāng)下廣泛使用的一種研究手段。在數(shù)值模擬中，連續(xù)的空間被劃分成離散的網(wǎng)格，進而計算每一個網(wǎng)格上的物理量。如何選取合適的網(wǎng)格大小(也稱空間分辨率)是所有數(shù)值模擬中的一個重要問題。因此，在上行先導(dǎo)數(shù)值模擬研究中，選取合理的上行先導(dǎo)始發(fā)判據(jù)和合適的空間分辨率都是至關(guān)重要的。

基于20世紀70年代Les Renardieres Group(“狐貍”團隊)開展的一系列正極性長空氣間隙放電實驗，分析提出了上行先導(dǎo)的起始可分為電暈起始、正先導(dǎo)始發(fā)和傳播三個階段。基于此，學(xué)者們提出了不同的閃電物理概念模型，這些模型利用實驗室長火花的定量結(jié)果來預(yù)測下行和上行先導(dǎo)的發(fā)展[6-11]。近些年，也有學(xué)者提出了更簡便的3種先導(dǎo)起始判據(jù)[12-14]。分別是平均電場判據(jù)、臨界長度判據(jù)和初始流注動態(tài)臨界長度判據(jù)。這些判據(jù)主要是用宏觀參量代替了先導(dǎo)發(fā)展的復(fù)雜計算過程。但是在這些模型中，被廣泛沿用至今的上行正先導(dǎo)始發(fā)判據(jù)主要有以下4種：介質(zhì)擊穿模型、臨界半徑準則、臨界長度法和臨界電場法。

介質(zhì)擊穿模型(dielectric breakdown model,DBM)是Niemeyer等[15]和Wiesmann等[16]開發(fā)的用于模擬雷暴云中放電通道的宏觀性質(zhì)和地面物體尖端的放電通道。先導(dǎo)的產(chǎn)生和傳播判據(jù)是導(dǎo)體或尖端周圍的電場是否超過臨界電場閾值(Ecrit)。但是臨界電場閾值在不同的模型中是不同的，例如，在云中臨界電場閾值Ecrit=100 kV/m[17]、Ecrit=150 kV/m[18]，以及在物體尖端臨界電場值Ecrit=500 kV/m[10]。此外，這個模型不能夠模擬微觀物理擊穿過程。

臨界半徑法與DBM相似，當(dāng)尖端電場等于電暈初始電場(通常Ecrit=3 000 kV/m)時，就可以達到穩(wěn)定的先導(dǎo)初始條件。這一標準是基于文獻[19]的實驗結(jié)果，并在文獻[6-8]的模型中得到應(yīng)用。由于這個概念是通過高壓實驗室長間隙放電實驗得到的，所以是否適用于自然雷電在更大維度范圍內(nèi)產(chǎn)生的上行先導(dǎo)還不清楚[20]。

臨界長度法是由Petrov等[10]提出的根據(jù)上行先導(dǎo)流光區(qū)電場強度以及長度提出的判據(jù)。該方法認為當(dāng)流光區(qū)長度發(fā)展達到了0.7 m，并且同時這個范圍里的電場強度超過了正極性流光區(qū)的平均電場值(約為500 kV/m)，則流光帶轉(zhuǎn)化為先導(dǎo)。

然而，以上3個標準對實驗結(jié)果無法體現(xiàn)先導(dǎo)發(fā)展的具體物理過程。因此，Becerra等[21-23]開發(fā)了一個自持式上行先導(dǎo)始發(fā)模型，該模型遵循了Goellian等[24]類似的方法。此標準中，當(dāng)流光區(qū)能產(chǎn)生大于或等于1 μC電暈電荷時[25-26]，方可形成不穩(wěn)定先導(dǎo)，而穩(wěn)定先導(dǎo)持續(xù)傳播的條件是不穩(wěn)定的先導(dǎo)長度達到2 m。盡管最后這個標準被認為是這4個標準里面最合理的，它仍舊不能應(yīng)用于所有的數(shù)值模擬中。例如，在模擬先導(dǎo)分支結(jié)構(gòu)或模擬區(qū)域太大而無法使用精細的空間分辨率時(當(dāng)最小網(wǎng)格的尺寸大于2 m)。因此，另外3種標準仍被廣泛使用。

但是，無論選擇哪種標準，首先都要計算避雷針尖端或物體尖端周圍的電勢或電場值。目前在大多數(shù)的先導(dǎo)模型中，主流模型采用的電位計算方法有有限差分法(finite difference method，F(xiàn)DM)、有限元法(finite element method，F(xiàn)EM)和電荷模擬法(charge simulation method，CSM)。由于可以模擬復(fù)雜邊界，F(xiàn)DM和FEM是兩種主要的方法。譚涌波等[27]在二維模型中采用FDM的方法，得出在相同建筑物形狀下尖端增強的電場實際上是相同的，但是在數(shù)值模擬中受網(wǎng)格大小的影響很大。網(wǎng)格尺寸越小，計算出的電場越大，尤其是在尖端處。因此，在上行先導(dǎo)模型中，網(wǎng)格的大小會影響對先導(dǎo)始發(fā)的判斷。然而以上4種判據(jù)都是基于實驗結(jié)果的，模擬中并沒有考慮網(wǎng)格大小對判據(jù)的影響。

此外，為了模擬雷暴云下的地面物體產(chǎn)生的上行先導(dǎo)，模擬域范圍至少要達到公里的量級，而建筑物接閃桿尖端通常僅有幾厘米。通常來說，空間分辨率越高模擬越精確，但實際上，特別是在三維(3D)閃電先導(dǎo)起始和傳播模型中，很難對大范圍(幾千米甚至幾十千米)進行非常小的網(wǎng)格(cm或mm)剖分。進而導(dǎo)致，學(xué)者們對數(shù)值模擬的定量計算結(jié)果一直持有懷疑態(tài)度。因此，如何選擇合適的網(wǎng)格大小是模擬上行先導(dǎo)始發(fā)的關(guān)鍵問題。換句話說，當(dāng)網(wǎng)格的最小尺寸固定時，應(yīng)該選擇怎樣的先導(dǎo)起始判據(jù)？但是現(xiàn)在對這一問題的研究討論較少。

為探討被廣泛沿用至今的4種先導(dǎo)起始判據(jù)，在不同空間分辨率模擬中的適用性，以提高上行先導(dǎo)數(shù)值模擬的可靠性。基于有限差分法，現(xiàn)建立一個計算建筑物周圍電場的3D變網(wǎng)格模型，分析網(wǎng)格大小與電場計算結(jié)果的關(guān)系，并取了3種不同建筑物高度，討論網(wǎng)格大小對4種先導(dǎo)起始判據(jù)的影響。最后，在此基礎(chǔ)上提出網(wǎng)格大小的選擇方法和不同始發(fā)判據(jù)在實際先導(dǎo)模擬中的適用性。

1 模型及計算方法介紹

1.1 模型建立

由于先導(dǎo)起始的判斷，均離不開精確的尖端電場計算。因此，為了更加精確地模擬建筑物周圍電場的分布特征，現(xiàn)基于有限差分法求解泊松方程，建立了三維變網(wǎng)格下的靜電場計算模型[28-32]。相對之前的研究，本文研究主要有兩點優(yōu)勢：三維模型克服了二維模型對實際三維空間模擬的不足；變網(wǎng)格技術(shù)克服了均勻網(wǎng)格對大模擬域中細小尖端的模擬不足。

使用變網(wǎng)格技術(shù)，在尖端處采用小網(wǎng)格，離尖端越遠網(wǎng)格按照一定比例逐漸增大，如圖1所示，以xy平面為例，尖端頂部范圍內(nèi)采用最小網(wǎng)格，后續(xù)按照比例增加。模擬域為地面上空500 m×500 m×500 m。僅考慮建筑物上尖端的特征而不考慮建筑物本身的幾何特征，因此將建筑物尖端模型簡化為一個接地的立方體形狀的理想導(dǎo)體，長、寬、高分別為20 cm、20 cm、200 m，位于模擬域的中間。此外，還建立了100 m和300 m高的建筑物模型。在雷暴云電場(Eb)下，近地面的電場可近似為均勻的垂直電場，同時不考慮空中自由電荷的存在。因此，在實際的電場計算時，建筑物尖端周圍的電場可通過求解已知邊界的拉普拉斯方程獲得，拉普拉斯方程為

r為避雷針邊緣半徑；dx、dy分別為建筑物在x、y方向上劃分的最小網(wǎng)格長度；m、n分別為建筑物在x、y方向上劃分的網(wǎng)格數(shù)量；f(dx)、f(dy)分別為關(guān)于dx與dy的函數(shù)，用來做變網(wǎng)格劃分

(1)

將模擬域的上邊界設(shè)置為第一類邊界，頂邊界的電位值φ(z= 500 m)為固定值，φ=zEb=-5 000 kV；底邊界為地面和建筑物尖端，也為第一類邊界，φ=0；其他4個側(cè)面均為第二類邊界，電位梯度在邊界上法向分量為0。由于，網(wǎng)格較多，采用超松弛迭代法求解變網(wǎng)格下的三維有限差分方程，從而獲得空間每個格點上的電位值，進而根據(jù)式(2)求得電場強度。

(2)

式(2)中：E為電場值，V/m。

1.2 模型驗證

為了驗證本文基于有限差分法建立的三維變網(wǎng)格電場計算模型的計算結(jié)果，利用基于有限元法的Comsol 軟件對上節(jié)所述的雷暴云電場下的建筑物尖端周圍的電場進行了計算。由于尖端處的網(wǎng)格尺寸大小，直接影響計算結(jié)果，而遠離尖端時，網(wǎng)格尺寸對尖端處的電場計算精度則影響較小。因此，在網(wǎng)格剖分時，則要保證尖端處的網(wǎng)格尺寸大小相當(dāng)。本文模型采用長方體對空間進行離散，在尖端處網(wǎng)格尺寸最小，為0.1 m，隨后逐步變大，最大網(wǎng)格尺寸為22.261 m；Comsol軟件中，采用四面體對空間進行離散，為了同本模型的結(jié)果對比，尖端處的最小網(wǎng)格尺寸也為0.1 m，隨后進行自由剖分，最大網(wǎng)格尺寸為20 m。以此網(wǎng)格剖分為例，利用上述兩種方法計算建筑物尖端周圍的電場大小。結(jié)果如圖2所示，當(dāng)環(huán)境電場為-10 kV/m 時，尖端上空電場隨高度的變化特征。不難發(fā)現(xiàn)，兩者計算結(jié)果基本相同。由此可見，本文建立的三維變網(wǎng)格電場計算模型的計算結(jié)果可靠。

圖2 背景電場為-10 kV/m時尖端上方的電場變化特征

2 空間分辨率對先導(dǎo)起始的判據(jù)的影響

2.1 介質(zhì)擊穿模型和臨界半徑準則

地面突出物(如樹木、房屋等)及地形的起伏都會對其周圍的大氣電場產(chǎn)生畸變作用，地面突出物周圍的電場值與環(huán)境電場值的比值稱為電場畸變系數(shù)[33-34]。基于本文模型的數(shù)據(jù)，不同空間分辨率下建筑物尖端拐角表面上方的電場畸變系數(shù)隨高度的變化曲線圖，如圖3所示。

由圖3可以看出，同一空間分辨率下尖端拐角表面上方的電場畸變系數(shù)值在1 m范圍內(nèi)隨高度迅速下降，且不同空間分辨率下尖端處電場畸變系數(shù)的值相差巨大，尖端網(wǎng)格劃分越細，尖端上方的電場畸變系數(shù)就越大。介質(zhì)擊穿模型和臨界半徑準則皆是以導(dǎo)體或尖端表面的電場是否超過臨界電場閾值來作為先導(dǎo)起始的判據(jù)。介質(zhì)擊穿模型在云中臨界電場閾值取Ecrit=100 kV/m的判據(jù)是根據(jù)實驗中在3個維度上網(wǎng)格間距皆取125 m得到的[11]。而云中臨界電場閾值Ecrit=150 kV/m的判據(jù)來自空間分辨率12.5 m[15]。臨界半徑準則的實驗室證明保持著7 m的恒定電極，并認為對于具有球形和圓柱對稱特征的電極臨界半徑在35 cm和10 cm左右趨于飽和。所以在使用這兩種判據(jù)進行數(shù)值模擬研究中對先導(dǎo)起始和傳播的研究時，會由于空間分辨率選取過大使得尖端表面電場計算值較小，從而導(dǎo)致判斷時先導(dǎo)難以始發(fā)，而空間分辨率選取過小時則會導(dǎo)致判斷先導(dǎo)過早產(chǎn)生。

ki為電場畸變系數(shù)；H為高度；hmin為空間分辨率

由于在數(shù)值模擬中，由一定空間分辨率所引起的系統(tǒng)誤差是一個固定的值，它與分辨率有關(guān)，而與建筑物結(jié)構(gòu)尺寸無關(guān)[35-36]。所以在使用這兩種判據(jù)時，可以假定0.01 m或更精細空間分辨率下計算出的尖端表面電場為標準實際電場值，將0.02、0.05、0.1 m等空間分辨率下尖端處的電場值與標準值相比進行修正，將得到的數(shù)個修正系數(shù)進行擬合，得到修正系數(shù)的擬合曲線圖。在之后采用過大空間分辨率的數(shù)值模擬研究中，可在擬合的修正系數(shù)曲線圖上得到修正系數(shù)，從而得到不同空間分辨率對應(yīng)的合適的判斷先導(dǎo)始發(fā)的尖端電場值，具體方法可參見文獻[37]。

2.2 Petrov 臨界長度法

Petrov根據(jù)確定垂直桅桿的雷擊距離實驗，提出了上行先導(dǎo)始發(fā)判據(jù)。認為當(dāng)流光區(qū)長度達到了0.7 m，并且同時流光區(qū)內(nèi)的電場強度超過了正極性流光區(qū)的平均電場值(約500 kV/m)，流光完成向先導(dǎo)狀態(tài)的轉(zhuǎn)化，認為上行先導(dǎo)起始。基于本文模型，更改不同空間分辨率下背景電場值，使?jié)M足0.7 m范圍內(nèi)電場強度皆超過500 kV/m。根據(jù)數(shù)據(jù)作出尖端拐角上方電場隨高度變化的曲線圖，如圖4所示，建筑物高100、200、300 m時當(dāng)最小空間分辨率分別取0.01、0.02、0.05、0.1 m時，根據(jù)Petrov臨界長度法均可判斷此時上行先導(dǎo)起始。

E為電場強度；rstreamer為流光區(qū)長度；hmin為空間分辨率；H為高度

按照Petrov的臨界長度法準則，用本文所建模型，將不同空間分辨率下在0.7 m處達到500 kV/m的電場強度所需要的背景電場值計算出來，如圖5所示。由圖5中曲線可知0.01 m的空間分辨率達到先導(dǎo)起始時所需要的背景電場值最大，隨著空間分辨率的增大，所需的背景電場值逐漸變小。100、200、300 m建筑物高度下0.01 m與0.1 m的空間分辨率所需先導(dǎo)起始的背景電場值相差分別約為4.96、2.4、1.46 kV/m。4個空間分辨率的平均背景電場值分別是22.31、11.7、7.93 kV/m。

Eb為背景電場值；hmin為空間分辨率

假設(shè)0.01 m分辨率下計算出的環(huán)境電場值最接近現(xiàn)實結(jié)果，將另外幾個分辨率下計算出的環(huán)境電場值除以0.01 m分辨率下的環(huán)境電場值，進行數(shù)據(jù)歸一化處理。圖6為3種不同建筑物高度在不同空間分辨率下計算出來的先導(dǎo)起始時環(huán)境電場結(jié)果的歸一化計算結(jié)果，由圖6中曲線可以看出在這3種建筑物高度模型中分辨率對Petrov的判據(jù)影響都不大。

Normalized Eb為電場歸一化后的值；hmin為空間分辨率

2.3 Becerra-Cooray上行先導(dǎo)自持始發(fā)模型

Gallimberti[25-26]建立了針對正極性先導(dǎo)的熱動力學(xué)模型，計算得到正極性先導(dǎo)起始所需的臨界電荷量約為1 μC。隨后，Petrov[10]在Gallimberti[25-26]建立的針對正極性先導(dǎo)的熱動力學(xué)模型的基礎(chǔ)上，考慮了二次起暈過程，對正先導(dǎo)熱力學(xué)模型進行了簡化和推廣。得出當(dāng)流光向不穩(wěn)定先導(dǎo)轉(zhuǎn)換發(fā)生在背景電勢分布高到足夠產(chǎn)生相當(dāng)于或高于1 μC的電暈電荷且當(dāng)不穩(wěn)定的先導(dǎo)持續(xù)傳播長度達到2 m時，產(chǎn)生穩(wěn)定的先導(dǎo)。

為了比較不同空間分辨率對該判據(jù)的影響，在基于有限差分法求解電場的程序基礎(chǔ)上，按照Beccera的自持上行先導(dǎo)始發(fā)模型建立了先導(dǎo)始發(fā)模型程序。

圖7為簡化的幾何近似法計算起始過程中電暈電荷量示意圖。

U為電勢；Utip為先導(dǎo)頭部電勢；U1為背景電勢；U2為初始電暈形成后的電勢分布；lL為先導(dǎo)長度；ls為流光長度；i為先導(dǎo)發(fā)展步伐

(3)

(4)

通過式(3)和式(4)可以得到電暈流光區(qū)前部的具體位置ls，m；且可以通過式(5)來表示：

(5)

同時，可以結(jié)合式(4)和式(5)得到二次起暈的電荷量ΔQ0，C/V。計算公式為

(6)

式(6)中：KQ為幾何參數(shù)，在簡化模型中采用3.5×10-11C/(V·m)的值。

(7)

(8)

式(8)中：E∞為穩(wěn)定狀態(tài)下正極性先導(dǎo)電位梯度，約為3×104V/m；x0為常量參數(shù)，由正極性先導(dǎo)速度v(約1.5×104m/s)和時間常量θ(約50×10-6s)相乘得到，m。在第i步時，先導(dǎo)頭部產(chǎn)生的電暈電荷量為

(9)

于是正極性先導(dǎo)以及頭部電暈區(qū)向前發(fā)展的距離為

(10)

(11)

按照背景電場值在10 s內(nèi)從0 kV/m發(fā)展到20 kV/m的增長率，計算出不同建筑物高度在不同空間分辨率下不穩(wěn)定先導(dǎo)始發(fā)和穩(wěn)定先導(dǎo)發(fā)展時所需的背景電場值，如圖8 所示。

圖8(a)為不穩(wěn)定先導(dǎo)起始時所需的背景電場值，圖8(b)為穩(wěn)定先導(dǎo)起始時所需的背景電場值。兩幅圖中不同建筑物高度起始所需環(huán)境電場值的曲線趨勢相同，0.01 m空間分辨率下先導(dǎo)起始所需的背景電場最小，隨著空間分辨率的增大先導(dǎo)起始所需的背景電場值逐漸增大。100、200、300 m建筑物高度模型中0.01 m和0.1 m的空間分辨率下不穩(wěn)定先導(dǎo)起始所需的背景電場值相差分別為0.282、0.18、0.148 kV/m，而穩(wěn)定先導(dǎo)起始所需背景電場值相差分別為1.888、1.392、1.204 kV/m。不同空間分辨率帶來的誤差不大，不穩(wěn)定先導(dǎo)起始平均背景電場分別為3.21、1.663 5、1.124 kV/m。穩(wěn)定先導(dǎo)起始平均電場分別為36.104、19.248、13.143 kV/m。

E為電場強度；hmin為空間分辨率

圖9為三種不同建筑物高度在不同空間分辨率下計算出來的先導(dǎo)起始時環(huán)境電場結(jié)果的歸一化計算結(jié)果。由圖9中曲線可以看出，在這三種建筑物高度模型中分辨率對becerra的判據(jù)影響都不大。

Normalized Eb為電場歸一化后的值；hmin為空間分辨率

通過研究可知，Petrov的臨界長度法和becerra的臨界準則法受分辨率影響較小。對于研究先導(dǎo)速度、電流大小等數(shù)值模擬研究使用becerra的判據(jù)，對于研究上行先導(dǎo)傳播分形等數(shù)值模擬研究使用Petrov的判據(jù)即可。

3 結(jié)論與討論

通過對建立的長、寬為0.2、0.2 m，高分別為100、200、300 m的建筑物尖端進行數(shù)值模擬研究，得出如下結(jié)論。

(1)針對以表面電場為判據(jù)的先導(dǎo)始發(fā)模型，空間分辨率帶來的影響較大，在使用此類判據(jù)時，可以通過擬合修正系數(shù)的方法來減小誤差。

(2)針對一定區(qū)域內(nèi)的電場或電勢作為判據(jù)時，只要空間分辨率小于判定范圍尺寸則不會對判定結(jié)果帶來較大影響。

綜上所述，當(dāng)僅以表面電場作為判據(jù)時，空間分辨率帶來的影響巨大；但當(dāng)以一定區(qū)域內(nèi)的電場或電勢作為判據(jù)，則發(fā)現(xiàn)空間分辨率所帶來的影響相對較小，幾乎可以忽略。因此，在實際的模擬中，針對不同的判據(jù)，需不同對待，不只是提高分辨率。為了節(jié)省計算機的計算耗時，根據(jù)所要研究的內(nèi)容，選取合適的先導(dǎo)始發(fā)判據(jù)，進而選擇能表征尖端大小和判據(jù)中出現(xiàn)的尺寸范圍的網(wǎng)格即可。此外，本文研究尚未考慮其他尺寸建筑物尖端始發(fā)上行先導(dǎo)時的情況，在后續(xù)工作中將對此進行研究。