曹鵬程， 廖紹凱,2， 張 研， 陳 達

(1.嘉興學院 建筑工程學院，嘉興 314001； 2.河海大學 力學與材料學院，南京 211100；3.河海大學 港口海岸與近海工程學院，南京 211100)

1 引 言

數(shù)值模擬鈍體湍流繞流問題主要有直接數(shù)值模擬法(DNS法)、大渦模擬法(LES法)和基于雷諾時均方程模擬法(RANS法)三種方法[1]。相關(guān)數(shù)值算法的研究已成為眾多學者的研究熱點。

DNS法能夠提供精確的解決方案，但是所需網(wǎng)格必須小于或等于流場中最小的漩渦結(jié)構(gòu)尺寸，對計算機容量要求極高，計算量巨大，限制了其廣泛應(yīng)用[2]。LES法的計算效率相比DNS法要提高很多，已成功應(yīng)用于求解各種形狀的湍流繞流問題[3]，但近壁區(qū)域仍需要非常精細的網(wǎng)格。RANS法可在較低的計算成本下實現(xiàn)高雷諾數(shù)下的湍流計算[4]。與二方程RNGk-ε模型、修正的k-ε模型[5]和k-ω模型[6]相比，一方程S -A模型[7]具有求解方程較少的優(yōu)點，可進一步降低計算成本。

近年來，基于一方程S -A模型的RANS法已成功求解了圓柱繞流和方柱繞流問題[8,9]，并大量運用在覆冰輸電線繞流數(shù)值計算中[10,11]。目前求解RANS方程主要采用基于特征線分裂法(CBS法)和迎風有限元法(SUPG法)，以克服方程中的非線性項帶來的數(shù)值振蕩問題。值得注意的是，這些算法在時間離散上采用一階離散格式[8,9]，因此，探索高階離散格式和減少剛度矩陣的更新次數(shù)以提高計算精度、效率和收斂性，有待于進一步研究。

本文將基于一方程 S -A 模型封閉RANS方程，采用沿均勻流線的坐標變換，消除其非線性項，采用沿均勻流線的三階Runge -Kutta法進行時間離散，采用Galerkin法進行空間離散，提出求解湍流模型的有限元算法。通過經(jīng)典的算例，進一步驗證算法的有效性、精度和收斂性。

2 基于S -A模型的有限元格式

2.1 控制方程

二維非定常不可壓縮流的無量綱RANS方程和S -A模型的標量方程為

(1)

cb 1=0.1355，σ=2/3,cb 2=0.622

cw 2=0.3,cw 3=2，cv 1=7.1

由于式(1)含有非線性對流項，為非自伴隨算子，當采用標準的Galerkin有限元方法進行空間離散時，得不到最優(yōu)解。因此，采用經(jīng)時間間隔Δt沿流線的坐標變換[12]，得到無對流項表達式為

(2)

2.2 時間離散

(3,4)

為了提高計算精度，基于沿均勻流線的三階Runge -Kutta法對式(3)進行時間離散。假定已知tn時刻流場的速度和壓力，則在給定的時間步長Δt=tn + 1-tn，tn + 1時刻中間輔助速度為

(5)

(6)

(7)

(8)

(9)

利用沿均勻流線的Taylor展開，并忽略其中的高階小量，式(5～7)的相關(guān)動坐標系下的量轉(zhuǎn)化為靜坐標系下的量為

(10)

(11)

(12)

(13)

將式(10～12) 代入式(5)，可得中間輔助速度為

(14)

基于式(4)，可得tn + 1時刻的速度為

(15)

顯然，求解式(15)前，需先求解tn + 1時刻的速度。為此，對式(15)兩邊取散度，并考慮tn + 1時刻的不可壓縮條件，可得壓力Poisson方程為

(16)

(17)

2.3 空間離散

(18)

權(quán)函數(shù)采用式(18)的單元插值函數(shù)，乘以式(13，14，17)的兩側(cè)，并在計算域中進行空間積分，可得其弱積分格式為

(19)

(20)

(21)

式中

類似地，對式(15,16)加權(quán)積分后的弱積分形式為

(22)

(23)

至此，獲得了基于 S -A 模型的沿均勻流線的三階 Runge -Kutta 法時間離散的有限元格式，其求解步驟歸納如下。

(3) 通過壓力Poisson方程(23)，求得tn + 1時刻的壓力pn +1。

(7) 返回步驟(1)繼續(xù)下一個時刻的計算。

3 算例應(yīng)用與驗證

3.1 方柱繞流

在流體動力學領(lǐng)域，方柱繞流是經(jīng)典的湍流數(shù)值計算案例。文獻[13-15]進行了實驗研究，并獲得了阻力系數(shù)、升力系數(shù)、Strouhal數(shù)和壓力系數(shù)。Shimada等[16]采用兩方程k-ε模型進行了相關(guān)的數(shù)值仿真模擬，文獻[8，9]在S -A模型的基礎(chǔ)上，提出了一階時間離散有限元格式并進行了數(shù)值驗算。下面將基于S -A模型的沿均勻流線的三階 Runge -Kutta 法時間離散的有限元格式來模擬上述方柱繞流問題，以進一步驗證該方法的有效性和計算精度。

