樓文娟，潘 晨，孫建平

(浙江大學(xué) 結(jié)構(gòu)工程研究所，杭州 310058)

0 引 言

架空輸電線路導(dǎo)線在表面覆冰后會(huì)形成非圓形的截面形狀，其風(fēng)致氣動(dòng)力可誘發(fā)導(dǎo)線舞動(dòng)，嚴(yán)重威脅輸電線路的安全[1-2]。

李萬(wàn)平[3-4]、樓文娟[5]、閻東[6]、李新民[7]等對(duì)新月形、D形覆冰導(dǎo)線進(jìn)行風(fēng)洞試驗(yàn)研究全風(fēng)攻角下氣動(dòng)力系數(shù)。Hisato[8]等對(duì)不同厚度的三角覆冰和圓角三角覆冰導(dǎo)線進(jìn)行表面測(cè)壓試驗(yàn)和天平測(cè)力試驗(yàn)，驗(yàn)證兩種試驗(yàn)結(jié)果的一致性，并獲得了覆冰導(dǎo)線氣動(dòng)力數(shù)據(jù)。

基于計(jì)算流體力學(xué)（CFD）的數(shù)值模擬方法能夠展現(xiàn)完備的流場(chǎng)信息，可與風(fēng)洞試驗(yàn)互補(bǔ)，便于更深層次地研究流場(chǎng)結(jié)構(gòu)的演化。李新民[9]用SA湍流模型對(duì)23 mm厚新月形覆冰導(dǎo)線進(jìn)行數(shù)值模擬，發(fā)現(xiàn)在風(fēng)攻角15°附近存在回流渦現(xiàn)象。Ishihara[10]利用大渦模擬（LES）對(duì)圓角三角形覆冰導(dǎo)線進(jìn)行模擬，模擬結(jié)果與風(fēng)洞試驗(yàn)數(shù)據(jù)相吻合。林巍[11]、呂翼[12]分別用k-ωSST湍流模型和SA湍流模型對(duì)覆冰6.7 mm（導(dǎo)線直徑26.82 mm）和覆冰18 mm （導(dǎo)線直徑32.76 mm）的新月形覆冰導(dǎo)線進(jìn)行數(shù)值模擬，并得到了較好的模擬結(jié)果。符玉珊[13]運(yùn)用格子玻爾茲曼方法對(duì)18 mm覆冰的新月形覆冰導(dǎo)線（導(dǎo)線直徑32.76 mm）進(jìn)行大渦模擬，在較低耗時(shí)情況下得到了較高精度的模擬結(jié)果。

現(xiàn)有對(duì)覆冰導(dǎo)線氣動(dòng)力特性的數(shù)值模擬研究中，覆冰厚度多數(shù)小于1D（D為導(dǎo)線直徑），而根據(jù)歷史觀測(cè)數(shù)據(jù)，覆冰厚度達(dá)到1D及以上的情況在覆冰嚴(yán)重地區(qū)經(jīng)常出現(xiàn)[14-15]。圖1根據(jù)文獻(xiàn)[5-6]的風(fēng)洞試驗(yàn)數(shù)據(jù)繪制，為不同厚度覆冰導(dǎo)線升力系數(shù)隨風(fēng)攻角變化情況的示意圖。當(dāng)覆冰厚度較大時(shí)，覆冰導(dǎo)線的升力系數(shù)曲線在風(fēng)攻角15°附近區(qū)間出現(xiàn)“尖峰”，并且隨著冰厚的增大更加明顯。升力系數(shù)的突變使得在該風(fēng)攻角區(qū)間易出現(xiàn)覆冰導(dǎo)線的風(fēng)致舞動(dòng)現(xiàn)象[6,16-17]，準(zhǔn)確模擬和理解這一升力系數(shù)的“尖峰”，對(duì)于覆冰導(dǎo)線風(fēng)致舞動(dòng)的驗(yàn)算和防治具有特別意義。

圖1 不同厚度新月形覆冰升力系數(shù)（I為湍流度）Fig. 1 Lift coefficients of the crescent-shape ice with different thickness (I indicates the turbulence intensity)

