新月形薄覆冰導線氣動力特性數值模擬研究

林 巍, 肖志斌, 樓文娟, 黃銘楓2,, 張躍龍

(1. 浙江大學 建筑設計研究院有限公司, 杭州 310027; 2. 浙江大學 平衡建筑研究中心, 杭州 310027;3. 浙江大學 結構工程研究所, 杭州 310058)

在特殊的氣候條件下,輸電導線可能會在其迎風面結冰,從而改變其截面形狀使受力變得復雜。在一定的風速和攻角條件下,導線容易產生舞動。舞動是一種低頻率(0.1～3 Hz)、大振幅(導線直徑的5倍～300倍)的自激振動。輸電線路的覆冰和積雪以及由此衍生的導線舞動已經嚴重威脅著電力及通信網絡的安全運行。

導線舞動的機理十分復雜,目前國際上比較公認的是Den Hartog[1]橫向舞動機理和Nigol[2]扭轉舞動機理。無論采用何種舞動機理,覆冰輸電導線的氣動力特性是導線舞動分析的重要參數。因此,研究覆冰輸電導線的氣動力特性對輸電導線舞動的研究及防治具有重要的意義[3]。

Masataka[4]對輸電線路124次舞動過程進行了實際觀測,歸納了不同覆冰形狀下對應的舞動次數,如表1所示,表中D為導線直徑。可以看出,在實際的舞動中導線新月形覆冰情況占據較大比重;對于新月形覆冰而言,0～0.5D的薄覆冰情況下導線舞動次數達總次數的92%。說明新月形薄覆冰對導線舞動起著關鍵作用。

表1 不同覆冰類型對應的舞動次數Tab.1 Number of galloping corresponding to different ice accretion types

風洞試驗和數值模擬是獲取覆冰導線氣動力特性的主要手段。Nigol在均勻流場風洞中測量了新月型覆冰導線的空氣動力參數;Chadha等[5]在風洞中測量了覆冰導線的三分力系數,考慮了湍流度對氣動力的影響; Chabart等[6]在對新月形覆冰截面進行氣動力特性試驗,得到了新月形覆冰導線的三分力系數及其擬合曲線;Ishihara[7]對 3種不同覆冰厚度的新月形覆冰的單導線及四分裂導線進行氣動力特性試驗,給出了均勻流場下的氣動力特性;張宏雁等[8]通過試驗研究新月形和扇形覆冰導線在不同風速下的氣動力特性;顧明等[9]通過節段模型高頻天平測力風洞試驗, 計算了準橢圓形和扇形兩種代表性覆冰導線的氣動力系數均值、均方根值;馬文勇等[10]也開展了準橢圓覆冰導線的氣動特性試驗研究。然而,由于覆冰的多樣性,目前有關覆冰導線的氣動力參數仍缺乏系統全面的數據。

姚育成[11]在高雷諾數下,在二維空間對不同冰厚的新月形覆冰導線的氣動力特性進行了數值模擬;呂翼等[12]通過對新月形和扇形覆冰的單導線、三分裂導線的氣動力特性進行了二維數值模擬;李新民等[13]對幾種典型覆冰導線的氣動力特性進行了二維數值模擬;韋遠武等[14]結合重疊網格技術實現覆冰導線靜態與動態氣動力模擬。李彭舉等[15]對新月形覆冰八分裂導線的氣動力特性進行了二維數值模擬。但已有的相關數值模擬研究大多把流域簡化為二維流場,其結果只能定性的反映氣動力隨攻角的變化趨勢,具體數值與試驗結果還有較大的差異。

圓柱繞流是空氣動力學中經典的繞流問題,國內外學者對此進行了大量的試驗研究和數值計算,已有結果包含的風速、風壓結果信息豐富,相對較為成熟。Breue[16]采用大渦模擬(LES)方法對孤立圓柱繞流問題(Re=140 000)進行了三維的數值模擬,通過不同離散格式和亞格子模型的比較,表明選擇合適的近壁處理和亞格子模型能和試驗結果吻合較好。 Kravchenkoa[17]采用動力亞格子模型對亞臨界狀態下圓柱繞流問題進行了三維數值模擬。輸電導線由于其直徑較小,其雷諾數基本處于亞臨界范圍,本文首先通過對圓柱繞流問題的LES數值模擬,考察近壁面網格對風速、風壓等結果的影響, 并同時驗證了采用k-ωSST 湍流模型的可行性。在此基礎上對新月形覆冰單導線的氣動力特性進行數值模擬,并與風洞試驗結果對比。為今后進行其他覆冰形狀氣動力數值計算提供參考。

1 試驗概況

