亞臨界區(qū)圓柱繞流相干結(jié)構(gòu)壁面?；旌蟁ANS/LES 模型

季夢(mèng) 尤云祥3)? 韓盼盼 邱小平 馬喬 吳凱健

1) (上海交通大學(xué),海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240)

2) (上海君昱信息科技有限公司,上海 201800)

3) (上海交通大學(xué)三亞崖州灣深?？萍佳芯吭?三亞 572000)

本文發(fā)展了一種具有壁面模化大渦模擬能力的雷諾平均納維-斯托克斯 (RANS)和大渦模擬(LES)方法的混合模型(簡(jiǎn)稱WM-HRL 模型),致力于對(duì)亞臨界區(qū)雷諾數(shù)鈍體繞流相干結(jié)構(gòu)這類復(fù)雜流動(dòng)現(xiàn)象進(jìn)行高置信度的CFD 解析模擬研究.該方法通過(guò)一個(gè)僅與當(dāng)?shù)鼐W(wǎng)格空間分布尺寸有關(guān)的湍動(dòng)能解析度指標(biāo)參數(shù)rk 即可實(shí)現(xiàn)從RANS 到LES 的無(wú)縫快速轉(zhuǎn)換,并且RANS/LES 混合轉(zhuǎn)換區(qū)的邊界位置及其各個(gè)分區(qū)(包括RANS區(qū)、LES 區(qū)及RANS/LES 混合轉(zhuǎn)換區(qū))對(duì)湍動(dòng)能的解析能力均可通過(guò)兩個(gè)指標(biāo)參數(shù) nrk1-Q和nrk2-Q 準(zhǔn)則進(jìn)行預(yù)先設(shè)定.通過(guò)對(duì)雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)的系列數(shù)值模擬研究,獲得了能夠高置信度解析并捕捉其繞流場(chǎng)中三維時(shí)空瞬態(tài)發(fā)展相干結(jié)構(gòu)特性的湍動(dòng)能解析度指標(biāo)參數(shù) nrk1-Q和nrk2-Q 準(zhǔn)則的組合條件.研究表明,該WM-HRL 模型不僅能夠準(zhǔn)確獲取圓柱繞流場(chǎng)中剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的精細(xì)譜結(jié)構(gòu),而且在同一套網(wǎng)格系統(tǒng)下通過(guò)變化湍動(dòng)能解析度指標(biāo)參數(shù) nrk2-Q和nrk1-Q 準(zhǔn)則的組合條件,還可以精細(xì)解析圓柱繞流場(chǎng)中兩類不同回流區(qū)的長(zhǎng)度結(jié)構(gòu)特征,及其對(duì)應(yīng)的圓柱尾部近壁面處V 和U 形兩個(gè)平均流向速度剖面的分支結(jié)構(gòu)特性.

1 引言

鈍體繞流問(wèn)題無(wú)論在理論上還是在工程實(shí)踐中,都有著重要的研究?jī)r(jià)值.所謂亞臨界區(qū)雷諾數(shù)鈍體繞流,即邊界層為層流狀態(tài),而尾流則為湍流狀態(tài)的一類流動(dòng)現(xiàn)象.在理論上亞臨界區(qū)雷諾數(shù)鈍體繞流場(chǎng)研究中,最具挑戰(zhàn)性的難題當(dāng)屬其復(fù)雜三維時(shí)空瞬態(tài)發(fā)展的相干結(jié)構(gòu)特性問(wèn)題,包括上游低強(qiáng)度湍流、自由剪切層中相干結(jié)構(gòu)的產(chǎn)生與發(fā)展、高強(qiáng)度湍流潰滅,以及隨后產(chǎn)生的渦泄現(xiàn)象等[1].

對(duì)圓柱繞流問(wèn)題,其亞臨界區(qū)雷諾數(shù)的范圍一般地定義為 400 ≤Re≤2.0×105.在這個(gè)雷諾數(shù)區(qū)圓柱繞流場(chǎng)中除了會(huì)出現(xiàn)大尺度渦泄這類相干結(jié)構(gòu)外,還會(huì)出現(xiàn)其他一些更為復(fù)雜的三維時(shí)空瞬態(tài)發(fā)展的相干結(jié)構(gòu)現(xiàn)象,包括剪切層小尺度Kelvin-Helmholtz 不穩(wěn)定性(K-H 不穩(wěn)定性)結(jié)構(gòu),在圓柱尾部會(huì)出現(xiàn)兩類不同長(zhǎng)度的回流區(qū)結(jié)構(gòu)現(xiàn)象,以及在圓柱尾部近壁面區(qū)的平均流向速度剖面會(huì)出現(xiàn)V 形和U 形兩個(gè)流動(dòng)分支等[2].

研究表明,在圓柱繞流雷諾數(shù)位于亞臨界區(qū)的情況,當(dāng)雷諾數(shù)Re從400 增大到約1200 時(shí),從圓柱表面分離的剪切層開(kāi)始出現(xiàn)不穩(wěn)定性現(xiàn)象[3],稱為K-H 不穩(wěn)定性.研究進(jìn)一步表明,圓柱繞流場(chǎng)中周期性渦泄(Karman 渦街)現(xiàn)象發(fā)生于雷諾數(shù)Re～190時(shí)[4],但當(dāng)雷諾數(shù)Re達(dá)到5000 附近時(shí),圓柱繞流渦泄出現(xiàn)突然轉(zhuǎn)變現(xiàn)象[3],其主要特征表現(xiàn)為泄渦結(jié)構(gòu)開(kāi)始變得不再具有周期性,這種現(xiàn)象可以維持到雷諾數(shù)Re=2.0×105.

所謂圓柱繞流場(chǎng)中的K-H 不穩(wěn)定性,是指一條速度不連續(xù)的切變線上產(chǎn)生渦度集中而導(dǎo)致的流動(dòng)不穩(wěn)定性現(xiàn)象.Karman 渦街屬于一類頻率相對(duì)較低(頻率記為fvs)的大尺度相干結(jié)構(gòu),而K-H不穩(wěn)定性則屬于一類頻率相對(duì)較高(頻率記為fkh)的小尺度相干結(jié)構(gòu),其主要特征表現(xiàn)為寬頻的信號(hào)特性,且其峰值頻率受雷諾數(shù)Re的影響而顯著變化.在亞臨界區(qū)雷諾數(shù)不大于5000 的情況,這兩種不穩(wěn)定性模式通常可以共存,且兩者的頻率近似滿足[3]fkh/fvs=0.023Re0.67.

在圓柱繞流泄渦為周期性大尺度Karman 渦街的情況,利用RANS 模型通常能夠獲取其較為準(zhǔn)確的渦泄頻率fvs的值,以及Karman 渦街的流態(tài)結(jié)構(gòu).由于RANS 模型只能提供圓柱繞流場(chǎng)的時(shí)均量信息,而不能獲得其三維空間瞬態(tài)發(fā)展的信息,因此對(duì)圓柱繞流大規(guī)模流動(dòng)分離及剪切層小尺度K-H 不穩(wěn)定性等這類復(fù)雜流動(dòng)問(wèn)題,RANS 模型并不適用[5].

直接數(shù)值模擬(DNS)[6-12]和大渦模擬 (LES)[6-8,13-23]以及部分平均N-S(PANS)[2,24-26]等尺度解析模擬(SRS)方法則可以彌補(bǔ)RANS 的這種缺陷,并已成為研究這類復(fù)雜常鈍體繞流問(wèn)題的主要手段.總體上,通過(guò)基于DNS,LES 及PANS 等大量CFD 數(shù)值模擬研究工作,并結(jié)合相關(guān)的模型實(shí)驗(yàn)[18,27]研究工作,對(duì)亞臨界區(qū)雷諾數(shù)下圓柱繞流場(chǎng)涉及的層流分離、層流-湍流轉(zhuǎn)捩、周期性渦脫落及剪切層不穩(wěn)定性等復(fù)雜流動(dòng)現(xiàn)象的形成機(jī)理及其特征等問(wèn)題,有了較為深入的認(rèn)識(shí).

在雷諾數(shù)Re=3900 下的圓柱繞流是典型的亞臨界區(qū)雷諾數(shù)流動(dòng),在其流動(dòng)中除了存在大尺度渦泄及剪切層小尺度K-H 不穩(wěn)定性這兩類相干結(jié)構(gòu)外,無(wú)論在模型實(shí)驗(yàn)中還是在CFD 數(shù)值模擬中,都發(fā)現(xiàn)在其流動(dòng)中還會(huì)出現(xiàn)兩類特殊的流動(dòng)結(jié)構(gòu)現(xiàn)象.其中,第一個(gè)特殊流動(dòng)現(xiàn)象為在其繞流場(chǎng)中會(huì)出現(xiàn)兩種不同長(zhǎng)度的回流區(qū)結(jié)構(gòu),而第二個(gè)特殊流動(dòng)現(xiàn)象是在圓柱尾部近壁面區(qū)的平均流向速度剖面會(huì)出現(xiàn)V和U 形兩個(gè)流動(dòng)分支結(jié)構(gòu).

Lourenco 和Shih[27]的實(shí)驗(yàn)發(fā)現(xiàn),該雷諾數(shù)下圓柱繞流的回流區(qū)長(zhǎng)度為L(zhǎng)r/D=1.18,并且在x1/D=1.06處的平均流向速度剖面呈V 形.然而,在同樣的雷諾下,Parnaudeau等[18]的實(shí)驗(yàn)發(fā)現(xiàn),圓柱繞流的回流區(qū)長(zhǎng)度為L(zhǎng)r/D=1.51,且在x1/D=1.06處的平均流向速度剖面呈U 形.同時(shí),兩個(gè)實(shí)驗(yàn)所測(cè)得的不同站位處平均流向和橫向速度及各向同性和異性雷諾應(yīng)力剖面特性等均不相同.為此,眾多學(xué)者采用CFD 數(shù)值模擬方法對(duì)此問(wèn)題進(jìn)行了長(zhǎng)時(shí)間的研究.

Tremblay[8]采用DNS 和LES,Breuer[15]及Wong 和Png[19]采用LES,Pereira等[2]及Luo等[24]采用PANS,均復(fù)現(xiàn)了Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果.Parnaudeau等[18]及Franke等[17]采用LES,Song等[11]及Ooi等[12]采用DNS,則均復(fù)現(xiàn)了Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果.Tremblay[8]認(rèn)為L(zhǎng)ourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果有誤,Afgan等[20]認(rèn)為L(zhǎng)ourenco 和Shih[27]的統(tǒng)計(jì)周期不夠進(jìn)而導(dǎo)致結(jié)果未收斂,Kravchenko等[16]認(rèn)為L(zhǎng)ourenco 和Shih[27]的實(shí)驗(yàn)由于受到了外部干擾而導(dǎo)致剪切層過(guò)早轉(zhuǎn)捩.

從目前的文獻(xiàn)資料看,許多學(xué)者認(rèn)為Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果更具有可信度,因此大部分學(xué)者都采用Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果作為CFD 數(shù)值模擬的基準(zhǔn)參考數(shù)據(jù),然而這種說(shuō)法很難令人信服.考慮到從實(shí)驗(yàn)、DNS、LES 及PANS等諸多方面均可獲得V 形和U 形結(jié)果,出現(xiàn)這兩種結(jié)果的背后必然有其更深層次的理論和物理方面的機(jī)制,如圓柱體長(zhǎng)徑比、展向網(wǎng)格分辨率及湍流模型對(duì)繞流場(chǎng)湍動(dòng)能的解析能力等.

Ma等[7]采用DNS 及LES,對(duì)圓柱體展向長(zhǎng)度的影響進(jìn)行了研究,發(fā)現(xiàn)當(dāng)展向長(zhǎng)度為2πD時(shí)產(chǎn)生V 形結(jié)果,而當(dāng)展向長(zhǎng)度為πD時(shí)產(chǎn)生U 形結(jié)果,并將產(chǎn)生兩種結(jié)果的原因歸于展向長(zhǎng)度設(shè)置問(wèn)題.Dong等[9]采用DNS 得出了類似結(jié)論.然而,Kravchenko等[16]基于LES 的研究發(fā)現(xiàn),當(dāng)展向網(wǎng)格分辨率為πD/8 時(shí),產(chǎn)生V 形結(jié)果,而當(dāng)展向網(wǎng)格分辨率為πD/48 時(shí),產(chǎn)生U 形結(jié)果,但在保持足夠密的展現(xiàn)網(wǎng)格分辨率時(shí),如果將展向長(zhǎng)度從πD增大至2πD時(shí)并不會(huì)產(chǎn)生明顯區(qū)別.

Xia等[6]采用展向很疏的網(wǎng)格分辨率開(kāi)展DNS 研究(在展向僅布置四層網(wǎng)格),卻驚奇地復(fù)現(xiàn)了V 形結(jié)果.Kim[13]基于LES 研究獲得了與Kravchenko等[16]一致的結(jié)論,即展向網(wǎng)格分辨率是導(dǎo)致U 形和V 形兩個(gè)分支結(jié)構(gòu)的主要原因.然而,Breuer[15]采用LES 的研究發(fā)現(xiàn),在設(shè)置展向πD/32 及πD/64 兩種網(wǎng)格分辨率的情況下均出現(xiàn)V 形結(jié)果.

對(duì)PANS 方法,其對(duì)鈍體繞流場(chǎng)湍動(dòng)能的解析能力通過(guò)一個(gè)濾波參數(shù)fk(=ku/k) 進(jìn)行控制.其中,k為總的湍動(dòng)能,ku為不可解湍動(dòng)能.Pereira等[2]及Ko?ínek等[26]采用PANS 對(duì)V 形和U 形問(wèn)題進(jìn)行了研究.研究發(fā)現(xiàn),當(dāng)fk大于某個(gè)值fk0時(shí),會(huì)產(chǎn)生V 形結(jié)果,而當(dāng)fk小于某個(gè)值fk0時(shí),會(huì)產(chǎn)生U 形結(jié)果,他們將此現(xiàn)象歸因于湍流模型對(duì)湍動(dòng)能的解析能,即具有高湍動(dòng)能解析能力的湍流模型產(chǎn)生U 形結(jié)果,而具有低湍動(dòng)能解析能力的湍流模型產(chǎn)生V 形結(jié)果.

為準(zhǔn)確獲取亞臨界區(qū)雷諾數(shù)圓柱繞流中剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)高頻成分的頻率fkh,PANS 模型中的濾波參數(shù)一般地需要滿足fk≤0.25,即PANS 對(duì)湍動(dòng)能的解析能力至少需要達(dá)到75%[2].這意味著,對(duì)PANS 而言,為準(zhǔn)確獲取亞臨界區(qū)雷諾數(shù)下圓柱繞流場(chǎng)中剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)高頻成分的頻率fkh,其所需網(wǎng)格量應(yīng)當(dāng)與LES 的網(wǎng)格量相當(dāng).