目前覆冰導(dǎo)線氣動(dòng)力系數(shù)的CFD模擬大多局限于無(wú)升力“尖峰”的較薄覆冰，通常選用雷諾時(shí)均法進(jìn)行氣動(dòng)力模擬。但對(duì)于厚覆冰導(dǎo)線的模擬結(jié)果并不理想，尤其是風(fēng)攻角15°附近升力系數(shù)“尖峰”模擬較為困難；另外，對(duì)于厚覆冰導(dǎo)線升力系數(shù)在風(fēng)攻角15°附近出現(xiàn)極大值的原因也尚不明確，鮮有討論。本文分別采用ANSYS Fluent軟件內(nèi)基于k-ωSST湍流模型的雷諾時(shí)均法和大渦模擬的數(shù)值方法，通過(guò)網(wǎng)格和時(shí)間步長(zhǎng)的優(yōu)化提高模擬精度，得到風(fēng)攻角10°～20°范圍內(nèi)1D新月形厚覆冰導(dǎo)線氣動(dòng)力參數(shù)和繞流風(fēng)場(chǎng)結(jié)果，并利用風(fēng)洞試驗(yàn)結(jié)果驗(yàn)證數(shù)值方法的適用性。通過(guò)對(duì)覆冰導(dǎo)線表面壓力分布和繞流風(fēng)場(chǎng)結(jié)構(gòu)分析，進(jìn)一步探究厚覆冰導(dǎo)線升力系數(shù)在風(fēng)攻角15°附近出現(xiàn)“尖峰”的原因。

1 數(shù)值方法和計(jì)算模型

1.1 大渦模擬

大渦模擬通過(guò)空間濾波將湍流分成大尺度渦和小尺度渦。其中，大尺度渦對(duì)湍流輸運(yùn)的動(dòng)量、質(zhì)量、能量等物理量起到?jīng)Q定性作用，因此對(duì)大尺度的渦進(jìn)行了不可壓縮Navier-Stokes方程直接模擬。小尺度的渦趨于各向同性，其對(duì)大尺度渦的影響通過(guò)亞格子模型模擬。

本文的大渦模擬選用局部渦黏度壁面自適應(yīng)（WALE）亞格子模型。WALE亞格子模型是在Smagorinsky-Lilly亞格子模型的基礎(chǔ)上發(fā)展過(guò)來(lái)的，它修正了Smagorinsky-Lilly亞格子模型在近壁面渦黏系數(shù)不為0的缺陷[18]。其亞格子應(yīng)力張量具體形式如下：

其中，μt為亞格子黏度系數(shù)，其定義式如下：

其中，

式中，Cw為WALE常數(shù)，取0.325；Δ為過(guò)濾之后的尺寸，為網(wǎng)格尺寸。

1.2 基于k-ω SST模型的雷諾時(shí)均法

雷諾時(shí)均法（RANS）將流場(chǎng)分為系綜平均量和湍流脈動(dòng)量，并通過(guò)湍流模型來(lái)封閉N-S方程。其中，基于湍流動(dòng)能k和比耗散率ω的k-ω模型[19]是一種較為常用的湍流模型。k-ω模型在湍流黏度的定義中考慮了湍流剪切應(yīng)力的傳遞，并通過(guò)設(shè)定混合函數(shù)將標(biāo)準(zhǔn)k-ε模型和k-ω模型進(jìn)行結(jié)合，在近壁面區(qū)域內(nèi)使用k-ω湍流模型求解，在邊界層外則使用k-ε模型。既能充分發(fā)揮k-ε模型在模擬遠(yuǎn)離壁面的充分發(fā)展湍流的優(yōu)勢(shì)，又能適應(yīng)于各種壓力梯度下的邊界層模擬，在模擬壓力梯度和分離流上具有較大優(yōu)勢(shì)，常被用于模擬鈍體繞流。

1.3 計(jì)算模型