導線覆冰情況極其復雜,具有很強的隨機性。 Masataka對124次舞動觀測中的覆冰形狀進行了歸類,發現在實際舞動的觀測中新月形覆冰占有很大的比重。新月形是輸電導線主要的覆冰形式且截面形狀易于量化,常作為氣動力特性的研究對象。本文試驗裸導線直徑采用LGJ-400/35(直徑26.8 mm),覆冰厚度為0.25D,D為裸導線直徑。試驗在浙江大學邊界層風洞(ZD-1)中完成,試驗來流風速取10 m/s,由截面的對稱性,試驗風攻角范圍取0°～180°,風向角間隔為5°。測力采用高頻測力天平,水平力量程為20 N,扭矩量程為 4 N·m,天平最高測力頻率1 000 Hz。圖2為模型的風洞試驗照片,試驗模型和設備詳見文獻[18]。

覆冰導線的氣動力可表示為升力FL、阻力FD和扭矩M,如圖1所示。并將氣動力無量綱化可得

圖1 氣動三分力方向及風向角定義Fig.1 Definition of aerodynamic force direction and wind angle

(1)

式中:CL、CD、CM分別為覆冰導線的升力系數、阻力系數和扭轉系數;U為試驗風速;ρ為試驗時空氣密度;D為裸導線的直徑;l為模型的有效長度。

2 圓柱繞流數值模擬

為驗證數值計算模型,本文將以圓柱為研究對象,分別進行二維和三維 LES 非定常繞流計算,并驗證采用k-ωSST 湍流模型的可行性。為與以往圓柱試驗文獻結果對比,雷諾數Re取3 900。計算流域及邊界條件,如圖3所示。流場區域足夠大,使尾流及渦脫落能充分發展。

圖2 新月形單導線模型試驗照片Fig.2 Crescent single conductor model test

圖3 圓柱繞流計算流域及邊界條件Fig.3 Calculation of basin and boundary conditions for flow around a cylinder

2.1 圓柱繞流計算基本參數

圓柱直徑D=2 cm,將整個流域劃分為9個子域,計算域網格劃分采用非均勻的結構網格,其中子域5為核心加密區,遠離圓柱表面的網格逐漸擴散。三維計算模型展向長度分別取3D和1D。為考察近壁面網格對計算結果的影響,綜合考慮y+值,x-y平面內網格采用了兩種不同密度的網格劃分方式:(1) mesh1:近壁面第一層網格尺度為0.003D,對應壁面y+<1,加密區網格以1.05的系數向外擴散,具體網格劃分見圖4(a)。(2) mesh2:近壁面第一層網格尺度為0.01D,對應壁面1

(a) mesh1(優質的)

(b) mesh2(粗糙的)圖4 圓柱繞流網格平面圖Fig.4 Grid plan of flow around a cylinder

數值計算的速度入口和初始條件如下:圓柱表面為無滑移邊界(u=v=0);速度入口為2.925 m/s均勻風速;上下邊界為對稱邊界(u=U,v=0);出口條件為壓力出口。

在非定常計算前,先采用基于RANS的Realizablek-ω湍流模型進行圓柱繞流定常計算,使其流場達到穩定的狀態,并以此作為非定常繞流計算的初始條件,這樣可以使圓柱尾流區較早進入充分發展階段。非穩態計算時間步長取0.000 5 s。

2.2 圓柱繞流計算工況

為考察亞格子模型和近壁面網格對數值結果的影響,分別選用標準Smagorinsky模型和Dynamic-Smagorinsky模型進行比較。通過計算可以發現,三維大渦模擬非定常繞流非常耗時。如果在現在普通8核微機上計算,要獲得一個統計意義上穩定的時程,計算時長往往在100 h以上甚至更長。對覆冰導線來說,需要計算不同風攻角下的氣動力系數,這樣長的計算時間是無法接受的。為此,同時選擇k-ωSST湍流模型進行比較。主要計算工況如表2所示。

表2 圓柱繞流計算工況表Tab.2 Calculation condition of flow around a cylinder

2.3 圓柱繞流計算結果

對表2中所有工況進行圓柱繞流CFD非定常計算,并將本文計算結果與以往風洞試驗結果[19-21]和數值模擬結果進行對比分析,主要針對平均風壓系數、平均速度場及斯托勞哈數(St)等。

定義測點的風壓系數如下

(2)

式中:Pi為測點的壓力;ρ為空氣密度;U∞為來流風速。

