張顯雄, 張志田, 張偉峰, 陳政清

(湖南大學(xué) 工程試驗(yàn)研究中心, 湖南 長(zhǎng)沙 410082)

0 引 言

隨著橋梁結(jié)構(gòu)趨向于細(xì)、長(zhǎng)、柔以及低阻尼方向發(fā)展，氣動(dòng)彈性問(wèn)題成為了大跨度橋梁的控制性因素，因此對(duì)各種橋梁抗風(fēng)問(wèn)題的研究越來(lái)越精細(xì)化。計(jì)算流體動(dòng)力學(xué)(Computational Fluid Dynamics, CFD)作為一種有可能替代風(fēng)洞試驗(yàn)的跨學(xué)科新型工具，已經(jīng)在橋梁抗風(fēng)領(lǐng)域展開(kāi)了初步應(yīng)用。相對(duì)航天航空、海洋勘探等領(lǐng)域而言，將CFD應(yīng)用于土木工程面臨的一大挑戰(zhàn)是：合理的求解結(jié)果需要對(duì)流體力學(xué)基本理論以及有限體積法這一數(shù)學(xué)工具有較好的掌握，而這正是目前國(guó)內(nèi)的土木工程師所缺乏的。利用CFD對(duì)橋梁抗風(fēng)問(wèn)題展開(kāi)研究時(shí)，橋梁工程師困惑于如何合理離散網(wǎng)格以及采用哪一種湍流模型展開(kāi)計(jì)算。

國(guó)內(nèi)外眾多關(guān)于采用CFD進(jìn)行橋梁抗風(fēng)研究的文獻(xiàn)對(duì)各種模擬結(jié)果的論述參差不一。就湍流渦粘模型為例，文獻(xiàn)[1]采用SSTk-ω湍流模型研究了二維箱梁斷面的靜氣動(dòng)力問(wèn)題，得到的氣動(dòng)力系數(shù)與大渦模擬值及試驗(yàn)值較一致，這一結(jié)果表明SSTk-ω湍流模型可以提供足夠的靜氣動(dòng)力模擬精度。文獻(xiàn)[2]同樣采用 SSTk-ω湍流模型對(duì)箱梁斷面的靜氣動(dòng)力問(wèn)題展開(kāi)了研究，得出SSTk-ω湍流模型用于二維分析時(shí)不能準(zhǔn)確預(yù)測(cè)Strouhal數(shù)的結(jié)論。對(duì)于氣動(dòng)力彈性問(wèn)題，上述狀況同樣存在，如文獻(xiàn)[3]采用RNGk-ε湍流模型對(duì)箱梁斷面計(jì)算的部分顫振導(dǎo)數(shù)與試驗(yàn)數(shù)據(jù)的相對(duì)誤差高達(dá)100%以上；而文獻(xiàn)[4]采用RNGk-ε湍流模型對(duì)箱梁斷面進(jìn)行顫振導(dǎo)數(shù)計(jì)算，卻得到了其計(jì)算結(jié)果能與試驗(yàn)結(jié)果吻合良好的結(jié)論。當(dāng)前，雖然已有學(xué)者提出了相關(guān)建議對(duì)CFD計(jì)算質(zhì)量進(jìn)行評(píng)判[5]，然而形成統(tǒng)一的檢驗(yàn)標(biāo)準(zhǔn)仍然有待發(fā)展。

鈍體橋梁截面繞流包含多種復(fù)雜的流動(dòng)現(xiàn)象，如剪切層的分離，漩渦脫落與再附，尾流區(qū)的發(fā)展等。矩形截面柱作為典型的鈍體構(gòu)件，因?yàn)榫邆鋷缀涡螤詈?jiǎn)單、分離點(diǎn)固定，流動(dòng)形態(tài)多樣等優(yōu)點(diǎn)，往往被用來(lái)驗(yàn)證網(wǎng)格離散方案的適應(yīng)性與湍流數(shù)值模擬的準(zhǔn)確性[6-8]。本文選取在橋梁風(fēng)工程領(lǐng)域應(yīng)用最廣泛的兩方程雷諾平均湍流渦粘模型計(jì)算方柱(長(zhǎng)寬比為1的矩形截面柱)繞流，綜合分析了采用五種兩方程湍流渦粘模型計(jì)算得到的鈍體繞流場(chǎng)特性，對(duì)各個(gè)湍流模型計(jì)算結(jié)果之間的差異展開(kāi)了討論。

1 兩方程湍流渦粘模型的基本理論

在橋梁風(fēng)工程領(lǐng)域，通常認(rèn)為氣體是恒溫不可壓縮的。通過(guò)雷諾平均分解[9](Reynolds Average Navier-Stokes, RANS)，可將廣義物理量的瞬時(shí)值分解為該物理量的時(shí)間平均值與脈動(dòng)值之和：

(1)

則可以建立如下的流體基本方程組。

質(zhì)量方程：

divV=0

(2)

動(dòng)量方程(Navier-Stokes方程)：

(4)

(5)

(6)

引入湍流動(dòng)力粘性系數(shù)μt后，6個(gè)二階雷諾應(yīng)力可以用1個(gè)湍動(dòng)能k表示，于是確定湍流動(dòng)力粘性系數(shù)μt的表達(dá)式就成為計(jì)算湍流流動(dòng)的關(guān)鍵。常用的處理法是將μt表征為兩個(gè)湍流脈動(dòng)特性參數(shù)的函數(shù)。所謂兩方程湍流渦粘模型是指在流體控制方程組(2)～(5)的基礎(chǔ)上引入這兩個(gè)參數(shù)的輸運(yùn)方程，形成封閉的方程組求解流場(chǎng)，其中最典型的是k-ε湍流模型與k-ω湍流模型。