DNS 需要解析邊界層內(nèi)所有尺度的湍流,對(duì)網(wǎng)格分辨率的要求特別高,需要的網(wǎng)格量特別巨大,不適合高雷諾數(shù)鈍體繞流CFD 計(jì)算.LES 直接解析大尺度的湍流,而小尺度湍流則用亞格子應(yīng)力?；?雖然與DNS 相比網(wǎng)格量要少很多,但對(duì)高雷諾數(shù)鈍體繞流的工程化應(yīng)用仍很遙遠(yuǎn)[28].由此可見(jiàn),對(duì)PANS 而言,其工程化應(yīng)用前景亦是如此.

在諸多湍流模型中,兼顧計(jì)算精度與資源消耗的混合RANS/LES(HRL)模型已成為當(dāng)今CFD領(lǐng)域研究與應(yīng)用的熱點(diǎn),包括脫體渦模擬(DES)、延遲脫體渦模擬(DDES)及增強(qiáng)版脫體渦模擬(IDDES)等[29].下文將這三類模型統(tǒng)稱為DES 類模型.

D'Alessandro等[30]基于DES,對(duì)不同網(wǎng)格分辨率及湍流模型的能力進(jìn)行了研究與評(píng)估,認(rèn)為V 形及U 形兩種結(jié)果與網(wǎng)格分辨率及相應(yīng)湍流模型性能關(guān)系密切.研究表明,標(biāo)準(zhǔn)SA-DES 模型在不同加密程度的網(wǎng)格下僅能預(yù)測(cè)V 形結(jié)果,而SA-IDDES 模型在很密網(wǎng)格下可預(yù)測(cè)U 形結(jié)果,在較疏網(wǎng)格下則可預(yù)測(cè)V 形結(jié)果.

綜上所述,目前的研究仍存在諸多未解問(wèn)題,主要為對(duì)V 形及U 形兩個(gè)平均流向速度剖面分支結(jié)構(gòu)產(chǎn)生的機(jī)理尚不能形成統(tǒng)一的認(rèn)識(shí).特別地,為何Ma等[7]的結(jié)果與其他相關(guān)文獻(xiàn)的結(jié)果矛盾? 第二,為何Breuer[15]所用展向網(wǎng)格分辨率與Kravchenko等[16]相同,但獲得的結(jié)果并不一致?第三,SA-DES 與SA-IDDES 模型所采用的RANS和LES 模型一樣,其區(qū)別主要在于RANS 與LES之間的轉(zhuǎn)換方式不同,為何網(wǎng)格分辨率對(duì)其結(jié)果卻有如此大的影響[30]?

另一方面,根據(jù)目前所能查閱到的文獻(xiàn)資料,利用混合RANS/LES 模型來(lái)研究亞臨界區(qū)雷諾數(shù)圓柱繞流的模型主要為DES 及DDES[24,30].雖然這兩類模型能夠較為準(zhǔn)確地獲得與實(shí)驗(yàn)結(jié)果一致的一階統(tǒng)計(jì)量(壓力系數(shù)、流向及橫向速度剖面分布等)的計(jì)算結(jié)果,但對(duì)二階統(tǒng)計(jì)量(各向同性及異性雷諾應(yīng)力剖面分布等)的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果仍有較大差異.同時(shí),由于這類模型對(duì)邊界層中小尺度流動(dòng)結(jié)構(gòu)的解析能力有限,因此難以準(zhǔn)確獲取圓柱繞流剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的精細(xì)信息,特別是其頻譜結(jié)構(gòu)及頻率fkh的準(zhǔn)確值.

有鑒于此,本文發(fā)展了一種壁面?；旌蟁ANS/LES 模型(WM-HRL),該方法也屬于一類HRL 方法.WM-HRL 模型與傳統(tǒng)DES 類模型的主要不同之處在于,可實(shí)現(xiàn)自剪切層小尺度K-H不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域即做具有至少80%湍動(dòng)能解析能力的完全LES 計(jì)算,不僅可有效地減少計(jì)算網(wǎng)格的數(shù)量,而且還可以有效解析剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)特征,并準(zhǔn)確地捕捉其頻譜結(jié)構(gòu)及特征頻率等信息.

在此基礎(chǔ)上,本文以亞臨界區(qū)雷諾數(shù)Re=3900 下的圓柱繞流問(wèn)題為對(duì)象,對(duì)該WM-HRL模型的能力進(jìn)行系列數(shù)值模擬和評(píng)估研究,包括對(duì)圓柱繞流場(chǎng)中剪切層小尺度K-H 不穩(wěn)定性和兩類不同長(zhǎng)度回流區(qū)精細(xì)結(jié)構(gòu)解析與捕捉能力的評(píng)估,以及對(duì)圓柱尾部x1/D=1.06 處平均流向速度剖面出現(xiàn)V 形和U 形兩個(gè)流動(dòng)分支結(jié)構(gòu)形成機(jī)理的進(jìn)一步認(rèn)識(shí)等.

2 理論模型

考慮不可壓縮流體的鈍體繞流問(wèn)題.設(shè)流體密度為ρ0,運(yùn)動(dòng)黏度系數(shù)為ν.建立直角坐標(biāo)系為o-x1x2x3,其中ox3軸垂直 向上為 正,u=(u1,u2,u3) 為流體運(yùn)動(dòng)的速度矢量.對(duì)各類混合RANS/LES 模型,雖然其在鈍體近壁面區(qū)取采用RANS進(jìn)行計(jì)算,而在遠(yuǎn)離鈍體近壁面的區(qū)域采用LES 進(jìn)行計(jì)算,但兩者均采用如下RANS 統(tǒng)一框架下的控制方程對(duì)流場(chǎng)進(jìn)行計(jì)算:

其中,頂標(biāo)“—”表示雷諾時(shí)均.τki為Reynolds 應(yīng)力,采用如下Boussinesq 近似進(jìn)行計(jì)算:

其中,k為湍動(dòng)能,ω為比耗率,ε=β*kω為耗散率,Sˉki為形變率張量,vt為渦黏系數(shù).

為封閉RANS 方程(1),需要引入相應(yīng)的湍流模型,本文采用SSTk-ω 模型.基于該湍流模型的混合RANS/LES 模型,可通過(guò)修改其k方程中的色散項(xiàng)而建立,具體如下[31]:

其中,α1為模型系數(shù),取值為0.31;F2為混合函數(shù)為形變率張量的模.

為保證迭代求解的穩(wěn)定性,Menter等[32]加入了湍流黏度的限制性條件,對(duì)湍流動(dòng)能生成項(xiàng)Pk進(jìn)行了如下的修正:

其中,β*為SST 的模型常數(shù),取值為0.09.

對(duì)DES 模型[29],定義為

其中,lRANS為RANS尺度,lLES為L(zhǎng)ES尺度,可分別表示為

其中,CDES為模型常數(shù),Δ為L(zhǎng)ES 濾波寬度.

在(8)式中,lRANS和lLES與當(dāng)?shù)鼐W(wǎng)格中湍流特征尺度相關(guān).當(dāng)lRANS＜lLES時(shí),?lhyb=lRANS,此時(shí)DES 為RANS 模式;當(dāng)lRANS＞lLES時(shí),?lhyb=lLES,此時(shí)DES 為L(zhǎng)ES 模式.由此可見(jiàn),DES 模型沒(méi)有明確的RANS/LES 混合轉(zhuǎn)換界面.該模型從RANS到LES 的轉(zhuǎn)換完全由RANS 和LES 尺度的相對(duì)大小決定,或者說(shuō)主要由網(wǎng)格尺度的空間分布決定.因此,在使用該模型時(shí),通常要求流向和展向網(wǎng)格尺度不能同時(shí)小于邊界層的厚度.

然而,當(dāng)邊界層內(nèi)流向網(wǎng)格和展向網(wǎng)格加密至某個(gè)中間情況時(shí),即當(dāng)網(wǎng)格尺度小于RANS 尺度而又遠(yuǎn)大于LES 計(jì)算壁湍流所需的網(wǎng)格尺度時(shí),此時(shí)DES 中的LES 尺度會(huì)提前進(jìn)入邊界層,導(dǎo)致邊界層內(nèi)?；睦字Z應(yīng)力不足(MSD),缺失的湍流脈動(dòng)又不足以被解析出來(lái),當(dāng)流向遇到一定的逆壓梯度時(shí),則產(chǎn)生非物理的流動(dòng)分離現(xiàn)象,即網(wǎng)格誘導(dǎo)非物理分離(GIS),甚至發(fā)生對(duì)數(shù)律層不匹配(LLM)現(xiàn)象[33].

有鑒于此,Spalart等[33]提出可在混合長(zhǎng)度尺度中引入邊界層的識(shí)別函數(shù)來(lái)延遲RANS 到LES 的轉(zhuǎn)換,從而防止LES 提前進(jìn)入邊界層,進(jìn)而得到DDES 方法,具體如下:

其中,fd為一個(gè)經(jīng)驗(yàn)性混合函數(shù),具體形式如下:

其中,dw為網(wǎng)格計(jì)算點(diǎn)到壁面的距離,κ為馮卡門常數(shù),取值為0.41.

函數(shù)fd的分布特點(diǎn)為,在近壁層的某個(gè)范圍內(nèi)(與網(wǎng)格及當(dāng)?shù)亓鲌?chǎng)有關(guān))等于0,即在此區(qū)域內(nèi)混合長(zhǎng)度=lRANS,此時(shí)DDES 模型還原為RANS模式.在此區(qū)域外的計(jì)算區(qū)域,混合長(zhǎng)度決定于兩個(gè)尺度lRANS和lLES的相對(duì)大小,此時(shí)DDES 與DES的表現(xiàn)是一樣的.

DDES 模型雖然解決了原始DES 模型存在的MSD 及GIS 等缺陷,但其仍繼承了DES 存在的其他問(wèn)題,如LLM 缺陷.為此,Shur等[34]提出了增強(qiáng)版的DDES 模型(IDDES),但在原始IDDES定義的混合長(zhǎng)度中,含有一個(gè)所謂的提升函數(shù)fe,這是一個(gè)純?nèi)斯?gòu)造的函數(shù).Gritskevich等[31]指出,該函數(shù)的作用僅為增加渦黏系數(shù),但這種人為增加渦黏系數(shù)的作用是否合理有待商榷,因此建議采用如下的簡(jiǎn)化版本:

hmax為計(jì)算單元的最大網(wǎng)格步長(zhǎng).

對(duì)DES,DDES 和IDDES 模型,當(dāng)其進(jìn)入LES后,在局部平衡流條件下,SST 模型k方程中的生成項(xiàng)與色散項(xiàng)相等[35],即

此外,由文獻(xiàn)[35]可得

由(14)式和(15)式可得

由(16)式可得

由(16)式和(17)式,可得:

(19)式正好與經(jīng)典Smagorinsky 亞格子渦黏模型一致.因此,DES,DDES 和IDDES 模型都是利用SST 模型k方程中產(chǎn)生項(xiàng)與色散項(xiàng)平衡這種極限情況下來(lái)間接等效LES 模式而實(shí)現(xiàn)的.

在(8)式中,濾波尺寸Δ控制LES 能否在Kolmogorov 能量譜尺度下解析盡可能多含能尺度的湍流場(chǎng).在經(jīng)典DES,DDES 及IDDES 模型(DES類模型)中,Δ一般取為最大網(wǎng)格尺寸Δmax.根據(jù)這3 類DES 類模型中“局部平衡流”的假定,即在邊界層外的區(qū)域,要實(shí)現(xiàn)DES 類模型從RANS 模式轉(zhuǎn)換為L(zhǎng)ES 模式,其中SST 湍流模型k方程中生成項(xiàng)與色散項(xiàng)需要達(dá)到平衡,此時(shí)DES 類模型相當(dāng)于經(jīng)典的Smagorinsky 型亞格子渦黏模型.

Breuer等[36]認(rèn)為,在DES 類模型中采用最大網(wǎng)格尺寸并不合適,進(jìn)而建議采用如下的體積立方根尺度:

其中,Δ1,Δ2和Δ3分別為3 個(gè)坐標(biāo)方向的網(wǎng)格尺寸.

Reddy等[37]針對(duì)DDES 模型建議了一個(gè)新的網(wǎng)格混合形式如下:

Shur等[34]則建議采用如下的定義:

其中,Cw=0.15,hwn為壁面法向網(wǎng)格步長(zhǎng).

對(duì)兩方程的SST 湍流模型,(8)式中CDES的取值可采用如下加權(quán)形式[31]:

在(23)式中,CDES,in為IDDES內(nèi)層RANS分支的系數(shù),CDES,out為IDDES 類模型外層LES 分支的系數(shù),一般地可按下式取值:

對(duì)IDDES 模型,在邊界層內(nèi)一般為完全RANS模式,而其完全LES 模式一般發(fā)生在邊界層外[38].由此可見(jiàn),IDDES 除了可以避免MSD 及GIS 的問(wèn)題外,也可以避免LLM 的問(wèn)題.本文研究的亞臨界雷諾數(shù)圓柱繞流問(wèn)題,其剪切層小尺度K-H不穩(wěn)定性發(fā)生在對(duì)數(shù)律層區(qū)內(nèi),由于RANS 模式難以準(zhǔn)確地捕捉到這類小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的三維時(shí)空瞬態(tài)發(fā)展流動(dòng)的精細(xì)結(jié)構(gòu),因此IDDES 同樣難以高置信度地解析這類非定常、非平衡流動(dòng)現(xiàn)象的精細(xì)結(jié)構(gòu)特性.

克服IDDES 模型缺陷的一種有效途徑是使其具有壁面?；鬁u模擬的能力.有鑒于此,本文構(gòu)造一種新的混合函數(shù)fnd,使新所構(gòu)造的混合RANS/LES 模型,除了具有延遲脫體渦模擬(DDES)的能力外,同時(shí)具有壁面?；鬁u模擬(WMLES)的能力,將其稱為WM-HRL 模型.