本文研究對(duì)象為1D新月形覆冰導(dǎo)線，裸導(dǎo)線直徑D為26.82 mm。計(jì)算域入口邊界至導(dǎo)線中心距離為30D，出口邊界距導(dǎo)線中心50D；計(jì)算域的橫風(fēng)向?qū)挾葹?0D，展向高度為2D，最大阻塞率為2.0%。計(jì)算域入口設(shè)為速度入口（入口風(fēng)速U= 10 m/s），出口為自由出流，覆冰導(dǎo)線表面為光滑無(wú)滑移壁面，展向邊界為周期性邊界，其余為對(duì)稱邊界。幾何模型和計(jì)算域如圖2所示。計(jì)算網(wǎng)格為使用Pointwise生成的結(jié)構(gòu)網(wǎng)格，近壁面第一層網(wǎng)格高度設(shè)為0.001D，以滿足LES計(jì)算時(shí)y+≈1要求，網(wǎng)格增長(zhǎng)率為1.05。圖3為計(jì)算網(wǎng)格局部示意圖。

圖2 幾何模型和計(jì)算域Fig. 2 Schematic of the model and the computational domain

圖3 近壁面網(wǎng)格Fig. 3 Near-wall grid

大渦模擬采用SIMPLEC算法求解壓力-速度耦合方程，壓強(qiáng)項(xiàng)為二階格式，動(dòng)量項(xiàng)為有界中心差分格式，時(shí)間離散方式為二階隱式格式。在使用基于k-ωSST模型的雷諾平均法在求解黏性子層時(shí)同樣需要滿足y+≈1的網(wǎng)格要求，選用與LES相同的計(jì)算網(wǎng)格。k-ωSST模型計(jì)算時(shí)湍動(dòng)能和耗散率項(xiàng)均設(shè)置為二階迎風(fēng)格式，其余和LES設(shè)置相同。

1.4 網(wǎng)格與時(shí)間步長(zhǎng)無(wú)關(guān)性驗(yàn)證

為分析網(wǎng)格數(shù)和時(shí)間步長(zhǎng)對(duì)大渦模擬計(jì)算結(jié)果的影響，分別計(jì)算了風(fēng)攻角15°下的4個(gè)算例，結(jié)果列于表1。

表1 不同計(jì)算工況設(shè)置Table 1 Parameters for different computation conditions

由表1可得，周向節(jié)點(diǎn)數(shù)對(duì)氣動(dòng)力模擬結(jié)果有較大影響。隨覆冰厚度增大，需要更多節(jié)點(diǎn)來(lái)捕捉覆冰導(dǎo)線壁面弧度細(xì)節(jié)。針對(duì)本文工況，600個(gè)周向節(jié)點(diǎn)能較好表現(xiàn)覆冰導(dǎo)線截面形狀特性。綜合考慮計(jì)算精度與耗時(shí)，本文LES選用Case2的網(wǎng)格與時(shí)間步長(zhǎng)設(shè)置；k-ωSST湍流模型選取同樣的網(wǎng)格，但由于本文算例中k-ωSST湍流模型對(duì)時(shí)間步設(shè)置不敏感，因此參考文獻(xiàn)[10]中的取值，無(wú)量綱時(shí)間步設(shè)置為0.186。

2 結(jié)果對(duì)比及分析

2.1 氣動(dòng)力系數(shù)

覆冰導(dǎo)線升力系數(shù)、阻力系數(shù)以及扭矩系數(shù)定義為：

式中，空氣密度ρ為1.225 kg/m3，導(dǎo)線外徑D為26.82 mm，展向長(zhǎng)度H為2D，U為試驗(yàn)風(fēng)速。升力FL、阻力FD、扭矩M的方向如圖4所示。

圖4 氣動(dòng)力方向及風(fēng)攻角定義Fig. 4 Definition of the aerodynamic force directions and wind attack angles

