基于蒙特卡洛方法的覆冰導(dǎo)線舞動跳閘概率計算

馬增泰，羅萍萍，孫天杰，林濟鏗，劉 輝

(1.上海電力大學(xué)電氣工程學(xué)院，上海 200090；2.同濟大學(xué)電子與信息工程學(xué)院，上海 201804；3.國網(wǎng)安徽省電力有限公司，安徽 合肥 230022)

架空導(dǎo)線的舞動是導(dǎo)線覆冰后，在橫向風(fēng)的作用下形成的低頻率、大振幅的自激振動，其對電力系統(tǒng)的安全運行構(gòu)成嚴(yán)重的威脅[1-4]。安徽、浙江、河南等省份在多次遭受冰凍天氣時，均出現(xiàn)了覆冰導(dǎo)線舞動導(dǎo)致的大范圍線路跳閘等惡性事件[5]。由于架空輸電線路的舞動受氣象條件影響較大，如何科學(xué)合理地根據(jù)氣象信息來計算架空線路的覆冰舞動跳閘概率，并對跳閘概率大的線路采取預(yù)防措施，對于提升電網(wǎng)應(yīng)對災(zāi)害天氣的能力，具有非常積極的理論和實際工程意義。

目前對覆冰架空導(dǎo)線舞動的研究主要集中在2個方面：舞動機理與模型的研究和基于線性化思想的舞動跳閘概率研究。舞動機理與模型的研究包括垂直方向舞動機理[6]、扭轉(zhuǎn)方向舞動機理[7]等。覆冰導(dǎo)線的舞動模型主要有單自由度舞動模型、雙自由度舞動模型及三自由度舞動模型[8]，并采用包括增量諧波法[9]、矩陣攝動法[10]、有限單元法[11]、牛頓法[12]等方法推導(dǎo)覆冰導(dǎo)線的舞動運動方程。文獻(xiàn)[12]使用牛頓法建立了單自由度舞動模型，并使用慢變參數(shù)法獲得舞動幅值的解析解，然而覆冰導(dǎo)線舞動是大幾何非線性運動，慢變參數(shù)法只適用于弱非線性運動方程；文獻(xiàn)[13]采用有限元的思想，使用索單元來建立覆冰導(dǎo)線的舞動模型，并使用振型疊加法求解該模型，得到了導(dǎo)線舞動的幅值；文獻(xiàn)[14]利用D’Alembert原理建立了覆冰導(dǎo)線的單自由度與雙自由度舞動模型，算例驗證了方法的正確性；文獻(xiàn)[15]使用Hamilton原理建立覆冰導(dǎo)線的三自由度舞動模型，并實現(xiàn)了對舞動形態(tài)的數(shù)值模擬，然而三自由度舞動模型不能得到舞動幅值的解析解。

基于線性化思想的舞動跳閘概率研究，其基本思想是先計算出覆冰導(dǎo)線舞動的幅值，然后與相間距離比較，進(jìn)而得到以導(dǎo)線舞動幅值為自變量的線性化的跳閘概率。文獻(xiàn)[16]提出了一種基于支持向量機(support vector machine,SVM)和自適應(yīng)增強(adaptive boosting,AdaBoost)分類器的覆冰架空線路舞動跳閘預(yù)警方法，通過挖掘易舞動區(qū)域的氣象數(shù)據(jù)來構(gòu)建基于SVM分類器的易舞氣象預(yù)報模型，從而得到舞動幅值的預(yù)測結(jié)果以及相應(yīng)的跳閘概率；文獻(xiàn)[17]使用拉力傳感器及覆冰導(dǎo)線實時水平荷載估算出導(dǎo)線舞動幅值，然后與導(dǎo)線相間距離比較得到舞動跳閘概率；文獻(xiàn)[18]基于現(xiàn)場舞動觀測點的溫度和濕度的信息對舞動幅值進(jìn)行實時預(yù)測；文獻(xiàn)[19]基于拉力與角度傳感器所得到的覆冰導(dǎo)線張力與風(fēng)攻角，推出覆冰導(dǎo)線舞動的幅值，并得到跳閘概率。基于線性化思想的舞動跳閘概率研究不足之處在于僅以舞動幅值與導(dǎo)線相間距離來計算跳閘概率，未考慮兩相導(dǎo)線的非同步舞動。