為此,首先引入湍動(dòng)能解析度指標(biāo)概念,即rk=ku/k,其中ku為未解湍動(dòng)能.當(dāng)rk=1 時(shí),可得ku=k,即湍動(dòng)能被完全?；?此時(shí)WM-HRL 為完全RANS 模式.當(dāng)rk=0 時(shí),即ku=0,即湍動(dòng)能被完全解析,此時(shí)WN-HRL 為完全DNS 模式.對(duì)LES 來(lái)說(shuō),為準(zhǔn)確地捕捉剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)這類特殊流動(dòng)現(xiàn)象,其對(duì)湍動(dòng)能的解析能力至少需要達(dá)到75%[2],即rk≤0.25 .

定義kc為截?cái)嗖〝?shù),它表示LES 的濾波寬度,可由當(dāng)?shù)鼐W(wǎng)格尺寸Δ*確定如下:

根據(jù)Kolmogorov的-5/3 譜定律,當(dāng)kc位于慣性亞區(qū)時(shí),ku可由下式進(jìn)行計(jì)算:

其中,Ck(≈1.5) 為Kolmogorov 常數(shù).由此可得:

其中,Lsgs為亞格 子濾波尺度,Ltur為湍流 含能尺度,它們可分別表示如下:

由(27)式和(28)式可得:

對(duì)當(dāng)?shù)鼐W(wǎng)格尺度Δ*,一種合理的選擇[39]為Δ*=hmax.此時(shí),rk可以改寫為

根據(jù)Nyquist-Shannon 樣本定理,為了能夠準(zhǔn)確解析相關(guān)尺度的湍流結(jié)構(gòu),湍流含能尺度Ltur需要滿足如下條件[39]:

其中,Nr為某個(gè)長(zhǎng)度尺度結(jié)構(gòu)能夠被解析的網(wǎng)格單元數(shù)量,它與單位波長(zhǎng)的波能夠被重構(gòu)所需的樣本點(diǎn)數(shù)量一致.

對(duì)RANS 模式,由于其不能對(duì)湍流結(jié)構(gòu)進(jìn)行直接解析,因此Ltur≤hmax.由(30)式可得,此時(shí)rk≥1.對(duì)LES 模式,Larsson等[39]指出,湍流大尺度脈動(dòng)場(chǎng)信息能夠被準(zhǔn)確捕捉的條件是Nr至少為2,即Ltur≥2hmax,將其代入(30)式,可得

有鑒于此,設(shè)rk1≤.本文將當(dāng)?shù)鼐W(wǎng)格滿足條件rk＞rk1的計(jì)算區(qū)域定義為WM-HRL 模型的完全RANS 區(qū)域,即

其中,P為計(jì)算區(qū)域Dcomp中的網(wǎng) 格單元,rk|P為網(wǎng)格單元P上的湍動(dòng)能解析度指標(biāo)值.

對(duì)WM-HRL 模型來(lái)講,其第2 個(gè)關(guān)鍵是如何合理地確定RANS/LES 混合轉(zhuǎn)化區(qū)域Dhyb.為此,設(shè)rk2為WM-HRL 模型進(jìn)入完全LES 模式時(shí)人們期望的大尺度湍流之解析度指標(biāo)值.據(jù)此,可以將WM-HRL 的RANS/LES 混合轉(zhuǎn)換區(qū)域Dhyb定義為

對(duì)PANS 模型,rk也是衡量其對(duì)鈍體繞流場(chǎng)大尺度湍流解析能力的關(guān)鍵性控制參數(shù).Pereira等[2]對(duì)Re=3900 下圓柱繞流問(wèn)題進(jìn)行了系列數(shù)值模擬,研究發(fā)現(xiàn),當(dāng)rk＞0.5 時(shí),PANS 模型尚不足以充分解析繞流場(chǎng)中大尺度渦結(jié)構(gòu)的信息.同時(shí),研究表明:只有當(dāng)rk≤0.5 時(shí),LES 才具有較好地解析大尺度渦結(jié)構(gòu)信息的能力.因此,把rk2取為小于0.5 的某個(gè)值將是一種合理的選擇.

當(dāng)rk2=0.5 時(shí),由(30)式可知,Ltur≥2.83hmax.再結(jié)合(31)式,可取Nr=3 .此時(shí),Ltur≥3hmax,將其代入(30)式可得rk≤0.48 .由此可見(jiàn),可將rˉk2的取值修正為rˉk2=0.48 .這意味著,當(dāng)WM-HRL模型進(jìn)入完全LES 模式后,其在區(qū)域Dhyb中對(duì)大尺度渦結(jié)構(gòu)的解析能力至少可達(dá)52%.

至此,對(duì)本文將構(gòu)造的一種新的WM-HRL 模型中如何合理地確定兩個(gè)關(guān)鍵區(qū)域DRANS和Dhyb進(jìn)行闡述,分別給出確定WM-HRL 模型進(jìn)行完全RANS 模式的區(qū)域DRANS之rk1準(zhǔn)則,以及確定WM-HRL 模型進(jìn)行RANS/LES 混合轉(zhuǎn)換模式的區(qū)域Dhyb之rk2準(zhǔn)則.其中,rk1≤rˉk1且rk2≤rˉk2.后文將其分別稱為rk1-Q準(zhǔn)則和rk2-Q準(zhǔn)則.

對(duì)如上 所述的rk1-Q和rk2-Q準(zhǔn)則,需要利 用(30)式進(jìn)行計(jì)算,其中RANS 含能尺度Ltur在進(jìn)行CFD 計(jì)算之前屬于未知量.因此,這兩個(gè)準(zhǔn)則無(wú)法用于在進(jìn)行CFD 計(jì)算網(wǎng)格設(shè)置時(shí)來(lái)具體確定兩個(gè)關(guān)鍵區(qū)域DRANS和Dhyb的邊界位置.

為此,下面繼續(xù)構(gòu)造rk的一種僅依賴于當(dāng)?shù)鼐W(wǎng)格尺寸的計(jì)算公式.Pope[40]指出,在對(duì)數(shù)律層區(qū),Ltur與dw成正比關(guān)系,即Ltur=Cwdw.在高雷諾數(shù)的情況,Han等[41]指出,Cw可近似取為Cw≈2.5 .由此可見(jiàn),(30)式可以改寫為:

對(duì)(35)式定義的rk,其值有可能會(huì)出現(xiàn)大于1 的情況.為避免這種情況,將 (35) 式修改為如下形式:

在(36)式的形式下,所定義的湍動(dòng)能解析度指標(biāo)參數(shù)rk將始終不會(huì)超過(guò)1.0,即rk≤1.0 .根據(jù)(36)式可得如下結(jié)論.

首先,當(dāng)rk≥rk1時(shí),dw/hmax≤0.4(rk1)-3/2.在rk≥時(shí),可得dw/hmax≤0.8,此時(shí)LES 將不能被真正激活.一般地,可將WM-HRL 模型為完全RANS 的區(qū)域定義為

其次,當(dāng)rk2＜rk＜rk1時(shí),0.4(rk1)-3/2＜dw/hmax＜0.4(rk2)-3/2.在rk2≤時(shí),可得dw/hmax≥1.2,此時(shí)LES 正好被激活.一般地,可將WM-HRL 模型為RANS/LES 混合轉(zhuǎn)換的區(qū)域定義為

將由(37)式確定的WM-HRL 模型為完全RANS 的區(qū)域稱為nrk1-Q準(zhǔn)則,而將由(38)式確定之WM-HRL 模型為RANS/LES 混合轉(zhuǎn)換的區(qū)域稱為nrk2-Q準(zhǔn)則.由此可知,與前面所述之原rk1-Q和rk2-Q準(zhǔn)則相比,新的nrk1-Q和nrk2-Q準(zhǔn)則僅與當(dāng)?shù)鼐W(wǎng)格尺寸相關(guān),因此可以很容易地根據(jù)這兩個(gè)準(zhǔn)則來(lái)進(jìn)行WM-HRL 模型中兩個(gè)關(guān)鍵區(qū)域DRANS和Dhyb的確定及其相應(yīng)的網(wǎng)格設(shè)置.

根據(jù)如上所述兩個(gè)nrk1-Q和nrk2-Q準(zhǔn)則,可構(gòu)造一種新的混合函數(shù)fnd如下:

其中,fs為待定函數(shù).同時(shí),將混合長(zhǎng)度修改為

由(39)式可知,在當(dāng)?shù)鼐W(wǎng)格滿足nrk1-Q準(zhǔn)則的情況,即rk≥rk1此時(shí)mk=rk1,可得fnd=1 .再結(jié)合(40)式可知,此時(shí)WM-HRL 模型為完全RANS模式.另一方面,當(dāng)rk＜rk1時(shí),由(39)式可知,mk=rk＜rrk1,此時(shí)fns=1-fs.

根據(jù)nrk2-Q準(zhǔn)則,函數(shù)fs需要滿 足如下條 件:第一,當(dāng)rk2＜mk＜rk1時(shí),fs需要滿足 0＜fs＜1 ;第二,當(dāng)mk≤rk2時(shí),fs需要滿足fs=1 ;第三,當(dāng)mk=rk1時(shí),fs需要滿足fs=0 .

能夠同時(shí)滿足上述3 個(gè)條件的函數(shù)可用一個(gè)關(guān)于mk的指數(shù)函數(shù)表示,具體如下:

其中,rk1＞rk2為自定義參數(shù),且取值在0 和1之間.

下面證明,在(42)式的條件下,由(41)式構(gòu)造的函數(shù)fs,滿足其所需的3 條要求:首先,由(42)式可知,當(dāng)mk=rk1時(shí),α=-0.75,由(41)式可知,此時(shí)fs=0 ;其次,當(dāng)rk2＜mk＜rk1時(shí),-0.75＜α＜-0.25,由(41)式可知,此時(shí) 0＜fs＜1 ;第三,當(dāng)mk≤rk2時(shí),α=-0.25,由(41)式可知,此時(shí)fs=1.

目前,常用的一類WMLES 模型為所謂的代數(shù)壁面?；鬁u模擬(Alg WMLES)[34],其基本思想是對(duì)Smagorinsky 型LES 的亞格子渦黏系數(shù)乘以一個(gè)阻尼系數(shù)Fr,即:

其中,vsgs為亞格子渦黏系數(shù);CSMAG為模型常數(shù),取值在0.1—0.18 之間;A為常數(shù),一般取為25;y+為無(wú)量綱壁面距離.

由(43)式可知,當(dāng)y+約大于60 時(shí),阻尼函數(shù)Fr趨于1,此時(shí)Alg WMLES 即為完全Smagorinsky 型亞格子渦黏模型.對(duì)y+≤60 的這個(gè)近壁區(qū)域,正好為黏性底層和過(guò)渡層,由于黏性底層很薄,其范圍約為y+≤10.0,此時(shí)Alg WLMES 的亞格子渦黏系數(shù)vsgs≈0.0,這相當(dāng)于DNS.在過(guò)渡層內(nèi),Fr從0 開(kāi)始增大到1,此時(shí)Alg WLMES 相當(dāng)于準(zhǔn)DNS.由此可見(jiàn),對(duì)高雷諾數(shù)鈍體繞流問(wèn)題,為了捕捉黏性底層及過(guò)渡層這兩個(gè)區(qū)域內(nèi)足夠精細(xì)的湍流信息,此時(shí)Alg WMLES 所需的網(wǎng)格量幾乎與DNS 的網(wǎng)格量相當(dāng).

另一類常用的WMLES 模型為所謂的壁面應(yīng)力?；鬁u模擬(WRMLES)[39],該模型根據(jù)湍流邊界層速度剖面的對(duì)數(shù)律來(lái)計(jì)算壁面剪切力,并輸入到LES 邊界網(wǎng)格作為邊界條件.在WRMLES中,網(wǎng)格間距(Δ)與局部邊界層厚度(δ)通常需要滿足條件δ/Δ≈20—30 .這種較粗的近壁網(wǎng)格有可能會(huì)導(dǎo)致缺乏湍流應(yīng)力的現(xiàn)象,同時(shí)基于最近鄰LES 速度對(duì)壁面應(yīng)力進(jìn)行建模的壁面應(yīng)力邊界條件會(huì)增大邊界層內(nèi)部區(qū)域的總應(yīng)力,因此可能會(huì)致壁面應(yīng)力被低估或高估,進(jìn)而發(fā)生LLM 問(wèn)題[40].

WM-HRL 模型與Alg WMLES 和WRMLES模型均不相同,該模型在黏性底層內(nèi)一般為完全RANS 解析模式,在過(guò)渡層內(nèi)的某個(gè)區(qū)域內(nèi)為RANS/LES 混合解析模式,在對(duì)數(shù)律層區(qū)域則為完全LES 模式.由于RANS 模式和RAN/LES 混合解析模式所需網(wǎng)格數(shù)量遠(yuǎn)小于DNS 和準(zhǔn)DNS模式的網(wǎng)格量,而且在混合函數(shù)(39)式—(42)式的雙重保護(hù)下,不僅可以避免MSD 及GIS 的問(wèn)題,同時(shí)也可避免LLM 的問(wèn)題.因此,該WM-HRL模型有望成為高雷諾數(shù)壁湍流三維時(shí)空瞬態(tài)發(fā)展湍流場(chǎng)高置信度解析的一種實(shí)用化的CFD 工具.

本文利用作者團(tuán)隊(duì)自研CFD 軟件NUWA:FLOWUV系統(tǒng),對(duì)如上所述WM-HRL 模型開(kāi)發(fā)了相應(yīng)的植入式程序.該自研CFD 軟件系統(tǒng),采用有限體積法(FVM)求解RANS 方程.其中,壓力速度求解采用PISO 算法,對(duì)流項(xiàng)及擴(kuò)散項(xiàng)的空間離散采用二階中心差分格式,時(shí)間離散格式為二階隱式格式.

3 計(jì)算模型與設(shè)置

本文對(duì)雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)問(wèn)題,采用WM-HRL 進(jìn)行數(shù)值模擬與分析.其中,Re=U∞D(zhuǎn)/v,U∞為來(lái)流速度,D為圓柱直徑.計(jì)算區(qū)域?yàn)榫匦螀^(qū)域,如圖1 所示.具體設(shè)置如下:圓柱底面中心位于坐標(biāo)原點(diǎn) (0,0,0),入口邊界位于x1/D=-10 ;出口邊界位于x1/D=15 ;前后兩個(gè)邊界分別位于x2/D=±4 ;展向高度為L(zhǎng)3/D=π .

圖1 計(jì)算區(qū)域設(shè)置Fig.1.Computational domain schematic.

