唐 聃 戴紹仕,2) 谷慧群 Bassam A.Younis

* (哈爾濱工程大學(xué)船舶工程學(xué)院,哈爾濱 150001)

? (加利福尼亞大學(xué)戴維斯分校,美國(guó)加利福尼亞州戴維斯 95616)

引言

邊界層分離所致的非定常、多尺度渦流結(jié)構(gòu)廣泛地存在于各類工程領(lǐng)域的湍流中[1-2],并影響著湍流的生成、耗散等過程,進(jìn)而直接影響著具有復(fù)雜幾何邊界結(jié)構(gòu)遭受的流體動(dòng)力.盡管大渦模擬(large eddy simulation,LES)方法已經(jīng)應(yīng)用于各類工程領(lǐng)域中湍流的預(yù)測(cè)問題上,但隨時(shí)間的推移,其計(jì)算性能的某些缺陷仍然存在[3-6].

目前,廣泛應(yīng)用的經(jīng)典常系數(shù)Smagorinsky 模型(Smagorinsky model,SM)[7]假設(shè)渦黏系數(shù)與湍流中的與局部、過濾的應(yīng)變率和長(zhǎng)度比尺成比例.這種假設(shè)導(dǎo)致了渦黏系數(shù)的過度預(yù)測(cè)[8-9],特別是對(duì)于以多尺度旋渦運(yùn)動(dòng)為主的湍流更加明顯.在實(shí)際的應(yīng)用研究中,眾多學(xué)者也發(fā)現(xiàn)了該模型在近壁區(qū)過度耗散的問題,以致它無(wú)法準(zhǔn)確地捕捉從層流到湍流的轉(zhuǎn)捩過程[3,8,10].為彌補(bǔ)這些缺陷,一些學(xué)者開始研究SM 模型中標(biāo)準(zhǔn)亞格子模型(SGS)的修正方法.Cheng 等[11]引入了基于壁面距離的van Driest阻尼函數(shù)來(lái)修正渦黏系數(shù),在靠近壁面處實(shí)現(xiàn)了亞格子應(yīng)力的強(qiáng)制衰減.但是在實(shí)際應(yīng)用中該模型對(duì)復(fù)雜邊界(尤其是大曲率邊界)流動(dòng)的計(jì)算并不理想,且無(wú)法準(zhǔn)確地預(yù)測(cè)雷諾應(yīng)力[12-13].隨后,Nicoud等[14]基于速度梯度張量的平方修正了渦流黏度的漸近性,進(jìn)而提出了壁面自適應(yīng)渦黏模型(wall adapting local eddy viscosity model,WALE).該模型很好地改進(jìn)了壁面邊界流的預(yù)測(cè)能力,但是對(duì)于高雷諾數(shù)時(shí)的流動(dòng)和強(qiáng)渦流系統(tǒng)的預(yù)報(bào),它過度耗散的行為仍然十分明顯[15-16].近些年,基于Germano恒等式[17]和Lilly 最小二乘法,具有動(dòng)態(tài)模型系數(shù)的Smagorinsky 模型(dynamic Smagorinsky model,DSM)得以發(fā)展,這種方法理論上要求在流場(chǎng)均勻的情況下計(jì)算平均模型系數(shù),但這種平均過程缺少物理上的合理解釋,二次過濾過程又顯著降低了計(jì)算效率,并且在不規(guī)則邊界處的計(jì)算還會(huì)出現(xiàn)錯(cuò)誤[18-19],此外對(duì)渦流場(chǎng)的求解仍呈現(xiàn)出過度耗散的行為[16,20].為了解決過度耗散的問題,Hughes 等[21]又提出了變分多尺度模型(variational multiscale model,VMM)來(lái)區(qū)分旋渦運(yùn)動(dòng)的大尺度和小尺度.隨后,一個(gè)正則化的版本(regularized varitional model,RVM) 由Jeanmart 等[22]初步提出,試圖提高渦流場(chǎng)的預(yù)報(bào)精度,進(jìn)而Bricteux 等[23]探究此方法對(duì)低Re數(shù)槽道流中水平壁面流的預(yù)測(cè)能力,發(fā)現(xiàn)此方法對(duì)復(fù)雜壁面近壁處漸近行為的預(yù)報(bào)并不正確[22-23].近期,為了克服過度耗散的問題,動(dòng)態(tài)模型系數(shù)的方法引起了很多學(xué)者的關(guān)注.如Montecchia 等[24]的顯式代數(shù)亞格子比尺模型,Tellez-Alvarez 等[25]的拉格朗日動(dòng)力學(xué)模型(LDM),Shukla 等[26]與Sircar 等[27]的動(dòng)態(tài)k方程亞格子模型,這些模型都因引入了額外的方程和濾波運(yùn)算[28],不僅對(duì)復(fù)雜壁面形狀黏性繞流的模擬效果產(chǎn)生了不利影響,而且還使計(jì)算效率顯著降低[29].