本文提出基于蒙特卡洛方法的導(dǎo)線舞動跳閘概率計算方法。由風(fēng)速與凍雨量等氣象信息，計算導(dǎo)線覆冰厚度并建立風(fēng)速的Weibull分布函數(shù)，然后根據(jù)牛頓定律建立覆冰導(dǎo)線單自由度舞動模型；使用Fluent軟件對不同冰厚的半橢圓型覆冰導(dǎo)線進(jìn)行氣動力仿真，得到氣動系數(shù)與覆冰厚度的關(guān)系式，進(jìn)而采用里茨—伽遼金方法求解舞動模型，獲得以風(fēng)速為自變量的舞動幅值解析表達(dá)式，并由風(fēng)速的分布函數(shù)獲得舞動幅值的概率分布函數(shù)；為了解決舞動跳閘數(shù)據(jù)收集過程的繁瑣且效率低下的問題，本文通過分析三相導(dǎo)線發(fā)生非同步舞動時的相位關(guān)系，得到某一舞動幅值下線路固有的跳閘概率計算模型；最后，使用蒙特卡洛方法模擬覆冰導(dǎo)線舞動跳閘的過程，獲得覆冰導(dǎo)線舞動越限跳閘概率的統(tǒng)計值。

1 舞動跳閘故障的概率模型

1.1 舞動幅值的概率模型的建立

導(dǎo)線的舞動是在覆冰的基礎(chǔ)上發(fā)生的，本文以極端天氣下凍雨時間里的風(fēng)速與凍雨量信息計算導(dǎo)線的覆冰厚度[20]，覆冰厚度為

rice=

(1)

式(1)中的rice為圓形覆冰冰厚，而實際導(dǎo)線的覆冰常是不規(guī)則的半橢圓型或新月型。本文結(jié)合實際情況[21]，將冰厚為rice的圓形覆冰截面等效為參數(shù)為rH、ry的半橢圓形覆冰截面，其中rH為橫向增長冰厚，ry為豎向增長冰厚，如圖1所示。假設(shè)2種覆冰形狀的覆冰截面面積相等和半橢圓形覆與冰橫豎向覆冰增長成比例，則rH、ry以及rice的關(guān)系為

圖1 覆冰導(dǎo)線截面Figure 1 Section of iced wire

(2)

式中r為裸導(dǎo)線導(dǎo)線半徑。對于新月形覆冰，可將導(dǎo)線ry方向的冰厚與rH方向的冰厚取相等值ra，令ra=2rice。

研究覆冰導(dǎo)線舞動的幅值與穩(wěn)定性的問題時，常用到單自由度舞動模型，即導(dǎo)線垂直舞動模型。采用多項式擬合覆冰截面的氣動系數(shù)后，單檔覆冰導(dǎo)線的單自由度舞動運動方程[12]可以表示為

(3)

式中m為覆冰導(dǎo)線質(zhì)量，m=πr2ρL+0.5πρice(rHry-r2)；ζ為阻尼比；ω為覆冰導(dǎo)線體系自振頻率；k為覆冰導(dǎo)線體系的剛度；ρair為空氣密度；D為覆冰導(dǎo)線直徑；y為覆冰導(dǎo)線在垂直方向的位移；c1、c2均為氣動參數(shù)。

由于覆冰導(dǎo)線舞動是小應(yīng)變、大位移的強幾何非線性運動，本文使用里茨—伽遼金法求解式(3)。設(shè)覆冰導(dǎo)線舞動在垂直方向的位移響應(yīng)為

y=ccos(ωt+φ)=acos(ωt)+bsin(ωt)

(4)

(5)

同時，對式(4)變分得：

δy=cos(ωt)δa+sin(ωt)δb

(6)

根據(jù)δa與δb的任意性，結(jié)合式(5)、(6)，在一個周期內(nèi)積分，整理可得覆冰導(dǎo)線舞動幅值為

(7)

由于覆冰導(dǎo)線覆冰厚度不同，其氣動系數(shù)也不同。本文使用Fluent仿真軟件對覆蓋不同冰厚的導(dǎo)線進(jìn)行氣動力仿真，獲取氣動系數(shù)三次多項式表達(dá)式，然后對該多項式的系數(shù)與冰厚rH建立函數(shù)關(guān)系式，從而得到式(7)中參數(shù)c1、c2的表達(dá)式，即

(8)

將式(8)改寫為

(9)

式中k1、k2為關(guān)于式(7)的系數(shù)。

由式(9)可知，在獲取導(dǎo)線覆冰厚度及相關(guān)線路參數(shù)后，覆冰導(dǎo)線舞動的幅值主要取決風(fēng)速。由于氣象觀測點與輸電線路地理位置、地形地貌以及海拔高度均不相同，故觀測站風(fēng)速需要經(jīng)過修正，即

(10)

式中α為導(dǎo)線所處位置地面粗糙度因子，取值在0.1～0.5。