邊界條件設(shè)置如下:在入口邊界,速度入口設(shè)置為自由來(lái)流條件,即=(U∞,0,0) ;湍動(dòng)能按湍流強(qiáng)度I=5% 設(shè)定,其中k=3(U∞I)2/2 ;比耗率按ω=k/vt設(shè)定,其中vt/v=10.0 .在出口邊界,設(shè)置為零壓梯度出口,即?pˉ=0 .在兩個(gè)垂直側(cè)面及上下邊界,均設(shè)置為對(duì)稱邊界條件.在圓球表面邊界,設(shè)置為無(wú)滑移條件,即==0,湍動(dòng)能設(shè)置為k=0,比耗率設(shè)置為ω=

采用ANSYS ICEM 軟件進(jìn)行分塊網(wǎng)格劃分,如圖2 所示.單元網(wǎng)格為六面體,圓柱體壁面處第一層網(wǎng)格位于y+=0.66 .沿圓柱體展向網(wǎng)格的尺寸為Δ3/D=π/64,水平方向的網(wǎng)格平均尺寸為Δ1(Δ2)/D=0.02.自第一層起,各層網(wǎng)格在徑向的增長(zhǎng)率為1.07,在邊界層內(nèi)總計(jì)布置了45 層網(wǎng)格,整個(gè)計(jì)算域的網(wǎng)格數(shù)量為143 萬(wàn).

圖2 計(jì)算網(wǎng)格剖面Fig.2.Computational grid configuration.

對(duì)SSTk-ω 模型,在鈍體近壁區(qū)通常采用所謂的壁面函數(shù)模型[42].該壁面函數(shù)模型對(duì)湍動(dòng)能k、比耗率ω 及壁面剪切應(yīng)力等,在黏性底層(y+＜5)均給出了明確的邊界條件.在本文所構(gòu)建的WM-HRL 模型中,在鈍體近壁區(qū)采用的是SSTk-ω 模型,為使所構(gòu)建的WM-HRL 模型在鈍體近壁區(qū)發(fā)揮出實(shí)際的壁面?；饔?第一層網(wǎng)格一般需要設(shè)置在黏性底層的y+＜1 內(nèi).

對(duì)雷諾數(shù)Re=3900 下的圓柱繞流問(wèn)題,本文通過(guò)系列數(shù)值模擬研究表明,將第一層網(wǎng)格設(shè)置在y+=0.66 處,可同時(shí)兼顧計(jì)算精度及網(wǎng)格總量控制等要求.另外,Lehmkuhl等[10]采用DNS 的研究表明,在湍流核心區(qū)Kolmogorov 尺度的平均值為D=0.02,為保證網(wǎng)格密度能夠捕捉到最小尺度之圓柱繞流的相干結(jié)構(gòu),網(wǎng)格尺寸需要滿足:=0.9.有鑒于此,本文將徑向的增長(zhǎng)率設(shè)置為1.07.

經(jīng)統(tǒng)計(jì),對(duì)雷諾數(shù)Re=3900 的情況,以圓柱體展向高度L3/D=π 為對(duì)象進(jìn)行CFD 計(jì)算的主要文獻(xiàn)如表1 所示.其中,Pereira等[2]采用的是L3/D=3.0 .表中還列出了展向網(wǎng)格分辨率Δ3/D以及所用總網(wǎng)格量等信息,并與本文所采用的相關(guān)網(wǎng)格信息進(jìn)行比較.

表1 在雷諾數(shù)Re=3900 下圓柱繞流文獻(xiàn)中所用計(jì)算模型與網(wǎng)格參數(shù)設(shè)置情況比較Table 1.Comparisons of computational models and grid parameters in references for flow around a circular cylinder at Reynolds number Re=3900.

由表1 可知,Lehmkuhl等[10]采用的展向網(wǎng)格分辨率最高,所用網(wǎng)格量也最大.Tremblay[8]和Breuer[15]采用的展向網(wǎng)格分辨率一致,但水平方向的網(wǎng)格分辨率不同.Luo等[24]采用的展向網(wǎng)格分辨率略低于前三篇文獻(xiàn)的情況,但水平方向的網(wǎng)格分辨率高于Breuer[15].Pereira等[2]和D'Alessandro等[30]采用的展向網(wǎng)格分辨率最低,兩者的總網(wǎng)格量接近.本文所采用的展向網(wǎng)格分辨率與Tremblay[8]和Breuer[15]的一致,但水平方向的網(wǎng)格分辨率均要低于這兩篇文獻(xiàn)的情況.總之,本文所用計(jì)算網(wǎng)格數(shù)量均少于DNS[10],LES[8,15],PANS[2,24]及DES類模型[24,30]的計(jì)算網(wǎng)格數(shù)量.

在數(shù)值計(jì)算中,無(wú)量綱化時(shí)間步長(zhǎng)Δt*(=ΔtU∞/D) 取值為 6.8×10-3,庫(kù)朗數(shù)CFL＜1,計(jì)算時(shí)長(zhǎng)為70 個(gè)泄渦周期,并對(duì)后45 個(gè)泄渦周期的數(shù)據(jù)做統(tǒng)計(jì)平均,以獲取圓柱繞流場(chǎng)統(tǒng)計(jì)量等信息,并與文獻(xiàn)中相關(guān)實(shí)驗(yàn)結(jié)果及CFD 數(shù)值模擬結(jié)果進(jìn)行比較分析,如表2 所示.

表2 文獻(xiàn)中雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)相關(guān)統(tǒng)計(jì)量的實(shí)驗(yàn)和數(shù)值結(jié)果Table 2.Experimental and numerical results for flow statistical characteristics from references for flow around a circular cylinder at Reynolds numbers Re=3900.

在表2 中,對(duì)PANS 模型,濾波控制函數(shù)fk定義為fk=ku/k.Lehmkuhl等[10]采用了兩種DNS模式.其中,“Mode H”表示“高能模態(tài)”,利用該模態(tài)的DNS 得到V 形結(jié)果.而“Mode L”表示“低能模態(tài)”,利用該模態(tài)的DNS 得到U 形結(jié)果.

由表2 可知,對(duì)Pereira等[2](PANS) (fk=0.25,0.5)的情況,與表中所列實(shí)驗(yàn)結(jié)果[18,27]相比,分離角的計(jì)算結(jié)果均明顯偏小.對(duì)Luo 等(PANS)[24](fk=0.5)的情況,分離角?s、阻力系數(shù)Cd及基線壓力系數(shù) -Cpb的計(jì)算結(jié)果與表中所列實(shí)驗(yàn)結(jié)果結(jié)果[27]相比均明顯偏大,相對(duì)誤差分別達(dá)到9.176%,37.76%和63.333%.對(duì)D'Alessandro 等(SA-DES)[30]的情況,阻力系數(shù)Cd的計(jì)算結(jié)果與表中所列實(shí)驗(yàn)結(jié)果[27]相比均明顯偏大,相對(duì)誤差達(dá)22.7%.

由表2 可進(jìn)一步發(fā)現(xiàn),在Parnaudeau等[18]和Lourenco 和Shih[27]的實(shí)驗(yàn)文章中未給出剪切層小尺度K-H 不穩(wěn)定性的頻譜特征及其無(wú)因次頻率的實(shí)驗(yàn)結(jié)果,而在Tremblay[8],Breuer[15],Luo等[24]及D'Alessandro等[30]的CFD 計(jì)算中,也未給出剪切層小尺度K-H 不穩(wěn)定性的頻譜特征及其無(wú)因次頻率的計(jì)算結(jié)果,但Lehmkuhl等[10]和Pereira等[2]給出了剪切層小尺度K-H 不穩(wěn)定性的頻譜特征及其無(wú)因次頻率的計(jì)算結(jié)果.

Parnaudeau等[18]的實(shí)驗(yàn)測(cè)得回流區(qū)長(zhǎng)度為L(zhǎng)r/D=1.51,并得到U 形流向速度剖面.Lourenco和Shih[27]的實(shí)驗(yàn)測(cè)得回流區(qū)長(zhǎng)度為L(zhǎng)r/D=1.18,并得到V 形流向速度剖面.Lehmkuhl等[10]利用DNS 之Mode H 計(jì)算得到的回流區(qū)長(zhǎng)度Lr/D=1.26,與Lourenco 和Shih[27]實(shí)驗(yàn)結(jié)果的相對(duì)誤差為6.78%.而Lehmkuhl等[10]利用DNS 之Mode L計(jì)算得到回流區(qū)長(zhǎng)度Lr/D=1.55,與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果的相對(duì)誤差為2.65%.

Tremblay[8]及Breuer[15]采用LES 均得到V形剖面,但是他們的CFD 計(jì)算所得回流區(qū)長(zhǎng)度與Lourenco 和Shih[27]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度有較大差異,相對(duì)誤差分別達(dá)到-11.86%和16.27%.Pereira等[2]采用PANS在fk=0.25,0.5 下分別得到U 形和V 形剖面.當(dāng)fk=0.5 時(shí),CFD 計(jì)算所得回流區(qū)長(zhǎng)度與Lourenco 和Shih[27]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度接近,相對(duì)誤差為-5.08%.當(dāng)fk=0.25時(shí),CFD 計(jì)算所得回流區(qū)長(zhǎng)度與Parnaudeau等[18]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度相比偏大,相對(duì)誤差達(dá)14.56%.此外,對(duì)Pereira等[2](PANS) (fk=0.5)的情況,與表中所列DNS[10]結(jié)果相比,K-H 不穩(wěn)定性頻率的計(jì)算結(jié)果明顯偏大,相對(duì)誤差達(dá)15.67%.

Luo等[24]采用PANS在fk=0.1 和0.5 下均得到V 形剖面.當(dāng)fk=0.1 時(shí),CFD 計(jì)算所得回流區(qū)長(zhǎng)度與Lourenco 和Shih[27]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度的相對(duì)誤差為7.63%.當(dāng)fk=0.5 時(shí),CFD計(jì)算所得回流區(qū)長(zhǎng)度與Lourenco等[27]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度相比明顯不符.Luo等[24]采用SSTDES 得到V 形剖面,但計(jì)算所得回流區(qū)長(zhǎng)度與Lourenco 和Shih[27]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度相比也明顯不符.

D'Alessandro等[30]采用SA-DES 也得到V 形剖面,但計(jì)算所得回流區(qū)長(zhǎng)度與Lourenco 和Shih[27]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度相比明顯不符.D'Alessandro等[30]采用SA-IDDES 及v2-f DES 均得到U 形剖面,采用SA-IDDES 計(jì)算得到的回流區(qū)長(zhǎng)度與Parnaudeau等[18]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度的相對(duì)誤差為-5.50%.但采用v2-f DES 計(jì)算得到的回流區(qū)長(zhǎng)度與Parnaudeau等[18]實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度相比偏大,相對(duì)誤差高達(dá)11.13%.

對(duì)鈍體繞流問(wèn)題,其邊界層包括黏性底層、過(guò)渡層及對(duì)數(shù)律層.對(duì)本文WM-HRL 模型來(lái)講,這三類邊界層范圍的信息將是至關(guān)重要的.為此,利用圖2 所示計(jì)算網(wǎng)格系統(tǒng),首先進(jìn)行RANS 數(shù)值模擬,結(jié)果表明:在雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)的邊界層厚度為y+≤180.0,黏性底層厚度為y+≤10.0,過(guò)渡層厚度為 10.0＜y+＜30.0 .

對(duì)雷諾數(shù)Re=3900 下的圓柱繞流場(chǎng)問(wèn)題,其相干結(jié)構(gòu)主要包含兩類結(jié)構(gòu).其中,第一類為與尾流中渦泄相關(guān)的大尺度不穩(wěn)定性結(jié)構(gòu),其無(wú)因次特征頻率為,而第二類為與分離剪切層脈動(dòng)相關(guān)的小尺度不穩(wěn)定性結(jié)構(gòu),稱為K-H 不穩(wěn)定性,其無(wú)因次特征頻率為.對(duì)本文WM-HRL 模型來(lái)講,圓柱繞流剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生的位置信息將是至關(guān)重要的.

剪切層的速度梯度較大,屬于一個(gè)位置狹小的窄帶,且對(duì)監(jiān)測(cè)點(diǎn)的分布位置敏感[23].為獲取圓柱繞流剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生的準(zhǔn)確邊界位置,結(jié)合相關(guān)文獻(xiàn)中DNS[10],PANS[2]及LES[23]等CFD 數(shù)值模擬結(jié)果,本文設(shè)置了如圖3所示的18 個(gè)監(jiān)測(cè)點(diǎn),其坐標(biāo)位置如表3 所示.其中,監(jiān)測(cè)點(diǎn)P1—P14 位于剪切層區(qū)域,監(jiān)測(cè)點(diǎn)P15 位于剪切層脫落區(qū),而監(jiān)測(cè)點(diǎn)P16—P18 位于尾流中心線上.表3 中所列的這些監(jiān)測(cè)點(diǎn)均位于圓柱體底部位置.同時(shí),對(duì)表3 中所列的每個(gè)監(jiān)測(cè)點(diǎn),沿圓柱體展向同時(shí)均勻布置了其他4 個(gè)相應(yīng)的監(jiān)測(cè)點(diǎn).因此,實(shí)際上總計(jì)布置了90 個(gè)監(jiān)測(cè)點(diǎn).

表3 監(jiān)測(cè)點(diǎn)坐標(biāo)信息Table 3.Coordinate information of the probes.

圖3 剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)監(jiān)測(cè)點(diǎn)分布Fig.3.Location configuration of the probes for the small scale K-H instability structure in the shear layer.

在此基礎(chǔ)上,利用本文WM-HRL 進(jìn)行系列數(shù)值模擬,對(duì)所得各監(jiān)測(cè)點(diǎn)處的橫向脈動(dòng)速度進(jìn)行Lomb 頻譜分析,詳細(xì)分析見(jiàn)4.3 節(jié).在做Lomb譜分析時(shí),用的是表3 中所列各監(jiān)測(cè)點(diǎn)相對(duì)應(yīng)的5 個(gè)監(jiān)測(cè)點(diǎn)處橫向脈動(dòng)速度做展向平均所得,其做法與Lehmkuhl等[10](DNS)的做法相同.

結(jié)果表明,在各個(gè)監(jiān)測(cè)點(diǎn)的Lomb 頻譜中,均出現(xiàn)一個(gè)頻率相對(duì)較低的譜峰,該譜峰所對(duì)應(yīng)頻率與渦泄頻率一致.同時(shí),對(duì)各種工況的系列CFD數(shù)值模擬結(jié)果中,在測(cè)點(diǎn)P4 的Lomb 頻譜中,均觀察到一個(gè)頻率相對(duì)較高的譜峰,該譜峰所對(duì)應(yīng)頻率與剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的頻率一致.由此可見(jiàn),雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)中剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生在對(duì)數(shù)律層中y+≥89.4的區(qū)域.