為了更準(zhǔn)確地捕捉大曲率邊界的邊界層分離和多尺度旋渦發(fā)放等流動(dòng)特征,本文以具有大曲率壁面、復(fù)雜物理性質(zhì)和廣泛工程應(yīng)用的圓柱為例,針對(duì)亞臨界雷諾數(shù)時(shí)圓柱黏性繞流的多尺度旋渦發(fā)放問題,采用正則化變分模型(RVM)與WALE 模型相結(jié)合的方法,推導(dǎo)濾波算子、建立改進(jìn)的大渦模型(MWALE),并采用C++語(yǔ)言開發(fā)計(jì)算代碼,在此基礎(chǔ)上研究改進(jìn)模型中濾波算子 (n) 對(duì)大曲率邊界、多尺度旋渦發(fā)放所致水動(dòng)力問題的計(jì)算能力和預(yù)報(bào)精度.

1 數(shù)值計(jì)算方法

1.1 控制方程

假設(shè)流體流動(dòng)是不可壓縮的,濾波后的連續(xù)性方程和Navier-Stokes (N-S)方程為

基于線性應(yīng)力?應(yīng)變關(guān)系的Smagorinsky 假設(shè)[7],用亞格子(SGS)模型來(lái)近似τij

1.2 WALE 模型的改進(jìn)

式中,Cw是WALE 模型參數(shù).

盡管WALE 模型改善了壁面附近的亞格子應(yīng)力分布,但是對(duì)大規(guī)模分離流的模擬,仍然產(chǎn)生了渦黏系數(shù)過高和能量過度耗散的問題[24].

將正則化變分模型(RVM)與LES 方法相結(jié)合后,明顯地,波數(shù)范圍從 [0,kmax/2] 平滑地過渡到[kmax/2,kmax].通過對(duì)LES 速度場(chǎng)中的進(jìn)行高通濾波(high-pass filtered,HPF)得到HPF 速度場(chǎng)中的(n),并將其應(yīng)用到SM 模型中渦黏系數(shù)和亞格子應(yīng)力的計(jì)算上,發(fā)現(xiàn)該模型對(duì)多尺度渦流運(yùn)動(dòng)的耗散問題有顯著的改善[25].采用LES 速度場(chǎng)的與低通濾波(low-pass filtered,LPF)速度場(chǎng)的之差求得

考慮在三維物理空間中,LPF 濾波場(chǎng)的速度記為(為濾波器,i=1,2,3).本文的濾波器使用的是僅依靠相鄰網(wǎng)格計(jì)算的緊湊型離散濾波器(stencil-3)[30-31].該濾波器是二階濾波器,經(jīng)過調(diào)諧,其傳遞函數(shù)G在LES 速度場(chǎng)中的截?cái)嗖〝?shù)(kh=π)處變?yōu)榱?

以n=1 為例,在一維空間中,將G應(yīng)用于函數(shù)f(x),濾波后的(x) 可表示為

式中,h是網(wǎng)格間距,δ2/4 是簡(jiǎn)化運(yùn)算符.hl和hr分別是非均勻網(wǎng)格左側(cè)和右側(cè)的間距.

擴(kuò)展到三維空間,是通過3 個(gè)方向的顯式離散濾波器(=??)求得.n=1,2 時(shí),LPF 速度場(chǎng)中速度的通式為

式中,I為單位算子;為不同方向的簡(jiǎn)化運(yùn)算符,如=(fi+1,j,k?2fi,j,k+fi?1,j,k)/fi,j,k.n是離散濾波器的濾波階數(shù);n=1 表示二階濾波器算子,n=2 表示4 階濾波器算子.

將式(10) 代入式(6),得到HPF 場(chǎng)中改進(jìn)的WALE 模型中的亞格子應(yīng)力

用式(11)來(lái)封閉控制方程,進(jìn)行數(shù)值求解.

1.3 程序?qū)崿F(xiàn)和數(shù)值求解

基于上述建立的MWALE 模型,采用C++語(yǔ)言在OpenFOAM 開源平臺(tái)上進(jìn)行程序編譯,實(shí)現(xiàn)了此程序的開發(fā)和數(shù)值求解.

程序主體在源文件(*.C)中執(zhí)行,源文件調(diào)用的各種庫(kù)和基本類在頭文件(*.H)中聲明.在make文件夾中,我們指定了編譯源文件的存儲(chǔ)位置和頭文件的聲明位置.然后,通過“wmake”命令生成庫(kù)文件(*.so),用于動(dòng)態(tài)調(diào)用.應(yīng)用MWALE 模型對(duì)多尺度湍流模擬的過程中,需要求解的偏微分方程可以采用OpenFOAM 中的“類”進(jìn)行表示.下面以標(biāo)準(zhǔn)的動(dòng)量方程為例,描述如下