圖5為風(fēng)攻角10°～20°范圍內(nèi)1D新月形覆冰導(dǎo)線的氣動(dòng)力系數(shù)。對(duì)比k-ωSST和LES的計(jì)算氣動(dòng)力系數(shù)可以發(fā)現(xiàn)，LES模擬結(jié)果與試驗(yàn)結(jié)果非常接近，在風(fēng)攻角15°下，升力系數(shù)、扭轉(zhuǎn)系數(shù)均出現(xiàn)“尖峰”，并且對(duì)角度一階求導(dǎo)后得到的Den Hartog系數(shù)和Nigol系數(shù)（圖6）在負(fù)斜率部分也和試驗(yàn)值非常接近；而k-ωSST模擬結(jié)果在風(fēng)攻角15°～20°范圍內(nèi)升力系數(shù)沒(méi)有降低，沒(méi)有出現(xiàn)明顯的“尖峰”。因此，LES模型比k-ωSST模型更能準(zhǔn)確模擬覆冰導(dǎo)線的氣動(dòng)力，能滿足覆冰導(dǎo)線舞動(dòng)分析的需求。

圖5 各風(fēng)攻角下1D新月形覆冰導(dǎo)線氣動(dòng)力系數(shù)Fig. 5 Aerodynamic coefficients of a 1D crescent-shape iced conductor at various wind attack angles

圖6 1D新月形覆冰導(dǎo)線Den Hartog系數(shù)和Nigol系數(shù)Fig. 6 Den Hartog and Nigol coefficients of a 1D crescent-shape iced conductor

2.2 平均流場(chǎng)特性

根據(jù)1D新月形覆冰導(dǎo)線表面的壓力分布特征，將其分為四個(gè)區(qū)域：覆冰端上側(cè)（0 ＜x/D＜ 1.5）、導(dǎo)線端上側(cè)（1.5 ＜x/D＜ 2）、導(dǎo)線端下側(cè)（2 ＜x/D＜ 2.5）和覆冰端下側(cè)（2.5 ＜x/D＜ 4），并分別將其記為Ⅰ區(qū)、Ⅱ區(qū)、Ⅲ區(qū)、Ⅳ區(qū)，如圖7所示。

圖7 覆冰導(dǎo)線分區(qū)及周向坐標(biāo)定義Fig. 7 Partitions of the iced conductor with the circumferential coordinates

圖8為兩種數(shù)值方法計(jì)算得到的平均流線圖。整體而言，兩種算法的平均流場(chǎng)在Ⅱ區(qū)和Ⅲ區(qū)均有較大尺度的回流，Ⅱ區(qū)逆時(shí)針大尺度渦緊貼著導(dǎo)線壁面，Ⅲ區(qū)渦為上側(cè)流體繞過(guò)Ⅱ區(qū)渦后再次附著在導(dǎo)線表面的二次渦，且Ⅱ區(qū)回流區(qū)域面積大于Ⅲ區(qū)回流區(qū)域面積。

k-ωSST模型計(jì)算的平均回流區(qū)長(zhǎng)度在風(fēng)攻角10°與12.5°較為接近，在15°附近最大，并在15°至20°范圍內(nèi)隨風(fēng)攻角增大而減小。LES計(jì)算得到的平均回流區(qū)長(zhǎng)度則在10°～15°下沒(méi)有明顯變化，在15°～20°下隨風(fēng)攻角增大而明顯增大。

相比較k-ωSST模型而言，LES能夠模擬出Ⅰ、Ⅱ區(qū)緊貼壁面逆時(shí)針的分離泡。隨著風(fēng)攻角的增大，這些分離泡的位置逐漸前移靠近導(dǎo)線覆冰端頂點(diǎn)，覆冰導(dǎo)線上側(cè)分離點(diǎn)也向覆冰頂端靠近，且在Ⅲ區(qū)存在小尺度渦結(jié)構(gòu)。Ⅲ區(qū)壁面附近渦結(jié)構(gòu)僅在15°風(fēng)攻角下為一對(duì)方向相反的小尺度渦，其余風(fēng)攻角下僅存在一個(gè)順時(shí)針渦。

2.3 表面風(fēng)壓系數(shù)

