劉 斐，陳方望

(1.中交第四公路工程局有限公司，北京 100022;2.武漢理工大學(xué)道路橋梁與結(jié)構(gòu)工程湖北省重點(diǎn)實驗室，武漢 430070)

日常生活中，大氣與海洋的流動、高速管流、葉輪機(jī)械等絕大多數(shù)流體運(yùn)動都屬于湍流模型。湍流特性在工程中是常見的，因此也是研究學(xué)者高度重視的。湍流流動的核心特征是尺度無窮性和非線性，理論分析、實驗研究抑或是計算機(jī)模擬來徹底認(rèn)識湍流都是十分困難的。計算機(jī)的發(fā)展使得通過數(shù)值模擬對湍流的研究具有更高的操作性，更高精度的數(shù)值模擬是一直被需要的[1]。

在實際工程中，后臺階流動是一種非常經(jīng)典的分離流動，其結(jié)構(gòu)形式簡單，但是能夠近似模擬大多數(shù)實際工程中流體運(yùn)動情況，如管道或河道突然變寬、風(fēng)吹過障礙物等。國內(nèi)外研究學(xué)者針對后臺階的流動已經(jīng)做了不少理論和實驗探究[2-5]，肖瀟[6]等人分析了后臺階流動的回流區(qū)特性和回流區(qū)長度變化規(guī)律；李通[7]等人分析了三維L型后臺階上半部分長高比值的差異對臺階表面流場特性的影響；王宇航[8]等人研究了三維后臺階模型在SA湍流模型、SST湍流模型和顯示代數(shù)雷諾應(yīng)力模型(EARSM)下的模擬結(jié)果，并對比實驗結(jié)果發(fā)現(xiàn)EARSM的模擬結(jié)果更加接近實際情況；Lee等[9]對后臺階分離和再附流動的主要結(jié)構(gòu)開展了試驗探究，對大尺度渦結(jié)構(gòu)的動力學(xué)特性進(jìn)行了詳細(xì)闡述；Otugen M V等[10]分析了再附位置對擴(kuò)張比的影響，發(fā)現(xiàn)擴(kuò)張比和再附點(diǎn)距臺階分離點(diǎn)的距離成正比關(guān)系。該文主要基于計算風(fēng)工程幾種常見的湍流模型，考慮后臺階中流體運(yùn)動作為分析對象，以FLUENT為計算平臺，對比分析了S-A模型、k-ε模型、RSM模型及LES模型等不同種湍流模型中模擬后臺階中的流體運(yùn)動。

1 湍流模型分類

湍流模型是基于一定的約束條件對湍流流體流動的數(shù)學(xué)描述，其發(fā)展起來的種類很多，常用的湍流模型包括單方程S-A(Spalart-Allmaras)模型[11]、雙方程k-ε模型、五方程雷諾力(RSM，Reynolds Stress Model)模型、大渦模擬(LES，Large Eddy Simulation)模型[12,13]等，以下對各個湍流模型原理做詳細(xì)介紹。

1.1 S-A模型

湍流時均的連續(xù)性方程與雷諾方程如下

(1)

(2)

(3)

式中，ρ為流體密度，kg/m3；u為流速，m/s；μ為流體動力黏度，Pa·s；φT為流體參數(shù)；S為源項。

單方程S-A模型是在湍流時均的連續(xù)性方程與雷諾方程的基礎(chǔ)之上，再建立了一個湍動能k的運(yùn)輸方程，將湍流黏度μt表示成k的函數(shù)，從而使方程組封閉。湍動能k的運(yùn)輸方程為

(4)

式中，σk、CD、Cμ為經(jīng)驗常數(shù)，σk=1,CD=0.09,Cμ=0.08～0.38。

(5)

1.2 k-ε模型

1)標(biāo)準(zhǔn)k-ε模型

在湍動能k方程的基礎(chǔ)上，引入了一個湍動耗散率ε的方程，形成k-ε雙方程模型，稱為標(biāo)準(zhǔn)k-ε模型。該模型由Launder和Spalding于1972年提出。模型中的ε定義為

(6)

于是，湍動黏度μt可表示為k與ε的函數(shù)

(7)

因此，標(biāo)準(zhǔn)k-ε模型的運(yùn)輸方程為

(8)

(9)

其中

(10)

式中，C1ε、C2ε、C3ε為經(jīng)驗常數(shù)，默認(rèn)值為C1ε=1.44、C2ε=1.90、C3ε=0.09；σk、σε分別為湍動能和湍動耗散率對應(yīng)的普朗特數(shù)，默認(rèn)值為σk=1.0、σε=1.3；Prt為湍動普朗特數(shù)，默認(rèn)值取Prt=0.85；gi為重力加速度在i方向上的分量；β為熱膨脹系數(shù)；Mt為湍動馬赫數(shù)；a為聲速。

2)重整化k-ε模型