基于有限體積法,該方程的離散化形式可表示為

針對(duì)黏性均勻流繞過孤立圓柱的問題,在數(shù)值模擬中時(shí)間項(xiàng)離散采用隱式的、二階精度向后差分(backward)算法、空間項(xiàng)離散也采用二階精度算法.鑒于在湍流研究中文獻(xiàn)[32-35]已證明: 在網(wǎng)格分辨率相同時(shí),對(duì)流項(xiàng)采用二階精度(或更高)的迎風(fēng)算法、二階中心差分算法在平均和脈動(dòng)的精度計(jì)算上無(wú)太大區(qū)別,都不會(huì)引起數(shù)值耗散,因此,本文對(duì)流項(xiàng)的離散采用線性迎風(fēng)穩(wěn)定輸運(yùn)(linear-upwind stabilized transport,LUST)算法,它是75%線性和25%線性迎風(fēng)的混合算法,以避免可能出現(xiàn)的短波振蕩問題[32-34].拉普拉斯項(xiàng)和梯度項(xiàng)的離散,分別采用高斯線性修正算法和高斯線性算法.代數(shù)方程的求解采用基于算子分裂的壓力隱式算法 (pressureimplicit with splitting of operators,PISO)進(jìn)行迭代求解.數(shù)值計(jì)算流程如圖1 所示.

圖1 數(shù)值計(jì)算流程圖Fig.1 Basic chart of numerical calculation

2 網(wǎng)格設(shè)計(jì)

具有大曲率邊界的圓柱黏性繞流問題,流動(dòng)中既有不固定的分離點(diǎn),又有邊界層分離后的多尺度旋渦發(fā)放,旋渦發(fā)放形態(tài)和尾流運(yùn)動(dòng)隨時(shí)空都有很大的變化,流動(dòng)現(xiàn)象非常復(fù)雜,實(shí)驗(yàn)資料也很豐富.考慮到Re=4×104時(shí)尾流場(chǎng)中有關(guān)脈動(dòng)特性的實(shí)驗(yàn)和數(shù)值資料尤為豐富,所以本文選擇黏性流繞過圓柱的流動(dòng)作為曲面邊界黏性繞流的典型代表,針對(duì)Re=4×104時(shí)圓柱流動(dòng)分離所致的旋渦發(fā)放問題,采用上述建立的MWALE 模型及其數(shù)值求解方法開展數(shù)值模擬和分析.在笛卡爾坐標(biāo)系中,計(jì)算域示意圖見圖2.圓柱底面的圓心位于(0,0,0)處,計(jì)算域的流向(x)、橫向(y)和展向(z)的尺度為29D×16D×πD,圓心距離計(jì)算域進(jìn)口和出口的距離分別為9D和20D,距離側(cè)面的距離為8D.

圖2 計(jì)算域示意圖Fig.2 Diagram of computational domain

為了更好地匹配圓柱體的曲面形狀,采用ANSYSICEM 軟件中的O-Block 技術(shù)進(jìn)行拓?fù)溆成?應(yīng)用六面體網(wǎng)格進(jìn)行計(jì)算域的空間離散(見圖3).為滿足大渦模擬的離散質(zhì)量要求,網(wǎng)格增長(zhǎng)率采用1.05.針對(duì)圓柱體近壁區(qū)的第一層網(wǎng)格(y+),本文設(shè)計(jì)了3 種方案:y+分別為5,2.5 和1.依據(jù)網(wǎng)格近密遠(yuǎn)疏的設(shè)計(jì)原則,沿柱體徑向相應(yīng)地布置95,115 和141 個(gè)種子;沿柱體周向和展向分別布置4×100 和11 個(gè)種子.對(duì)應(yīng)3 種設(shè)計(jì)方案的網(wǎng)格定義為: 粗網(wǎng)格(M1)、中網(wǎng)格(M2)和細(xì)網(wǎng)格(M3),離散單元總數(shù)見表1.

表1 圓柱水動(dòng)力系數(shù)對(duì)比表Table 1 Comparisons of hydrodynamic coefficients of the cylinder

圖3 流場(chǎng)計(jì)算域的網(wǎng)格分布和邊界條件示意圖Fig.3 Diagram of grid distribution and boundary conditions in the computational domain

