新月形薄覆冰導(dǎo)線氣動力特性數(shù)值模擬研究

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

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

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

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

Masataka[4]對輸電線路124次舞動過程進行了實際觀測,歸納了不同覆冰形狀下對應(yīng)的舞動次數(shù),如表1所示,表中D為導(dǎo)線直徑?？梢钥闯?在實際的舞動中導(dǎo)線新月形覆冰情況占據(jù)較大比重;對于新月形覆冰而言,0～0.5D的薄覆冰情況下導(dǎo)線舞動次數(shù)達總次數(shù)的92%。說明新月形薄覆冰對導(dǎo)線舞動起著關(guān)鍵作用。

表1 不同覆冰類型對應(yīng)的舞動次數(shù)Tab.1 Number of galloping corresponding to different ice accretion types

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

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

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

1 試驗概況

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

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

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

(1)

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

2 圓柱繞流數(shù)值模擬

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

圖2 新月形單導(dǎo)線模型試驗照片F(xiàn)ig.2 Crescent single conductor model test

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

2.1 圓柱繞流計算基本參數(shù)

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

(a) mesh1(優(yōu)質(zhì)的)

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

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

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

2.2 圓柱繞流計算工況

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

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

2.3 圓柱繞流計算結(jié)果

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

定義測點的風壓系數(shù)如下

(2)

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

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

(a) 二維計算結(jié)果

(b) 展向長度1D

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

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

(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可知,對圓柱流向中心線平均風速來說:① 在圓柱尾流近壁區(qū)(圖中x/D<1.0區(qū)域),不同工況下的計算結(jié)果與風洞試驗結(jié)果均比較吻合;② 在圓柱遠離壁面尾流區(qū)(圖中x/D>1.0區(qū)域),不同工況下的計算結(jié)果與風洞試驗結(jié)果吻合度不盡相同?；趍esh1網(wǎng)格的計算結(jié)果優(yōu)于mesh2網(wǎng)格的結(jié)果;亞格子模型對計算精度的影響較大。展向長度1D與3D時所得計算結(jié)果無明顯差別。

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

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

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

表3 圓柱繞流氣動力統(tǒng)計值Tab.3 Aerodynamic statistics of cylinder flow

(3)

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

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

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

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

3 新月形覆冰導(dǎo)線氣動力特性數(shù)值模擬

3.1 計算參數(shù)

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

(a) 計算域及邊界條件

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

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

3.2 氣動力特性結(jié)果

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

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

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

圖12 覆冰單導(dǎo)線氣動三分力系數(shù)計算結(jié)果Fig.12 Results of aerodynamic three component coefficient of iced single conductor

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

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

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

(a) α=0°

(b) α=10°

(c) α=20°

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

4 結(jié) 論

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

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

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

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