根據(jù)文獻(xiàn)[22]可知，Weibull分布能夠較好地描述風(fēng)速的分布，假設(shè)風(fēng)速序列(v1,v2,…,vn)服從Weibull分布，其概率分布函數(shù)為

(11)

式中σ為尺度參數(shù)；ξ為形狀參數(shù)；v為風(fēng)速。通過極大似然估計法等方法來確定式(11)中的參數(shù)σ、ξ，構(gòu)造對數(shù)似然函數(shù)為

L(ξ,σ)=

(12)

(13)

式(12)、(13)較為復(fù)雜，采用牛頓—拉夫遜迭代法求解其中參數(shù)，由式(13)可得：

(14)

其中，雅可比元素為

(15)

為求解式(14)，選取合適的初值，經(jīng)過反復(fù)迭代直到滿足max(Δσ，Δξ)<ε后，獲得Weibull分布的尺度參數(shù)σ與形狀參數(shù)ξ。獲取兩參數(shù)σ與ξ后，可以建立風(fēng)速的Weibull概率分布。由式(9)可知，舞動的幅值與風(fēng)速有關(guān)，經(jīng)過式(9)、(11)及概率論知識可將風(fēng)速的概率分布變換為舞動幅值c的概率分布，其表示為

(16)

由式(16)可知，設(shè)舞動幅值超過安全間距的概率為

Po=1-F(cse)

(17)

式中cse為覆冰導(dǎo)線舞動至相間最小安全間距的幅值，垂直方向上cse=Lx-Lxmin；Po為當(dāng)前風(fēng)速及覆冰條件下線路舞動幅值超過cse的概率。

1.2 覆冰導(dǎo)線非同步舞動單一幅值固有跳閘概率模型的建立

覆冰導(dǎo)線在某一綜合條件下的舞動幅值超過相間電氣安全限度時，存在發(fā)生相間閃絡(luò)跳閘故障的可能。由于目前缺少導(dǎo)線相間非同步舞動的跳閘模型，本文以垂直布線的三相導(dǎo)線為例，建立導(dǎo)線非同步舞動的某一幅值的固有跳閘概率模型，如圖2所示。

圖2 導(dǎo)線垂直布線Figure 2 Schematic diagram of vertical arrangement of wires

假設(shè)覆冰導(dǎo)線的A相與C相發(fā)生舞動，將舞動時的兩相導(dǎo)線中點位置分別為

(18)

式中ω為yA導(dǎo)線舞動時的角速度；φ為以C相導(dǎo)線為參照，A相導(dǎo)線舞動的相位；Lx為導(dǎo)線靜態(tài)時的垂直相間距離。其中φ為隨機值，體現(xiàn)了A相導(dǎo)線相對C相發(fā)生舞動的隨機性。

當(dāng)φ為0時，A相與C相導(dǎo)線同向舞動，由于兩相導(dǎo)線始終保持安全間距Lx，因而不會發(fā)生跳閘事故[23]；當(dāng)φ為π時，A相與C相導(dǎo)線反向運動，如果舞動幅值超過一定限度，兩相導(dǎo)線將因電氣間距過小或碰線而導(dǎo)致閃絡(luò)跳閘。

基于以上分析，當(dāng)導(dǎo)線舞動幅值越限后，A相導(dǎo)線在相位φ=π左右舞動時存在相間閃絡(luò)跳閘區(qū)域，如圖3所示。

由圖3可知，如果兩相導(dǎo)線舞動的幅值超過安全間距且其值為c，當(dāng)A相導(dǎo)線舞動的相位φ=φx+π在區(qū)間[φ1+π,φ2+π]內(nèi)時，A、C兩相導(dǎo)線將發(fā)生相間跳閘，此時φx滿足：

圖3 導(dǎo)線跳閘Figure 3 Schematic diagram of line trip

csin(ωt+φx+π)+Lx-csin(ωt)≤Lxmin

(19)

取ωt=π/2，并將式(19)取等號可得邊界跳閘相位，其表示為

(20)

由于相位為φ1+π與φ2+π的舞動波形關(guān)于ωt=π/2對稱，所以在舞動幅值c超過安全間距的條件下，發(fā)生跳閘的相位區(qū)間的大小為|2φ1|。因此，舞動幅值為c時的固有跳閘概率為

(21)

2 舞動跳閘故障過程的蒙特卡洛模擬