在圖2 所述網(wǎng)格系統(tǒng)設(shè)置下,利用(36)式進(jìn)行計(jì)算,結(jié)果表明:在對(duì)數(shù)律層的y+≥67.4 的區(qū)域中,rk均不大于0.2,即rk≤0.2 .由于圓柱繞流剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生在y+≥89.4的區(qū)域,因此本文WM-HRL 模型對(duì)該剪切層小尺度K-H 不穩(wěn)定結(jié)構(gòu)為完全LES 解析模式,且其對(duì)湍動(dòng)能的解析度至少為80%.

4 數(shù)值結(jié)果與分析

當(dāng)利用WM-HRL 進(jìn)行數(shù)值模擬時(shí),影響圓柱繞流場(chǎng)計(jì)算置信度的主要因素包括兩個(gè)邊界位置和3 個(gè)區(qū)域的湍動(dòng)能解析度.其中,兩個(gè)邊界分別為RANS 結(jié)束邊界(記為ΓRANS)和LES 啟動(dòng) 邊界(記為ΓLES),3 個(gè)區(qū)域分別為RANS 區(qū)、RANS/LES 混合轉(zhuǎn)換區(qū)和LES 區(qū).為下文陳述簡(jiǎn)便計(jì),記為RANS 結(jié)束邊界ΓRANS的位置,而記為L(zhǎng)ES 啟動(dòng)邊界ΓLES的位置.

4.1 WM-HRL 模型參數(shù)影響規(guī)律

首先討論LES 啟動(dòng)邊界ΓLES位于剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域的情況,并考慮兩種情況:第一種情況為=105.8,rk2=0.1556 ;第 二種情況為=113.9,rk2=0.1484 .在每種情況下,均設(shè)置11 種RANS 結(jié)束邊界ΓRANS的位置及其相應(yīng)的rk1,如表4 所示.其中,=7.9位于黏性底層,=13.3,14.9,18.4,20.4,27.1 位于過(guò)渡層,=38.4,41.7,49.2,72.7 位于對(duì)數(shù)律層但在剪切層K-H 不穩(wěn)定區(qū)外,而=91.2則位于剪切層K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域內(nèi).

表4 當(dāng) ΓLES 位于剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域內(nèi)時(shí),相關(guān)流場(chǎng)統(tǒng)計(jì)量的數(shù)值結(jié)果Table 4.Numerical results for flow statistic characteristics when ΓLES is located in the K-H instability region of the shear layer.

下面討論LES 啟動(dòng)邊界ΓLES位于剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域外,但仍位于對(duì)數(shù)律層區(qū)的情況,并考慮3 種情況:第一種情況為=49.2,rk2=0.2546 ;第二種 情況為=72.4,rk2=0.1983 ;第三種情況為=84.6,rk2=0.1713.在每種情況下,均設(shè)置若干種RANS 結(jié)束邊界ΓRANS的位置及其相應(yīng)的rk1,其設(shè)置原則與表4 一致,如表5 所示.在各種及rk1和及rk2的組合下,利用WM-HRL 進(jìn)行數(shù)值模擬所得相關(guān)流場(chǎng)統(tǒng)計(jì)量的結(jié)果如表5 所示.

表5 當(dāng) ΓLES 位于剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域外且在對(duì)數(shù)律層內(nèi)時(shí),相關(guān)流場(chǎng)統(tǒng)計(jì)量的數(shù)值結(jié)果Table 5.Numerical results for flow statistic characteristics when ΓLES is located in the log-law layer and outside the K-H instability region of the shear layer.

由表5 可知,只有當(dāng)RANS 結(jié)束邊界ΓRANS位于過(guò)渡層,且rk1≤=0.63 時(shí),才能同時(shí)準(zhǔn)確獲得和Lr/D的計(jì)算結(jié)果.在此條件下,當(dāng)計(jì)算得到V 形速度剖面時(shí),Lr/D的計(jì)算結(jié)果與Lourenco和Shih[27]的實(shí)驗(yàn)結(jié)果相比,相對(duì)誤差最大為-9.32%,而的計(jì)算結(jié)果與Lehmkuhl等[10]的DNS 之Mode H 計(jì)算結(jié)果相比,相對(duì)誤 差最大為12.68%.當(dāng)計(jì)算得到U 形速度剖面時(shí),Lr/D的計(jì)算結(jié)果與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果相比,相對(duì)誤差最大為-12.58%,而的計(jì)算結(jié)果與Lehmkuhl等[10]的DNS 之Mode H 計(jì)算結(jié)果相比,相對(duì)誤差最大為-20.15%.

最后討論LES 啟動(dòng)邊界ΓLES位于過(guò)渡層的情況,并考慮5 種情況,分別為=29.6 .在每種情況下,均 設(shè)置若干種RANS 結(jié)束邊界ΓRANS的位置及其相應(yīng)的rk1.在各種的組合下,利用WM-HRL 進(jìn)行數(shù)值模擬所得相關(guān)流場(chǎng)統(tǒng)計(jì)量的結(jié)果如表6 所示.

表6 當(dāng) ΓLES 位于過(guò)渡層時(shí),相關(guān)流場(chǎng)統(tǒng)計(jì)量的數(shù)值結(jié)果Table 6.Numerical results for flow statistic characteristics when ΓLES is located in the buffer layer.

需要指出的是,在上文統(tǒng)計(jì)中,對(duì)?s之CFD計(jì)算值,已去掉了Pereira等[2](PANS) (fk=0.25和0.5)的異常計(jì)算結(jié)果,以及Luo等[24](PANS)(fk=0.5)的異常計(jì)算結(jié)果.對(duì)Cd之CFD 計(jì)算值,已去掉Luo等[24](PANS)(fk=0.5)的異常計(jì)算結(jié)果,以及D'Alessandro等[30](SA-DES)的異常計(jì)算結(jié)果.對(duì)-Cpb之CFD 計(jì)算值,已去掉Luo等[24](PANS)(fk=0.5)的異常計(jì)算結(jié)果.

由于剪切層小尺度K-H 不穩(wěn)定性發(fā)生的位置是未知的,但位于對(duì)數(shù)律層區(qū)中.系列數(shù)值模擬結(jié)果表明,對(duì)本文所構(gòu)造的WM-HRL 模型,為了能夠準(zhǔn)確解析并捕捉剪切層小尺度K-H 不穩(wěn)定性的結(jié)構(gòu)特征及其頻譜特性,可將其LES 模式的啟動(dòng)邊界均設(shè)置在過(guò)渡層內(nèi),而RANS 模式的結(jié)束邊界則既可設(shè)置在過(guò)渡層內(nèi)也可設(shè)置在黏性底層內(nèi),并且使rk2≤0.2,即在對(duì)數(shù)律層區(qū)的網(wǎng)格具有至少80%的湍動(dòng)能解析度.

在圓柱體展向長(zhǎng)度及其網(wǎng)格系統(tǒng)已經(jīng)確定的條件下,對(duì)DNS[6-12]及LES[6-8,13-23]來(lái)講,通過(guò)CFD計(jì)算只能得到回流區(qū)長(zhǎng)度及流向速度剖面形狀的其中一種分支結(jié)構(gòu).對(duì)DES 類模型[24,30]來(lái)講,不同RANS 湍流模型的類型及RANS/LES 之間不同的混合轉(zhuǎn)換方式等,均可能會(huì)影響其對(duì)回流區(qū)長(zhǎng)度及流向速度剖面形狀分支結(jié)構(gòu)的計(jì)算結(jié)果.

對(duì)PANS 模型[2,24-26]而言,其核心思想是通過(guò)引入一個(gè)fk參數(shù)來(lái)調(diào)控流場(chǎng)可解/不可解湍流量的比例而實(shí)現(xiàn).對(duì)基于SSTk-ω 的PANS 模型,fk參數(shù)可調(diào)控其ω 方程中耗散項(xiàng)之值的變化.一般地,fk值越小,耗散項(xiàng)的值也越小,從而使求解得到的比耗率ω 增大,進(jìn)而使湍流渦黏系數(shù)vt減小,可解尺度就釋放的越多[25].

對(duì)PANS 模型,Pereira等[2]指出:為準(zhǔn)確獲取圓柱繞流回流區(qū)結(jié)構(gòu)及其長(zhǎng)度信息,fk的值不能大于0.5,即fk≤0.5 .同時(shí),當(dāng)fk從0.5 減小到某個(gè)值(文中為0.25)后,圓柱尾部近壁面處流向速度剖面的形狀從V 形轉(zhuǎn)變?yōu)閁 形.但Luo等[24]在fk=0.5 時(shí)通過(guò)CFD 計(jì)算并沒(méi)有得到準(zhǔn)確的回流區(qū)長(zhǎng)度,而且在fk=0.1 時(shí),得到的是V 形剖面.導(dǎo)致這種相矛盾結(jié)果的原因之一,可能與兩位作者所使用網(wǎng)格結(jié)構(gòu)及其空間分辨率分布不同有關(guān)(具體見(jiàn)表2).

本文采用WM-HRL 模型的系列數(shù)值模擬結(jié)果表明,在同一套網(wǎng)格系統(tǒng)下,通過(guò)改變WM-HRL模型中兩個(gè)邊界位置及3 個(gè)區(qū)域的湍動(dòng)能解析度的組合,既可以通過(guò)數(shù)值模擬獲得與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果一致的回流區(qū)長(zhǎng)度及圓柱近壁面處平均流向速度的U 形剖面,也可以通過(guò)數(shù)值模擬獲得與Lourenco 和Shih[27]實(shí)驗(yàn)結(jié)果一致的回流區(qū)長(zhǎng)度及圓柱近壁面處平均流向速度的V 形剖面.

這表明對(duì)本文所建立的WM-HRL 模型,可以在同一套網(wǎng)格系統(tǒng)下通過(guò)變化湍動(dòng)能解析度指標(biāo)參數(shù)nrk1-Q和nrk2-Q準(zhǔn)則的組合條件,精細(xì)地解析圓柱繞流場(chǎng)中兩類不同回流區(qū)長(zhǎng)度結(jié)構(gòu)特征,及其對(duì)應(yīng)的圓柱尾部近壁面處U 形和V 形兩個(gè)流向速度剖面的分支結(jié)構(gòu).

4.2 一階和二階統(tǒng)計(jì)量特性

選取表6 中相關(guān)數(shù)值模擬結(jié)果案列,討論雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)的一階和二階統(tǒng)計(jì)量特性的計(jì)算結(jié)果,并與Parnaudeau等[18]及Lourenco和Shih[27]的相關(guān)實(shí)驗(yàn)結(jié)果進(jìn)行比較分析.對(duì)U形流向速度剖面的情況,選取4 個(gè)組合工況如下.Case AU:=18.4;Case BU:=7.9,=29.6.其中,對(duì)Case AU 工況,其RANS 結(jié)束邊界位于黏性底層和過(guò)渡層交界處,而LES 啟動(dòng)邊界位于過(guò)渡層內(nèi);對(duì)Case BU 工況,其RANS 結(jié)束邊界和LES 啟動(dòng)邊界均位于過(guò)度層內(nèi);對(duì)Case CU 工況,其RANS 結(jié)束邊界位于黏性底層內(nèi),而LES啟動(dòng)邊界位于過(guò)渡層內(nèi);對(duì)Case DU 工況,其RANS 結(jié)束邊界位于過(guò)渡層內(nèi),而LES 啟動(dòng)邊界位于過(guò)渡層與對(duì)數(shù)律層交界處.

對(duì)V 形流向速度剖面的情況,選取4 個(gè)組合工況如下.Case AV:=13.3;Case BV:=14.9,=29.6.其中,對(duì)Case AV 工況,其RANS 結(jié)束邊界位于黏性底層內(nèi),而LES 啟動(dòng)邊界位于黏性底層與過(guò)渡層交界處;對(duì)Case BV 工況,其RANS 結(jié)束邊界位于黏性底層和過(guò)渡層交界處,而LES 啟動(dòng)邊界位于過(guò)渡層內(nèi);對(duì)Case CV 工況,其RANS結(jié)束邊界和LES 啟動(dòng)邊界均位于過(guò)渡層內(nèi);對(duì)Case DV 工況,其RANS 結(jié)束邊界位于過(guò)渡層內(nèi),而LES 啟動(dòng)邊界位于過(guò)渡層與對(duì)數(shù)律層交界處.

在針對(duì)一階和二階統(tǒng)計(jì)量的分析中,對(duì)Case AU—Case DU 這4 種工況,除了將其與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果進(jìn)行比較外,同時(shí)還與Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果進(jìn)行了比較.對(duì)Case AV—Case DV 這4 種工況,除了將其與Lourenco 和Shih[27]的實(shí)驗(yàn)進(jìn)行比較外,同時(shí)還與Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算結(jié)果進(jìn)行了比較.

4.2.1 一階統(tǒng)計(jì)量特性

在圖4 中,分別給出了上述8 種工況下圓柱表面周向系數(shù)Cp分布特性的數(shù)值結(jié)果.由于Parnaudeau等[18]和Lourenco 和Shih[27]的實(shí)驗(yàn)沒(méi)有給出圓柱表面壓力系數(shù)的數(shù)據(jù),此處采用Norberg[43]的實(shí)驗(yàn)結(jié)果進(jìn)行比較分析.

圖4 圓柱表面周向壓力系數(shù) Cp 分布特性Fig.4.Azimuthal distribution characteristics for pressure coefficient along the circular cylinder surface.

由圖4 可知,無(wú)論是對(duì)Case AU—Case DU這4 種工況,還是對(duì)Case AV—Case DV 這4 種工況,利用本文WM-HRL 模型所得沿圓柱表面周向壓力系數(shù)Cp分布的數(shù)值結(jié)果之間的差異均很小.對(duì)圖4(a)的情況,本文計(jì)算結(jié)果與Norberg[43]實(shí)驗(yàn)結(jié)果相比的相對(duì)誤差最大為10.3%,而Lehmkuhl等[10]利用DNS(Mode L)計(jì)算結(jié)果與Norberg[43]實(shí)驗(yàn)結(jié)果相比的相對(duì)誤差最大為3.1%.對(duì)圖4(b)的情況,本文計(jì)算結(jié)果與Norberg[43]實(shí)驗(yàn)結(jié)果相比的相對(duì)誤差最大為14.0%,而Lehmkuhl等[10]利用DNS(Mode H)計(jì)算結(jié)果Norberg[43]實(shí)驗(yàn)結(jié)果相比的相對(duì)誤差最大為0.6%.