重整化k-ε模型(RNGk-ε模型)是由Yakhot和Orzag提出的，它是對瞬時N-S方程用重整化的數(shù)學(xué)方法推導(dǎo)出來的模型。模型中的常數(shù)與標(biāo)準(zhǔn)k-ε模型不同，而且方程中也出現(xiàn)了新的函數(shù)或項。所得到的湍動能和耗散率方程與標(biāo)準(zhǔn)k-ε模型相似，為

(11)

(12)

(13)

3)可實現(xiàn)k-ε模型

標(biāo)準(zhǔn)k-ε模型對時均應(yīng)變率有特別大的情性，有可能導(dǎo)致負(fù)的正應(yīng)力。為使流動符合湍流的物理定律，需要對正應(yīng)力進(jìn)行某種數(shù)學(xué)約束?？蓪崿F(xiàn)k-ε模型正是基于此提出，其湍動能及耗散率輸運(yùn)方程為

(14)

(15)

(16)

式中，C1ε、C2ε、C3ε、C2、A0為經(jīng)驗常數(shù)，默認(rèn)值為C1ε=1.44、C2ε=1.92、C3ε=0.09、C2=1.9、A0=4.0；σk,σε分別為湍動能和湍動耗散率對應(yīng)的普朗特數(shù)，默認(rèn)值為σk=1.0σε=1.3。

1.3 RSM模型

上述介紹的各種雙方程模型都采用各向同性的湍流黏度來計算湍流應(yīng)力，這些模型很難模擬出旋轉(zhuǎn)流動及流動方向表面曲率變化的影響。為克服這些缺點(diǎn)，有必要直接對雷諾方程中的湍流脈動應(yīng)力建立微分方程式并進(jìn)行求解。

雷諾應(yīng)力模型就是求解雷諾應(yīng)力張量的各個分量的輸運(yùn)方程，具體形式為

(17)

(18)

1.4 LES模型

大渦模型的控制方程是對N-S方程在波數(shù)空間或物理空間進(jìn)行過濾得到的。過濾的過程是去掉比過濾寬度或給定物理寬度小的渦旋，從而得到大渦旋的控制方程為

(19)

(20)

2 不同湍流模型的對比計算與分析

模擬在后臺階裝置中流體流動的現(xiàn)象。如圖1為幾何模型尺寸圖，其中模型左側(cè)設(shè)置為速度進(jìn)口，尺寸為0.1 m，水流以0.2 m/s流入；右側(cè)為自由出口，尺寸為0.2 m；上側(cè)為對稱邊界，尺寸為2 m；下側(cè)為壁面邊界，為臺階狀，高度為0.1 m，上下臺階長度均為1 m。

依據(jù)所述的后臺階裝置的相關(guān)參數(shù)在ANSYS自帶ICEM CFD建立有限元模型如圖2所示，單元總數(shù)為12 730，節(jié)點(diǎn)個數(shù)為12 251，單元尺寸為0.005 m。對模型采用不同的湍流模型進(jìn)行計算，包括k-epsilon(2eqn)標(biāo)準(zhǔn)模型、k-epsilon(2eqn)重整化模型、k-epsilon(2eqn)可實現(xiàn)模型、k-omega(2eqn)標(biāo)準(zhǔn)模型、k-omega(2eqn)BSL模型、k-omega(2eqn)SST模型等，計算結(jié)果如圖3～8所示，包括壓力云圖和速度云圖。

比較分析圖3～圖8可以發(fā)現(xiàn)：

1)流體在截面小的區(qū)域流動時速度要比在截面大的區(qū)域流動時大得多，且在截面變化處會存在速度趨近于零的流體，在流體出口處，流體速度并不會均勻或接近均勻流出，而是出現(xiàn)了很明顯的分層，上部速度最大，底層速度趨近于零。

2)流體在流動過程中與壁面接觸的流體速度會小于沒有與壁面接觸的流體速度。

3)在截面發(fā)生突變的位置，流體壓力會急劇降低稱為負(fù)值，且入口處的壓力小于出口處的壓力。

4)湍流模型選取的不同，對計算的結(jié)果有一定的影響；其中不同湍流模型得到的壓力云圖在總體分布上有較大差異，且最大值相差也很明顯，k-epsilon(2eqn)的三種模型在壓力云圖最大值的差別上要小于k-omega(2eqn)的三種模型計算得到的壓力云圖最大值之間的差異。

5)不同種模型計算得到的速度云圖在速度分布上差異不是很明顯，且速度最大值也相差不大，但還是有微量的差異。

3 結(jié) 語

針對已有學(xué)者提出的不同種湍流模型，對比分析了不同種湍流模型在同一個后臺階裝置中應(yīng)用時，后臺階內(nèi)部流體流動的情況，包括壓力云圖和速度云圖。結(jié)果表明：湍流模型選取的不同，對計算的結(jié)果有一定的影響；其中不同湍流模型得到的壓力云圖在總體分布上有較大差異，且最大值相差也很明顯，k-epsilon(2eqn)的三種模型在壓力云圖最大值的差別上要小于k-omega(2eqn)的三種模型計算得到的壓力云圖最大值之間的差異；不同種湍流模型計算得到的速度云圖在速度分布上差異不是很明顯，且速度最大值也相差不大，僅存在有微量的差異。