計(jì)算域的邊界條件和初始條件如圖3 所示.進(jìn)口邊界采用狄利克雷條件,即u=U;出口邊界采用黎曼條件,即 ?u/?n=0.圓柱體壁面采用無(wú)滑移邊界條件.計(jì)算域的兩側(cè)和上下邊界均采用無(wú)窮遠(yuǎn)邊界條件.為保證計(jì)算收斂,時(shí)間步長(zhǎng)(Δt)滿足庫(kù)朗數(shù)C=U·Δt/Δx≤1 (Δx為網(wǎng)格間隔),Δt取5×10?3s.

表2 給出了3 種網(wǎng)格設(shè)計(jì)方案(M1,M2 和M3),n=1,2 時(shí)MWALE 模型預(yù)報(bào)的水動(dòng)力系數(shù)對(duì)比.由表可見,n=1 時(shí)采用M1 和M2 預(yù)報(bào)的結(jié)果差別很小,但與Bouak 等[36]、West 等[37]和Schewe 等[38]實(shí)驗(yàn)測(cè)試的結(jié)果相比偏離很大,尤其是和的預(yù)報(bào)誤差大于10%,然而采用M3 預(yù)報(bào)的,,和S t分別為1.26,0.15,0.42 和0.19,與實(shí)驗(yàn)結(jié)果[36-38]相比誤差都小于10%,明顯更吻合.n=2 時(shí),3 種網(wǎng)格方案預(yù)報(bào)的結(jié)果與實(shí)驗(yàn)結(jié)果更接近.與n=1相比,顯然M3 網(wǎng)格的預(yù)報(bào)精度更高.這些充足的證據(jù)證明了第3 種網(wǎng)格設(shè)計(jì)方案(M3)具有更好數(shù)值計(jì)算精度,因此后續(xù)開展的一系列研究都采用M3方案.

3 結(jié)果討論與分析

為了獲悉MWALE 模型對(duì)曲面結(jié)構(gòu)黏性繞流問題的計(jì)算能力,本文以圓柱為代表,詳盡地研究了MWALE 模型中兩種濾波算子(n=1,2)對(duì)圓柱黏性繞流多尺度旋渦發(fā)放的流動(dòng)特性和力學(xué)特性的預(yù)報(bào)精度的影響.

圓柱體的近壁尾流區(qū)是旋渦生成和發(fā)展的重要區(qū)域.在此區(qū)域中,渦量和速度的分布特性直接影響圓柱體的受力特性.圖4 首先給出了基于兩種濾波算子(n=1,2)的MWALE 模型所預(yù)測(cè)的平均渦量(＜?＞)沿圓柱體表面分布的對(duì)比.由圖可見,在圓柱體迎流區(qū)(0°≤θ ≤90°)內(nèi),兩種算子預(yù)測(cè)的平均渦量都隨 θ 的增加先增大后減小,在 θ=45°時(shí)都捕捉到了渦量的峰值;并且預(yù)測(cè)的平均渦量分布與Son 等[39]的實(shí)驗(yàn)測(cè)試結(jié)果吻合得都很好,但是WALE模型預(yù)測(cè)結(jié)果明顯大于實(shí)驗(yàn)值.在圓柱體的背壓區(qū)(90°＜θ ≤180°),也是多尺度旋渦發(fā)放的區(qū)域,WALE模型和n=1 的MWALE 模型所預(yù)測(cè)的平均渦量都出現(xiàn)顯著的波動(dòng),且預(yù)測(cè)的渦量明顯偏大,而n=2的MWALE 模型預(yù)測(cè)的平均渦量卻非常平穩(wěn),且與實(shí)驗(yàn)結(jié)果[39]吻合得很好.

圖4 圓柱表面上平均渦量分布的對(duì)比Fig.4 Distributions of mean vorticity on the cylinder surface

圖5 進(jìn)一步給出了MWALE 模型中兩種濾波算子所預(yù)測(cè)的一個(gè)旋渦脫落周期內(nèi)、典型時(shí)刻渦量等值面的分布和SGS 渦黏系數(shù)分布云圖.由圖5(a)可以明顯地看出,圓柱壁面發(fā)生流動(dòng)分離后形成了多尺度的旋渦結(jié)構(gòu),在隨主流下泄的過程中發(fā)放的旋渦結(jié)構(gòu)都出現(xiàn)了破碎現(xiàn)象,但與WALE 模型和n=1 的MWLAE 模型相比,n=2 的MWALE 模型預(yù)測(cè)的尾流渦量波動(dòng)和破碎現(xiàn)象明顯地減弱,并且圓柱壁面上渦量的分布也更均勻、更光滑,遠(yuǎn)離圓柱壁面區(qū)的渦量耗散也明顯減小.由圖5(b)給出在不同相位時(shí)SGS 渦黏系數(shù)分布云圖的對(duì)比,顯然可見,3 種計(jì)算模型中預(yù)測(cè)的 νSGS都呈現(xiàn)擴(kuò)散現(xiàn)象,但是WALE 模型預(yù)測(cè)的 νSGS明顯最大,并且在擴(kuò)散中耗散現(xiàn)象最嚴(yán)重(尤其在尾流區(qū)).相比之下,n=1的MWALE 模型預(yù)測(cè)的 νSGS有所降低,耗散現(xiàn)象也有所削弱;而n=2 的MWALE 模型預(yù)測(cè)的 νSGS擴(kuò)散性明顯最小,黏性耗散也最輕.這說(shuō)明n=2 的濾波算子更好地改善了WALE 模型對(duì)SGS 渦黏系數(shù)預(yù)測(cè)過高和耗散過大的問題.這也正好解釋并揭示了下文中旋渦流的脈動(dòng)特性和雷諾應(yīng)力偏大的實(shí)質(zhì)內(nèi)因.