在圖5 中,分別給出了前述8 種工況下沿尾流中心線平均流向速度剖面特性的數(shù)值結(jié)果.由圖5可知,在圓柱表面正后方中心線上的點(diǎn)P(0.5D,0)處,平均流向速度的值為0,在達(dá)到一個(gè)負(fù)的最小值后開(kāi)始增大,并在中心線上的點(diǎn)Q(Lr+0.5D,0)處又變?yōu)?.其中,Lr/D為回流區(qū)長(zhǎng)度.Parnaudeau等[18]的實(shí)驗(yàn)測(cè)量獲得Lr/D=1.51,Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算得到Lr/D=1.55,如圖5(a)所示.而Lourenco 和Shih[27]的實(shí)驗(yàn)測(cè)量獲得Lr/D=1.18,Lehmkuhl等[10]利用Mode H之DNS 計(jì)算得到Lr/D=1.26,如圖5(b)所示.此即在雷諾數(shù)Re=3900 下圓柱繞流兩個(gè)實(shí)驗(yàn)及DNS 計(jì)算中出現(xiàn)兩類不同回流區(qū)分支結(jié)構(gòu)的情況.

圖5 沿尾流中心線平均流向速度剖面特性Fig.5.Distribution characteristics of mean stream-wise velocities along the wake centerline.

由圖5 進(jìn)一步可發(fā)現(xiàn),無(wú)論是對(duì)Case AU—Case DU 這4 種工況,還是對(duì)Case AV—Case DV這4 種工況,利用本文WM-HRL 模型所得沿尾流中心線平均流向速度分布的數(shù)值結(jié)果之間的差異均很小.對(duì)圖5(a)的情況,在回流區(qū)外,本文數(shù)值結(jié)果與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果之間的相對(duì)誤差最大為5.5%,而與Lehmkuhl等[10]利用Mode L 之DNS 的計(jì)算結(jié)果之間的相對(duì)誤差最大為3.2%.

對(duì)圖5(b)的情況,本文數(shù)值結(jié)果與Lehmkuhl等[10]利用Mode H 之DNS 的計(jì)算均吻合良好,兩者之間的相對(duì)誤差最大為4.2%.但無(wú)論是本文計(jì)算結(jié)果還是Lehmkuhl等[10]利用Mode H 之DNS的計(jì)算結(jié)果,在x1/D∈[2.36,3.26]的范圍內(nèi),它們與Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果均有一定的差異.

圖6 中,分別給出了前述8 種工況下在3 個(gè)不同站位(x1/D=1.06,1.54,2.02)處平均流向速度剖面特性的數(shù)值結(jié)果.站位x1/D=1.06 位于剪切層K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域,同時(shí)也位于回流區(qū)內(nèi).Parnaudeau等[18]的實(shí)驗(yàn)測(cè)量獲得U形速度剖面,而Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算也得到U 形剖面,如圖6(a)所示.而Lourenco 和Shih[27]的實(shí)驗(yàn)測(cè)量獲得V 形速度剖面,而Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算也得到V 形剖面,如圖6(b)所示.此即在雷諾數(shù)Re=3900 下,圓柱繞流兩個(gè)試驗(yàn)中所測(cè)量得到的平均流向速度剖面出現(xiàn)U 形和V 形兩個(gè)流動(dòng)分支結(jié)構(gòu)的情況.

圖6 圓柱后方不同站位處平均流向速度剖面特性Fig.6.Distribution characteristics of mean stream-wise velocities at different locations in the backside of the circular cylinder.

由圖6 可知,無(wú)論是對(duì)Case AU—Case DU 這4 種工況,還是對(duì)Case AV—Case DV 這4 種工況,利用本文WM-HRL 模型所得該站位處平均流向速度分布的數(shù)值結(jié)果之間的差異均很小.在位于剪切層K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域的站位x1/D=1.06 處,對(duì)Case AU—Case DU 這4 種工況,本文計(jì)算所得平均流向速度剖面均為U 形,與Parnaudeau等[18]的實(shí)驗(yàn)及Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算所得速度剖面形狀均一致.對(duì)Case AV—Case DV 這4 種工況,本文計(jì)算所得平均流向速度剖面均為V 形,與Lourenco 和Shih[27]的實(shí)驗(yàn)及Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算所得速度剖面形狀均一致.

站位x1/D=1.54 和2.02 位于回流區(qū)中.由圖6可見(jiàn),Parnaudeau等[18]和Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果均獲得V 形流向速度剖面,而Lehmkuhl等[10]利用Mode L 和H 之DNS 計(jì)算也均獲得V形流向速度剖面.同時(shí),由圖6 還可觀察到,對(duì)Case AU-Case DU 這4 種工況,利用本文WMHRL 模型計(jì)算所得平均流向速度剖面均為V 形,與Parnaudeau等[18]的實(shí)驗(yàn)及Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算所得速度剖面形狀一致.對(duì)Case AV—Case DV 這4 種工況,利用本文WMHRL 模型計(jì)算所得平均流向速度剖面也均為V 形,與Lourenco 和Shih[27]的實(shí)驗(yàn)及Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算所得速度剖面形狀均一致.

圖7 分別給出了前述8 種工況下在3 個(gè)不同站位(x1/D=1.06,1.54,2.02)處平均橫向速度剖面特性的數(shù)值結(jié)果.

圖7 圓柱后方不同站位處平均橫向速度剖面特性Fig.7.Distribution characteristics of mean cross-flow velocities at different locations in the backside of the circular cylinder.

由圖7 可見(jiàn),無(wú)論是對(duì)Case AU—Case DU 這4 種工況,還是對(duì)Case AV—Case DV 這4 種工況,利用本文WM-HRL 模型所得圓柱后方不同站位處平均橫向速度分布的數(shù)值結(jié)果差異均很小.

對(duì)圖7(a)的情況,在x1/D=1.06 站位處,本文計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果之間的相對(duì)誤差最大為3.5%,而與Lehmkuhl等[10]利用Mode L 的DNS 計(jì)算結(jié)果之間的相對(duì)誤差最大為3.8%.在x1/D=1.54 站位處,本文計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果之間的相對(duì)誤差最大為9.9%,而與Lehmkuhl等[10]利用Mode L 的DNS計(jì)算結(jié)果之間的相對(duì)誤差最大為10.3%.在x1/D=2.02 站位處,本文計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果之間的相對(duì)誤差最大為15%,而與Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算結(jié)果之間的相對(duì)誤差最大為16.1%.

對(duì)圖7(b)的情況,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果一致,但兩者與Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果均有一定的差異.在x1/D=1.06 站位處,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 的DNS 計(jì)算結(jié)果之間的相對(duì)誤差最大為7.8%.在x1/D=1.54 站位處,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果之間的相對(duì)誤差最大為5.4%.在x1/D=2.02 站位處,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果之間的相對(duì)誤差最大為10.0%.

Tremblay[8]采用DNS 和LES 均復(fù)現(xiàn)了Lourenco 和Shih[27]對(duì)雷諾數(shù)Re=3900 下,圓柱繞流在其后方不同站位處平均流向及橫向速度剖面特性的實(shí)驗(yàn)結(jié)果.其中,LES 亞格子模型分別采用了經(jīng)典Smagorinsky 模型(Smag)及動(dòng)態(tài)模型(Dyn).結(jié)果表明,DNS的結(jié)果要優(yōu) 于LES的結(jié)果,而Dyn 亞格子模型的結(jié)果則要優(yōu)于Smag 亞格子模型的結(jié)果.

對(duì)PANS 模型,其對(duì)圓柱后方各站位處平均流向及橫向速度剖面的計(jì)算精度依賴于濾波控制參數(shù)fk的取值方式.對(duì)DES 類模型,其對(duì)圓柱后方各站位處平均流向及橫向速度剖面的計(jì)算精度則依賴于所用RANS 模型的類型及RANS/LES的混合轉(zhuǎn)換方法等.總體上,在網(wǎng)格空間分辨率分布及相關(guān)的模型參數(shù)設(shè)置合理的條件下,這兩類模型對(duì)圓柱后方各站位處平均流向及橫向速度剖面的計(jì)算精度,一般地可以達(dá)到與LES 相當(dāng)?shù)挠?jì)算精度[2,24-26,30].

本文采用WM-HRL 模型的系列數(shù)值模擬結(jié)果表明,在同一套網(wǎng)格系統(tǒng)下,通過(guò)改變WM-HRL模型中兩個(gè)邊界位置及3 個(gè)區(qū)域的湍動(dòng)能解析度的組合,既可以通過(guò)數(shù)值模擬獲得與Parnaudeau等[18]實(shí)驗(yàn)所得一致的圓柱后方各站位處的平均流向及橫向速度剖面,也可以通過(guò)數(shù)值模擬獲得與Lourenco 和Shih[27]實(shí)驗(yàn)所得一致的圓柱后方各站位處的平均流向及橫向速度剖面.

4.2.2 二階統(tǒng)計(jì)量特性

圖8 分別給出了前述8 種工況下在3 個(gè)不同站位(x1/D=1.06,1.54,2.02)處各向同性流向雷諾應(yīng)力剖面特性的數(shù)值結(jié)果.

圖8 圓柱后方不同站位處各向同性流向雷諾應(yīng)力剖面特性Fig.8.Distribution characteristics of isotropic stream-wise Reynolds stresses at different locations in the backside of the circular cylinder.

由圖8 可見(jiàn),對(duì)Case AU—Case DU 的情況,在位于剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域的站位x1/D=1.06 處,本文計(jì)算結(jié)果與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果及Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算結(jié)果相比,主要差異在兩個(gè)“峰”處.在位于回流區(qū)的站位x1/D=1.54,2.02處,本文計(jì)算結(jié)果與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果及Lehmkuhl等[10]利用Mode L 的DNS 計(jì)算結(jié)果相比,在兩個(gè)“峰”處的值小于相應(yīng)的實(shí)驗(yàn)值及DNS 計(jì)算值.

對(duì)Case AV—Case DV 的情況,本文計(jì)算結(jié)果均與Lourenco 和Shih[27]實(shí)驗(yàn)結(jié)果及Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果具有良好的一致性.在位于剪切層小尺度K-H 不穩(wěn)定區(qū)的站位x1/D=1.06 處,本文對(duì)BV 和DV 工況的計(jì)算結(jié)果與Lourenco 和Shih[27]實(shí)驗(yàn)結(jié)果之間的最大相對(duì)誤差分別為1.7%和5%,而與Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果之間的最大誤差分別為0.3%和3.7%.在x1/D=1.54 站位處,本文計(jì)算結(jié)果與Lourenco 和Shih[27]實(shí)驗(yàn)結(jié)果之間的最大相對(duì)誤差分別為16.1%,3.4%,3.8%,6.4%.在x1/D=2.02 站位處,本文計(jì)算結(jié)果與Lourenco和Shih[27]實(shí)驗(yàn)結(jié)果之間的最大相對(duì)誤差分別為4.6%,11%,3.5%,6.7%.

圖9 分別給出了前述8 種工況下在3 個(gè)不同站位(x1/D=1.06,1.54,2.02)處各向同性橫向雷諾應(yīng)力剖面特性的數(shù)值結(jié)果.

圖9 圓柱后方不同站位處各向同性橫向雷諾應(yīng)力剖面特性Fig.9.Distribution characteristics of isotropic cross-flow Reynolds stresses at different locations in the backside of the circular cylinder.

由圖9 可見(jiàn),對(duì)Case AU—Case DU 的情況,在位于剪切層小尺度K-H 不穩(wěn)定區(qū)的站位x1/D=1.06 處,本文計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果相比的最大相對(duì)誤差分別為23.6%,7.7%,15.7%,9.1%,而Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果相比的最大相對(duì)誤差為14.5%.

在位于回流區(qū)的站位x1/D=1.54 處,本文計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果相比的最大相對(duì)誤差分別為13.7%,25.3%,6.4%,13.6%,而Lehmkuhl等[10]利用Mode L 之DNS 計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果相比的最大相對(duì)誤差為15%.

在位于回流區(qū)的站位x1/D=2.02 處,本文計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果相比的最大相對(duì)誤差分別為3.3%,7.2%,12.9%,0.5%,而Lehmkuhl等[10]利用Mode L 的DNS 計(jì)算結(jié)果與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果相比的最大相對(duì)誤差為7%.

對(duì)Case AV—Case DV 的情況,無(wú)論是本文計(jì)算結(jié)果還是Lehmkuhl等[10]利用Mode H 的DNS 計(jì)算結(jié)果一致,兩者均與Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果有一定的差異.在位于剪切層小尺度K-H 不穩(wěn)定區(qū)的站位x1/D=1.06 處,本文對(duì)Case BV 和Case DV 工況的計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 的DNS 計(jì)算結(jié)果之間的相對(duì)誤差較小,分別為6.8%和3.7%,而對(duì)Case AV 和CaseCV工況,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H之DNS 計(jì)算結(jié)果的相對(duì)誤差較大,分別為29.5%和29.3%.

在位于回流區(qū)的站位x1/D=1.54 處,對(duì)Case AV—Case DV 工況,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 的DNS 計(jì)算結(jié)果的相對(duì)誤差分別為7.75%,15.2%,14%和1%.在位于回流區(qū)的站位x1/D=2.02 處,對(duì)Case AV—Case DV工況,本文計(jì)算結(jié)果與Lehmkuhl等[10]利用Mode H 之DNS 計(jì)算結(jié)果的相對(duì)誤差分別為1.7%,1.2%,0.5%,2.2%.

圖10 分別給出了上述8 種工況下在3 個(gè)不同站位(x1/D=1.06,1.54,2.02)處各向異性雷諾應(yīng)力剖面特性的數(shù)值結(jié)果.對(duì)此情況,Lehmkuhl等[10]沒(méi)有給出相關(guān)的DNS 之計(jì)算結(jié)果.

圖10 圓柱后方不同站位處各向異性雷諾應(yīng)力剖面特性Fig.10.Distribution characteristics of anisotropy cross-flow Reynolds stresses at different locations in the backside of the circular cylinder.

由圖10 可見(jiàn),對(duì)圖10(a)的情況,在位于剪切層小尺度K-H 不穩(wěn)定性發(fā)生區(qū)域的站位x1/D=1.06 處,本文在Case AU—Case DU 這4 種工況下的計(jì)算結(jié)果均與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果具有良好的一致性.在位于回流區(qū)的站位x1/D=1.54 處,本文對(duì)Case AU 和Case BU 的計(jì)算結(jié)果與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果也具有良好的一致性,但在 |x2/D|＜0.48 的范圍內(nèi)本文對(duì)Case CU和Case DU 的計(jì)算結(jié)果,與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果有一定的差異,最大相對(duì)誤差分別為38.7%和14.8%.在位于回流區(qū)的站位x1/D=2.02處,本文對(duì)Case BU—Case DU 的計(jì)算 結(jié)果 與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果具有良好的一致性,但在 |x2/D|＜0.48 的范圍內(nèi)本文對(duì)Case AU 的計(jì)算結(jié)果,與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果有一定的差異,最大相對(duì)誤差為30.8%.

