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

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

1 引 言

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

DNS法能夠提供精確的解決方案，但是所需網格必須小于或等于流場中最小的漩渦結構尺寸，對計算機容量要求極高，計算量巨大，限制了其廣泛應用[2]。LES法的計算效率相比DNS法要提高很多，已成功應用于求解各種形狀的湍流繞流問題[3]，但近壁區(qū)域仍需要非常精細的網格。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法進行空間離散，提出求解湍流模型的有限元算法。通過經典的算例，進一步驗證算法的有效性、精度和收斂性。

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)解。因此，采用經時間間隔Δ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)的相關動坐標系下的量轉化為靜坐標系下的量為

(10)

(11)

(12)

(13)

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

(14)

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

(15)

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

(16)

(17)

2.3 空間離散

(18)

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

(19)

(20)

(21)

式中

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

(22)

(23)

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

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

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

3 算例應用與驗證

3.1 方柱繞流

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

典型的方柱無量綱模型和邊界條件如圖1所示，其中以方柱的邊長取為特征尺寸和以左側的來流速度為特征速度進行無量綱化。本次模型的計算域取為Ω=[-10,25]×[-10,10]，采用四邊形四節(jié)點線性結構單元劃分網格，網格數(shù)為41461，其中方柱近壁區(qū)域采用細網格進行劃分，方柱表面布置了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)合作用下，覆冰輸電線受到風激勵作用后產生流致效應，當滿足特定條件后，覆冰輸電線會發(fā)生一種低頻率、大振幅的自激振動，嚴重威脅電力供應系統(tǒng)的安全運行。為此，眾多學者致力于輸電線覆冰情況下的氣動系數(shù)及其演化規(guī)律的研究，可為舞動的判定提供基礎性數(shù)據，具有重要的科學意義。

參照國標(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]風洞試驗的工況一致。通過以輸電線圓心旋轉覆冰輸電線來表征風攻角，如圖6所示。

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

圖5 新月形覆冰模型

圖6 風攻角

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

圖8 覆冰輸電線繞流網格劃分

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

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

4 結 論

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

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