1.1 k-ε湍流渦粘模型

k-ε湍流渦粘模型是一種半經(jīng)驗(yàn)半理論性的高Re數(shù)湍流模型，只適應(yīng)充分發(fā)展的湍流。不同k-ε湍流模型的主要區(qū)別在于：1) 湍流動(dòng)力渦粘系數(shù)μt的計(jì)算方法；2) 控制湍動(dòng)能k與湍動(dòng)能耗散率ε的湍流Prandtl數(shù)σk與σε取值；3) 湍動(dòng)能耗散率ε的生成項(xiàng)與耗散項(xiàng)構(gòu)成形式。

Standardk-ε湍流模型[10]是最早提出的k-ε湍流模型，其輸運(yùn)方程如下：

Gk-ρε

(7)

(8)

式中，Gk表示根據(jù)平均速度梯度計(jì)算的湍動(dòng)能生成項(xiàng)；C1ε與C2ε為常數(shù)。湍流動(dòng)力粘性系數(shù)μt的表達(dá)式為：

(9)

Cμ為常數(shù)，取值為0.09。

由Yakhot與Orzag[11]提出的RNGk-ε湍流模型具有與Standardk-ε湍流模型相似的輸運(yùn)方程：

Gk-ρε

(10)

(11)

Realizablek-ε湍流模型[12]是比RNGk-ε湍流模型應(yīng)用更廣的一種k-ε湍流模型，Shih在文獻(xiàn)[8]中提出湍流動(dòng)力粘性系數(shù)μt的計(jì)算與應(yīng)變率有直接聯(lián)系，其輸運(yùn)方程發(fā)生了明顯變化：

Gk-ρε

(12)

(13)

式中，湍流動(dòng)力粘性系數(shù)μt的表達(dá)式依然與式(9)相同，但此時(shí)Cμ不再為常數(shù)，而是關(guān)于k和ε的函數(shù)。

1.2 k-ω湍流渦粘模型

基于Wilcox[13]方法改進(jìn)提出的Standardk-ω湍流模型也是一種半經(jīng)驗(yàn)湍流模型，ω被定義為特定耗散率。輸運(yùn)方程如下：

(15)

式中，Yk與Yω分別表示因湍流引起的湍動(dòng)能k以及耗散率ω的耗散項(xiàng)；k與ω的有效擴(kuò)散系數(shù)分別用Γk和Γω表示；Gω代表耗散率ω的生成項(xiàng)，湍流動(dòng)力粘性系數(shù)μt的表達(dá)式為：

(16)

α*為降低湍流渦粘性的低雷諾數(shù)修正系數(shù)，是關(guān)于k和ω的函數(shù)。

Menter[14]在Standardk-ω湍流模型的基礎(chǔ)上進(jìn)行修正，提出了的分區(qū)剪應(yīng)力輸運(yùn)k-ω湍流模型(SSTk-ω)，即在近壁面區(qū)域與遠(yuǎn)壁面區(qū)域分別采用Standardk-ε湍流模型與Standardk-ω湍流模型。其輸運(yùn)方程為：

(17)

(18)

Dω為交叉擴(kuò)散項(xiàng)，起協(xié)調(diào)Standardk-ε湍流模型與standardk-ω湍流模型交界區(qū)域的作用。修正后湍流模型的湍流動(dòng)力粘性系數(shù)μt的計(jì)算公式為：

(19)

2 模型控制與數(shù)值計(jì)算

2.1 湍流渦粘模型選取

方柱繞流做為一種典型鈍體構(gòu)件繞流，高Re數(shù)下其流動(dòng)具有固定分離點(diǎn)、尾流區(qū)渦旋運(yùn)動(dòng)不規(guī)則等特性。流體繞經(jīng)方柱時(shí)，流線彎曲，近壁面處流場(chǎng)變化劇烈，方柱表面壓力變化迅速；流動(dòng)在迎風(fēng)面的角點(diǎn)上分離后，產(chǎn)生旋渦脫落，在尾流區(qū)形成自由剪切湍流。方柱繞流因流動(dòng)現(xiàn)象復(fù)雜，物理試驗(yàn)易于實(shí)現(xiàn)而被常用來(lái)檢驗(yàn)湍流模型計(jì)算結(jié)果的準(zhǔn)確程度；另一方面兩方程湍流渦粘模型具有形式簡(jiǎn)單，計(jì)算快捷等優(yōu)點(diǎn)常用來(lái)模擬鈍體繞流。

本文選取了Standardk-ε、RNGk-ε、Realizablek-ε、Standardk-ω以及SSTk-ω這5個(gè)兩方程湍流渦粘模型分別模擬方柱繞流場(chǎng)，并結(jié)合以往數(shù)據(jù)，從三分力系數(shù)、Strouhal數(shù)、回流區(qū)長(zhǎng)度、脈動(dòng)力時(shí)程特性及尾流區(qū)速度場(chǎng)等方面進(jìn)行綜合對(duì)比，分析了各個(gè)湍流渦粘模型的優(yōu)劣與差異。

2.2 計(jì)算域布置