對(duì)圖10(b)的情況,在位于剪切層小尺度KH 不穩(wěn)定性發(fā)生區(qū)域的站位x1/D=1.06 處,以及在位于回流區(qū)的站位x1/D=1.54 處,本文對(duì)Case AV—Case DV 這4 種工況下的計(jì)算結(jié)果均與Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果具有良好的一致性.在位于回流區(qū)的站位x1/D=2.02 處,本文對(duì)Case AV—Case DV 的計(jì)算結(jié)果與Lourenco和Shih[27]的實(shí)驗(yàn)結(jié)果相比,在“峰”和“谷”處有一定的差異.其中,在“峰”處的相對(duì)誤差分別為30.2%,29.5%,38%,43.8%,而在“谷”處的相對(duì)誤差分別為32.7%,31.8%,41.1%,40.8%.

Tremblay[8]采用DNS 和LES 同時(shí)也均復(fù)現(xiàn)了Lourenco 和Shih[27]對(duì)雷諾數(shù)Re=3900 下,圓柱繞流后方不同站位處雷諾應(yīng)力剖面特性的實(shí)驗(yàn)結(jié)果.其中,LES 亞格子模型分別包括經(jīng)典Smagorinsky 模型(Smago)及動(dòng)態(tài)模型(Dyn).結(jié)果表明,DNS 的結(jié)果要優(yōu)于LES 的結(jié)果,而Dyn 亞格子模型的結(jié)果要優(yōu)于Smago 亞格子模型的結(jié)果.

對(duì)PANS 模型,其對(duì)圓柱后方各站位處雷諾應(yīng)力剖面的計(jì)算精度同樣地依賴于濾波控制參數(shù)fk的取值方式.對(duì)DES 類模型,其對(duì)圓柱后方各站位處雷諾應(yīng)力剖面的計(jì)算精度則依賴于所選用RANS 模型的類型及RANS/LES 之間的混合轉(zhuǎn)換方法等.

由于這兩類模型均采用了相關(guān)的RANS 模型,而RANS 模型對(duì)湍流尺度的解析能力是基于對(duì)渦黏系數(shù)的調(diào)控而實(shí)現(xiàn),但這種調(diào)控與計(jì)算網(wǎng)格的空間分布特性對(duì)湍流尺度的系統(tǒng)性分辨率直接相關(guān).通常情況下,雖然可以通過(guò)加密網(wǎng)格來(lái)調(diào)控這兩類模型對(duì)湍流尺度的解析能力,但在近壁面網(wǎng)格較密的條件下,會(huì)出現(xiàn)對(duì)網(wǎng)格空間分布特性過(guò)于敏感的缺陷,導(dǎo)致近壁面RANS 區(qū)域被破壞,而網(wǎng)格空間分布尺度還沒(méi)有小到足夠進(jìn)行大渦模擬的程度[25].

本文采用WM-HRL 模型的系列數(shù)值模擬結(jié)果表明,在同一套網(wǎng)格系統(tǒng)下,通過(guò)改變WM-HRL模型中“兩個(gè)邊界位置”及“三個(gè)區(qū)域的湍動(dòng)能解析度”的組合,既可通過(guò)數(shù)值模擬獲得與Parnaudeau等[18]實(shí)驗(yàn)結(jié)果一致的圓柱后方各站位處的雷諾應(yīng)力剖面,也可以通過(guò)數(shù)值模擬獲得與Lourenco 和Shih[27]實(shí)驗(yàn)結(jié)果一致的圓柱后方各站位處的雷諾應(yīng)力剖面.

研究發(fā)現(xiàn),對(duì)本文所構(gòu)建的WM-HRL 模型,在RANS/LES 混合轉(zhuǎn)換區(qū)邊界位置及其兩個(gè)指標(biāo)參數(shù)nrk1-Q和nrk2-Q準(zhǔn)則取值的某些組合下,計(jì)算所得二階統(tǒng)計(jì)量與相關(guān)文獻(xiàn)的實(shí)驗(yàn)結(jié)果相比,存在較大的誤差,其主要原因如下.

第一,在本文所構(gòu)建的WM-HRL 模型中,其RANS 模型采用的是SSTk-ω模型.該湍流模型在邊界層內(nèi)為標(biāo)準(zhǔn)的Wilcoxk-ω兩方程模型(簡(jiǎn)稱SKW 模型).對(duì)WM-HRL 模型,其RANS 模型必須直接作用于圓柱近壁面處,但對(duì)SKW 模型,會(huì)存在如下缺陷:包括可能會(huì)過(guò)大地預(yù)測(cè)回流區(qū)的長(zhǎng)度,還可能會(huì)過(guò)大地預(yù)測(cè)壁面剪切應(yīng)力的值及過(guò)小地預(yù)測(cè)壁面湍動(dòng)能的值[44].

第二,在本文所構(gòu)建的WM-HRL 模型中,LES模型采用的是經(jīng)典的Smargorinsky 型亞格子模型.該亞格子模型雖然概念簡(jiǎn)單且數(shù)值求解穩(wěn)定性高,但不僅會(huì)在近壁面產(chǎn)生過(guò)大的數(shù)值耗散,而且定常的Smagorinsky 常數(shù)難以描述復(fù)雜時(shí)空變化下的湍流場(chǎng).

為克服上述兩個(gè)缺陷,對(duì)RANS 模型可改用Peng等[44]提出的低雷諾數(shù)Wilcoxk-ω兩方程模型的修正版本.Peng等[44]的研究表明,該修正版本可有效地提高其對(duì)近壁面流的解析能力,包括可減小近壁面處ω的值及增大近壁面處k的值等,并具有更好地模擬分離、回流及再附等復(fù)雜流動(dòng)現(xiàn)象的能力.

另一方面,對(duì)LES 模型可改用動(dòng)態(tài)(Dyn)亞格子模型.Tremblay[8]的研究表明,對(duì)雷諾數(shù)Re=3900 下的圓柱繞流場(chǎng)問(wèn)題,與基于經(jīng)典Smagorinsky 型亞格子渦黏系數(shù)的LES 模型相比,采用基于動(dòng)態(tài)(Dny)亞格子渦黏系數(shù)的LES 模型,可以更好地獲得圓柱繞流場(chǎng)相關(guān)統(tǒng)計(jì)量的計(jì)算結(jié)果.這些正是本文作者團(tuán)隊(duì)將在下一階段重點(diǎn)開(kāi)展的相關(guān)研究工作.

4.3 相干結(jié)構(gòu)頻譜及流場(chǎng)特征

首先研究剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生區(qū)域的范圍問(wèn)題.為此,選取Case CU 和Case CV 兩種工況進(jìn)行分析.研究對(duì)象為這兩種工況下,在圖3 中P1—P12 這12 個(gè)監(jiān)測(cè)點(diǎn)處橫向脈動(dòng)速度的Lomb 譜特征,結(jié)果如圖11 和圖12 所示.

圖11 在Case CU 工況下,在P1—P12 監(jiān)測(cè)點(diǎn)處橫向脈動(dòng)速度的Lomb譜Fig.11.Lomb spectrums of the cross-stream fluctuation velocities at different probes P1-P12 for the Case CU.

圖12 在Case CV 工況下,在P1—P12 監(jiān)測(cè)點(diǎn)處橫向脈動(dòng)速度的Lomb譜Fig.12.Lomb spectrums of the cross-stream fluctuation velocities at different probes P1-P12 for the Case CV.

由圖11 和圖12 可見(jiàn),在12 個(gè)監(jiān)測(cè)點(diǎn)P1—P12處橫向脈動(dòng)速度的Lomb 譜中,均出現(xiàn)一個(gè)頻率相對(duì)較低的譜峰,該峰值所對(duì)應(yīng)頻率與渦泄頻率一致,其無(wú)因次頻率≈0.22.同時(shí),在12個(gè)監(jiān)測(cè)點(diǎn)P1—P12處橫向脈動(dòng)速度的Lomb 譜中,還會(huì)出現(xiàn)一個(gè)與該頻率峰值相對(duì)應(yīng)的倍頻()成分.前者為圓柱體升力脈動(dòng)的主頻成分,而后者為圓柱體阻力脈動(dòng)的主頻成分.

對(duì)Case CU 工況,在測(cè)點(diǎn)P4—P10 處橫向脈動(dòng)速度的Lomb 譜中(見(jiàn)圖11),還可以觀察到一個(gè)頻率更高的譜峰,該峰值所對(duì)應(yīng)頻率與剪切層小尺度K-H不穩(wěn)定性的頻率一致,其無(wú)因次頻率在1.33—1.47之間.進(jìn)一步的研究表明,在此工況下圓柱繞流場(chǎng)中剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生在對(duì)數(shù)律層中 89.4 ≤y+≤253.5 的區(qū)域.

對(duì)Case CV 工況,在測(cè)點(diǎn)P4—P7 處橫向脈動(dòng)速度的Lomb 譜中(見(jiàn)圖12),也可以觀察到一個(gè)頻率更高的譜峰,該峰值所對(duì)應(yīng)頻率也與剪切層小尺度K-H 不穩(wěn)定性的頻率一致,其無(wú)因次頻率也在1.43—1.44 之間.進(jìn)一步的研究表明,在此工況下圓柱繞流場(chǎng)中剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)發(fā)生在對(duì)數(shù)律層中 89.4 ≤y+≤167.4 的區(qū)域.

討論其他6 個(gè)監(jiān)測(cè)點(diǎn)P13—P18 處流向和橫向脈動(dòng)速度的Lomb 譜特征,計(jì)算工況選取Case CU 和Case CV 兩種情況,結(jié)果如圖13 所示.其中,監(jiān)測(cè)點(diǎn)P13 和P14 分別位于剪切層偏上及偏下的位置,監(jiān)測(cè)點(diǎn)P15 位于剪切層脫落區(qū),監(jiān)測(cè)點(diǎn)P16 位于回流區(qū)內(nèi)的尾流中線上,而監(jiān)測(cè)點(diǎn)P17 和P18 位于回流區(qū)外的尾流中線上.

圖13 在Case CU (第1 和第2 列)和Case CV (第3 和第4 列)工況下,在P13—P18 監(jiān)測(cè)點(diǎn)處流向(第1 和第3 列)及橫向(第2 和第4 列)脈動(dòng)速度的Lomb譜Fig.13.Lomb spectrums of the stream-wise (from the first to third rows) and cross-stream (from the second to fourth rows) fluctuation velocities at different probes P13-P18 for the Case CU (from the first to second rows) and the Case CV (from the third to fourth rows).

由圖13 可見(jiàn),在Case CU 和Case CV 兩種工況下,對(duì)P13—P18 這6 個(gè)監(jiān)測(cè)點(diǎn),無(wú)論是在其流向脈動(dòng)速度的Lomb 譜中,還是在其橫向脈動(dòng)速度的Lomb 譜中,均清晰地出現(xiàn)兩個(gè)頻率較低的譜峰,一個(gè)譜峰對(duì)應(yīng)渦泄頻率(),另一個(gè)譜峰對(duì)應(yīng)渦泄頻率的倍頻().

對(duì)位于剪切層脫落區(qū)的監(jiān)測(cè)點(diǎn)P15 的情況,在Case CU 和Case CV 兩種工況下,無(wú)論是在其流向脈動(dòng)速度的Lomb 譜中,還是在其橫向脈動(dòng)速度的Lomb 譜中,除了均清晰地出現(xiàn)單頻頻率及其倍頻的譜峰外,還可觀察到一個(gè)三倍頻率的譜峰.

對(duì)位于其他區(qū)域的監(jiān)測(cè)點(diǎn)P16—P18 的情況,在Case CU 和Case CV 兩種工況下,從其流向脈動(dòng)速度的Lomb 譜中只能清晰地看到單頻頻率及其倍頻的譜峰,而從其橫向瞬態(tài)速度的Lomb 頻譜圖中則只能清晰地看到單頻頻率及其三倍頻率的譜峰,但沒(méi)有出現(xiàn)明顯的二倍頻率的譜峰.

由圖13 進(jìn)一步可見(jiàn),在Case CU 和Case CV兩種工況下,對(duì)位于剪切層偏上方的監(jiān)測(cè)點(diǎn)P13,無(wú)論是在其流向脈動(dòng)速度的Lomb 譜中,還是在其橫向脈動(dòng)速度的Lomb 譜中,除了均清晰地出現(xiàn)與渦泄相關(guān)的單頻頻率及倍頻頻率的譜峰外,還可清晰地觀察到一個(gè)頻率更高的譜峰,其對(duì)應(yīng)頻率與剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的頻率一致,其無(wú)因次頻率在1.43—1.44 之間.

在Case CU 和Case CV 兩種工況下,對(duì)位于剪切層偏下方的監(jiān)測(cè)點(diǎn)P14,在其流向脈動(dòng)速度的Lomb 譜圖中,只清晰地出現(xiàn)了與渦泄相關(guān)的單頻頻率及倍頻頻率的譜峰.然而,在其橫向脈動(dòng)速度的Lomb 譜圖中,除了均清晰地出現(xiàn)與渦泄相關(guān)的單頻頻率及倍頻頻率的譜峰外,還可清晰地觀察到一個(gè)與剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)頻率一致的譜峰,其無(wú)因次頻率也在1.43—1.44 之間.

最后,從圖13 可見(jiàn),在Case CU 和Case CV 兩種工況下,對(duì)監(jiān)測(cè)點(diǎn)P13 和P14,無(wú)論是在其流向脈動(dòng)速度的Lomb 頻譜圖中,還是在其橫向脈動(dòng)速度的Lomb 譜中,均沒(méi)有出現(xiàn)明顯的三倍頻率的譜峰.

對(duì)雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)中出現(xiàn)三倍頻率的現(xiàn)象,在Parnaudeau等[18]的DNS 及Pereira[2]等的PANS 計(jì)算中都有發(fā)現(xiàn),而在Lehmkuhl等[10]采用DNS 及Tremblay[8]和Breuer[15]采用LES 的計(jì)算中則沒(méi)有發(fā)現(xiàn).Ong 和Wallace[45]指出,流體行為中的峰值現(xiàn)象可以通過(guò)線性穩(wěn)定性理論進(jìn)行預(yù)測(cè),用于研究各種不穩(wěn)定性問(wèn)題.然而,這個(gè)理論是否能夠適用于本文所觀察到的這種復(fù)雜現(xiàn)象,或者本文算例所出現(xiàn)的這類復(fù)雜頻譜現(xiàn)象是否還存在其他影響因素,值得在后續(xù)工作中做進(jìn)一步的研究.