從圖5可看到，k-ωSST模型和LES模擬的升力系數(shù)在15°～20°有明顯區(qū)別，圖9為該范圍覆冰導(dǎo)線壁面平均風(fēng)壓系數(shù)分布圖。由圖9可得，1D新月形覆冰導(dǎo)線表面平均風(fēng)壓系數(shù)在Ⅱ、Ⅲ區(qū)均保持在0附近，靠近導(dǎo)線覆冰端頂點(diǎn)的平均壓力系數(shù)存在大幅度突變。平均壓力系數(shù)均在覆冰頂端的上表面和下表面附近（約x/D= 0.1與x/D= 3.9處）達(dá)到最小值和最大值。k-ωSST模型和LES計(jì)算結(jié)果差異主要體現(xiàn)在Ⅰ區(qū)。

圖8 不同風(fēng)攻角下平均流線圖Fig. 8 Averaged streamline patterns under various wind attack angles

圖9 風(fēng)攻角15° ～ 20°下的平均風(fēng)壓系數(shù)Fig. 9 Averaged wind pressure coefficient at wind attack angles ranging from 15° to 20°

圖10為各風(fēng)攻角的平均速度和壓力云圖?？梢钥吹?，LES結(jié)果中最大平均速度都在15 m/s左右，k-ωSST模型計(jì)算得到的最大平均速度和LES基本一致，僅在20°風(fēng)攻角下突然增大至17 m/s。在各風(fēng)攻角下均存在兩個(gè)高流速區(qū)域（x/D= 0.1和x/D= 2.5附近），此處流場(chǎng)對(duì)覆冰導(dǎo)線的吸力最大，對(duì)應(yīng)壓力曲線上的兩個(gè)極小值。

在覆冰端高速區(qū)域內(nèi)，k-ωSST模型的高速區(qū)域面積更大，呈半圓形；LES模型的高速區(qū)域呈長(zhǎng)條狀，對(duì)應(yīng)于圖9所示的LES模型在覆冰端上側(cè)的平均壓力分布明顯小于k-ωSST模型。LES方法捕捉到的小尺度渦脫引起覆冰導(dǎo)線上側(cè)出現(xiàn)風(fēng)速介于最大速度和0之間的中速區(qū)，高速區(qū)向上抬升，減小了導(dǎo)線在此范圍內(nèi)的吸力。

k-ωSST和LES模擬結(jié)果存在區(qū)別，其主要原因是LES能夠模擬出緊貼覆冰導(dǎo)線上側(cè)壁面的多個(gè)小尺度渦，而這些小尺度渦對(duì)覆冰導(dǎo)線壁面壓力有著不可忽視的影響；k-ωSST模型僅模擬出了Ⅰ區(qū)和Ⅱ區(qū)的大尺度回流渦，無(wú)法模擬出Ⅰ區(qū)壁面上的小尺度渦，夸大了Ⅰ區(qū)壁面的風(fēng)吸力。并且隨著覆冰厚度的增大，覆冰端上側(cè)的小尺度渦結(jié)構(gòu)會(huì)更復(fù)雜，對(duì)覆冰導(dǎo)線的影響更大。因此對(duì)于較厚覆冰導(dǎo)線，LES是數(shù)值模擬的首選。

圖10 各風(fēng)攻角下的平均速度云圖和平均壓力系數(shù)Fig. 10 Averaged velocity contours and averaged wind pressure coefficients at various wind attack angles

3 升力系數(shù)突變?cè)?/h2>

由圖8可知，風(fēng)攻角10°至20°是覆冰導(dǎo)線尾流處的渦不斷擴(kuò)大并附著于覆冰導(dǎo)線上壁面的過(guò)程；并隨著覆冰厚度增大，Ⅰ、Ⅱ區(qū)的小尺度渦結(jié)構(gòu)更加復(fù)雜。在該過(guò)程中，覆冰導(dǎo)線周圍流場(chǎng)變化致使導(dǎo)線氣動(dòng)力出現(xiàn)較大變化。由圖5可知，LES方法能更有效地模擬1D新月形覆冰導(dǎo)線氣動(dòng)力系數(shù)隨風(fēng)攻角的變化情況，其表面壓力分布結(jié)果更準(zhǔn)確。下面根據(jù)LES方法模擬所得平均流場(chǎng)結(jié)果，分析升力系數(shù)在風(fēng)攻角15°出現(xiàn)“尖峰”的流體動(dòng)力學(xué)機(jī)制。