流體為不可壓縮空氣，方柱邊長(zhǎng)D=0.065 m，雷諾數(shù)Re=22 000。因?yàn)闃蛄簽榇罂绺呗柦Y(jié)構(gòu)，展向不同截面處的流動(dòng)具有一定的相似性，而且展向的流動(dòng)尺度遠(yuǎn)小于其余兩個(gè)方向的流動(dòng)尺度，因此本文采用二維分析。參照Tamura等的觀點(diǎn)[15]將計(jì)算域設(shè)置為圓形并在全計(jì)算域采用對(duì)稱離散方案，避免網(wǎng)格不均勻加劇數(shù)值誤差。方柱中心到計(jì)算域周邊的距離為50D，阻塞率為2%，滿足風(fēng)工程研究的要求[16]。其中計(jì)算域入口邊界條件為來(lái)流速度入口，水平來(lái)流速度U=4.95 m/s，來(lái)流湍流度0.5%；出口為壓力出口邊界條件，本文指定出口邊界上的壓力恒定為0；方柱表面為固壁邊界條件。CFD計(jì)算域布置如圖1所示。

圖1 計(jì)算域布置方式Fig.1 Arrangement of computational domain

采用結(jié)構(gòu)化網(wǎng)格離散計(jì)算域。所有湍流模型在近壁面處采用增強(qiáng)壁面處理方式(Enhanced Wall Treatment,EWT)[17]考慮邊界層效應(yīng)，保證不同的湍流模型可以采用相同的離散方案。本文分別采用了5種離散方案，各離散方案的網(wǎng)格數(shù)呈倍數(shù)增加，詳細(xì)的離散方案見(jiàn)表1。

以離散方案三為例，方柱近壁處網(wǎng)格劃分圖2所示。

表1 網(wǎng)格離散參數(shù)Table 1 Parameters of discretization

圖2 方柱近壁處網(wǎng)格Fig.2 Grid near the square cylinder

本文在計(jì)算域內(nèi)布置了226個(gè)測(cè)點(diǎn)用來(lái)監(jiān)測(cè)不同湍流模型下的方柱繞流特性。在尾流區(qū)離方柱0.5D、1D、2D、3D、5D、8D的六個(gè)剖面位置處各自沿Y方向均勻布置了17個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)編號(hào)為1～102；在尾流區(qū)流向中心線上，均勻布置76個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)編號(hào)為103～178，用來(lái)監(jiān)測(cè)尾流區(qū)的速度場(chǎng)。為獲得方柱表面的非定常壓力脈動(dòng)，本文在方柱各邊均勻布置48個(gè)壓力監(jiān)測(cè)點(diǎn)，測(cè)點(diǎn)布置與編號(hào)如圖3所示。根據(jù)有限體積法(Finite Volume Method，F(xiàn)VM)計(jì)算原理，文章在四個(gè)角點(diǎn)位置處(P179，P191，P203，P215)設(shè)置測(cè)點(diǎn)的計(jì)算結(jié)果并不具備數(shù)學(xué)意義，只是用來(lái)標(biāo)明各側(cè)面的長(zhǎng)度界限。

圖3 方柱表面壓力測(cè)點(diǎn)布置圖Fig.3 Gaging point of static pressure

2.3 數(shù)值方法控制

對(duì)于高Re數(shù)下的方柱繞流，本文采用二階隱式格式對(duì)流體控制方程的瞬態(tài)項(xiàng)進(jìn)行離散，空間離散采用二階迎風(fēng)格式；壓力與速度耦合方程采用扭曲修正率為1的SIMPLEC隱式方法進(jìn)行解耦；根據(jù)以往報(bào)道數(shù)據(jù)，雷諾數(shù)Re=22000時(shí)方柱渦脫頻率約為9.5 Hz～10.5 Hz，計(jì)算抽樣頻率取渦脫頻率的200倍，符合采樣定理的要求；無(wú)量綱時(shí)間步長(zhǎng)T=UΔt/D設(shè)置為0.0381；每一時(shí)間步長(zhǎng)內(nèi)迭代20次或者當(dāng)流體控制方程的各項(xiàng)殘差小于10-6時(shí)，認(rèn)為當(dāng)前時(shí)間步內(nèi)的數(shù)值計(jì)算已經(jīng)收斂。

3 模擬結(jié)果分析

在同等計(jì)算資源配置與網(wǎng)格離散方案狀態(tài)下，不同的湍流模型在本文的計(jì)算過(guò)程中體現(xiàn)了不同的耗時(shí)表現(xiàn)。RNGk-ε湍流模型計(jì)算效率最高，其次是Standardk-ω湍流模型，再次是Standardk-ε湍流模型，Realizablek-ε與SSTk-ω湍流模型的計(jì)算效率最低，最高效率與最低效率約相差1.5倍機(jī)時(shí)。

為保證各個(gè)湍流模型的計(jì)算結(jié)果具有可對(duì)比性，開(kāi)始記錄數(shù)據(jù)時(shí)要求各個(gè)湍流模型的計(jì)算結(jié)果均已進(jìn)入統(tǒng)計(jì)定常階段且至少達(dá)到20個(gè)以上渦脫周期。本文記錄了20 000個(gè)時(shí)間步長(zhǎng)的流場(chǎng)參數(shù)，對(duì)應(yīng)的物理時(shí)間為10 s，獲得的方柱旋渦脫落在100個(gè)周期以上。

3.1 整體氣動(dòng)力特性