覆冰導(dǎo)線舞動的閃絡(luò)跳閘故障的模擬可以認(rèn)為是有多個可能性的抽樣[24]，其中包括確定導(dǎo)線舞動幅值是否越限、舞動的越限幅值及是否跳閘。以舞動的幅值是否越限為例，將舞動幅值超過安全間距的概率Po置于[0,1]區(qū)間內(nèi)，并產(chǎn)生一個在[0,1]區(qū)間均勻分布的隨機數(shù)R，R的位置表明舞動幅值是否越限，如圖4所示；風(fēng)速抽樣與舞動是否跳閘抽樣分別如圖5、6所示。

圖4 舞動幅值是否越限Figure 4 Galloping amplitude over limit or not

圖5中，vmax為通過氣象觀測獲取的最大風(fēng)速，vo為導(dǎo)線剛好舞動至相間電氣安全間距時的幅值co為所對應(yīng)的風(fēng)速，由式(9)可得：

圖5 風(fēng)速抽樣Figure 5 Sampling wind speed

(22)

使用蒙特卡洛方法可以計算一條線路第i檔導(dǎo)線的舞動幅值越限跳閘概率Pi的統(tǒng)計值，根據(jù)該方法，第j次模擬過程如圖7所示。

圖6 舞動是否跳閘抽樣Figure 6 Galloping trip or not

圖7 蒙特卡洛方法計算舞動跳閘概率Figure 7 Flow chart of the calculation of the galloping tripping probability using Monte Carlo method

如果設(shè)置模擬次數(shù)為N，則每次模擬都會產(chǎn)生一個值為0或1的Tj。由于大數(shù)法則保證在模擬次數(shù)足夠多即獲得的樣本足夠多之后，蒙特卡洛法獲取的估計值將收斂于真值，則第i檔覆冰導(dǎo)線舞動越限跳閘概率的統(tǒng)計值為

(23)

對于整條輸電線路，其整體舞動跳閘概率由各檔導(dǎo)線舞動跳閘概率邏輯串聯(lián)進(jìn)行計算，即

(24)

式中ws為第s條線路檔數(shù)。

3 求解步驟

依據(jù)線路所跨區(qū)域的大小劃分氣象分區(qū)，收集各分區(qū)的風(fēng)速及風(fēng)向、凍雨量、溫度、相對濕度等信息后，基于蒙特卡洛方法的覆冰輸電線路的舞動跳閘概率計算步驟如下：

1)計算各氣象分區(qū)覆冰厚度，利用Fluent軟件獲取各冰厚條件下的氣動參數(shù)；

2)利用最大似然方法估計式(16)中風(fēng)速Weibull分布的參數(shù)σ與ξ；

3)依據(jù)覆冰導(dǎo)線舞動條件[21]，如表1所示，判斷當(dāng)前線路第i檔導(dǎo)線是否發(fā)生舞動；

表1 導(dǎo)線起舞判斷因素Table 1 Conductor galloping judgment factors

4)如果發(fā)生舞動，則轉(zhuǎn)向步驟4，否則Pi=0，i=i+1，返回步驟3；

5)獲取線路各檔導(dǎo)線參數(shù)，由式(16)、(17)計算舞動幅值越限概率Po；

6)設(shè)置蒙特卡洛模擬次數(shù)為N，獲取每次模擬所得的Tj并進(jìn)模擬，由式(23)計算當(dāng)前線路的第i檔舞動跳閘概率Pi；

7)判斷當(dāng)前線路檔數(shù)i是否小于ws，如果是，則i=i+1，返回步驟3，否則轉(zhuǎn)到步驟8；

8)由式(24)計算當(dāng)前整條輸電線路的舞動跳閘概率PL。

本文方法計算流程如圖8所示。

圖8 計算舞動跳閘概率流程Figure 8 Flow chart of the calculation of the galloping tripping probability

4 算例分析

本文使用Fluent軟件對2～25 mm冰厚的半橢圓形型覆冰導(dǎo)線進(jìn)行氣動系數(shù)模擬。Fluent軟件通常用來模擬可壓縮或不可壓縮氣體的流動，由于Fluent已經(jīng)集合在ANSYS軟件中，所以本文直接使用ANSYS軟件中的Geometry模塊建立半橢圓形覆冰導(dǎo)線截面模型，并使用ICEM CFD模塊對覆冰導(dǎo)線截面進(jìn)行網(wǎng)格劃分，如圖9所示。

圖9 覆冰導(dǎo)線截面網(wǎng)格Figure 9 Gridding of iced wire cross section