圖5為本文LES數值模擬所得的平均風壓系數與以往風洞試驗結果的比較,圖中“expRe=3 900”表示Re=3 900時文獻[19]中的試驗結果。在平均風壓系數方面:① 如圖5(a)所示,基于兩種網格及不同亞格子模型下,二維LES計算結果與風洞試驗結果差別都較大。主要在圓柱側面的分離區和背后的尾流區差別十分嚴重。幾種工況下結果差別不大;② 由圖5(b)、(c)可見,表面的壓力系數在正對來流處為1,隨著向周圍擴展,壓力系數值減小,達到最小值后很快增大到一個較為穩定的值并在圓柱背后形成平坦的穩態壓力分布,變化趨勢與試驗結果吻合。mesh1網格下計算結果優于mesh2網格,Dynamic-Smagorinsky亞格子模型計算結果優于Smagorinsky下計算結果。展向長度取1倍特征長度時已顯示很好的計算精度,展向長度取3倍特征長度對計算精度改善不明顯。

(a) 二維計算結果

(b) 展向長度1D

(c) 展向長度3D圖5 平均風壓系數比較Fig.5 Comparison of average wind pressure coefficient

LES計算所得圓柱流向(x)中心線(y=0,z=H/2)的平均速度場和橫風向(y)近壁區(x/D=1.54,z=H/2)脈動速度場的比較結果如圖6 ～圖8所示,U0為入口來流速度,u′,v′分別為x,y方向的脈動速度。從風壓比較圖中知,二維計算結果與試驗結果差別較大,因此對二維速度場計算結果此處不再列出。

(a) 展向長度1D

(b) 展向長度3D圖6 x流向平均速度比較Fig.6 Comparison of mean velocity in x direction

(a) 展向長度1D

(b) 展向長度3D圖7 x向脈動速度比較(x/D=1.54)Fig.7 Comparison of pulsating velocity in x direction(x/D=1.54)

(a) 展向長度1D

(b) 展向長度3D圖8 y向脈動速度比較(x/D=1.54)Fig.8 Comparison of pulsating velocity in y direction(x/D=1.54)

由圖6可知,對圓柱流向中心線平均風速來說:① 在圓柱尾流近壁區(圖中x/D<1.0區域),不同工況下的計算結果與風洞試驗結果均比較吻合;② 在圓柱遠離壁面尾流區(圖中x/D>1.0區域),不同工況下的計算結果與風洞試驗結果吻合度不盡相同。基于mesh1網格的計算結果優于mesh2網格的結果;亞格子模型對計算精度的影響較大。展向長度1D與3D時所得計算結果無明顯差別。

由圖7可知,對于圓柱x/D=1.54處x向脈動風速來說,除1D_mesh2_dsma工況外,其他各工況下計算結果差別不大,與風洞試驗結果吻合較好。