非定常計(jì)算得到的流場(chǎng)變量都表現(xiàn)出脈動(dòng)特性，氣動(dòng)力參數(shù)需要采用統(tǒng)計(jì)方法獲得。本文采用Lyn的水洞試驗(yàn)結(jié)果(Re=21 400)[18-19]，F(xiàn)ranke的三維RANS模擬結(jié)果(Re=22 000)[20]以及Grigoriadis的三維大渦模擬結(jié)果(Re=22 000)[21]作為整體氣動(dòng)力的參考值。

各湍流模型在不同離散方案下計(jì)算的Strouhal數(shù)如圖4所示。Realizablek-ε湍流模型的計(jì)算結(jié)果偏大，但受網(wǎng)格變化的影響較小，與試驗(yàn)值的相對(duì)誤差最小為4.6%，最大為5.3%；SSTk-ω湍流模型的計(jì)算結(jié)果受網(wǎng)格數(shù)量變化的影響最大，與試驗(yàn)值的相對(duì)誤差最小為5.3%，最大為15.9%；Standardk-ω湍流模型的計(jì)算結(jié)果與Franke的三維模擬結(jié)果吻合較好，與試驗(yàn)值的相對(duì)誤差最小為0.8%，最大為2.3%；RNGk-ε湍流模型的計(jì)算結(jié)果與Standardk-ω湍流模型較接近，與試驗(yàn)值的相對(duì)誤差最小為0.8%，最大為4.5%；Standardk-ε湍流模型的計(jì)算結(jié)果與試驗(yàn)值相差最大，相對(duì)誤差最大為15.9%，最小為13.6%。

Strouhal數(shù)作為表征渦脫頻率的參數(shù)之一，從圖4可以看出Standardk-ε湍流模型模擬的渦脫頻率最小，Realizablek-ε湍流模型計(jì)算的渦脫頻率普遍較大；對(duì)于二維RANS模擬，Standardk-ω湍流模型對(duì)方柱繞流渦脫的捕捉能力明顯優(yōu)于其余湍流模型。

圖4 不同工況下的Strouhal對(duì)比圖Fig.4 Compariosn of Strouhal under various cases

阻力系數(shù)是方柱繞流數(shù)值模擬的難點(diǎn)之一。圖5顯示了各個(gè)湍流模型統(tǒng)計(jì)的阻力系數(shù)平均值隨網(wǎng)格數(shù)增加而變化的趨勢(shì)。SSTk-ω湍流模型的統(tǒng)計(jì)結(jié)果與試驗(yàn)值的相對(duì)誤差最小為0.9%，最大為14.8%；RNGk-ε湍流模型的統(tǒng)計(jì)結(jié)果受網(wǎng)格變化的影響程度與SSTk-ω湍流模型相似，與試驗(yàn)值的相對(duì)誤差最小為2.7%，最大為14.3%；Realizablek-ε湍流模型的統(tǒng)計(jì)結(jié)果分布于試驗(yàn)值與大渦模擬值之間；Standardk-ω湍流模型的統(tǒng)計(jì)結(jié)果與試驗(yàn)值的吻合程度優(yōu)于其余湍流模型，其與試驗(yàn)值的相對(duì)誤差最小為1.8%，最大為3.4%；Standardk-ε湍流模型的計(jì)算結(jié)果最小。

圖5 不同工況下的阻力系數(shù)平均值對(duì)比圖Fig.5 Comparison of under various cases

回流區(qū)長(zhǎng)度Lr是數(shù)值模擬的另一難點(diǎn)。對(duì)比圖6所示各個(gè)湍流模型統(tǒng)計(jì)的回流區(qū)長(zhǎng)度可知，除Standardk-ε湍流模型外，其余湍流模型的計(jì)算結(jié)果普遍低于Lyn的試驗(yàn)值，在RNGk-ε與Realizablek-ε湍流模型算得最小試驗(yàn)相對(duì)誤差為7.9%以及SSTk-ω湍流模型算得最大試驗(yàn)相對(duì)誤差為21%之間隨網(wǎng)格數(shù)變化而上下浮動(dòng)；Standardk-ε與RNGk-ε湍流模型統(tǒng)計(jì)的回流區(qū)長(zhǎng)度受網(wǎng)格變化的影響程度大致相似。

圖6 不同工況下的回流區(qū)長(zhǎng)度對(duì)比圖Fig.6 Comparison of Lr under various cases

升力系數(shù)脈動(dòng)均方根分布如圖7所示。Standardk-ε湍流模型的計(jì)算結(jié)果受網(wǎng)格變化的影響最?。籗tandardk-ω與Realizablek-ε湍流模型湍流模型的統(tǒng)計(jì)結(jié)果隨網(wǎng)格數(shù)量增加的變化趨勢(shì)基本相同；RNGk-ε湍流模型的統(tǒng)計(jì)結(jié)果隨網(wǎng)格數(shù)增加而降低；SSTk-ω湍流模型的計(jì)算結(jié)果整體大于其余湍流模型的計(jì)算結(jié)果。

圖7 不同工況下的升力系數(shù)脈動(dòng)均方根對(duì)比圖Fig.7 Comparison of under various cases

阻力系數(shù)脈動(dòng)均方根分布如圖8所示。SSTk-ω湍流模型的計(jì)算結(jié)果整體大于其余湍流模型的計(jì)算結(jié)果且受網(wǎng)格變化的影響最明顯；其余湍流模型計(jì)算的阻力系數(shù)脈動(dòng)均方根小于Grigoriadis的大渦模擬結(jié)果，Standardk-ε湍流模型的計(jì)算結(jié)果最小。

