周超 ,孟凡澤 ,黃鋼 ,姬昆鵬
(1.華北電力大學 電站能量傳遞轉化與系統教育部重點實驗室,北京 102206;2.華北電力大學 河北省電力機械裝備健康維護與失效預防重點實驗室,保定 071003;3.中國電力科學研究院,北京 100192)
輸電導線重覆冰導致的輸電線路倒塔、斷線及舞動等次生危害時有發生,對跨區域輸電造成嚴重威脅[1-3].由于覆冰環境因素、導線自身運動特性的復雜性,現有實驗及預測模型僅限于準靜態導線,較少涉及運動導線覆冰形狀數值模擬計算[4-5].因此,亟需建立考慮扭轉因素的輸電導線覆冰預測模型.
覆冰模型主要分為三類:經驗模型、理論模型和數值計算模型[6].Langmuir 等[7]研究了水滴在無限長圓柱體上的運動軌跡,提出水滴碰撞系數概念.Makkonen 等[8]通過研究覆冰表面熱平衡過程建立覆冰時變物理模型,結合導線直徑、風速、溫度等環境參數確定了導線覆冰質量與時間的函數關系.Myers等[9]提出基于潤滑理論的水膜流動模型,研究了冰層內部的傳熱效應。Fu 等[10]利用邊界單元法建立簡化二維圓柱導線覆冰仿真計算模型,并計算了導線局部水滴碰撞系數。蔣興良等[11]提出覆冰預測最優時間步長模型,研究了風速、溫度等氣象參數對霧凇覆冰質量和冰形的影響.侯碩等[12]基于潤滑理論建立了任意二維截面上積冰數值模型,并與Messinger 模型預測結果進行了對比.劉春城等[13]建立雨凇覆冰模型,研究了覆冰質量與水滴中值直徑、風速及導線直徑等參數隨時間的變化規律.梁曦東等[14]利用邊界單元法建立了導線時變覆冰預測模型,計算時可以針對覆冰過程中環境參數的變化進行更新.何青等[15]提出考慮碰撞系數、焦耳熱和導線表面溫度的凍結系數計算方法,并指出考慮碰撞系數和導線表面溫度使得凍結系數變大.周超等[16]建立了三維時變導線覆冰預測模型,基于多步積冰算法計算得到的覆冰預測精度較單步積冰算法提高8%.
覆冰過程中,導線受覆冰偏心和風載荷等的影響,容易發生扭轉位移.導線扭轉后,由于冰形的影響,等效直徑增大,可捕獲更多水滴,覆冰冰形由單側覆冰逐漸轉變為近似圓形[17].韓興波等[18]研究發現,相同冰期內導線中心點扭轉角度最大,兩端最小且呈對稱分布.樊社新等[19]對不同覆冰厚度導線扭轉剛度進行實驗測量,使用最小二乘法得到了導線扭轉剛度隨冰厚變化的公式.李嘉祥等[20]提出了一種考慮覆冰偏心的分裂導線扭轉剛度計算方法,認為不均勻的覆冰會對導線扭轉剛度有顯著影響.傅觀君等[21]實驗研究發現,分裂導線覆冰具有強不均勻性的特點,且隨著覆冰厚度的增加,不均勻性更加明顯.胡琴等[22]對比了分裂導線和單導線的覆冰扭轉特性,發現分裂導線扭轉角度僅為單導線的10.2%.
準確預測導線覆冰扭轉后的冰形,對預防線路災害具有重要意義,現有研究較少考慮導線扭轉覆冰時冰形的演變過程.因此,本文通過考慮覆冰偏心、風速等對導線的耦合扭轉特性以及冰形對導線扭轉過程中剛度的影響,建立導線動態扭轉覆冰預測模型,并與實驗數據對比.基于新建立的預測模型,研究導線不同傾斜角度下冰形的變化及溫度、風速等環境參數對導線扭轉速度的影響.
1)空氣流場計算
空氣流場計算采用基于RNG(重整化群)方法的k~ε湍流模型,該模型對覆冰導線空氣流場計算具有較高的精度[23].將導線周圍的流場看作黏性不可壓縮的流體,使用N-S 方程計算空氣流場的速度分布[24],即
2)過冷水滴碰撞軌跡計算
在流場中,過冷水滴在碰撞導線的過程中受到多種力的作用[25],可表示為
式中:FN為水滴所受拖拽力,FG為水滴自身重力,Ff為水滴所受氣流的浮力,FM為其它作用力.使用歐拉法求解氣液兩相流模型,得到水滴在流場中的速度分布,進而求得導線各位置的局部水滴碰撞系數.局部水滴碰撞系數[26]定義為
式中:dY是流場中相鄰兩條水滴流徑的垂直距離,dL是捕獲的相鄰水滴在導線表面相隔的長度.實際導線覆冰的環境中,水滴直徑大小不唯一,根據文獻[10]數值模擬時使用水滴中值直徑(MVD)進行計算,水滴分布模型為Monodisperse.
3)導線表面結冰計算
過冷水滴在導線表面經過對流、蒸發、升華和凝固等熱量傳遞后凍結成冰,整個過程中能量守恒.基于導線表面的熱交換過程,建立能量守恒方程[27],表示為
式中:LWC 為液態水含量;hf為水膜高度,cf和cs為水和冰的比熱容,Levap和Lfusion為水的蒸發潛熱和冰的融化潛熱,E∞和Ef為流場溫度和水膜溫度,為液滴運動速度和壁面局部液滴撞擊速度,T∞和Tf為環境溫度和覆冰表面溫度,Q為防冰熱通量.
覆冰密度與環境溫度、風速、液態水含量等環境參數有關,根據文獻[24]采用Makkonen-Stallabrass密度模型,計算式為
其中RM為計算參數,其表達式為
覆冰開始后,導線在耦合力矩的作用下發生扭轉.對于長度為N的導線,將其等分為2n+1 個節點,中間節點為N0.相鄰節點間長度為x,如圖1所示.