圖5 圓柱尾流渦量等值面和圓柱中截面SGS 渦黏系數(shù)云圖分布Fig.5 Distributions of vorticity isosurfaces in the wake of the cylinder and SGS viscosity contour in the middle section of the cylinder

邊界層發(fā)生分離后,在圓柱背壓區(qū)附近的旋渦運(yùn)動(dòng)得以發(fā)展,會(huì)出現(xiàn)明顯的回流現(xiàn)象,此區(qū)定義為回流區(qū).回流區(qū)長(zhǎng)度與再循環(huán)旋渦長(zhǎng)度(Lr,為沿尾流中心線從圓心到平均流向速度符號(hào)變化處的距離)相等.由于回流區(qū)內(nèi)流體的剪切應(yīng)力較大,并且剪切力與平均流的相互作用非常強(qiáng)烈,從而多尺度旋渦在此生成、發(fā)展,形成湍流.回流區(qū)的長(zhǎng)度被認(rèn)為是評(píng)估圓柱繞流特性的重要指標(biāo)之一.為了評(píng)估MWALE 模型中兩種濾波算子對(duì)回流區(qū)內(nèi)回流的捕捉能力和時(shí)間平均的回流區(qū)長(zhǎng)度的預(yù)測(cè)精度,圖6 給出了沿圓柱尾流中心線(y/D=0)流向無(wú)量綱平均速度分布的對(duì)比.由圖可見: WALE 和MWALE模型預(yù)測(cè)的＜＞/U速度形線在 0.5 ≤x/D≤1.5 都成拋物線形分布,且＜＞/U都出現(xiàn)負(fù)值.MWALE 模型(n=1,2) 預(yù)測(cè)的Lr/D(用D無(wú)量綱) 分別為1.31 和1.40,與Unal 等[40]實(shí)驗(yàn)測(cè)試(Re=4.13×104)的結(jié)果1.45 相比,n=2 預(yù)測(cè)的誤差僅為3.5%,顯然誤差更小;與WALE 模型和DTKE 模型[41]的預(yù)測(cè)值1.33 和1.13 相比,n=2 的MWALE 模型預(yù)測(cè)的精度明顯更高.關(guān)于回流區(qū)內(nèi)的速度形線,相比于n=1 的WALE 模型和DTKE 模型的預(yù)測(cè)結(jié)果[41],n=2 的MWALE 模型與實(shí)驗(yàn)測(cè)試結(jié)果[40]吻合的更好,尤其是預(yù)測(cè)的＜＞/U最小值(?0.28)與實(shí)驗(yàn)值(?0.26)非常接近.流過回流區(qū)后,因旋渦發(fā)放下泄速度的增大,＜＞/U隨x/D＞1.5 增大而增大,并逐漸逼近來(lái)流速度(U),顯然,n=2 的預(yù)測(cè)＜＞/U隨x/D分布與實(shí)驗(yàn)結(jié)果吻合得更好.

圖6 沿尾流中線的平均流向速度Fig.6 Mean stream-wise velocity along the wake centerline