從上述對(duì)比分析可知：湍流渦粘模型的各項(xiàng)計(jì)算結(jié)果隨網(wǎng)格數(shù)量增加呈現(xiàn)出不同的變化規(guī)律。一種湍流模型采用不同網(wǎng)格離散方案導(dǎo)致不同的計(jì)算結(jié)果；不同的湍流模型采用同一種網(wǎng)格離散方案也導(dǎo)致不同的計(jì)算結(jié)果。對(duì)于Strouhal數(shù)，離散方案三計(jì)算的結(jié)果優(yōu)于其余方案；對(duì)于阻力系數(shù)平均值，離散方案三與方案四的結(jié)果比較好；對(duì)于回流區(qū)長(zhǎng)度，各個(gè)方案的計(jì)算結(jié)果大致相當(dāng)；對(duì)于升力系數(shù)脈動(dòng)均方根，離散方案三與方案四的結(jié)果優(yōu)于其余方案；對(duì)于阻力系數(shù)脈動(dòng)均方根，離散方案三的結(jié)果最優(yōu)。后續(xù)分析采用離散方案三的計(jì)算結(jié)果，表2展示了采用離散方案三時(shí)各個(gè)湍流模型計(jì)算得到的Y+極值。

圖8 不同工況下的阻力系數(shù)脈動(dòng)均方根對(duì)比圖Fig.8 Comparison of under various cases

模型Standardk-εRNGk-εRealizablek-εStandardk-ωSSTk-ω極小值0.0260.0570.0570.0210.040極大值0.6090.6290.6450.6180.689

3.2 方柱表面壓力分布

方柱所受的氣動(dòng)阻力大小與方柱表面的壓力分布有直接聯(lián)系。定義平均壓力系數(shù)的表達(dá)式為：

(20)

其中p為監(jiān)測(cè)點(diǎn)壓力，p0為參考?jí)毫?，本文參考?jí)毫榱?。監(jiān)測(cè)點(diǎn)壓力系數(shù)平均值在方柱表面的分布如圖9所示。Chen的試驗(yàn)研究[22]對(duì)應(yīng)的來(lái)流風(fēng)湍流度小于0.5%，Re=16 000，與本文Re=22 000接近。圖中各條曲線表明：方柱的來(lái)流迎風(fēng)面(179～191號(hào)測(cè)點(diǎn))為正壓區(qū)，各個(gè)湍流模型的計(jì)算得到的壓力系數(shù)平均值與試驗(yàn)值吻合良好；在方柱的背風(fēng)面(203～215號(hào)測(cè)點(diǎn))及上下兩側(cè)面(192～203號(hào)測(cè)點(diǎn)，216～226號(hào)測(cè)點(diǎn))為負(fù)壓區(qū)，因?yàn)樵诜街吔俏恢冒l(fā)生了流動(dòng)分離，導(dǎo)致分離層內(nèi)的壓力顯著降低。除Standardk-ε湍流模型獲得的負(fù)壓力系數(shù)平均值明顯較大以外，RNGk-ε、Realizablek-ε以及Standardk-ω湍流模型三者的計(jì)算結(jié)果較接近且與試驗(yàn)值具有相同的走勢(shì)，SSTk-ω湍流模型計(jì)算得到的統(tǒng)計(jì)值最小。這與圖5描述的平均阻力系數(shù)分布相吻合，表明Standardk-ε湍流模型預(yù)測(cè)流動(dòng)分離存在明顯的劣勢(shì)。方柱繞流阻力主要來(lái)自于迎風(fēng)面與背風(fēng)面的壓力差，Standardk-ε與SSTk-ω湍流模型分別具有最小與最大的壓力差，代表最大的阻力結(jié)果與最小的阻力結(jié)果。

圖9 方柱表面壓力系數(shù)平均值 around the square cylinder

方柱表面的脈動(dòng)壓力主要受渦旋運(yùn)動(dòng)的影響，渦旋運(yùn)動(dòng)越劇烈，脈動(dòng)越明顯。因?yàn)檩斶\(yùn)方程的差異，不同湍流模型計(jì)算得到不同程度的渦旋運(yùn)動(dòng)。圖10描述了方柱表面壓力非定常特性的壓力系數(shù)脈動(dòng)均方根值，反映出流場(chǎng)漩渦做非定常運(yùn)動(dòng)的強(qiáng)弱。在迎風(fēng)面中點(diǎn)位置處(185號(hào)測(cè)點(diǎn))和背風(fēng)面中點(diǎn)處(209號(hào)測(cè)點(diǎn))分別出現(xiàn)壓力系數(shù)脈動(dòng)均方根極小值，表明此處壓強(qiáng)受渦脫的影響最小；角點(diǎn)179與角點(diǎn)191附近的壓力系數(shù)脈動(dòng)均方根出現(xiàn)劇烈變化，表明此兩處壓強(qiáng)受渦脫的影響最大，即渦脫發(fā)生位置。圖中各曲線的排列表明：Standardk-ε湍流模型計(jì)算的渦旋運(yùn)動(dòng)最弱，SSTk-ω湍流模型計(jì)算的渦旋運(yùn)動(dòng)最強(qiáng)。

圖10 方柱表面壓力系數(shù)均方根 around the square cylindert