圖1 導線扭轉示意圖Fig.1 Schematic diagram of conductors torsion
在對稱的情況下只需計算其中一半節點的扭轉角度,控制公式為[28]
式中:k為時間步,i為節點坐標,j為力矩坐,G為導線剪切強度,為導線覆冰后轉動慣量,ΔTn為不同時間步作用于導線的扭轉力矩增量.各參數計算式為
根據文獻[10],覆冰計算時將鋼芯鋁絞線簡化為圓柱體,數值模擬使用Fluent 及Fluent Icing 完成.流體域為大小20D×20D×40D(D為輸電導線直徑)的三維計算域,導線距速度入口及上下邊界的距離均為10D.為提高覆冰計算的準確度,在導線表面使用適于流體計算的膨脹網格,并對速度入口處網格進行了加密.網格劃分如圖2 所示,計算域邊界面參數見表1.

表1 流場邊界參數

圖2 流體模型網格設置Fig.2 Mesh settings of the fluid model
圖3 為導線覆冰扭轉計算流程.與靜態導線相比,導線扭轉后其周圍空氣流場和水滴碰撞導線的軌跡隨之發生變化.

圖3 導線覆冰扭轉計算流程圖Fig.3 Flow chart of conductor icing torsion calculation
時間步長ts的選擇影響到覆冰計算的精度和效率,根據文獻[11]選定單步步長為360 s,此時覆冰冰形準確度與實驗吻合度較好.由圖3 可知,計算時,首先對模型的初始邊界條件及風速、溫度、液態水滴濃度及MVD 等外部環境參數進行設定.由設定的環境參數對導線外部的空氣流場和水滴碰撞軌跡進行計算,獲得導線壁面的水滴收集效率.根據所得的水滴收集效率,通過式(5)、式(6)求解得到導線各節點的冰形和覆冰質量,從而得到導線不同位置處覆冰質心的坐標.將上述計算所得覆冰質量、冰質心和風速等相關參數代入式(9)、式(10),得到導線覆冰后各節點所受耦合扭轉力矩,從而根據式(8)獲得導線各節點扭轉角度并加以修正.總覆冰時長范圍內,每完成1 個時間步的計算后,需對扭轉后的覆冰導線重新劃分網格,并對模型邊界條件及外部環境參數進行更新.重復迭代上述計算過程,直至到達計算所設定的總覆冰時長.
為驗證本文導線覆冰計算方法的有效性,對文獻[17]所做導線覆冰扭轉風洞實驗進行驗證.實驗記錄了導線在不同時刻的扭轉角度,并分別在第125 min、245 min、363 min和444 min 4個時間節點提取了導線覆冰后的冰形,導線的扭轉角度分別為140°、223°、263°和291°.文獻[28]采用時不變參數對實驗進行了模擬,實驗及文獻相關參數見表2.