表2為各風(fēng)攻角下壁面四個(gè)區(qū)域的升力系數(shù)分量，從中可以看到覆冰導(dǎo)線升力主要來(lái)源于Ⅳ區(qū)正壓以及Ⅰ區(qū)負(fù)壓。對(duì)Ⅳ區(qū)，在10°至20°范圍內(nèi)隨著風(fēng)攻角的增大，Ⅳ區(qū)前沿駐點(diǎn)的位置遠(yuǎn)離覆冰頂點(diǎn)，駐點(diǎn)處的壁面曲率減小，來(lái)流與壁面夾角增大，根據(jù)動(dòng)量定律，壁面壓力增大；同時(shí)，隨著風(fēng)攻角的增大，Ⅳ區(qū)正壓沿升力方向的分量減小，沿阻力方向的分量增大。如表2所示，Ⅳ區(qū)正壓提供的升力在15°風(fēng)攻角下達(dá)到最大值。對(duì)Ⅰ區(qū)，在10°～20°范圍，隨著風(fēng)攻角增大，分離點(diǎn)向來(lái)流方向前移，并出現(xiàn)緊貼壁面的小尺度渦，其與尾流區(qū)域渦旋結(jié)構(gòu)的相互作用致使Ⅰ區(qū)壁面吸力先減后增。同時(shí)，在風(fēng)攻角15°下，覆冰導(dǎo)線Ⅲ區(qū)附近的一對(duì)方向相反的小尺度渦增大了沿壁面法向的壓力梯度，使Ⅲ區(qū)出現(xiàn)較大的正壓。綜合覆冰導(dǎo)線上側(cè)吸力略微減小、下側(cè)（Ⅲ、Ⅳ區(qū)）的壁面壓力因駐點(diǎn)所在處的壁面曲率、來(lái)流夾角和壁面切線方向等因素影響達(dá)到最大值，致使導(dǎo)線上下表面壓力差明顯大于其余風(fēng)攻角，1D厚覆冰導(dǎo)線的升力系數(shù)曲線在15°風(fēng)攻角附近出現(xiàn)“尖峰”。

表2 各風(fēng)向角下各壁面分區(qū)的升力系數(shù)分量Table 2 Components of the lift coefficient in each surface partition under various wind attack angles

4 結(jié) 論

本文分別采用基于k-ωSST湍流模型的雷諾時(shí)均法和大渦模擬的數(shù)值方法模擬1D新月形覆冰導(dǎo)線在風(fēng)攻角10°～20°范圍內(nèi)的氣動(dòng)力特性和繞流風(fēng)場(chǎng)，并結(jié)合風(fēng)洞試驗(yàn)數(shù)據(jù)對(duì)結(jié)果進(jìn)行對(duì)比分析。主要結(jié)論如下：

1）相比于覆冰導(dǎo)線氣動(dòng)力模擬中廣泛應(yīng)用的雷諾時(shí)均法，大渦模擬能夠準(zhǔn)確捕捉壁面附近的小尺度渦結(jié)構(gòu)，有較準(zhǔn)確的模擬結(jié)果。對(duì)于較厚覆冰導(dǎo)線的氣動(dòng)力數(shù)值模擬，應(yīng)首選LES。

2）1D新月形覆冰導(dǎo)線在上側(cè)壁面處的渦結(jié)構(gòu)影響整體流場(chǎng)，并在下側(cè)壁面曲率、來(lái)流夾角和壁面切線方向共同作用下導(dǎo)致風(fēng)攻角15°時(shí)的升力系數(shù)突變。

根據(jù)本文的模擬結(jié)果以及出現(xiàn)升力“尖峰”的原因闡釋，得到了較厚新月形覆冰的流場(chǎng)特征，但具體的防舞措施還需進(jìn)一步探究嘗試。