圖7 又進(jìn)一步給出了兩種算子預(yù)測(cè)的在尾流中心線上流向和橫向脈動(dòng)速度(＜＞/U2,＜＞/U2)沿流向的分布.由圖7(a)可見,在近壁尾流區(qū)(0.5 ≤x/D≤2.5),MWALE 模型和WALE 模型所預(yù)測(cè)的＜＞/U2都隨x/D的增加先迅速增加,后迅速降低,并且都在x/D=1.1 處出現(xiàn)峰值,相應(yīng)的峰值分別為0.165,0.132 和0.222.與Szepessy 等[42]、Lim 等[43]和Unal 等[44]的實(shí)驗(yàn)值(0.125～0.142)相比,顯然基于n=1 的MWALE 和WALE模型預(yù)測(cè)的值都明顯偏大,而n=2 的預(yù)測(cè)值在實(shí)驗(yàn)值范圍內(nèi);在尾流遠(yuǎn)場(chǎng)區(qū)(x/D≥2.5),因不同尺度旋渦的充分發(fā)展,隨x/D的增加＜＞/U2的變化逐漸趨于平緩,兩個(gè)算子預(yù)測(cè)的＜＞/U2線形也逐漸趨于一致,并且與Lim 等[43]實(shí)驗(yàn)值十分接近,相比之下,WALE 模型的預(yù)測(cè)值明顯小于實(shí)驗(yàn)值.圖7(b)又給出兩種算子預(yù)測(cè)的橫向脈動(dòng)速度(＜＞/U2)沿流向的分布.明顯地可以看出: 兩種算子預(yù)測(cè)的＜＞/U2峰值都出現(xiàn)在x/D=1.3 附近,但是n=1 的MWALE 模型和WALE 模型所預(yù)測(cè)的＜＞/U2峰值都遠(yuǎn)大于Unal 等[44]的實(shí)驗(yàn)值,并且它們所預(yù)測(cè)的＜＞/U2分布線形明顯地偏離了實(shí)驗(yàn)測(cè)試曲線;而n=2 預(yù)測(cè)的＜＞/U2峰值(0.46)與實(shí)驗(yàn)值[44](0.47)吻合得更好.此外,n=2 預(yù)測(cè)的＜＞/U2分布的線形也與實(shí)驗(yàn)測(cè)試結(jié)果更一致.由此可見,n=2的濾波算子對(duì)WALE 模型的改進(jìn)效果更好,對(duì)渦流場(chǎng)中脈動(dòng)速度的預(yù)測(cè)能力更好,能夠更準(zhǔn)確地捕捉曲面邊界發(fā)放的多尺度旋渦結(jié)構(gòu)引起的流體脈動(dòng)特性.

為進(jìn)一步評(píng)估MWALE 中濾波算子的階數(shù)對(duì)尾流場(chǎng)中速度分布預(yù)測(cè)能力的影響,截取約40 個(gè)旋渦脫落周期(t*=200～400,t*=Ut/D))進(jìn)行數(shù)據(jù)分析.圖8 給出了4 個(gè)典型剖面(x/D=0.6,1.5,2.0,2.5)處無(wú)量綱平均速度分量(＜＞/U,＜＞/U) 的橫向分布.由圖8(a)可見: 在x/D=0.6 時(shí),MWALE 模型(n=1,2)和WALE 模型預(yù)測(cè)的＜＞/U沿尾流中心線都呈現(xiàn)橫向放置的“U”形分布,雖然與Unal 等[40]的實(shí)驗(yàn)形線和DTKE 模型預(yù)測(cè)的行線[45]都很吻合,但是WALE 模型預(yù)測(cè)的結(jié)果在尾流中線附近出現(xiàn)明顯的波動(dòng);隨x/D的增加,3 種模型(MWALE,WALE 和DTKE[45])預(yù)測(cè)的＜＞/U的形線都由倒向的“U”形變?yōu)榈瓜虻摹癡”形分布,并且速度變化幅度都逐漸減小,與實(shí)驗(yàn)測(cè)試[40]的趨勢(shì)符合的也很好(見圖8(b)～圖8(d)).盡管n=1 預(yù)測(cè)的＜＞/U的最小值略大,但是WALE 模型和濾波算子對(duì)＜＞/U預(yù)測(cè)的差異可以忽略.可見,3 種模型對(duì)＜＞/U的預(yù)測(cè)能力相當(dāng).從圖8(a)～圖8(d)中＜＞/U分布可見: 兩種濾波算子和WALE 模型預(yù)測(cè)的尾流中心線上的橫向速度都呈現(xiàn)反對(duì)稱分布,但與WALE 模型、DTKE 模型[45]預(yù)測(cè)結(jié)果相比,基于兩種濾波算子的MWALE 模型預(yù)報(bào)的＜＞/U分布與實(shí)驗(yàn)結(jié)果[40]吻合的更好;隨x/D的增加,因多尺度旋渦的充分發(fā)展,＜＞/U變化梯度明顯增大,與實(shí)驗(yàn)測(cè)試結(jié)果[40]相比,n=2 的MWALE 模型預(yù)測(cè)的誤差明顯最小(見圖8(a)～圖8(d)),這說(shuō)明n=1 的濾波算子對(duì)WALE模型耗散問題的改進(jìn)不佳,導(dǎo)所預(yù)測(cè)的多尺度旋渦發(fā)展中的耗散性偏大.這些證據(jù)再次說(shuō)明n=2 的MWALE 模型預(yù)測(cè)尾流平均速度的能力相對(duì)更好.