表2 覆冰扭轉實驗及文獻[28]相關參數
圖4 為上述4 個扭轉角度下數值模擬導線中間節點(0.22 m 長度處)冰形與文獻[17]、文獻[28]冰形的對比.對比發現,扭轉角度為140°時數值模擬覆冰面積相比文獻[17]中導線覆冰面積大16.8%,其余3 種角度下數值模擬結果比文獻[17]小13.75%、9.56%、7.54%.同文獻[28]相比,4種扭轉角度下數值模擬計算得到的覆冰截面積精度分別提高36.60%、42.21%、-7.3%和6.0%.

圖4 不同扭轉角度下覆冰形狀對比Fig.4 Comparison of ice cover shapes at different torsion angles
分析圖4 可知,導線覆冰扭轉角度為140°時,數值模擬所得覆冰在導線表面覆蓋面積較文獻[17]大6.45%,與文獻[28]相比,準確度提高11.35%。數值模擬所得覆冰面積大于文獻[17]所得,原因是覆冰早期冰質量較小,風載荷力矩對導線扭轉影響較大,模擬中風載荷力矩小于實驗中的實際值.隨著覆冰質量增加,風載荷力矩對導線扭轉的影響減弱,數值模擬所得冰形也更接近實驗結果.扭轉角度為223°時,數值模擬所得冰形在導線迎風面與實驗冰形符合較好,迎風面下側有突起,背風面的厚度略小于實驗冰形.扭轉角度為263°時,實驗中導線迎風面的覆冰厚度均勻,與數值模擬所得結果相比,在迎風面下側的覆冰厚度接近,迎風面上側覆冰厚度較大.扭轉角度為291°時,數值模擬所預測的冰形與實驗得到的冰形整體擬合程度良好.覆冰過程中,數值模擬所預測的冰形與實驗得到的冰形相比更為光滑.原因可解釋為:實驗中水滴直徑及風速等環境參數的變化使得水滴與導線表面碰撞軌跡發生改變,導致實驗中冰形表面凹凸程度超過數值模型.
以導線直徑為19.6 mm、檔距100 m 為例進行覆冰扭轉模擬,相關參數如表3 所示.提取8 個不同時間點的導線中間節點(50 m 長度處)冰形進行分析,結果如圖5所示.

表3 覆冰模擬參數

圖5 導線覆冰扭轉冰形變化Fig.5 Conductor ice covering and twisting ice shape change
分析圖5 可知,導線受到覆冰偏心的作用逐漸扭轉.隨著數值模擬時間的增加,覆冰向背風側移動,冰形由單側覆冰變為近似圓形.
在役運行輸電線路中,受弧垂、高差等因素影響,輸電導線必然存在一定程度的傾斜.短時間內導線傾角不會對覆冰特性產生明顯的影響.然而,隨著覆冰時間的增加,傾斜導線下游則可能出現更多的覆冰.以表2 條件為工況參數,設置水平傾角α為30°、45°和60°典型導線傾角工況,研究導線傾角對覆冰扭轉的影響,導線傾角定義見圖6.導線長度為0.44 m,模擬時長為245 min,提取導線中間節點冰形進行分析.

圖6 導線傾角定義Fig.6 Inclination angle definition for transmission conductor
2.3.1 計算結果分析
3種工況計算結果如圖7~圖9所示.