最后討論圓柱繞流場(chǎng)中兩類相干結(jié)構(gòu)的流場(chǎng)特征,計(jì)算工況選取Case AU—DU 及Case AV—DV 兩類情況,結(jié)果如圖14 和圖15 所示.其中,圖14 的左列為x3/D=π/2 斷面上展向渦量ω3=?u2/?x1-?u1/?x2的云圖,右列為流向速度云圖,且圖中白色垂直虛線為回流區(qū)末端位置(即u1/U∞=0位置).圖15 為沿圓柱體長(zhǎng)度的三維展向渦量云圖.

圖14 在Case AU—DU(前4 行)和Case AV—DV(后4 行)工況下,圓柱繞流渦量(左)及流向速度(右)云圖Fig.14.Contours of the span-wise vorticity (left) and stream-wise velocity (right) for both Case AU-DU (the first four lines) and Case AV-DV (the last four lines).

圖15 在Case AU—DU (左)和Case AV—DV (右)工況下,圓柱繞流展向三維渦量云圖Fig.15.Contours of the three-dimensional span-wise vorticities both Case AU-DU (left) and Case AV-DV (right).

在Case AU—DU 及Case AV—DV 兩類工況下,從圖14 (左)均清晰地可見(jiàn)附著于圓柱上下側(cè)壁的兩個(gè)層流分離剪切層結(jié)構(gòu),表現(xiàn)為兩條位置狹小的窄帶結(jié)構(gòu).同時(shí),從圖14 (右)則可清晰地觀察到圓柱尾部因流動(dòng)分離而形成的兩個(gè)回流區(qū)分支結(jié)構(gòu).具體來(lái)說(shuō),第1 個(gè)分支結(jié)構(gòu)發(fā)生在x1/D=1.06 處平均流向速度剖面為U 形的情況,回流區(qū)長(zhǎng)度較長(zhǎng);而第2 個(gè)分支結(jié)構(gòu)發(fā)生在x1/D=1.06 處平均流向速度剖面為V 形的情況,回流區(qū)長(zhǎng)度較短.

結(jié)合圖11 和圖12 的結(jié)果進(jìn)一步可知,在Case AU—DU 及Case AV—DV 兩類工況下,兩個(gè)剪切層中小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的觸發(fā)位置都基本相同,均約位于y+=89.4 處,但兩個(gè)剪切層中小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的消逝位置不同.在Case AU—DU 工況下,其消逝位置約為y+=253.5.而在Case AV—DV 工況下,其消逝位置約為y+=167.4 .

這意味著,當(dāng)x1/D=1.06 處平均流向速度剖面為U 形時(shí),剪切層中小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的長(zhǎng)度較大,而且相應(yīng)的回流區(qū)長(zhǎng)度也較長(zhǎng).而當(dāng)x1/D=1.06 處平均流向速度剖面為V 形時(shí),剪切層中小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的長(zhǎng)度較小,而且相應(yīng)的回流區(qū)長(zhǎng)度也較短.因此,對(duì)雷諾數(shù)Re=3900 下圓柱繞流中出現(xiàn)兩類回流區(qū)分支結(jié)構(gòu),以及在圓柱尾部近壁面處平均流向速度出現(xiàn)U 形和V 形兩類剖面結(jié)構(gòu)形狀的形成機(jī)理,可能與WMHRL 模型中兩個(gè)邊界位置及3 個(gè)區(qū)域的湍動(dòng)能解析度發(fā)生變化時(shí),該湍流模型對(duì)圓柱繞流剪切層中小尺度層流分離流動(dòng)的解析能力不同有關(guān).

對(duì)DNS 模型,在同一套網(wǎng)格系統(tǒng)下,其對(duì)圓柱繞流剪切層中小尺度層流分離流動(dòng)的解析能力是確定的,因此只能獲得某一種回流區(qū)分支結(jié)構(gòu)及與之對(duì)應(yīng)的圓柱尾部近壁面平均流向速度剖面的某一種形狀.但對(duì)不同的網(wǎng)格系統(tǒng),由于其對(duì)圓柱繞流剪切層中小尺度層流分離流動(dòng)的解析能力相應(yīng)地會(huì)不同,進(jìn)而可能會(huì)出現(xiàn)在不同網(wǎng)格系統(tǒng)下獲得不同的回流區(qū)分支結(jié)構(gòu)及與之對(duì)應(yīng)的圓柱尾部近壁面流向速度剖面的U 形或V 形結(jié)構(gòu).

對(duì)PANS 及DES 類模型,它們可選用不同的RANS 湍流模型,也可以改變湍流模型中的相關(guān)模型參數(shù).對(duì)DES 類模型,其執(zhí)行RANS/LES 之間轉(zhuǎn)換過(guò)渡的方式還可以不同.另一方面,對(duì)LES 模型,其可以選用經(jīng)典Smagorinsky 及Dyn等[9]不同的亞格子渦黏模型.這些復(fù)雜的因素都可能會(huì)導(dǎo)致即使在同一套網(wǎng)格系統(tǒng)下,使得這些湍流模型對(duì)圓柱繞流剪切層中小尺度層流分離流動(dòng)的解析能力不同,進(jìn)而可能會(huì)獲得不同的回流區(qū)分支結(jié)構(gòu)及與之對(duì)應(yīng)的圓柱尾部近壁面平均流向速度剖面的U 形或V 形結(jié)構(gòu).

亞臨界雷諾數(shù)區(qū)圓柱繞流中主要有兩類相干結(jié)構(gòu).其中,一類為剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu),而另一類為大尺度的Karman 渦街結(jié)構(gòu).為進(jìn)一步闡述它們的形成機(jī)理及其特征,在圖15中,給出了 在Case AU—DU 和Case AV—DV兩類工況下,圓柱繞流展向三維渦量云圖的數(shù)值結(jié)果.

由圖15 并結(jié)合圖14 可見(jiàn),由于剪切層兩側(cè)的速度差誘發(fā)K-H 不穩(wěn)定性,使剪切層失穩(wěn)并出現(xiàn)滾卷現(xiàn)象.一方面,這種現(xiàn)象會(huì)在剪切層中導(dǎo)致大量小尺度的旋渦產(chǎn)生于此區(qū)域,在橫向瞬態(tài)速度的Lomb 譜中表現(xiàn)為寬頻信號(hào)的凸起(見(jiàn)圖11 和圖12).另一方面,這種現(xiàn)象還會(huì)在圓柱尾跡區(qū)形成另一類大尺度的旋渦結(jié)構(gòu),即Karman 渦街結(jié)構(gòu),在橫向瞬態(tài)速度的Lomb 譜中表現(xiàn)為頻率相對(duì)較低的窄頻信號(hào)的凸起(見(jiàn)圖11—圖13).

由圖15 還可以清晰地觀察到,附著于圓柱壁面的剪切層伴隨著其失穩(wěn)、滾卷等復(fù)雜現(xiàn)象在圓柱尾流區(qū)形成流向與展向相嵌套的三維渦體結(jié)構(gòu).由于本文所構(gòu)建的WM-HRL 模型具有能夠精細(xì)解析湍流譜結(jié)構(gòu)的能力(見(jiàn)圖11—圖13),進(jìn)而不僅能夠精細(xì)地解析各種小尺度渦結(jié)構(gòu)(如K-H 不穩(wěn)定性結(jié)構(gòu)等),又能夠精細(xì)地解析各種大尺度渦結(jié)構(gòu)(如Karman 渦街等).因此,該WM-HRL 模型表現(xiàn)出預(yù)期的相當(dāng)于LES 甚至DNS 對(duì)各種小尺度及大尺度運(yùn)動(dòng)解析的能力.

5 結(jié)論

本文通過(guò)修改SST-k-ω湍流模型中k方程的色散項(xiàng),構(gòu)造了一種新的具有壁面模化能力的鈍體繞流混合RANS/LES 模型,即WM-HRL 模型.該模型在鈍體近壁面的某個(gè)計(jì)算區(qū)域內(nèi)采用RANS 模式,在經(jīng)歷一個(gè)RANS/LES 混合轉(zhuǎn)換過(guò)渡區(qū)后,在其余計(jì)算區(qū)域則采用LES 模式,且該LES 模式等效為經(jīng)典Smagorinsky 亞格子模型.

通過(guò)構(gòu)造一個(gè)僅與當(dāng)?shù)鼐W(wǎng)格空間分布尺寸相關(guān)的湍動(dòng)能解析度指標(biāo)參數(shù)rk,提出了確定該WM-HRL 模型兩個(gè)邊界位置和3 個(gè)區(qū)域湍動(dòng)能解析度的nrk1-Q和nrk2-Q準(zhǔn)則.其中,nrk1-Q準(zhǔn)則用于確定WM-HRL 模型中RANS 模式結(jié)束邊界位置,及其在RANS 區(qū)湍動(dòng)能的最大解析度(即rk1),而nrk2-Q準(zhǔn)則用于確定WM-HRL 模型中LES 模式啟動(dòng)邊界位置及其在LES 區(qū)對(duì)湍動(dòng)能的最小解析度(即rk2).

在此基礎(chǔ)上,本文構(gòu)造了一個(gè)新的具有雙重保護(hù)功能的混合函數(shù)fnd.第一重保護(hù)通過(guò)nrk1-Q準(zhǔn)則實(shí)現(xiàn),保證在rk≥rk1時(shí)WM-HRL 為完全RANS模式.通過(guò)第一重保護(hù)層的設(shè)置,可使WM-HRL既具有延遲脫體渦模擬(DDES)的能力,又具有壁面?；?WM)的能力.在RANS 區(qū)最大湍動(dòng)能解析度指標(biāo)rk1的合理設(shè)置下,可避免MSD 和GIS的問(wèn)題.

第二重 保護(hù)通過(guò)nrk2-Q準(zhǔn)則實(shí) 現(xiàn),可保證在rk≤rk2時(shí)WM-HRL 為完全LES 模式.在此第二重保護(hù)層設(shè)置的條件下,可通過(guò)rk2值的合理設(shè)置,一方面在進(jìn)行計(jì)算網(wǎng)格設(shè)置時(shí)即可實(shí)現(xiàn)對(duì)LES 啟動(dòng)邊界位置的預(yù)先設(shè)定,另一方面還可自定義LES 區(qū)中WM-HRL 對(duì)鈍體繞流各種復(fù)雜湍流場(chǎng)的解析能力.

在合理設(shè)置nrk1-Q和nrk2-Q準(zhǔn)則的條件下,既可保證WM-LES 模型在RANS/LES 混合轉(zhuǎn)換區(qū)即可激活LES 模式,同時(shí)還可保證在完全LES 區(qū)即可完全激活LES 模式,進(jìn)而避免LLM 的問(wèn)題.

基于該WM-HRL 模型,本文以雷諾數(shù)Re=3900 下圓柱繞流為對(duì)象,開(kāi)展了系列數(shù)值模擬分析與評(píng)估,得到如下主要結(jié)論.

1)即使當(dāng)LES 的啟動(dòng)邊界位于圓柱繞流剪切層小尺度K-H 不穩(wěn)定性區(qū)內(nèi)時(shí),在合理設(shè)置的nrk2-Q準(zhǔn)則(一般需要rk2≤0.5)的條件下,該WMHRL 模型也能獲取準(zhǔn)確的渦泄頻率、分離角、阻力系數(shù)及基準(zhǔn)線壓力系數(shù)等統(tǒng)計(jì)量信息.

2)為使WM-HRL 模型能夠準(zhǔn)確獲取圓柱繞流剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)特征及其頻譜特性,以及回流區(qū)兩個(gè)分支結(jié)構(gòu)的長(zhǎng)度信息,可將LES 模式的啟動(dòng)邊界位置均設(shè)置在過(guò)渡層內(nèi),而RANS 模式的結(jié)束邊界位置則既可在黏性底層也可在過(guò)渡層內(nèi),并且使rk2≤0.2,即在對(duì)數(shù)律層區(qū)LES 對(duì)湍動(dòng)能具有至少具有80%的解析能力.

3)在上述2)的條件下,WM-HRL 模型能夠獲取與相關(guān)實(shí)驗(yàn)精度相當(dāng)?shù)膱A柱繞流一階和二階統(tǒng)計(jì)量的信息,以及剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的無(wú)因次頻率fˉkh及回流區(qū)兩個(gè)分支結(jié)構(gòu)的長(zhǎng)度信息,并且當(dāng)回流區(qū)長(zhǎng)度的數(shù)值結(jié)果與Parnaudeau等[18]的實(shí)驗(yàn)結(jié)果一致,在x1/D=1.06 處的平均流向速度剖面為U 形,而當(dāng)回流區(qū)長(zhǎng)度的數(shù)值結(jié)果與Lourenco 和Shih[27]的實(shí)驗(yàn)結(jié)果一致,則為V 形.

特別地,在本文的系列數(shù)值模擬中,通過(guò)兩個(gè)邊界位置及3 個(gè)區(qū)域的湍動(dòng)能解析度的各種不同組合,在僅用一套網(wǎng)格系統(tǒng)的條件下,利用新構(gòu)造的WM-HRL 模型,針對(duì)雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)問(wèn)題,獲得了兩類長(zhǎng)度不同的剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)特性,具體如下.

第一類剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的長(zhǎng)度較長(zhǎng),此時(shí)回流區(qū)長(zhǎng)度也較長(zhǎng),與之相應(yīng)的圓柱尾部近壁面處的平均流向速度剖面為U 形.第二類剪切層小尺度K-H 不穩(wěn)定性結(jié)構(gòu)的長(zhǎng)度較短,此時(shí)回流區(qū)長(zhǎng)度也較短,與之相應(yīng)的圓柱尾部近壁面處的平均流向速度剖面為V 形.

關(guān)于雷諾數(shù)Re=3900 下圓柱繞流場(chǎng)中回流區(qū)的兩個(gè)分支結(jié)構(gòu),以及與之相應(yīng)的圓柱尾部近壁面處U 形和V 形兩類平均流向速度剖面形狀等問(wèn)題,已經(jīng)成為困擾學(xué)界三十余年的難題.本項(xiàng)研究為后續(xù)更深層次認(rèn)識(shí)這一難題提供了一種新的手段.