渦脫是控制流場(chǎng)變化的主要因素之一。通過(guò)對(duì)各個(gè)湍流模型監(jiān)測(cè)的測(cè)點(diǎn)脈動(dòng)壓力時(shí)程作功率譜密度(PSD)分析發(fā)現(xiàn)：除SSTk-ω湍流模型外，方柱各測(cè)點(diǎn)脈動(dòng)壓力的卓越頻率與各自對(duì)應(yīng)的湍流模型計(jì)算得到的渦脫頻率相等，沒(méi)有出現(xiàn)低于主渦脫頻率的運(yùn)動(dòng)。這不符合Norberg[23]關(guān)于渦脫應(yīng)該集中在一個(gè)頻帶上的觀點(diǎn)，說(shuō)明兩方程雷諾平均渦粘模型進(jìn)行二維模擬具有局限性，不能準(zhǔn)確反映湍流的三維基本流動(dòng)特征。

3.3 尾流中心線上的湍動(dòng)能與速度

定義統(tǒng)計(jì)意義下的無(wú)量綱湍動(dòng)能K=k/U2，本文算得的尾流區(qū)中心線上的湍動(dòng)能如圖11所示。從各湍流模型的計(jì)算結(jié)果與Lyn的試驗(yàn)結(jié)果和Bouris的大渦模擬結(jié)果(Re=22 000)[24]對(duì)比圖可以看出：本文計(jì)算結(jié)果的峰值相對(duì)參考值要小。因?yàn)橥牧鞅举|(zhì)上是三維的，二維計(jì)算使得尾流區(qū)湍流模擬不充分是導(dǎo)致本文計(jì)算的湍動(dòng)能小于參考值的原因之一。SSTk-ω湍流模型計(jì)算的湍動(dòng)能值高于其余湍流模型，但其峰值仍比試驗(yàn)結(jié)果要小29%，比LES模擬結(jié)果要小20%，但三者的下降速率大體一致；RNGk-ε與Realizablek-ε、Standardk-ω這三個(gè)湍流模型的峰值基本相同，但下降速度慢于SSTk-ω湍流模型。

圖11 尾流中心線湍動(dòng)能分布Fig.11 Turbulence kinetic energy along thecenter line in the wake

方柱尾流區(qū)中心線上的平均流向速度Uavg分布如圖12所示。SSTk-ω湍流模型與Lyn的試驗(yàn)結(jié)果吻合良好；RNGk-ε湍流模型、Realizablek-ε湍流模型以及Standardk-ω湍流模型的計(jì)算結(jié)果相互接近且與Grigoriadis的大渦模擬結(jié)果吻合較好；Standardk-ε湍流模型的遠(yuǎn)場(chǎng)流域結(jié)果靠近RNGk-ε湍流模型。

除Standardk-ε湍流模型外，在X<2的范圍內(nèi)，所有算例的恢復(fù)速度明顯上升；在X>2后，SSTk-ω湍流模型的計(jì)算結(jié)果增長(zhǎng)速度變慢，增長(zhǎng)至0.58U后基本保持平穩(wěn)；其余湍流模型約在X>6后計(jì)算結(jié)果保持平穩(wěn)，遠(yuǎn)場(chǎng)速度恢復(fù)至0.8U附近。

尾流區(qū)中心線上速度的脈動(dòng)情況是表征湍流特性的重要參數(shù)。圖13與圖14分別表示方柱尾流區(qū)中心線上的流向速度脈動(dòng)均方根urms與橫向速度脈動(dòng)均方根vrms的分布。SSTk-ω湍流模型的計(jì)算結(jié)果雖然與參考值存在一定的偏差，但與試驗(yàn)值的變化趨勢(shì)一致；而對(duì)于橫向速度的脈動(dòng)均方根，除Standardk-ε湍流模型外，其余湍流模型在達(dá)到峰值前均與參考值吻合良好，約在X=2D位置處達(dá)到峰值；之后SSTk-ω湍流模型與參考值的吻合程度明顯優(yōu)于其余湍流模型。

圖12 尾流中心線流向速度平均值Fig.12 Uavg along the center line in the wake

圖13 尾流中心線流向速度脈動(dòng)均方根Fig.13 urms along the center line in the wake

圖14 尾流中心線橫向速度脈動(dòng)均方根Fig.14 vrms along the center line in the wake

對(duì)比圖13與圖14可知：本文計(jì)算的橫向速度脈動(dòng)均方根峰值遠(yuǎn)大于流向速度脈動(dòng)均方根峰值，表明方柱尾流渦的橫向運(yùn)動(dòng)比流向運(yùn)動(dòng)更劇烈。

3.4 尾流區(qū)的平均速度剖面

本文模擬得到的尾流區(qū)流向平均速度剖面如圖15所示。因?yàn)榉街把卮嬖诹鲃?dòng)分離現(xiàn)象，渦旋運(yùn)動(dòng)的動(dòng)量虧損導(dǎo)致尾流區(qū)內(nèi)的速度降低。流體隨著剖面位置向下游移動(dòng)，相同湍流模型在同一剖面上不同位置處的速度差顯著減小，且與參考值的吻合程度明顯降低。不同剖面上的速度分布曲線變化體現(xiàn)了尾流區(qū)主剪切層的鈍化規(guī)律，各個(gè)湍流模型統(tǒng)計(jì)結(jié)果的差異也主要體現(xiàn)在這一區(qū)域內(nèi)。從圖15可以看出流向速度剖面沿中心線呈對(duì)稱分布，本文模擬的主剪切層厚度約在中心線附近-1.5