圖7 傾斜30°導線與水平導線冰形對比Fig.7 Comparison of tilted 30° conductor and horizontal conductor
傾斜30°導線:覆冰125 min 時覆冰扭轉角度為143.6°,覆冰面積較水平導線計算結果大8.4%.覆冰245 min 時覆冰扭轉角度為230.2°,覆冰面積較水平導線計算結果大15.16%.覆冰早期,30°傾斜導線與水平導線冰形相比沒有明顯差異.隨著計算時間增加,傾斜導線迎風面的覆冰厚度較水平導線不均勻性增大.原因可解釋為:未凍結的水滴在重力和氣流的作用下,由導線上游向斜下方流動,更多的水滴在導線迎風面的下側發生凍結,并累積.
傾斜45°導線:覆冰125 min 時扭轉角度為146.2°,覆冰面積較水平導線計算結果大15.31%。覆冰245 min 時扭轉角度為234.4°,覆冰面積相較水平導線大20.5%.與導線傾角為30°時相比,導線傾角為45°時,覆冰在迎風面下側積累更多,覆冰在導線表面的覆蓋面積更大,但整體覆冰形貌類似.
傾斜60°導線:覆冰125 min 時扭轉角度為150.4°,覆冰面積較水平導線大18.26%.覆冰245 min時扭轉角度為239.2°,覆冰面積較水平導線大26.3%.傾斜導線迎風面的冰形呈“溝壑”狀,覆冰在迎風面上側的厚度略小于水平導線,下側的厚度則出現了明顯的增加.
對比圖7、圖8 和圖9 可知,覆冰時間較短時,傾斜導線和水平導線的覆冰計算結果不會出現顯著差異.隨著時間的增加,覆冰在導線迎風面厚度出現分化,上側較小,下側明顯增加.導線傾斜角度較小時,冰形表面較為光滑,傾斜角度增大后,覆冰形狀逐漸粗糙,并在導線的迎風面出現了“溝壑狀”分層.原因可解釋為:導線傾斜角度的增加使未凍結水滴受到重力、離心力和氣流的影響增大,并隨著水膜更快地向導線的斜下方流動,迎風面上側的水滴減少,下側增多,覆冰出現了分層,導致“溝壑”的出現.

圖8 傾斜45°導線與水平導線冰形對比Fig.8 Comparison of tilted 45° conductor and horizontal conductor

圖9 傾斜60°導線與水平導線冰形對比Fig.9 Comparison of tilted 60° conductor and horizontal conductor
由中間節點向兩側各取50 mm 長導線(60°傾角)展示三維扭轉覆冰形貌[29],結果如圖10 所示.覆冰125 min 時導線上游覆冰厚度較為均勻,下游迎風面下側的覆冰厚度大于上側.覆冰245 min 時導線覆冰厚度的不均勻性加劇,上下游覆冰均出現了分層.原因可解釋為:覆冰過程中,多數未凍結水滴隨著水膜的流動在導線下游凍結,使導線下游覆冰厚度更早出現分層.剩余部分未凍結水滴在流動過程中凍結,因此隨著時間的增加,導線上游覆冰也逐漸呈“溝壑”狀.圖11 和圖12 分別為圖10 所示覆冰形貌左右兩端冰形剖面圖.

圖10 傾斜60°導線三維覆冰形貌Fig.10 Three-dimensional ice topography of tilted 60°conductor

圖11 覆冰125 min左右兩端冰形剖面圖Fig.11 Ice profile at both ends of the ice cover for about 125 min

圖12 覆冰245 min左右兩端冰形剖面圖Fig.12 Ice profile at both ends of the ice cover for about 245 min
圖13 為導線覆冰扭轉角度隨時間變化曲線.由圖中可以看出,覆冰初期水平導線及各傾斜角度導線的扭轉角度差別不大.隨著時間的增加,大傾角導線扭轉角度相較小傾角導線更大.

圖13 α=0°、30°、45°、60°導線扭轉角度隨時間變化曲線Fig.13 α=0°,30°,45°,60° conductor torsion angle change curve with time
導線覆冰的結果受到液態水含量(LWC)、液態水滴中值直徑(MVD)、溫度和風速等因素的影響[30].環境參數的差異改變了導線的覆冰形狀和質量,導致扭轉過程中扭轉剛度及扭轉力矩發生變化,最終影響導線覆冰預測精度.以直徑34.9 mm、檔距100 m導線中間節點(50 m長度處)為對象,對比分析不同溫度、風速、MVD影響下的導線覆冰扭轉特性,不同工況下的環境參數如表4所示,結果如圖14~圖16所示.

表4 覆冰工況條件參數

