基于S -A湍流模型的Runge -Kutta有限元算法

曹鵬程， 廖紹凱,2， 張 研， 陳 達(dá)

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

1 引 言

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

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

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

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

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

2.1 控制方程

二維非定常不可壓縮流的無(wú)量綱RANS方程和S -A模型的標(biāo)量方程為

(1)

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

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

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

(2)

2.2 時(shí)間離散

(3,4)

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

(5)

(6)

(7)

(8)

(9)

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

(10)

(11)

(12)

(13)

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

(14)

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

(15)

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

(16)

(17)

2.3 空間離散

(18)

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

(19)

(20)

(21)

式中

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

(22)

(23)

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

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

(7) 返回步驟(1)繼續(xù)下一個(gè)時(shí)刻的計(jì)算。

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

3.1 方柱繞流

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

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

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

圖2 方柱繞流有限元模型

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

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

圖4 阻力Cd和升力Cl時(shí)程曲線

3.2 覆冰輸電線繞流

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

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

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

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

圖5 新月形覆冰模型

圖6 風(fēng)攻角

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

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

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

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

4 結(jié) 論

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

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