五種湍流渦粘模型在二維方柱繞流數(shù)值模擬中的對比研究

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

(湖南大學(xué) 工程試驗研究中心, 湖南 長沙 410082)

0 引 言

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

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

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

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

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

(1)

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

質(zhì)量方程：

divV=0

(2)

動量方程(Navier-Stokes方程)：

(4)

(5)

(6)

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

1.1 k-ε湍流渦粘模型

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

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

Gk-ρε

(7)

(8)

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

(9)

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

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

Gk-ρε

(10)

(11)

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

Gk-ρε

(12)

(13)

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

1.2 k-ω湍流渦粘模型

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

(15)

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

(16)

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

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

(17)

(18)

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

(19)

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

2.1 湍流渦粘模型選取

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

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

2.2 計算域布置

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

圖1 計算域布置方式Fig.1 Arrangement of computational domain

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

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

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

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

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

圖3 方柱表面壓力測點布置圖Fig.3 Gaging point of static pressure

2.3 數(shù)值方法控制

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

3 模擬結(jié)果分析

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

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

3.1 整體氣動力特性

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

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

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

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

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

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

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

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

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

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

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

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

圖8 不同工況下的阻力系數(shù)脈動均方根對比圖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 方柱表面壓力分布

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

(20)

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

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

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

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

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

3.3 尾流中心線上的湍動能與速度

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

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

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

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

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

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

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

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

對比圖13與圖14可知：本文計算的橫向速度脈動均方根峰值遠大于流向速度脈動均方根峰值，表明方柱尾流渦的橫向運動比流向運動更劇烈。

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

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

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

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

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

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

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

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

(21)

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

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

4 結(jié) 論

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

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

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

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

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

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

參考文獻：

[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].中國公路學(xué)報, 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.