圖14 溫度對覆冰質量及扭轉角度的影響Fig.14 Effect of temperature on ice cover mass and torsion angle
分析圖14 可知,導線覆冰質量和扭轉角度隨著溫度降低而增大.原因可解釋為:溫度降低使得過冷水滴在導線表面的凍結系數增大,單位時間內的覆冰質量增加,冰形外擴,覆冰重心外移.由式(10)可知覆冰質量增加及覆冰重心外移使得冰載荷力矩增加,導線扭轉速度加快.但隨著溫度的降低,溫度對凍結系數的影響效果減弱,覆冰質量增速放緩,導線扭轉增速隨之下降.
分析圖15 可知,導線覆冰質量和扭轉角度隨著風速增大而增大.原因可解釋為:風速增大使得相同時間內導線表面捕獲的水滴數量增多,覆冰質量隨之增加,導線表面覆冰由駐點位置向兩側生長,覆冰厚度增加.根據式(10)可知冰載荷力矩和風載荷力矩的增加使得導線扭轉的速度加快。覆冰早期,風載荷力矩作為導線扭轉的主要因素,風速增加對導線扭轉速度影響更為明顯.

圖15 風速對覆冰質量及扭轉角度的影響Fig.15 Effect of wind speed on ice cover mass and torsion angle
分析圖16 可知,導線覆冰質量和扭轉角度隨著水滴MVD 增大而增加.原因可解釋為:隨著水滴MVD 增大,水滴受氣流擾動效應減弱,導線表面的水滴碰撞系數增加,覆冰質量和面積增大,導線的扭轉速度加快.但隨著水滴MVD 的增長,導線駐點位置覆冰厚度增加較少.與溫度、風速等環境參數相比,對導線扭轉速度的影響較弱.

圖16 水滴MVD對覆冰質量及扭轉角度的影響Fig.16 Effect of MVD of raindrops on ice cover mass and torsion angle
圖17 為導線覆冰扭轉后不同時間收集效率對比.分析可知,導線扭轉后覆冰形狀的變化,使導線及覆冰表面的局部水滴收集效率降低.但扭轉后的導線,等效直徑增大,獲得更大的水滴捕獲面積,整體水滴收集效率提高.工況2 條件下,隨著導線扭轉角度的增加,覆冰60 min 和覆冰90 min 時導線表面水滴收集效率較覆冰30 min 時分別增加6.38%和11.4%.導線扭轉使單位時間內覆冰質量增加,結果如圖18 所示.覆冰初期,靜態覆冰和扭轉覆冰的覆冰質量沒有明顯差異,隨著覆冰時間的增加,扭轉導線的覆冰質量逐漸大于靜態導線的覆冰質量.

圖17 不同時間導線表面液滴收集效率Fig.17 Droplet collection efficiency on the surface of the wire at different times

圖18 扭轉覆冰與未扭轉覆冰質量隨時間變化曲線Fig.18 Curve of torsional icing and untorsional icing mass with time
本文通過考慮覆冰偏心、風速等對導線的耦合扭轉特性及覆冰形狀對導線扭轉剛度的影響,建立了導線動態扭轉覆冰模型.與文獻[17]實驗數據對比,驗證了模型的有效性.在此基礎上研究了導線傾角和溫度、風速等環境參數對冰形和扭轉的影響.得到的主要結論如下:
1)扭轉角度為140°時,數值模擬得到的覆冰面積較文獻[17]大16.8%;扭轉角度為291°時,覆冰面積小7.54%。覆冰早期,導線扭轉受風載荷力矩影響較大,隨著覆冰質量提升,冰載荷力矩成為導線扭轉的主要因素.
2)受導線傾斜角度的影響,未凍結水滴在重力、氣流和離心力的作用下向導線下游移動.覆冰 245 min 時,傾斜導線覆冰面積較水平導線(α=0°)增大15.16%(傾角30°)、20.50%(傾角45°),26.30%(傾角60°).隨傾斜角度的增加,冰形呈“溝壑狀”.
3)溫度降低、風速及MVD 增大提高了導線表面水滴收集效率,單位時間內的覆冰質量增加,導線的扭轉速度加快.風速增大還進一步提高了風載荷力矩對導線扭轉的影響.
4)覆冰扭轉后的導線等效直徑增大,則水滴捕獲面積增大.隨著導線的扭轉,壁面水滴收集效率相比覆冰30 min 時分別增加6.38%(覆冰60 min)和11.4%(覆冰90 min).扭轉導線單位長度覆冰質量逐漸大于靜態導線.