圖8 尾流不同位置處平均速度分量的分布Fig.8 Distribution of mean velocity component at different locations in the wake

圓柱黏性繞流因邊界層分離發(fā)生的多尺度渦流的脈動(dòng)特性,至今鮮有關(guān)于大渦模擬模型對(duì)其預(yù)測(cè)性能評(píng)估的詳細(xì)報(bào)道(尤其是Re≥4.0×104).為了進(jìn)一步評(píng)估MWALE 模型對(duì)多尺度旋渦所致雷諾應(yīng)力分布的計(jì)算能力,針對(duì)Re=4.0×104時(shí)圓柱繞流旋渦發(fā)放所致雷諾應(yīng)力分布預(yù)測(cè)能力問題,本文進(jìn)一步討論Re=4.0×104時(shí)MWALE 模型所預(yù)測(cè)的圓柱近壁處雷諾應(yīng)力分布.圖9 給出了MWALE 模型(n=1,2)和WALE 模型預(yù)測(cè)的穩(wěn)態(tài)時(shí)、流向和橫向雷諾應(yīng)力分布云圖的對(duì)比.由圖9(a)可見,流向雷諾應(yīng)力(＜＞/U2)沿尾流中心線都呈現(xiàn)對(duì)稱的漿片狀分布,并都有明顯的鞍點(diǎn)特征,但是n=1 的MWALE模型和WALE 模型預(yù)測(cè)的最大＜＞/U2范圍明顯比實(shí)驗(yàn)結(jié)果[44](Re=4.13×104)大,而n=2 預(yù)測(cè)的與實(shí)驗(yàn)結(jié)果十分接近;實(shí)驗(yàn)測(cè)得的橫向雷諾應(yīng)力(＜＞/U2)沿尾流中心線的分布呈現(xiàn)對(duì)稱的橢圓狀[44](見圖9(b)).顯然,WALE 模型和n=1 的MWALE 模型過度地預(yù)測(cè)了＜＞/U2值,使得圓柱背壓區(qū)處＜＞/U2值明顯偏大,并且尾流區(qū)＜＞/U2的分布形狀也由橢圓變成了近似的矩形.然而,n=2 預(yù)測(cè)的＜＞/U2分布與文獻(xiàn)[44]的測(cè)試結(jié)果十分相近.

圖9 圓柱尾流雷諾應(yīng)力云圖對(duì)比(Re=4.0×104)Fig.9 Comparisons of contour of the fluctuating velocity component in the wake (Re=4.0×104)

為了定量的評(píng)估濾波算子對(duì)雷諾應(yīng)力預(yù)測(cè)精度的影響,圖10 又給出了Re=3.9×103時(shí)在圓柱尾流兩個(gè)典型位置處的雷諾應(yīng)力分量(＜＞/U2,＜＞/U2)分布的對(duì)比.在x/D=1.06 時(shí),與Lourenco等[46]的實(shí)驗(yàn)結(jié)果相比,MWALE 模型、WALE 模型和DTKE 模型[47]預(yù)測(cè)雷諾應(yīng)力分布趨勢(shì)與實(shí)驗(yàn)測(cè)試的趨勢(shì)一致,但是在近尾流中心線上WALE 模型和n=1 的MWALE 模型預(yù)測(cè)的雷諾應(yīng)力橫向(y/D)分布都明顯偏大(見圖10(a)～圖10(b)),WALE 模型尤為明顯,而DTKE 模型[47]預(yù)測(cè)的流向和橫向雷諾應(yīng)力又非常小,幾乎為0.隨著旋渦運(yùn)動(dòng)順流的發(fā)展,沿流向x/D=2.02,1D尾流寬度范圍內(nèi),與實(shí)驗(yàn)值相比,DTKE 模型[47]和n=1 的MWALE 模型預(yù)測(cè)的＜＞/U2峰值都明顯比實(shí)驗(yàn)值大,WALE 模型的預(yù)測(cè)值又略小,而n=2 測(cè)預(yù)測(cè)的＜＞/U2分布和峰值都與實(shí)驗(yàn)測(cè)試[46]結(jié)果非常吻合.對(duì)于＜＞/U2的預(yù)測(cè),WALE 模型和DTKE 模型的預(yù)測(cè)值明顯小于實(shí)驗(yàn)值,而基于n=1 的MWALE 模型預(yù)測(cè)的又明顯偏大,但是n=2 預(yù)測(cè)的＜＞/U2分布與實(shí)驗(yàn)結(jié)果很好地吻合在一起.n=2 的MWALE模型對(duì)雷諾應(yīng)力良好的預(yù)測(cè),進(jìn)一步充分地證明了n=2 的MWALE 模型對(duì)多尺度渦流脈動(dòng)特性的預(yù)測(cè)具有更高精度.