設(shè)置圖9中的右邊界為氣體速度入口的邊界條件，左邊界為氣體流出的邊界條件。完成導(dǎo)線覆冰截面的網(wǎng)格劃分后，將網(wǎng)格文件導(dǎo)入Fluent模塊中，使用適用于圓柱體繞流的K-ω湍流模型，并在導(dǎo)線的垂直方向上設(shè)置升力監(jiān)測器以獲取當(dāng)前覆冰厚度與風(fēng)攻角的截面升力系數(shù)。

通過旋轉(zhuǎn)覆冰導(dǎo)線截面可改變迎面風(fēng)攻角，本文對導(dǎo)線在不同覆冰厚度時的風(fēng)攻角取值范圍均設(shè)定為[-60°,+60°]，且以10°為增量。在完成全部覆冰厚度及風(fēng)攻角的升力系數(shù)的測定后，先建立的各覆冰厚度時的風(fēng)攻角與升力系數(shù)的函數(shù)關(guān)系式，部分覆冰厚度的氣動系數(shù)曲線如圖10所示。

圖10 不同覆冰厚度下的氣動系數(shù)曲線Figure 10 Aerodynamic coefficient curves under different icing thicknesses

然后，使用三次多項式曲線擬合方法建立氣動參數(shù)c1與氣動參數(shù)c2關(guān)于rH的表達(dá)式，即

(25)

式(25)可代入到式(7)中以代替氣動參數(shù)。

本文對2019年2月遭受冰凍天氣影響的安徽宣城地區(qū)的500 kV輸電線路覆冰舞動短路跳閘概率進(jìn)行了計算。由于輸電線路所跨區(qū)域較大，將線路所經(jīng)區(qū)域劃分為4個氣象分區(qū)，分別獲取各分區(qū)的風(fēng)速、凍雨量、溫度以及濕度等信息并對風(fēng)速的Weibull分布進(jìn)行擬合，分布函數(shù)如圖11所示。

圖11 風(fēng)速的Weibull概率分布函數(shù)Figure 11 Weibull probability distribution function of wind speed

2月10日到12日，凍雨持續(xù)時間6 h，結(jié)合各分區(qū)氣象信息判定輸電線路發(fā)生舞動。根據(jù)本文方法，獲取各檔導(dǎo)線的參數(shù)及其氣象分區(qū)下的Weibull分布函數(shù)，然后使用蒙卡羅方法由式(23)計算得各線路及各檔舞動跳閘概率，其中山瀝線單檔導(dǎo)線舞動跳閘率如表2所示。

表2 山瀝線單檔導(dǎo)線舞動跳閘率Table 2 Single-speed wire galloping trip rate of Shanli line

對于不同線路，在各氣象分區(qū)下，經(jīng)式(23)計算得到各檔線路Pi，由式(24)計算得到整條線路舞動越限跳閘概率PL。各條線路的舞動越限跳閘率與實際跳閘記錄的對比如表3所示。

表3 線路舞動跳閘率Table 3 Line galloping trip rate

由表3可知，線路舞動跳閘率超過0.6的一些線路均發(fā)生舞動跳閘事件。計算獲得的山橋線舞動跳閘概率值較小，實際中也沒有發(fā)生舞動跳閘事件，分析認(rèn)為山橋線處于覆冰區(qū)域邊緣，且部分檔導(dǎo)線與風(fēng)向較為一致，從而破壞了起舞條件。由此可知由本文方法所獲得覆冰導(dǎo)線舞動跳閘概率與實際故障存在關(guān)聯(lián)性。從表3還可以看到舞動跳閘事件多集中發(fā)生在連續(xù)檔上，考慮到連續(xù)檔導(dǎo)線線路參數(shù)及覆冰、風(fēng)速狀況等較為一致，由此可知導(dǎo)線的舞動及其相應(yīng)的跳閘事件往往集中發(fā)生在易舞動的同一氣象條件下(包括風(fēng)速、降雨、溫度等)。

5 結(jié)語

結(jié)合安徽的舞動跳閘概率分析，本文方法具有以下特點：

1)對氣動系數(shù)與覆冰厚度的函數(shù)擬合解決了覆冰舞動時冰厚的不確定性對氣動系數(shù)的影響；

2)基于風(fēng)速威布爾分布、里茨—伽遼金法幅值表達(dá)式而獲得的舞動幅值概率分布函數(shù)能夠描述舞動越限概率；

3)通過分析三相導(dǎo)線舞動時的相位關(guān)系，提出舞動幅值固有跳閘概率計算公式，該公式解決了三相覆冰導(dǎo)線舞動跳閘概率因非同期舞動的隨機性而缺少理論計算方法的問題；

4)基于蒙特卡洛方法獲得的覆冰導(dǎo)線舞動越限跳閘概率結(jié)果準(zhǔn)確可靠。