在X=0.5D剖面，各個(gè)湍流模型的模擬結(jié)果在主剪切層內(nèi)與參考值吻合良好，在剪切層邊緣處的模擬結(jié)果均小于參考值；RNGk-ε與Realizablek-ε湍流模型的流向速度梯度最大，SSTk-ω湍流模型最??；在X=1.5D剖面，Standardk-ε湍流模型在主剪切層內(nèi)的流向速度梯度最大，其在主剪切層外緣的模擬速度與大渦模擬結(jié)果吻合良好，其余湍流模型在主剪切層內(nèi)的速度與大渦模擬結(jié)果吻合較好，RNGk-ε、Realizablek-ε與Standardk-ω湍流模型的流向速度梯度基本一致，SSTk-ω湍流模型的流向速度梯度依然最小；當(dāng)尾流擴(kuò)展至X=3.0剖面位置，SSTk-ω湍流模型的模擬結(jié)果全剖面與試驗(yàn)結(jié)果吻合良好，Standardk-ε湍流模型在主剪切層內(nèi)的流向速度梯度依然最大，其次是SSTk-ω湍流模型；在X=5.0D剖面，主剪切層厚度略微擴(kuò)大，Standardk-ε湍流模型在主剪切層內(nèi)縮小了與其余湍流模型的流向速度梯度差異，SSTk-ω湍流模型擴(kuò)大了與RNGk-ε、Realizablek-ε以及Standardk-ω湍流模型在主剪切層內(nèi)的流向速度梯度的差異；在X=8.0D剖面，主剪切層厚度已經(jīng)不明顯，SSTk-ω湍流模型在全剖面的速度梯度最大，Standardk-ω湍流模型則進(jìn)一步縮小了與其余湍流模型在速度梯度上差異。

各個(gè)湍流模型在尾流區(qū)不同剖面處的橫向平均速度分布如圖15所示。所有湍流模型的計(jì)算結(jié)果均表明：橫向速度在中心線附近-1.0

綜合對(duì)比分析圖16可知：各個(gè)湍流模型計(jì)算結(jié)果之間的差異隨著剖面位置向下游移動(dòng)而明顯減小，而且同一剖面上不同位置處的橫向速度差也隨剖面位置向下游移動(dòng)而明顯降低；Standardk-ε湍流模型的計(jì)算結(jié)果沿流場(chǎng)變化最小，SSTk-ω湍流模型的計(jì)算結(jié)果沿流場(chǎng)變化最大，RNGk-ε、Realizablek-ε以及Standardk-ω湍流模型的計(jì)算結(jié)果始終相互接近。圖中各曲線顯示本文的模擬結(jié)果更接近大渦模擬結(jié)果。

圖15 尾流區(qū)不同截面處流向速度平均值Fig.15 Comparison of Uavg in the wake

圖16 尾流區(qū)不同截面處橫向速度平均值Fig.16 Comparison of Vavg in the wake

方柱尾流區(qū)的流動(dòng)為自由剪切流。尾流特征半厚度作為一種表征尾流區(qū)特征長(zhǎng)度的物理量，可以反映出尾流區(qū)自由剪切流的擴(kuò)散率。定義尾流區(qū)半厚度為Z1/2，在該厚度邊緣處的流向速度為：

(21)

式中，u0(x)為尾流區(qū)中心線上的流向平均速度。本文計(jì)算的各個(gè)湍流模型的尾流特征半厚度沿程分布如圖17所示。

圖17 尾流特征半厚度流程變化Fig.17 Half self-similar profile of the wake

4 結(jié) 論

1) 網(wǎng)格數(shù)量變化對(duì)具體湍流模型各項(xiàng)計(jì)算結(jié)果精度的影響不盡相同。在確定的計(jì)算域內(nèi)，網(wǎng)格劃分?jǐn)?shù)量并非越密集越好，本文的結(jié)果表明：二維模型網(wǎng)格數(shù)量約為Re數(shù)的3倍時(shí)， RANS模擬可取得較優(yōu)的結(jié)果。

2) 在渦脫頻率方面(Strouhal數(shù))，SSTk-ω湍流模型的計(jì)算值受網(wǎng)格數(shù)量的影響十分明顯；其它幾類(lèi)模型均對(duì)網(wǎng)格數(shù)量不敏感。與試驗(yàn)值相比， Standardk-ε模型幾類(lèi)網(wǎng)格方案的計(jì)算值明顯偏小，差值在13%左右。

3) 在平均阻力系數(shù)方面，SSTk-ω湍流模型與RNGk-ε模型的計(jì)算結(jié)果受網(wǎng)格數(shù)量影響顯著，且隨著網(wǎng)格加密，阻力系數(shù)單調(diào)減小，這兩類(lèi)模型明顯存在一個(gè)最優(yōu)網(wǎng)格劃分方案；其它幾類(lèi)模型均對(duì)網(wǎng)格數(shù)量不敏感。但同樣的，與渦脫頻率類(lèi)似，Standardk-ε模型的阻力系數(shù)計(jì)算值明顯低于試驗(yàn)值，相差達(dá)28%左右，劣于其它模型。

4) 在方柱表面壓力分布方面，各類(lèi)模型在迎風(fēng)面上的平均風(fēng)壓基本一致，明顯的差異體現(xiàn)在側(cè)風(fēng)與背風(fēng)面的負(fù)壓區(qū)。Standardk-ε模型負(fù)壓絕對(duì)值十分顯著地低于其它模型，合理解釋了該模型阻力系數(shù)偏低的原因。此外，在脈動(dòng)風(fēng)壓方面， SSTk-ω湍流模型模擬的脈動(dòng)量最高，RNGk-ε、Realizablek-ε的模擬結(jié)果大致相當(dāng)，Standardk-ε湍流模型的壓力脈動(dòng)量最低。