典型的方柱無量綱模型和邊界條件如圖1所示，其中以方柱的邊長取為特征尺寸和以左側(cè)的來流速度為特征速度進行無量綱化。本次模型的計算域取為Ω=[-10,25]×[-10,10]，采用四邊形四節(jié)點線性結(jié)構(gòu)單元劃分網(wǎng)格，網(wǎng)格數(shù)為41461，其中方柱近壁區(qū)域采用細網(wǎng)格進行劃分，方柱表面布置了320個節(jié)點，表層單元厚度為0.002，方柱有限元模型如圖2所示。

圖1 方柱繞流模型和邊界條件

圖2 方柱繞流有限元模型

圖3 平均壓力系數(shù)分布圖比較

圖4給出了本文算法和文獻[8]的阻力Cd和升力Cl時程曲線，可以看出，本文算法在無量綱時間30處已獲得穩(wěn)定的收斂解，而文獻[8]的無量綱收斂時間為130，說明本文算法有助于提高算法的收斂性，具有較低的計算成本。

圖4 阻力Cd和升力Cl時程曲線

3.2 覆冰輸電線繞流

在冬季風、雨和雪冰凍等惡劣氣象聯(lián)合作用下，覆冰輸電線受到風激勵作用后產(chǎn)生流致效應(yīng)，當滿足特定條件后，覆冰輸電線會發(fā)生一種低頻率、大振幅的自激振動，嚴重威脅電力供應(yīng)系統(tǒng)的安全運行。為此，眾多學者致力于輸電線覆冰情況下的氣動系數(shù)及其演化規(guī)律的研究，可為舞動的判定提供基礎(chǔ)性數(shù)據(jù)，具有重要的科學意義。

參照國標(LGJ GB1179-83)的型號規(guī)定，輸電線為LGJ-240(直徑D=26.8 mm)型號，并采用文獻[17,18]以固定尺寸的小圓描述覆冰外形的頂端，小圓和輸電線輪廓的大圓之間用切線連接描述覆冰的外形，通過移動小圓來控制覆冰厚度的新月形模型，本文考察覆冰厚度H=D的情況，如圖5所示。

采用均勻來流風場，空氣密度為1.25 kg/m3，動力粘度系數(shù)為1.72×10-5Pa·s，風速為10 m/s與文獻[17]風洞試驗的工況一致。通過以輸電線圓心旋轉(zhuǎn)覆冰輸電線來表征風攻角，如圖6所示。

本文數(shù)值模擬風攻角α從0°～40°(Δα=10°)的情況，并在氣動力系數(shù)改變劇烈的地方進行加密計算。0°風攻角時覆冰輸電線繞流無量綱化后的邊界條件和網(wǎng)格如圖7和圖8所示。

圖5 新月形覆冰模型

圖6 風攻角

圖7 覆冰輸電線繞流邊界條件

圖8 覆冰輸電線繞流網(wǎng)格劃分

圖9為采用本文算法得到的1.0D新月形覆冰輸電線的平均阻力系數(shù)和平均升力系數(shù)隨風攻角的變化關(guān)系，并將計算結(jié)果與已有的風洞試驗和數(shù)值模擬結(jié)果對比分析。文獻[17，19]通過風洞試驗進行了相同風速下1.0D覆冰厚度的輸電線氣動力系數(shù)的測定，Shinichi等[18]基于LES模型對此進行了數(shù)值模擬。可以看出，本文的計算結(jié)果與試驗和文獻[18]數(shù)值模擬的結(jié)果基本保持一致，并很好地展現(xiàn)了升力系數(shù)尖峰突跳現(xiàn)象，這種尖峰突跳現(xiàn)象將給輸電線路的穩(wěn)定性產(chǎn)生影響，嚴重時將產(chǎn)生舞動現(xiàn)象。

圖9 氣動力系數(shù)隨風攻角的變化曲線比較

4 結(jié) 論

本文基于一方程 S -A 湍流模型封閉RANS方程求解湍流繞流問題。通過坐標變換得到了無對流項RANS方程，以使Galerkin法空間離散時獲得最優(yōu)解，時間離散時采用沿均勻流線的三階 Runge -Kutta 進行離散，并利用Taylor展開將動作標系下的量轉(zhuǎn)換成靜坐標系下的量，建立了求解RANS方程的有限元算法。

基于本文有限元算法，數(shù)值模擬了經(jīng)典的方柱繞流和覆冰輸電線繞流問題。數(shù)值結(jié)果表明，得到的平均阻力系數(shù)和均方根升力系數(shù)與實驗結(jié)果基本一致，且優(yōu)于一階時間離散算法的結(jié)果，表明該算法可以有效地預(yù)測湍流繞流問題，且有利于提高數(shù)值解的計算精度和收斂性。