劉 斐,陳方望
(1.中交第四公路工程局有限公司,北京 100022;2.武漢理工大學(xué)道路橋梁與結(jié)構(gòu)工程湖北省重點(diǎn)實(shí)驗(yàn)室,武漢 430070)
日常生活中,大氣與海洋的流動(dòng)、高速管流、葉輪機(jī)械等絕大多數(shù)流體運(yùn)動(dòng)都屬于湍流模型。湍流特性在工程中是常見(jiàn)的,因此也是研究學(xué)者高度重視的。湍流流動(dòng)的核心特征是尺度無(wú)窮性和非線性,理論分析、實(shí)驗(yàn)研究抑或是計(jì)算機(jī)模擬來(lái)徹底認(rèn)識(shí)湍流都是十分困難的。計(jì)算機(jī)的發(fā)展使得通過(guò)數(shù)值模擬對(duì)湍流的研究具有更高的操作性,更高精度的數(shù)值模擬是一直被需要的[1]。
在實(shí)際工程中,后臺(tái)階流動(dòng)是一種非常經(jīng)典的分離流動(dòng),其結(jié)構(gòu)形式簡(jiǎn)單,但是能夠近似模擬大多數(shù)實(shí)際工程中流體運(yùn)動(dòng)情況,如管道或河道突然變寬、風(fēng)吹過(guò)障礙物等。國(guó)內(nèi)外研究學(xué)者針對(duì)后臺(tái)階的流動(dòng)已經(jīng)做了不少理論和實(shí)驗(yàn)探究[2-5],肖瀟[6]等人分析了后臺(tái)階流動(dòng)的回流區(qū)特性和回流區(qū)長(zhǎng)度變化規(guī)律;李通[7]等人分析了三維L型后臺(tái)階上半部分長(zhǎng)高比值的差異對(duì)臺(tái)階表面流場(chǎng)特性的影響;王宇航[8]等人研究了三維后臺(tái)階模型在SA湍流模型、SST湍流模型和顯示代數(shù)雷諾應(yīng)力模型(EARSM)下的模擬結(jié)果,并對(duì)比實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)EARSM的模擬結(jié)果更加接近實(shí)際情況;Lee等[9]對(duì)后臺(tái)階分離和再附流動(dòng)的主要結(jié)構(gòu)開(kāi)展了試驗(yàn)探究,對(duì)大尺度渦結(jié)構(gòu)的動(dòng)力學(xué)特性進(jìn)行了詳細(xì)闡述;Otugen M V等[10]分析了再附位置對(duì)擴(kuò)張比的影響,發(fā)現(xiàn)擴(kuò)張比和再附點(diǎn)距臺(tái)階分離點(diǎn)的距離成正比關(guān)系。該文主要基于計(jì)算風(fēng)工程幾種常見(jiàn)的湍流模型,考慮后臺(tái)階中流體運(yùn)動(dòng)作為分析對(duì)象,以FLUENT為計(jì)算平臺(tái),對(duì)比分析了S-A模型、k-ε模型、RSM模型及LES模型等不同種湍流模型中模擬后臺(tái)階中的流體運(yùn)動(dòng)。
湍流模型是基于一定的約束條件對(duì)湍流流體流動(dòng)的數(shù)學(xué)描述,其發(fā)展起來(lái)的種類很多,常用的湍流模型包括單方程S-A(Spalart-Allmaras)模型[11]、雙方程k-ε模型、五方程雷諾力(RSM,Reynolds Stress Model)模型、大渦模擬(LES,Large Eddy Simulation)模型[12,13]等,以下對(duì)各個(gè)湍流模型原理做詳細(xì)介紹。
湍流時(shí)均的連續(xù)性方程與雷諾方程如下
(1)
(2)
(3)
式中,ρ為流體密度,kg/m3;u為流速,m/s;μ為流體動(dòng)力黏度,Pa·s;φT為流體參數(shù);S為源項(xiàng)。
單方程S-A模型是在湍流時(shí)均的連續(xù)性方程與雷諾方程的基礎(chǔ)之上,再建立了一個(gè)湍動(dòng)能k的運(yùn)輸方程,將湍流黏度μt表示成k的函數(shù),從而使方程組封閉。湍動(dòng)能k的運(yùn)輸方程為
(4)
式中,σk、CD、Cμ為經(jīng)驗(yàn)常數(shù),σk=1,CD=0.09,Cμ=0.08~0.38。
(5)
1)標(biāo)準(zhǔn)k-ε模型
在湍動(dòng)能k方程的基礎(chǔ)上,引入了一個(gè)湍動(dòng)耗散率ε的方程,形成k-ε雙方程模型,稱為標(biāo)準(zhǔn)k-ε模型。該模型由Launder和Spalding于1972年提出。模型中的ε定義為
(6)
于是,湍動(dòng)黏度μt可表示為k與ε的函數(shù)
(7)
因此,標(biāo)準(zhǔn)k-ε模型的運(yùn)輸方程為
(8)
(9)
其中
(10)
式中,C1ε、C2ε、C3ε為經(jīng)驗(yàn)常數(shù),默認(rèn)值為C1ε=1.44、C2ε=1.90、C3ε=0.09;σk、σε分別為湍動(dòng)能和湍動(dòng)耗散率對(duì)應(yīng)的普朗特?cái)?shù),默認(rèn)值為σk=1.0、σε=1.3;Prt為湍動(dòng)普朗特?cái)?shù),默認(rèn)值取Prt=0.85;gi為重力加速度在i方向上的分量;β為熱膨脹系數(shù);Mt為湍動(dòng)馬赫數(shù);a為聲速。
2)重整化k-ε模型
重整化k-ε模型(RNGk-ε模型)是由Yakhot和Orzag提出的,它是對(duì)瞬時(shí)N-S方程用重整化的數(shù)學(xué)方法推導(dǎo)出來(lái)的模型。模型中的常數(shù)與標(biāo)準(zhǔn)k-ε模型不同,而且方程中也出現(xiàn)了新的函數(shù)或項(xiàng)。所得到的湍動(dòng)能和耗散率方程與標(biāo)準(zhǔn)k-ε模型相似,為
(11)
(12)
(13)
3)可實(shí)現(xiàn)k-ε模型
標(biāo)準(zhǔn)k-ε模型對(duì)時(shí)均應(yīng)變率有特別大的情性,有可能導(dǎo)致負(fù)的正應(yīng)力。為使流動(dòng)符合湍流的物理定律,需要對(duì)正應(yīng)力進(jìn)行某種數(shù)學(xué)約束??蓪?shí)現(xiàn)k-ε模型正是基于此提出,其湍動(dòng)能及耗散率輸運(yùn)方程為
(14)
(15)
(16)
式中,C1ε、C2ε、C3ε、C2、A0為經(jīng)驗(yàn)常數(shù),默認(rèn)值為C1ε=1.44、C2ε=1.92、C3ε=0.09、C2=1.9、A0=4.0;σk,σε分別為湍動(dòng)能和湍動(dòng)耗散率對(duì)應(yīng)的普朗特?cái)?shù),默認(rèn)值為σk=1.0σε=1.3。
上述介紹的各種雙方程模型都采用各向同性的湍流黏度來(lái)計(jì)算湍流應(yīng)力,這些模型很難模擬出旋轉(zhuǎn)流動(dòng)及流動(dòng)方向表面曲率變化的影響。為克服這些缺點(diǎn),有必要直接對(duì)雷諾方程中的湍流脈動(dòng)應(yīng)力建立微分方程式并進(jìn)行求解。
雷諾應(yīng)力模型就是求解雷諾應(yīng)力張量的各個(gè)分量的輸運(yùn)方程,具體形式為
(17)
(18)
大渦模型的控制方程是對(duì)N-S方程在波數(shù)空間或物理空間進(jìn)行過(guò)濾得到的。過(guò)濾的過(guò)程是去掉比過(guò)濾寬度或給定物理寬度小的渦旋,從而得到大渦旋的控制方程為
(19)
(20)
模擬在后臺(tái)階裝置中流體流動(dòng)的現(xiàn)象。如圖1為幾何模型尺寸圖,其中模型左側(cè)設(shè)置為速度進(jìn)口,尺寸為0.1 m,水流以0.2 m/s流入;右側(cè)為自由出口,尺寸為0.2 m;上側(cè)為對(duì)稱邊界,尺寸為2 m;下側(cè)為壁面邊界,為臺(tái)階狀,高度為0.1 m,上下臺(tái)階長(zhǎng)度均為1 m。
依據(jù)所述的后臺(tái)階裝置的相關(guān)參數(shù)在ANSYS自帶ICEM CFD建立有限元模型如圖2所示,單元總數(shù)為12 730,節(jié)點(diǎn)個(gè)數(shù)為12 251,單元尺寸為0.005 m。對(duì)模型采用不同的湍流模型進(jìn)行計(jì)算,包括k-epsilon(2eqn)標(biāo)準(zhǔn)模型、k-epsilon(2eqn)重整化模型、k-epsilon(2eqn)可實(shí)現(xiàn)模型、k-omega(2eqn)標(biāo)準(zhǔn)模型、k-omega(2eqn)BSL模型、k-omega(2eqn)SST模型等,計(jì)算結(jié)果如圖3~8所示,包括壓力云圖和速度云圖。
比較分析圖3~圖8可以發(fā)現(xiàn):
1)流體在截面小的區(qū)域流動(dòng)時(shí)速度要比在截面大的區(qū)域流動(dòng)時(shí)大得多,且在截面變化處會(huì)存在速度趨近于零的流體,在流體出口處,流體速度并不會(huì)均勻或接近均勻流出,而是出現(xiàn)了很明顯的分層,上部速度最大,底層速度趨近于零。
2)流體在流動(dòng)過(guò)程中與壁面接觸的流體速度會(huì)小于沒(méi)有與壁面接觸的流體速度。
3)在截面發(fā)生突變的位置,流體壓力會(huì)急劇降低稱為負(fù)值,且入口處的壓力小于出口處的壓力。
4)湍流模型選取的不同,對(duì)計(jì)算的結(jié)果有一定的影響;其中不同湍流模型得到的壓力云圖在總體分布上有較大差異,且最大值相差也很明顯,k-epsilon(2eqn)的三種模型在壓力云圖最大值的差別上要小于k-omega(2eqn)的三種模型計(jì)算得到的壓力云圖最大值之間的差異。
5)不同種模型計(jì)算得到的速度云圖在速度分布上差異不是很明顯,且速度最大值也相差不大,但還是有微量的差異。
針對(duì)已有學(xué)者提出的不同種湍流模型,對(duì)比分析了不同種湍流模型在同一個(gè)后臺(tái)階裝置中應(yīng)用時(shí),后臺(tái)階內(nèi)部流體流動(dòng)的情況,包括壓力云圖和速度云圖。結(jié)果表明:湍流模型選取的不同,對(duì)計(jì)算的結(jié)果有一定的影響;其中不同湍流模型得到的壓力云圖在總體分布上有較大差異,且最大值相差也很明顯,k-epsilon(2eqn)的三種模型在壓力云圖最大值的差別上要小于k-omega(2eqn)的三種模型計(jì)算得到的壓力云圖最大值之間的差異;不同種湍流模型計(jì)算得到的速度云圖在速度分布上差異不是很明顯,且速度最大值也相差不大,僅存在有微量的差異。