由圖8可知,對于圓柱x/D=1.54處y向脈動風速來說,在圓柱特征長度以外(圖中y/D<-0.5和y/D>-0.5),各工況所得計算結果無明顯差別,與風洞試驗結果吻合較好。但對圓柱特征長度范圍內(圖中-0.5

各工況下圓柱氣動力系數統計值與風洞試驗和以往數值模擬的比較如表3所示。表中St為斯托勞哈數,定義為

表3 圓柱繞流氣動力統計值Tab.3 Aerodynamic statistics of cylinder flow

(3)

式中,f為渦脫頻率,Hz。升力系數和阻力系數均用來流平均風速無量綱化。CD,mean表示阻力系數均值。表中Cp,back為圓柱尾部點處風壓系數。

如圖9所示是計算時間100D/U0(約22個渦脫周期) 后1D_mesh1_dsma工況和k-ωSST工況下升、阻力系數的時程曲線。由時程曲線上可知,除大渦的周期脫落外還包括小渦的脫落(在圖中表現為高頻脈動部分),這是與k-ωSST、k-ω等時均湍流模型結果只含有低頻部分不同之處。因此,LES模型需要更長的計算時間來獲得如表3中具有統計意義的時均量。

(b) SST模型圖9 圓柱升、阻力系數時程Fig.9 Lift and drag coefficient time history of the cylinder

通過對風壓、速度、氣動力及尾流渦脫頻率等的比較,對于圓柱繞流來說,近壁面網格的疏密對LES計算結果有較大影響,近壁網格越密計算結果越精確。Dynamic-Smagorinsky亞格子模型因為CS不再固定具有更好的精度;對平均氣動力而言,二維網格過高的估計了圓柱的平均阻力系數,這與圖5(a)中圓柱負壓區風壓系數數值偏小有關;但展向長度取1倍特征長度的三維模型所得計算結果已經足夠精確。k-ωSST模型對此類鈍體繞流問題也具有較高的精度,考慮到計算效率問題,遇到工況較多時,采用k-ωSST湍流模型是可行的。

3 新月形覆冰導線氣動力特性數值模擬

3.1 計算參數

覆冰導線截面尺寸同風洞試驗,覆冰0.25D。為再次考察展向長度對導線氣動力的影響,對比二維模型與三維模型的區別。計算時取二維模型和展向長度分別為1倍導線直徑(1D)和3倍導線直徑(3D)的三維模型。計算域大小取為35D×20D×1D,如圖10(a)所示。采用區域分塊非均勻結構化網格,附面層網格與外部網格間為加密的非結構網格過渡區,展向劃分為10份。為保證壁面附近y+≈1,最小網格尺寸為0.000 8D。風攻角定義如圖1所示,0°風攻角下網格劃分平面圖如圖10(b)、(c)所示,總網格數約23萬。

(a) 計算域及邊界條件

(c) 局部網格圖10 邊界條件、計算域和離散網格平面示意圖Fig.10 Schematic diagram of boundary conditions, calculation basin and grid plan

采用速度入口邊界條件,來流風速為10 m/s,湍流度為5%;出口為壓力出口邊界條件;y向和z向側面為對稱邊界條件;導線表面為無滑移壁面,如圖10(a)所示。在進行非穩態計算前,先進行穩態計算,然后將穩態計算的流場作為k-ωSST非定常計算的初始流場。非穩態時間步長取0.000 5 s。

3.2 氣動力特性結果

覆冰導線的三分力系數定義如前文所述,取計算所得時程穩定后的若干個周期的平均值,得到覆冰導線平均三分力系數。20°攻角下模型展向長度取3D的氣動三分力系數時程,如圖11所示。

圖11 20°攻角下氣動三分力系數時程Fig.11 Time history of aerodynamic three component force coefficient at 20 °angle of attack

圖12給出了新月形覆冰單導線氣動三分力系數數值模擬與文獻[18]風洞試驗結果的對比。

圖12 覆冰單導線氣動三分力系數計算結果Fig.12 Results of aerodynamic three component coefficient of iced single conductor

由圖12可以看出,三維模型計算所得的阻力系數的數值模擬結果與風洞試驗結果在各攻角下符合較好。升力系數曲線在兩個尖峰對應的攻角處存在一定的誤差,展向長度3D優于1D結果,其他攻角下與試驗結果較為一致。阻力系數曲線展向長度取3D與1D相差不大。扭轉系數的數值模擬結果與風洞試驗結果存在較大差異。主要原因在于,該薄覆冰模型的扭矩量級很小,小量級下本次測力天平已不能精確反應扭轉系數的大小,因此,隨著攻角的變化,測力結果已不能反應扭轉系數值的變化。而數值模擬恰好能彌補這一缺點。

總的來說,三維計算模型的計算結果與測力結果是比較一致的,雖然個別攻角下存在一定誤差,但計算結果仍能較好的反應氣動力的隨攻角的變化。二維計算模型與試驗結果差別較大,采用三維模型是必要的。展向長度取3D略優于取1D的計算結果。

展向長度取3D時計算的平均風壓分布結果如圖13所示。由圖可知,風壓分布對風攻角的變化非常敏感,當α=0°時,覆冰導線截面上面對稱,表面風壓也呈對稱分布,升力系數平均值等于零;當α=10°和α=20°時,上下表面風壓不再對稱,上表面負壓大于下表面,這也是升力系數大于零的原因。

(a) α=0°

(b) α=10°

(c) α=20°

(d) 平均風壓系數分布圖13 單導線表面平均風壓分布結果Fig.13 Average wind pressure distribution on the surface of single conductor

4 結 論

本文通過CFD數值模擬方法對雷諾數為3 900的經典圓柱繞流問題進行了非定常計算,通過不同計算工況下平均風壓、平均風速、脈動風速、斯托勞哈數(St)等參數與以往試驗結果的對比,研究了近壁面網格及展向長度對結果的影響。在此基礎上,通過k-ωSST湍流模型計算了新月形薄覆冰單導線二維和三維網格模型的氣動三分力系數。

(1) 近壁面網格的疏密對風壓系數、氣動力和流場的計算結果有較大影響,近壁網格越密計算結果越精確。通常要求壁面y+值在1左右。

(2) 圓柱繞流計算表明,一個完整的大渦脫落周期內將出現一次升力的最大值和兩次阻力的最大值。從精度上講,采用壁面y+<1的網格并配合使用Dynamic-Smagorinsky亞格子模型可以得到很好的鈍體繞流大渦模擬計算結果,但若僅以平均氣動力為目標,采用k-ωSST計算是可行且高效的。

(3) 薄覆冰單導線的計算結果與試驗結果對比表明,二維與三維模型存在較大區別,采用三維計算模型是必要的。展向長度取3倍導線直徑能取得較好的計算結果。