5) 在尾流區(qū)中心線上，SSTk-ω湍流模型模擬的流向速度平均值與橫向速度脈動(dòng)均方根與試驗(yàn)數(shù)據(jù)吻合良好，優(yōu)于其它RANS模型；RNGk-ε、Realizablek-ε以及Standardk-ω湍流模型的計(jì)算結(jié)果大致相當(dāng)；對(duì)于湍動(dòng)能分布，二維RANS湍流模型的計(jì)算結(jié)果明顯低于試驗(yàn)值。

(6) 隨著無(wú)量綱距離X/D的增加，各模型模擬得到的尾流區(qū)剪切層厚度相差明顯。Standardk-ω模型的剪切層厚度最大。Standardk-ε模型的剪切層厚度最小，最小的尾跡寬度再次表明了其明顯偏小的阻力系數(shù)。

參考文獻(xiàn)：

[1]Miranda S D, Patruno L, Ricci M, et al.Numerical study of a twin box bridge deck with increasing gap ration by using RANS and LES approaches[J].Engineering Structrues, 2015, 99: 546-558.

[2]Zhu Zhiwen.Feasibility investigation on prediction of vortex shedding of flat box girders based on 2D RANS models[J].Chinese Journal of Highway and Transport, 2015, 28(6): 24-33.(in Chinese)祝志文.基于二維RANS模型計(jì)算扁平箱梁漩渦脫落的可行性分析[J].中國(guó)公路學(xué)報(bào), 2015, 28(6): 24-33.

[3]Giuseppe V.A numerical model for wind loads simulation on long-span bridges[J].Simulation Modelling Practice and Theory, 2003, 11: 315-351.

[4]Lin Huang, Haili Liao, Bin Wang, et al.Numerical simulation for aerodynamic derivatives of bridge deck[J].Simulation Modelling Practice and Theory, 2009, 17: 719-729.

[5]Tominaga Y, Mochida A, Yoshie R, et al.AIJ guiderlindes for practical applications of CFD to pedestrian wind environment around buildings[J].Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96(10/11): 1749-1761.

[6]Sohankar A.Large eddy simulation of flow past rectangular-section cylinders: Side ratio effects[J].Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96: 640-655.

[7]Yu D H, Butler K, Kareem A, et al.Simulation of the influence of aspect rations on the aerodynamics of rectangular prisms[J].Journal of Engineering Mechanics, 2013, 139(4): 429-438.

[8]Bruno L, Salvetti M V, Ricciardelli F.Benchmark on the aerodynamics of a rectangular 5∶1 cylinder: An overview after the first four years of activiey[J].Journal of Wind Engineering and Industrial Aerodynamic, 2014, 126: 87-106.

[9]Markatos N C.The mathematical modelling of turbulent flows[J].Applied Mathematical Modelling, 1986, 10(3): 190-220.

[10]Launder B E, Spalding D B.Lectures in mathematical models of turbulence[M].London: Academic Press, 1972.

[11]Yakhot V, Orszag S A.Renormalization group analysis of turbulence.I.Basic theory[J].Journal of Scientific Computing, 1986, 1(1): 3-51.

[12]Shih T H, Liou W W, Shabbir A, et al.A newk-εeddy viscosity model for high reynolds number turbulent flows[J].Computers & Fluids, 1995, 24(3): 227-238.

[13]Wilcox D C.Turbulence modeling for CFD[M].California: DCW Industries, Inc., 1998.

[14]Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA, 1994, 32(8): 1598-1605.

[15]Tamura T, Itoh Y, Wada A, et al.Numerical study of pressure fluctuations on a rectangular cylinder in aerodynamic oscillation[J].Journal of Wind Engineering and Industrial Aerodynamics, 1995, 54: 239-250.

[16]Tominaga Y, Mochida A, Yoshie R, et al.AIJ guidelines for

practical applications of CFD to pedestrian wind environment around buildings[J].Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96(10): 1749-1761.

[17]Kader B.Temperature and concentration profiles in fully turbulent boundary layers[J].International Journal of Heat Mass Transfer, 1981, 24(9): 1541-1544.

[18]Lyn D A, Einav S, Rodi W, et al.A laser-doppler velocimetry study of ensemble-averaged characteristics of the turbulent near wake of a square cylinder[J].Journal of Fluid Mechanics, 1995, 304: 285-319.

[19]Lyn D A, Rodi W.The flapping shear layer formed by flow separation from the forward corner of a square cylinder[J].Journal of fluid Mechanics, 1994, 267: 353-376.

[20]Franke R, Rodi W.Calculation of vortex shedding past a square cylinder with various turbulence models[M]//Turbulent Shear Flows 8.Springer Berlin Heidelberg, 1993: 189-204.

[21]Grigoriadis D G E, Bartzis J G, Goulas A.LES of the flow past a rectangular cylinder using the immersed boundary concept[J].International Journal for Numerical Methods in Fluids, 2003, 41(6): 615-632.

[22]Chen J M, Liu C H.Vortex shedding and surface pressures on a square cylinder at incidence to a uniform air stream[J].International Journal of Heat and Fluid Flow, 1999, 20(6): 592-597.

[23]Norberg C.Interaction between freestream turbulence and vortex shedding for a single tube in cross-flow[J].Journal of Wind Engineering and Industrial Aerodynamics, 1986, 23: 501-514.

[24]Bouris D, Bergeles G.2D LES of vortex shedding from a square cylinder[J].Journal of Wind Engineering and Industrial Aerodynamics, 1999, 80(1): 31-46.