在多尺度渦流的作用下,圓柱體會(huì)受到升、阻力.濾波算子對(duì)圓柱體的無(wú)量綱壓力(中截面處)、升力和阻力的預(yù)報(bào)精度影響見圖11～圖12.由圖11(a)可見,n=2 的MWALE 模型預(yù)測(cè)的穩(wěn)態(tài)壓力系數(shù)(Cp)周向(θ)分布與實(shí)驗(yàn)結(jié)果[48-50]吻合地很好,而n=1 的MAWLE 模型和WALE 模型預(yù)測(cè)的Cp值明顯小于實(shí)驗(yàn)結(jié)果[48-50];對(duì)于脈動(dòng)壓力系數(shù)()周向的預(yù)測(cè)見圖11(b).顯然,當(dāng) θ＜60°時(shí),n=1 或2 預(yù)測(cè)的值都與實(shí)驗(yàn)值[37](Re=4.0×104) 相吻合,但是WALE 模型的預(yù)測(cè)結(jié)果卻略大;當(dāng)θ＞60°時(shí),n=2 預(yù)測(cè)的值也與實(shí)驗(yàn)值吻合的非常好,但是n=1 和WALE 模型預(yù)測(cè)的值波動(dòng)較大,且大于實(shí)驗(yàn)值[37],WALE 模型預(yù)測(cè)值偏大更顯著,這是因耗散偏大所致.

圖12 圓柱的升力和阻力系數(shù)時(shí)歷曲線Fig.12 Time histories of lift and drag coefficients

圖12 給出了兩種濾波算子預(yù)報(bào)的無(wú)量綱升、阻力(Cl,Cd)的時(shí)歷曲線.由圖可見,WALE 模型和n=1 的MWALE 模型預(yù)測(cè)的Cd值波動(dòng)都很大.兩種模型預(yù)測(cè)的穩(wěn)態(tài)和脈動(dòng)阻力系數(shù)(和)分別為(1.35,0.16)和(1.26,0.15),與文獻(xiàn)[36-38]的實(shí)驗(yàn)結(jié)果(1.13～1.22,0.11～0.14)相比都偏大,但是n=2預(yù)報(bào)的和為1.17 和0.12,都在實(shí)驗(yàn)結(jié)果范圍內(nèi);對(duì)于升力,WALE 模型預(yù)測(cè)的值為0.49,而MWALE 模型(n=1,2)預(yù)測(cè)的值分別為0.42,041,都在實(shí)驗(yàn)[36-38]范圍內(nèi)(0.38～0.45).因此,整體來(lái)看,n=2 的MWALE 模型對(duì)圓柱遭受的渦激勵(lì)力的預(yù)報(bào)精度更高.

4 結(jié)論

針對(duì)曲面邊界結(jié)構(gòu)因邊界層分離所致的非定常、多尺度渦流的大渦模擬過度耗散問題,本文采用正則化變分模型(RVM)與WALE 模型相結(jié)合的方法,獲得了改進(jìn)的大渦模型(MWALE),開發(fā)了數(shù)值計(jì)算程序.在此基礎(chǔ)上,以經(jīng)典的圓柱繞流作為曲面邊界黏性繞流的代表,研究了MWALE 模型中濾波算子(n=1,2)對(duì)亞臨界雷諾數(shù)時(shí)圓柱黏性繞流的瞬時(shí)、多尺度渦流的流動(dòng)特性和力學(xué)特性預(yù)測(cè)能力的影響.研究得出以下主要結(jié)論.

(1)改進(jìn)的WALE 模型采用4 階濾波算子(n=2)能更好地反映出圓柱壁面處平均渦量及尾流中瞬時(shí)渦量分布特性,對(duì)圓柱體背壓區(qū)處再循環(huán)旋渦長(zhǎng)度(Lr)的預(yù)測(cè)也更準(zhǔn)確.

(2)與WALE 模型、DTKE 模型和二階濾波算子(n=1)的MWALE 模型相比,4 階濾波算子(n=2)的MWALE 模型能更準(zhǔn)確地預(yù)測(cè)圓柱黏性繞流尾流的流向、橫向脈動(dòng)速度的分布特性和雷諾應(yīng)力的分布特性.

(3) 在圓柱體背壓區(qū),穩(wěn)態(tài)壓力和脈動(dòng)壓力對(duì)MWALE 模型中濾波算子階數(shù)很敏感.因?yàn)? 階濾波算子很好地解決了WALE 模型過渡耗散問題,所以它對(duì)多尺度旋渦發(fā)放所致的脈動(dòng)壓力及升、阻力的預(yù)測(cè)精度更高.